Variant 分析阶段小结2- 变异寻找碎碎念

生信技能树

发布于 2018-07-27 12:57:44

3.8K0

发布于 2018-07-27 12:57:44

写在前面：『思考问题的熊』专栏上次更新还要追溯到4月19号的 Variant 分析阶段小结1-基础碎碎念，过去接近一个月的时间里我分别经历了两次长途出差和电脑无法连接特定网络的持续尴尬。特定网络正是所有以 qq.com 结尾的网站，当然包括微信公众平台，所以文章都编辑不了。身体和心理的双重袭击让我只能选择围笑:)

前几天经过网络驱动重装升级又降级，DNS修改修改又修改，交换机重启重启再重启等操作，终于能继续展示真正的技术了。今天是 variant 分析的第二部分小节，三步寻找突变。写这些文章的时候我还在用GATK3的流程，后面整理好新的内容再做补充。

题图：葫芦娃之二娃

8800字，约19分钟，思考问题的熊专栏10

整体流程

1.Pre-processing: raw DNA(RNA)seq sequence reads (FASTQ files) to analysis-ready reads (BAM files)

原始数据预处理比对参考基因组原始 bam 文件进行处理

2.Variant discovery: reads (BAM files) to variants (VCF files)

利用GATK, bcftools 或 freebayes 等软件找到 raw Variants利用各种策略对 raw Variants 进行筛选

3.Refinement and evaluation 根据需要进行后续分析和统计

统计 variant 的各种分布情况和基因型信息对数据进行需要的注释下游个性化分析

mapping

对于DNA-seq数据，常用软件有 bwa 或 bowtie；对于RNAseq数据可以使用STAR(后续文章再具体说 RNA-seq 找 SNP 的问题)

首先对参考基因组进行index，然后可以使用bwa mem 进行比对。这里需要说明的是如果在分析过程中但凡要涉及到使用 GATK 相关的流程，比对后产生的 bam 文件必须包含@RG tag 信息，如果没有的在后续分析中会各种报错。

使用 samtools and bcftools

如果使用samtools 和 bcftools 找snp，官方网站给了一些建议的步骤，大致如下：

bwa index <reference.fa>
bwa mem -t 10 -R '@RG\tID:group1\tSM:sample1\tPL:illumina\tLB:library1' \
<reference.fa> <read1.fa> <read1.fa> > aligned.sam

# 官方Workflows建议步骤，删掉一些bwa产生的不常用且可能会对下游分析有影响的flags
samtools fixmate -O bam <aligned.sam> <aligned_fixmate.bam>

# To sort them from name order into coordinate order
samtools sort -@ 5 -O bam -o <aligned_sorted.bam> -T </tmp/aligned_temp> <aligned_fixmate.bam>

# remove duplicate
samtools rmdup <aligned_sorted.bam> <aligned_rmdup.bam>
# 会直接取出掉那些重复的reads

# samtools index
samtools index <aligned_rmdup.bam>

GATK and picard

bwa mem -M -t 10 -R '@RG\tID:group1\tSM:sample1\tPL:illumina\tLB:library1' \
<reference.fa> <read1.fa> <read1.fa> > aligned.sam
#-M flag causes BWA to mark shorter split hits
# as secondary (essential for Picard compatibility)

# picard sort
java -Xmx10g -jar /home/zf/software/picard/build/libs/picard.jar SortSam \
INPUT=aligned.sam OUTPUT=aligned_sort.bam SORT_ORDER=coordinate

# picard rmeomv dupilcate
java -Xmx10g -jar /home/zf/software/picard/build/libs/picard.jar \
MarkDuplicates INPUT=aligned_sort.bam \
OUTPUT=aligned_rmdup.bam METRICS_FILE=metrics.txt
# MarkDuplicates处理完成不会删除这些reads而是加上一个标签
#带有这些标签的reads在后续处理的时候会被跳过

# picard index
java -Xmx10g -jar /home/zf/software/picard/build/libs/picard.jar \
BuildBamIndex INPUT=aligned_rmdup.bam

Variants calling

samtools http://samtools.sourceforge.net/
bcftools http://www.htslib.org/doc/bcftools.html

samtools and bcftools

