二.变异结果注释与统计
书接上回https://www.jianshu.com/p/ab6a35502786
如果是使用转录组重测序的话需要修改过滤条件，会发现大多数变异都发生在基因区，其他的差别都不大。比对参考基因组之后进行变异检测，同样可以用GATK的流程来做。会多一个cigar字符串的处理步骤。因为在基因组上会比对出跨区域的reads，所以要把这样一整条reads（中间是非基因区）当做两条来处理。
三代测序由于错误率高，本身是不适合用于变异检测的，尤其是SNP和InDel。
我们要统计的SNP变异类型包括：
1.变异位点是否存在于基因区or上下2kbp范围内，如果在2k之内的话可能存在一些转录元件。太远的话就可以视为基因间区了
2.如果在基因区的话，就要区分是否在intron区。在UTR区的变异影响不大，而在CDS区的话又要进一步进行分析。
3.如果变异影响了氨基酸类型的话，就会对基因功能（蛋白）造成影响了。
4.更加严重的变异类型有：对起始密码子的突变，突变后形成终止密码子（提前终止），干扰可变剪切位点。
在InDel变异中，造成的非3倍数的移码突变会非常严重。即使是3倍数的移码突变也是较严重的突变。
annovar 注释
• 文件准备：
基因组的位置文件gtf
变异位点文件vcf
gtf第八列是相位（phase）信息，意为处在CDS交界处的碱基偏移
SNP 和 INDEL
# 格式转换
$ gtfToGenePred -genePredExt \
genome.gtf \
genome_refGene.txt#下划线后的东西不能变
# 提取转录本（来自annovar软件）
$ retrieve_seq_from_fasta.pl --format refGene \
--seqfile genome.fasta \
--out genome_refGeneMrna.fa \
genome_refGene.txt
输出得到fa文件之后进行注释。annovar可以基于基因位置进行注释，也可以基于指定区域注释。
#格式转换
$ convert2annovar.pl \
-format vcf4 \ #输入文件格式
-allsample \ #不加的话只有第一个样本的信息
-withfreq \ #输出alt频率
all.filtered.snp.vcf \ #输入文件，来自GATK的输出
>all.snp.vcf.annovar.input #输出文件
#注释
$ annotate_variation.pl \
-geneanno \ #基于基因进行注释
--neargene 2000 \ #需要考虑的gene上下游范围
-buildver genome \ #基因组数据库名称
-dbtype refGene \ #数据库类型
-outfile all.anno.snp \ #输出文件前缀
-exonsort \ # 结果按exon进行排序
all.snp.vcf.annovar.input ./#最后一项是数据库所在的目录
#统计基因组各区域中snp的个数
$ less -N all.anno.variant_function |cut -f 1|sort | uniq -c
分别会得到全基因组变异文件和外显子变异文件，后者详细记载有变异类型。InDel变异只需要改变文件的名字就可以了。最后得到的输出文件中，碱基的突变会不止一个，而变异类型的种类也会更多。
当一个变异可能属于多种变异类型的时候，annovar会按变异造成的后果严重性进行优先级排序，只输出最严重的。如果需要所有类型的输出结果可以用SnpEff软件或者加参数。
SV与CNV
可以继续用上一步构建好的数据库，命令几乎一样。
$ convert2annovar.pl -format vcf4 -allsample -withfreq all_SV.vcf >SV_annovar.input
$ annotate_variation.pl -geneanno --neargene 2000 -buildver genome -dbtype refGene -outfile SV_annovar.output -exonsort SV_annovar.input ./
三.群体结构分析
包括structure，PCA，进化树分析和LD decay等
01.snp数据过滤
需要过滤maf次等位基因频率，missing，还有LD过滤。前二者用于去掉低质量snp位点，后者用于生成符合哈温平衡的群体。
$ plink --vcf all.vcf \ # 输入文件
--geno 0.1 \ # 设置缺失率阈值
--maf 0.01 \ # 设置maf阈值
--biallelic-only strict \ # 去掉三等位位点
--out all.missing_maf \ # 输出文件前缀
--recode vcf-iid \ # 输出vcf格式
--allow-extra-chr \ # 允许其他类型的染色体
--set-missing-var-ids @:# \ # 给没有id的snp取一个id
--keep-allele-order #让plink不要给ref和alt换位置
用plink过滤maf和missing，也可以用vcftools来完成。
$ plink --vcf all.missing_maf.vcf \
--indep-pairwise 50 10 0.2 \ # LD过滤阈值
--out tmp.ld \#前缀
--allow-extra-chr \
--set-missing-var-ids @:#
$ plink --vcf all.missing_maf.vcf \
--extract tmp.ld.prune.in \ # 保留的SNP列表
--out all.LDfilter \
--recode vcf-iid \#输出文件的格式
--keep-allele-order \
--allow-extra-chr \
--set-missing-var-ids @:#
第一步输出是可以保留下来的的位点的id信息，所以上一步过滤的时候最好可以给它加一个id。第二步再利用位点id信息来提取对应的vcf文件。
显然，LD过滤会越来越少。如果不需要那么多snp位点去进行分析的话可以把阈值卡得更严格，从而节省计算资源。
02.进化树构建
• 文件准备
vcf 文件：all.LDfilter.vcf
$ run_pipeline.pl -Xms1G -Xmx5G \
-SortGenotypeFilePlugin \
-inputFile all.LDfilter.vcf \
-outputFile all.LDfilter.sort.vcf \
-fileType VCF
$ run_pipeline.pl -Xms1G -Xmx5G \
-importGuess all.LDfilter.sort.vcf \
-ExportPlugin \
-saveAs sequences.phy \
-format Phylip_Inter
phylip格式转换之前需要先排序。虽然表面上是一个pl程序，但其中有调用java的部分，所以仍然需要设置一些java参数
dnadist可以用于nj法交互式建树，不过由于群体遗传的数据量往往非常大，脚本式利用配置文件在后台运行会更加合适。
$ echo -e "sequences.phy\nY" > dnadist.cfg #-e打入换行符
$ dnadist < dnadist.cfg >dnadist.log ##距离矩阵
$ mv outfile infile.dist
$ echo -e "infile.dist\nY" > neighbor.cfg
$ neighbor < neighbor.cfg >nj.log
#距离矩阵
$ less infile.dist | tr '\n' '|'| sed 's/| / /g' | tr '|' '\n' >infile.dist.table
#去掉换行符
$ less outtree | tr '\n' ' '|sed 's/ //g' > outtree.nwk
dnadist发现outfile会报错而非覆盖，所以最好把它名字改掉。
这样进化树文件就做好了，再改用其他软件来作图。听说ggtree包很好看。
03.structure分析
# vcf格式转bed格式
$ plink --vcf ../00.filter/all.LDfilter.vcf \
--make-bed --out all --allow-extra-chr \
--keep-allele-order --set-missing-var-ids @:#
# 生成k=2~4的脚本
$ seq 2 4 | \
awk '{print "admixture --cv -j2 all.bed "$1" 1>admix."$1".log 2>&1"}'\
> admixture.sh
$ cat admixture.sh
admixture --cv -j2 all.bed 2 1>admix.2.log 2>&1
admixture --cv -j2 all.bed 3 1>admix.3.log 2>&1
admixture --cv -j2 all.bed 4 1>admix.4.log 2>&1
$ nohup sh admixture.sh &
#绘图
$ mkdir result
$ cp ./*.Q result/
$ less outtree.nwk |sed 's/:/\n/g'|awk -F '[,(]' '{print $NF}'|awk '$1 !~")"' >sample.pop.order
$ Rscript draw_admixture.R result all.fam sample.pop.order structure
在log文件里面看各K值的CV error
Rscript第一个参数是Q矩阵；第二个参数来自格式转换时生成的fam文件，也可能用别的含样品顺序的文件；第三个参数是输出的样品顺序文件；最后是输出的前缀。可以基于树文件中给出的大概分类情况来调整顺序。
其他的看这篇文章吧https://www.jianshu.com/p/21c2a683c6fb
04.PCA分析
• 文件准备：
vcf 文件;all.LDfilter.vcf
分组信息表：sample.pop
plink --vcf all.LDfilter.vcf \
--pca 10 \ # 主成分个数
--out PCA_out \
--allow-extra-chr --set-missing-var-ids @:#
利用输出的PCA_plink.out.eigenvec和分组信息，在R中进行绘图。如果要分析三个主成分的话可以输出两张图，一张3dPCA的效果并不好。
05.LDdecay和LDBlock
随着两个site在基因组上的距离变远，连锁也会变弱，二者间的R²会变小。将所有SNP的LD系数按距离求平均作为纵轴，将距离作为横轴
LD衰减距离取决于群体类型、世代间隔和染色体相对位置.
• 文件准备
亚群样品列表文件：
pop.SC.table
pop.YZR.table
$ head -n 3 pop.SC.table
GA0008
GA0035
GA0038
#两个亚群
$ PopLDdecay -InVCF all.vcf \ # 输入vcf文件
-SubPop pop.SC.table \ # 指定要分析的亚群
-MaxDist 500 \ # SNP对距离上限500kb
-OutStat pop.SC.stat
$ PopLDdecay -InVCF ../data/all.vcf -SubPop ../data/pop.YZR.table -MaxDist 500 -OutStat pop.YZR.stat
#多个亚群共同绘图
### 准备配置文件，可以从文件名中awk出来
$ ls pop.*.stat.gz |awk -F"." '{ print $0"\t"$2 }' > ld_stat.list
$ Plot_MultiPop.pl -inList ld_stat.list \
-output ld_stat.multi \
-keepR#用于生成R脚本
LDBlock用于展示染色体局部的相关性，所以需要自己先有一个关注的预选区域
$ LDBlockShow -InVCF all.missing_maf.vcf \
-OutPut LDBlockShow \ # 输出图片前缀
-Region 1:300000-500000 \ # 绘图范围
-SeleVar 2 \ # 绘制Rˆ2
-OutPdf \ # 输出pdf图片
-OutPng # 输出png图片
LDBlock展示了连锁遗传信息，一般会和曼哈顿图一起用
四.群体选择分析
01.π值，Fst，Tajima'D检验
• 文件准备
亚群样品列表文件 pop.SC.table pop.YZR.table 群体变异检测结果文件：all.vcf
#π值计算
$ vcftools --vcf all.vcf \
--window-pi 100000\ #设置窗口大小
--window-pi-step 10000\ #设置步长
--keep pop.SC.table \ #亚群样品列表
--out ./Pi.SC #输出文件前缀
###同理指定其他亚群
$ vcftools --vcf all.vcf --window-pi 100000 --window-pi-step 10000\ --keep pop.YZR.table --out ./Pi.YZR
#π值沿染色体的分布
$ ls Pi.*pi | awk -F "." '{print $2"\t" $0}' > Pi.list
接下来就可以利用Pi.list来进行CMplot绘图。染色体中间π值很低的区域往往是着丝粒区域。同样，ROD的曼哈顿图也可以用CMplot绘制。
#TajimaD
$ vcftools --vcf all.vcf \
--TajimaD 100000 \
--keep ../../data/pop.SC.table \
--out TajimaD.SC#前缀名称
$ vcftools --vcf all.vcf --TajimaD 100000 --keep pop.YZR.table --out TajimaD.YZR
$ ls TajimaD.*.Tajima.D |awk -F"." '{print $2"\t"$0}' > TajimaD.list
#Fst
$ vcftools --vcf all.vcf \
--fst-window-size $window \ # 窗口大小
--fst-window-step $step \ # 设置step距离
--weir-fst-pop pop.SC.table \
--weir-fst-pop pop.YZR.table \
--out ./Fst.pop1.pop2
一系列的统计值计算都是用vcftools完成的，在log中有Fst的平均值和加权值的统计，可以用加权平均数来进行CMplot绘图。输出文件为.windowed.weir.fst，两个样本分化越大，这个值也会越大。
接下来可以进行Fst和ROD的联合分析，以及Fst，D和π的联合绘图。
02.XP-CLR
这里使用后发表的Python版本的XPCLR软件
$ xpclr -I ../data/all.vcf \
-O Chr01.xpclr-python.out \
-C Chr01 \
-Sa pop.SC.table \#两个亚群信息
-Sb pop.YZR.table \
--rrate 1e-8 \#单位碱基对应的遗传距离
--ld 0.95 \
--maxsnps 200 \
--size 100000 \#滑窗
--step 20000
$ cut -f 2,3,4,12 Chr01.xpclr-python.out | awk 'NF == 4 ' > Chr01.xpclr-python.table
#提取两个chr中xpclr的结果(第12列）
$ xpclr -I ../data/all.vcf -O Chr02.xpclr-python.out -C Chr02 -Sa pop.SC.table -Sb pop.YZR.table --rrate 1e-8 --ld 0.95 --maxsnps 200 --size 100000 --step 20000
$ cut -f 2,3,4,12 Chr02.xpclr-python.out | awk 'NF == 4 ' > Chr02.xpclr-python.table
#如果是第一行，或者最后一列是xpclr，而且没有空值的列输出
$ cat Chr*.xpclr-python.table | awk '{if((NR==1 || $NF != "xpclr") && NF == 4 ){print $0 }}' > all.xpclr-python.table
最后的结果也可以在R中绘图，受到选择越强的区域xpclr值会越大。
