Pre-release!

在上一篇中，我们通过Sentieon得到了vcf文件，接下来要对vcf进行过滤，以适用不同的分析要求。首先是用于构建系统发生树的过滤流程

进化树

King 亲缘关系判断

这一步在大多数植物群体分析的文章中并没有应用到，只是在一篇文章中偶然看到，抱着试试的心态看看结果，后来发现结果还不错，就保留了这个步骤。所以这一步是可选的。

King https://www.kingrelatedness.com/ 主要是在GWAS中用于推断亲缘关系的。在查找相关教程的时候看到有的说法是由于亲缘关系较近的个体在后续的计算中会导致误差，例如群体历史推算和模拟等，因此需要删除一些亲缘关系较近的个体 https://github.com/shangshanzhizhe/Work_flow_of_population_genetics/blob/master/Work_flows/king.md 因此决定尝试一下。

king查看亲缘关系,要先把vcf转成plink的格式

King的输入文件要求是PLINK的二进制文件，因此需要先将vcf转成plink的格式，主要有两个步骤，第一个是通过vcftools把vcf转成map和ped格式，第二个是通过plink生成bed bim fam 和 nosex

vcftools --gzvcf ../../00.vcf/406cannabis.jc.vcf.gz --plink --chrom-map chrom.txt --out 406cannabis.jc

命令中通过–chrom-map参数确定染色体编号，chrom.txt的格式参考vcftools https://vcftools.github.io/man_latest.html 官网，第一列为vcf文件中的染色体编号，第二列为染色体的序号，这里的染色体序号要跟给的格式一致，不能自己写chrxx，后面会识别不出来。两列之间用制表符分隔。

#chrom.txt
Chr01	1
Chr02	2
Chr03	3
Chr04	4
Chr05	5

这一步生成map文件和ped文件。map有四列，参考 https://zzz.bwh.harvard.edu/plink/data.shtml#map ped文件格式参考 https://zzz.bwh.harvard.edu/plink/data.shtml#ped。

plink --file 406cannabis.jc --make-bed --noweb --out 406cannabis.jc

这一步生成bed bim fam nosex四个文件。https://zzz.bwh.harvard.edu/plink/data.shtml#bed bed是二进制文件，bed会将谱系/表型信息存储在fam文件中(fam其实就是前六列的ped)，还会生成一个包含等位基因名称的bim文件

使用king计算样本的亲缘关系

king主要有两种算法，一种是kinship,一种是ibdseg,这里主要介绍使用kingship计算亲缘关系，然后根据kinship的值区分关系等级

king -b 406cannabis.jc.bed --kinship --prefix kinship --cpus 16

根据kinship.kin0文件中中最后一列kinship的值划分亲缘关系。kinship >0.354, [0.177, 0.354], [0.0884, 0.177], [0.0442, 0.0884] 分别对应duplicate/MZ twin, 1st-degree, 2nd-degree 和 3rd-degree。我们取2级以上的关系认为是亲缘关系很近的个体，摒弃。

鉴定单拷贝基因

找单拷贝基因,必须用CDS

这里采取的是用参考基因组CDS序列自己比对自己的方式鉴定单拷贝基因。以CDS序列构建blast数据库，然后自己比对回自己，在结果中挑选只能比对上自己的基因，就认为是单拷贝基因。

构建blast库

makeblastdb -in modified.cds.fa -dbtype nucl -parse_seqids -out CsLZ.cds

自己比对上自己

blastn -query modified.cds.fa -db CsLZ.cds -outfmt 6 -num_threads 128 -out blast.result

用脚本提取单拷贝基因的bed格式

get_SCG_bed.py

vcftools过滤

本文过滤SNP的时候主要有几个参考标准：等位基因的数量、位点深度、质量、次等位基因频率MAF和缺失率（位点和个体）。由于原始vcf文件过大，在统计这些参数的时候都要跑很久或者占用大量内存，因此可以先过滤几个固定的参数，例如个体缺失率，MAF和等位基因的数量，再根据统计的结果过滤其余的参数。vcftools提供了–geno-depth –freq –counts –site-quality –missing-site –missing-indv等参数统计这些信息，这些参数不能一次性在一个命令中同时添加，需要分别提交，例如计算每个位点在每个样本中的深度

vcftools --gzvcf ../../00.vcf/406cannabis.ht.vcf.gz --geno-depth

–freq 生成frq文件, N_ALLELES是等位基因的数量，是后续的–min-alleles和–max-alleles的判断标准。最后面是每个等位基因的频率。alleles的指定实际是看AC的数量（不是数值），当在ALT列有两种及以上的突变情况（即逗号分隔）时，alleles被认为大于2，这可能会被过滤掉

–geno-depth生成gdepth文件，是每个位点在每个样本中的深度，与miss结合是–minDP和–maxDP的判断标准。对于深度的过滤，–minDP必须搭配–max-missing进行过滤。具体的操作逻辑是这样的，–max-missing指定允许出现的缺失（./.）比例，在非缺失样本中，每个样本的深度必须符合DP的取值范围。当–max-missing设置为0时，不会进行过滤，当不指定–max-missing时，–minDP将无效。–meanDP是不参与判断每个样本各自的深度的，只计算平均深度。假设有100个样品，你希望平均每个样品在2～3x那么应该设置–min-meanDP 200 –max-meanDP 300

–missing-indv生成imiss文件，是每个样本在全部位点的缺失率，结合–remove过滤高缺失的样本得到的out.imiss里面，找最后一列大于0.5的，剔除掉 awk ‘$5 > 0.5’ out.imiss | cut -f1 > lowDP.indv #注意外类群要保留着,把亲缘关系近的也加进来 –missing-site生成lmiss文件，是每个位点在全部样本的缺失率，是–max-missing的判断标准

–site-quality生成lqual文件，是每个位点的质量，对应vcf文件中的QUAL列，是–minQ的判断标准

统计好各项指标之后可以对他们在snp中的分布情况绘制直方图，这样可以方便的确定参数的选择，之后便可以加上单靠本基因的bed文件进行过滤

vcftools --gzvcf ../../00.vcf/406cannabis.ht.vcf.gz --remove lowDP.indv --bed ../../01.SCG/CsLZ.SCG.bed --max-alleles 2 --min-alleles 2 --maxDP 50 --minDP 7 --minQ 30 --max-missing 0.8 --maf 0.05 --recode --recode-INFO-all --stdout | bgzip -c > 406cannabis.ht.DP7.mis8.vcf.gz

vcf转phy

使用vcf2phylip.py把vcf转换成phy文件

python3 vcf2phylip.py --input ../02.filtsnp/haplotyper/406cannabis.ht.DP7.mis8.vcf.gz -f -n -r --output-folder ht/

Pre-release!#

进化树#

King 亲缘关系判断#

鉴定单拷贝基因#

vcftools过滤#

vcf转phy#

Pre-release!

进化树

King 亲缘关系判断

鉴定单拷贝基因

vcftools过滤

vcf转phy