WGS分析笔记（4）—— Call SNVs/indels
最近一直忙期末汇报和脚本编写，没来的及接着往下写文章，年前把这一块写了，然后再往下的分析流程就比较特异性了，做一步写一步那种，相对于这种已经常规化的流程（除了一些细节上的差异，别的都是大同小异的了），再往下的可能就不是很常规的分析手段了，不同的实验室有不同的分析方法，希望能够有大佬提点意见，多多交流。
“突变”区别
对于WGS的数据，在处理完bam文件以后，就是call variations了，之后所有的工作其实都是针对variations进行分析，那么说到变异呢，其实中文的“变异”在英文里对应多个单词，经常看到傻傻分不清，这里稍微区别一下。
mutation：核苷酸序列的永久性改变（来源于ACMG），在人群中小于1%或5%（来源于北大生物信息学公开课）
polymorphism：在人群中频率超过1%（来源于ACMG），或超过5%（来源于北大生物信息学公开课）
variation/variant：以上两个的总和（来源于北大生物信息学公开课）
变异类型
我们一般所说的variations主要有四大类：SNV，INDEL，SV，CNV。
SNV即单核苷酸位点变异（single nucleotide variants）
INDEL即小片段插入缺失（insertion and deletion）
SV即结构变异（structural variation）
CNV即拷贝数变异（copy number variation）
SNV与SNP
说到SNV，可能大家也经常看到SNP（single nucleotide polymorphism），同样是混淆使用，这俩其实还是有一些区别的，看到上面对变异的区分，其实也能看出这俩的区别来：
一般SNP是二态的，SNV没有这样的限制，如果在一个物种中该单碱基变异的频率达到一定水平（1%）就叫SNP，而频率未知（比如仅仅在一个个体中发现）就叫SNV，SNV包含SNP。
变异简述
人基因组通常有4.1~5.0M的变异，但是99.9%都是由SNV和short indel造成。
通常，一个人全基因组内会有约 3.6~4.4 M 个 SNVs，绝大数（大于 95%）的高频（群体中等位基因频率大于 5%）的 SNP 在 dbSNP中有记录，高频的SNP一般都不是致病的主要突变位点。
通常，一个人全基因组内会有约 600K 的 Indel（<50bp的插入缺失为small indel）。
编码区或剪接位点处发生的插入缺失都可能会改变蛋白的翻译。
SV
结构变异指的是在基因组上一些大的结构性的变异，比如大片段丢失（deletion）、大片段插入（insertion）、大片段重复（duplication）、拷贝数变异（copy number variants）、倒位（inversion）、易位（translocation）。一般来说，结构变异涉及的序列长度在1kb到3Mb之间。结构变异普遍存在于人类基因组中，是个人差异和一些疾病易感性的来源。结构变异还可能导致融合基因的发生，一些癌症已经证实和结构变异导致的基因融合事件有关。
CNV
拷贝数变异指的是基因组上大片段序列拷贝数的增加或者减少，可分为缺失（deletion）和重复(duplication)两种类型，是一种重要的分子机制。CNV能够导致孟德尔遗传病与罕见疾病, 同时与包括癌症在内的复杂疾病相关，因此对于染色体水平的缺失、扩增的研究已经成为疾病研究热点。
以上数据在不同的数据库或文献上可能有所差异，但相去不远，具体的变异分类以及分布就不详述，可参考一些文献，下面说说怎么进行其中的SNV和indel的检测。
以下我将提供两种分析方式的脚本（DeepVariant以及bcftools）和流程，为啥没有GATK？因为我觉得黄老师的这篇GATK分析流程写得已经很好了，大家可以参考一下。
Bcftools
与旧版的samtools+bcftools不同，作者为了避免bcftools和samtools的版本不同导致的不兼容，新版的bcftools可以自己完成call snv/indel的工作。
使用bcftools进行变异检测，一般分为三部曲，分别为三个模块mpileup、call、filter，当然，我们一般也不会分三步进行操作，而是使用管道（pipeline）进行编写脚本，这样能减少产生一些不必要的过程文件，同时提高自动化和效率，下面是我的实际使用脚本。
$ bcftools mpileup --threads 12 -q 20 -Q 20 -Ou -f /your/path/of/reference /your/bamfile/after/sorted.merged.markdup | bcftools call --threads 12 -vm -Ov | bcftools filter --threads 12 -s FILTER -g 10 -G 10 -i "%QUAL>20 && DP>6 && MQ>50 && (DP4[2]+DP4[3])>4" > raw.tmp.vcf
## bcftools mpileup检测变异;
# --threads线程数
# -q表示reads比对质量选择，MAPQ，默认0;
# -Q表示reads碱基对质量选择，默认13
# -O表示输出格式,u表示未压缩bcf格式;
# -f参考序列位置
## bcftools call参数
# --threads线程;
# -v只输出变异位点;
# m为克服-c调用模型中已知的局限性(与-c冲突)而设计的多等位和罕见变异调用的替代模型;
# -O输出文件格式v未压缩vcf;
## bcftools filter筛选变异;-i只保留后面条件的;-s对不符合的变异打上标签
$ awk -F "\t" '{if($1~/#/){print}else if($7~/PASS/){print}}' raw.tmp.vcf > var.flt.vcf
# 将标为低质量的变异去掉
有几个点值得讨论一下，首先是mpileup的-q参数，这个和之前提到的samtools view的-q是一样的，前面的文章有大篇幅说过这个MAPQ的数值，20翻译过来的意思其实就是比对正确率99%。
filter中的-s是软过滤的意思，就是把不符合后面条件的variations打上标签，但不过滤掉；-g，-G这对参数是说，indel附近的indel或snp是不准确的，大多是假阳性，过滤掉，这里我设的10bp，这个值还是比较合理的；-i是保留后面符合条件的变异，刚好和-e相反，两者选其一，我这里用的-i编写过滤表达式
QUAL：基于Phred格式的表示ALT的质量，也可以理解为可靠性；可以理解为所call出来的变异位点的质量值。Q=-10lgP，Q表示质量值；P表示这个位点发生错误的概率。因此，如果想把错误率从控制在90%以上，P的阈值就是1/10，那lg（1/10）=-1，Q=（-10）*（-1）=10。同理，当Q=20时，错误率就控制在了0.01。这个参数其实和mpileup的-Q是重复的。
DP是碱基的覆盖深度，一般很多公司和课题组会选择10，但我这里选择的是6，本着宽进严出的原则，保留更多的阳性variations，10也是没有关系的。
MQ
参考上一篇里的MAPQ分布可以看到，使用bwa mem后，其实大多数的reads的MAPQ都在60，这样的话，MQ最好也就是60，但是在40有一个小突起，对于MQ的筛选，大家的选择可以酌情而定，50是一个大家使用较多的阈值点。
MAPQ分布图
MAPQ累积分布曲线
DP4分别是正反链上REF和ALT的深度，我用“DP4[2]+DP4[3]>4”筛选ALT的深度至少是5的variations。
那么，以上的脚本其实是符合单个样本进行分析的，但如果是家系样本进行分析，其实是有问题的，因为在call这一步的时候用-v参数只输出变异位点，在获得了三个人的vcf文件进行merge的时候，对于一些变异（这些变异只存在于三个人的某一个或某两个），你就不知道不存在的那个人身上，是因为无变异还是没有覆盖到reads。这不利于做trio分析。那么要克服这个问题只需要去掉-v参数即可。
然后进行merge，是在完成了一个家系的call snp/indel以后。
$ bgzip -c -f -@ 12 var.flt.vcf > var.vcf.gz
# 进行文件压缩;-c不改变内容; -f强制输出，存在就覆盖;-@线程数
$ bcftools index -t var.vcf.gz
# 建立索引，merge需要; -t建立tbi格式索引
$ bcftools merge -Ov --force-samples -l file.list -o merge.var.vcf
# 进行合成，-O输出文件格式,v表示vcf格式文件;
# --force-samples，对于重名样本强制合成;-o输出文件
# -l包含文件名的文件，一行一个文件名，将先证者放在第一位
之后可以用bcftools对结果做一下统计处理
## 统计结果plot-vcfstats
$ bcftools stats -F/your/path/of/reference -s - merge.var.vcf > merge.var.vcf.stats && \
$ plot-vcfstats merge.var.vcf.stats -p vars_output
# plot-vcfstats程序在bcftools下的misc目录中
DeepVariant
谷歌提供了分别适用于WGS和WES的脚本，可供大家参考。我亲测了一下这个软件，速度有点慢……完全没有谷歌自己描述的那么快，做为尝鲜吧，把当时的脚本放上来，这里要特别感谢师姐的指导！虽然最后我也没打算用这个软件完成我的课题吧。
https://github.com/google/deepvariant/tree/r0.7/scripts
Deepvariant无需安装，直接拉docker下来就行，不然要是安装，这个环境配置怕是要折磨死人的。
# docker安装，仅针对ubuntu用户
$ sudo apt-get install apt-transport-https ca-certificates curl software-properties-common
$ curl -fsSL https://download.docker.com/linux/ubuntu/gpg | sudo apt-key add -
$ sudo apt-key fingerprint 0EBFCD88
$ sudo add-apt-repository "deb [arch=amd64] https://download.docker.com/linux/ubuntu $(lsb_release -cs) stable"
$ sudo apt-get update
$ sudo apt-get install docker-ce
$ sudo docker run hello-world
$ sudo usermod -aG docker $USER
拉取一下docker就可以使用deepvariant了
$ sudo docker pull dajunluo/deepvariant
$ sudo docker images
$ sudo docker run -it dajunluo/deepvariant:latest
$ sudo docker run -it --name deepvariant -v /you/path/of/refrence/dir/:/home/ref_hg19 -v /home/biowork/:/home/biowork dajunluo/deepvariant
# -v是把我本地的数据对接到docker环境里，这个请自行根据需要对接
$ nohup python /home/bin/make_examples.zip --mode calling --ref /you/path/of/refrence --reads //you/path/of/bam --examples example.gz > 1.log &
$ nohup python /home/bin/call_variants.zip --outfile call_variants_output.gz --examples example.gz --checkpoint /home/models/model.ckpt > 2.log &
$ nohup python /home/bin/postprocess_variants.zip --ref /you/path/of/refrence --infile call_variants_output.gz --outfile output.vcf.gz > 3.log &
之后便可以用bcftools或者gatk进行merge，对于这一步可以参考上面的gatk链接或者bcftools步骤。
到了这一步就可以完了么？显然不是，你用上述两个方法或者参考黄老师的脚本用gatk得到的结果，如果你仔细看，就会发现。
vcf示例
这是什么鬼，这又是什么鬼，为什么有那么多的多等位位点，你要是去统计，会发现这样的多等位位点还挺多的，所以你还不能直接过滤。对于这些变异我们一般是不会考虑嵌合的，因为平均30X的WGS是没法检测出嵌合体的。那么对于这样的位点怎么去考虑分析呢？
比较便捷的方法就是用bcftools里面的norm工具了。
$ bgzip -c -f -@ 12 merge.var.vcf > merge.var.vcf.gz
# 进行文件压缩;-c不改变内容; -f强制输出，存在就覆盖;-@线程数
$ bcftools index -t merge.var.vcf.gz
# 建立索引, -t建立tbi格式索引
$ bcftools norm -Ov -m-any -f /you/path/of/refrence merge.var.vcf.gz > norm.vcf
这样，多等位位点就变成了二等位位点了，便于后续的分析。今天的内容就到这里了，下一篇就是注释以及各种过滤了。
水平有限，要是存在什么错误请评论指出！请大家多多批评指正，相互交流，共同成长，谢谢！！！