samtools mpileup -d 1000 -Bugf reference.fa aligned_rmdup_samtools.bam | \
bcftools call --threads 8 -vmO z -o raw_samtools.vcf.gz

关于这几个参数的含义，B是指 Disable probabilistic realignment for the computation of base alignment quality (BAQ). BAQ is the Phred-scaled probability of a read base being misaligned.

官方说法是试用这个参数可以有效的降低由比对错误(misalignments)带来的假阳性SNP，但是实际效果如何，或者会不会矫枉过正大家可以自行测试。另外，在 bcftools 中的 M 是指 alternative modelfor multiallelic and rare-variant calling.

结果展示

# samtools snps
CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  test
Chr1    22321   .       A       G       8.99921 .       DP=1;SGB=-0.379885;MQ0F=0;AC=2;AN=2;DP4=0,0,1,0;MQ=60   GT:PL   1/1:38,3,0

# samtools indels
CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  test
Chr1    31696   .       TGGG    TGG     5.04449 .       INDEL;IDV=1;IMF=1;DP=1;SGB=-0.379885;MQ0F=0;AC=2;AN=2;DP4=0,0,1,0;MQ=60 GT:PL   1/1:33,3,0

GATK and picard

说明此文章针对的是GATK 3.0 版本。且分析物种为植物。更新内容后续会发布在博客。

GATK call variant 的工具非常之多，但是如果观察官方提供的最佳实践步骤的话多数都是使用HaplotypeCaller(HC)，这厮在前任的基础上引入了实时de novo算法，能够通过对活跃区域（变异热点区域）进行局部重组装的同时寻找 SNP 和 INDEL。HC 可以处理非二倍体物种和混池实验数据，也可以处理可变剪切。

因为它的存在，使得用Unified Genotyper 时的 indel realign, BQSR步骤能够省略（建议保留，不过植物里面就省了心有余而力不走），但是HC速度慢相对慢，某些情况下3.x 的版本跑多线程还会报bug。

更具体的内容可以参考该软件官方介绍地址，和它的工作原理。

工作流程如下

寻找突变热点区域
进行局部重组装，确定单倍型（haplotypes）
利用成对隐马尔科夫模型，计算给定数据的单倍型可能性。
确定样本基因型

关于这个软件具体怎么用，可以参考这两篇文章：

Call variants with HaplotypeCaller

Call germline SNPs and indels via local re-assembly of haplotypes

java -Xmx60g -jar GenomeAnalysisTK.jar \
-T HaplotypeCaller \
-R all.con.fa  \
-I aligned_rmdup.bam \
--genotyping_mode DISCOVERY \
-stand_call_conf 30 # for DNA-seq(20 RNA-seq)
-o raw_gatk.vcf

结果展示

#gatk snps
CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  test
Chr1    55398   .       A       G       56.74   .       AC=2;AF=1.00;AN=2;DP=2;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=28.37;SOR=2.303 GT:AD:DP:GQ:PL  1/1:0,2:2:6:84,6,0

# gatk indels
CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  test
Chr1    217501  .       G       GAT     52.70   .       AC=2;AF=1.00;AN=2;DP=2;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=26.35;SOR=2.303 GT:AD:DP:GQ:PL  1/1:0,2:2:6:89,6,0

freebayes

freebayes 的用法比GATK简单很多，比bcftools的参数还要简单一些，也是一个常用的call variant 软件，但是软件本身不支持多线程，相对较慢（不知道目前有没有更新版本）。其输入文件为和samtoosl 一样，也是去重数据且要有组信息。

freebayes -f <reference.fa> <aligned_rmdup.bam> >raw_freebayes.vcf

结果展示

Chr1    248217  .       G       T       60.1096 .       AB=0;ABP=0;AC=2;AF=1;AN=2;AO=2;CIGAR=1X;DP=2;DPB=2;DPRA=0;EPP=7.35324;EPPR=0;GTI=
0;LEN=1;MEANALT=1;MQM=60;MQMR=0;NS=1;NUMALT=1;ODDS=7.37776;PAIRED=1;PAIREDR=0;PAO=0;PQA=0;PQR=0;PRO=0;QA=76;QR=0;RO=0;RPL=1;RPP=3.0103;RP
PR=0;RPR=1;RUN=1;SAF=1;SAP=3.0103;SAR=1;SRF=0;SRP=0;SRR=0;TYPE=snp;technology.illumina=1      GT:DP:AD:RO:QR:AO:QA:GL 1/1:2:0,2:0:0:2:76:
-7.21481,-0.60206,0

和上面两个略有不同的是，在freebayes中，SAF 和 SAR 展示了alles的strand信息，RPL和RPR 展示了reads的方向信息。这些都可以做为后续的筛选依据。

get highest-quality call set

得到的原始vcf文件可能包含有非常多的假阳性结果，这个时候对原始结果进行各种各样的filter就变得非常关键。如果使用的GATK，可以采用variant recalibration（VQSR）或者 hard-filtering 进行矫正。如果是做非人物种的研究，VQSR这步基本可以放弃，无论是bcftools 出来的结果还是 GATK 出来的结果，我们都可以使用hard-filtering。

进一步问 hard-filter的标准是什么？答案似乎是没有标准。

关于数据过滤的问题，在GATK 论坛上有几篇非常好的讨论文章。

Understanding and adapting the generic hard-filtering recommendations
Using JEXL to apply hard filters or select variants based on annotation values
(howto) Apply hard filters to a call set

GATK filter

使用工具 VariantFiltration和 SelectVariants。

其中VariantFiltration的功能是Filter variant calls based on INFO and FORMAT annotations，而SelectVariants 的功能是Select a subset of variants from a larger callset。这两个工具的官方主页有不少例子可供参考。

Extract the SNPs and Indels from the call set

一般情况我们会将SNP和indel数据分开筛选和后续分析，所以首先将snp和indel数据分开。

java -Xmx8g -jar GenomeAnalysisTK.jar \
-T SelectVariants \
-R all.con.fa \
-V raw_gatk.vcf \
-selectType SNP \
-o raw_snps_gatk.vcf

java -Xmx8g -jar GenomeAnalysisTK.jar \
-T SelectVariants \
-R all.con.fa \
-V raw_gatk.vcf \
-selectType INDEL \
-o raw_indels_gatk.vcf

官方推荐标准

For SNPs:

QD<2.0
MQ<40.0
FS>60.0
SOR>3.0
MQRankSum<-12.5
ReadPosRankSum<-8.0

For indels:

QD<2.0
ReadPosRankSum<-20.0
InbreedingCoeff<-0.8
FS>200.0
SOR>10.0

Apply the filter to the SNP call set

java -Xmx8g -jar GenomeAnalysisTK.jar \
-T VariantFiltration \
-R all.con.fa \
-V raw_snps_gatk.vcf \
--filterExpression "( vc.hasAttribute('QD') && QD<2.0 ) || FS > 60.0 \
|| MQ < 40.0 || ( vc.hasAttribute('ReadPosRankSum' ) && ReadPosRankSum < -8.0 ) \
|| ( vc.hasAttribute('MQRankSum') && MQRankSum < -12.5 ) \
|| ( vc.hasAttribute('SOR') && SOR > 10.0 ) || QUAL < 30.0 " \
--filterName "my_snp_filter" \
-o filter_snps_gatk.vcf

java -Xmx8g -jar GenomeAnalysisTK.jar \
-T VariantFiltration \
-R all.con.fa \
-V raw_indels_gatk.vcf \
--filterExpression " ( vc.hasAttribute('QD') && QD<2.0 ) \|| FS > 200.0 \
|| ( vc.hasAttribute ('ReadPosRankSum') && ReadPosRankSum < -20.0 ) \
|| ( vc.hasAttribute('SOR') && SOR > 10.0 )|| QUAL < 30.0" \
--filterName "my_indel_filter" \
-o filter_indels_gatk.vcf

到这一步需要注意的是，不同于bcftools的筛选，这里的VariantFiltration只是把不符合的内容标记出来，在新生成的文件中，符合筛选标准的variant在FILTER一列会显示pass ，如果没有通过会显示上面命令中的 filterName

结果展示

# 通过
Chr1    55398   .   A   G   56.74   PASS    AC=2;AF=1.00;AN=2;DP=2;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=60.00;QD=28.37;SOR=2.303 GT:AD:DP:GQ:PL  1/1:0,2:2:6:84,6,0

# 未通过
Chr1    1343566 .       G       T       62.74   my_snp_filter   AC=2;AF=1.00;AN=2;DP=2;ExcessHet=3.0103;FS=0.000;MLEAC=2;MLEAF=1.00;MQ=38.21;QD=31.37;SOR=2.303 GT:AD:DP:GQ:PL  1/1:0,2:2:6:90,6,0

提取通过filter的结果

java -Xmx8g -jar GenomeAnalysisTK.jar \
-T SelectVariants \
-R all.con.fa \
-V filter_snps_gatk.vcf \
-ef -o filtered_snps_gatk.vcf

bcftools filter

bcftools filter的方式有两种，一种是使用 bcftools filter，另一种是直接使用 bcftools view

# 首先对结果进行index
tabix -p vcf raw_samtools.vcf.gz
# 注意tabix 使用时输入文件一定是使用bgzip压缩的vcf文件

# 分别提取snp和indel结果
bcftools view -i "type='snps'" raw_samtools.vcf.gz -O z\
 -o raw_samtools_snps.vcf.gz && tabix \
 -p vcf raw_samtools_snps.vcf.gz

bcftools view -i "type='indels'" raw_samtools.vcf.gz -O z \
-o raw_samtools_indels.vcf.gz && tabix \
-p vcf raw_samtools_indels.vcf.gz

#
bcftools filter -e "MQ < 40 || QUAL < 30" -s LOWQUAL \
raw_samtools_snps.vcf.gz |bcftools view \
-f PASS - > filtered_samtools_snps.vcf.gz

freeabyes filter

可以和bcftools 的过滤方式保持一致也可以加入一下独有的筛选标准。比如 SAF >0 & SAR>0; PRR>1 & RPL>1 ; QUAL/AO>10 等等。

关于筛选用到的QUAL和MQ 解释如下：

MQ is the mapping quality, which is the fifth column in SAM record. QUAL, meanwhile, is the base quality score, which is derived from the 11th column in SAM record. MQ is typically an indication of how unique the region's sequence is, the higher the MQ, the more unique the sequence. QUAL, is the sequencing quality, which can be platform biased, e.g. Ion seemed to have lower QUAL compared to Illumina.

Software comparison

仅仅是某一次测试数据且和参数有很大的关系

用时：freebayes 和 GATK 的时间要多于bcftools，还长不少。

raw variant 数量

原始数据中，gatk和freebayes产生的结果相对保守，samtools 的snp明显更多，如果进一步统计可以发现主要原因是后者的QUAL没有提前进行过滤造成，且在低QUAL位点找到了大量的snp，过滤后的数据整体基本保持一致。freebayes 找到了更多的MNPs。indel方面，freebayes 的长度最长可以找到40bp，gatk则可能找到很长的indel，大于60bp的有40个，而bcftools则找到了50bp的indel。

软件上手难以程度： bcftools 和 freebayes 要明显简单与gatk，一方面是因为gatk的功能实在是太多，另一方面是因为步骤比价复杂。但是上手难度大的另一个意思是更加值得深入学习。

售后支持： gatk的官方网站和社区实在是强大，你碰到的或者想到的问题在它的社区应该都可以找到答案。

综合考虑：gatk 的支持度最高，处理流程和方法最完善。bcftools 处理vcf文件方便耗时少。顺手也可以把freebayes 做一下。如果时间允许可以使用三个软件共同操作，当然，到了这里变异相关的分析不是结束而是刚刚开始……

Variant 分析阶段小结1-基础碎碎念
谁来拯救你我的屁屁踢
RNA-seq 从原理到应用
生物统计学与R极简手册
用 Excel 怎么了，你咬我啊？
高效使用云笔记只需“一五一十”

作者博客地址：http://kaopubear.top

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2018-05-17，如有侵权请联系 cloudcommunity@tencent.com 删除

其他

本文分享自生信技能树微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体分享计划，欢迎热爱写作的你一起参与！

其他

登录后参与评论

0 条评论

热度