【www.shanpow.com--热门范文】
gatk篇1:GATK使用(一)
GATK使用(一)
(2013-01-18 09:05:42)转载▼
标签：
gatk
分类：
Bioinformatics
首先说说GATK可以做什么。它主要用于从sequencing 数据中进行variant
calling，包括SNP、INDEL。比如现在风行的exome
sequencing找variant，一般通过BWA+GATK的pipeline进行数据分析。
要run
GATK，首先得了解它的网站(http://www.broadinstitute.org/gatk/)。
一.准备工作
(1)程序下载
http://www.broadinstitute.org/gatk/download
务必下载最新版，注意不是那个GATK-lite，而是GATK2.（目前是2，不知何时会更新到3.更新很快啊有木有，所以一定要用新版）解压到某个目录即可，无需安装。打开解压缩后的目录，会发现GenomeAnalysisTK.jar这个文件。我们要用的各种工具都在这个包里。
(2)Reference
sequence和annotation下载
务必使用GATK
resource bundle上各种序列和注释来run GATK。
GATK resource
bundle介绍：
http://gatkforums.broadinstitute.org/discussion/1213/whats-in-the-resource-bundle-and-how-can-i-get-it
GATK resource bundle
FTP地址：
http://gatkforums.broadinstitute.org/discussion/1215/how-can-i-access-the-gsa-public-ftp-server
例如呢，以hg19(human genome build
19.此外，b37也很常用)为参考序列,
输入如下命令：
lftp ftp.broadinstitute.org -u
gsapubftp-anonymous
回车后键入空格，便可进入resource
bundle。进入其中名为bundle的目录，找到最新版的hg19目录。
要下载的文件包括：
ucsc.hg19.fasta
ucsc.hg19.fasta.fai
1000G_omni2.5.hg19.vcf
1000G_omni2.5.hg19.vcf.idx
1000G_phase1.indels.hg19.vcf
1000G_phase1.indels.hg19.vcf.idx
dbsnp_137.hg19.vcf
dbsnp_137.hg19.vcf.idx
hapmap_3.3.hg19.vcf
hapmap_3.3.hg19.vcf.idx
Mills_and_1000G_gold_standard.indels.hg19.vcf
Mills_and_1000G_gold_standard.indels.hg19.vcf.idx
(3)R
设置
GATK某些步骤需要使用R画一些图，此时会调用一些R
packages。如果你使用的R没有装这些packages，就无法生成pdf文件。因此得安装一些必要的R
packages，包括（但不限于）如下：
bitops
caTools
colorspace
gdata
ggplot2
gplots
gtools
RColorBrewer
reshape
gsalib
如何check这些packages有没有安装？以ggplot2为例，进入R，输入library(ggplot2)，如果没有error信息，就表明安装了ggplot2；如果没有安装，输入命令install.packages("ggplot2")，按照提示操作便可安装ggplot2。如果package安装不成功，很可能是权限问题：默认使用的R是系统管理员安装的，而你在/usr等目录下没有write权限。这是可以选择直接在自己目录下安装了R，然后install各种packages。也可以在 install.packages时指定将packages安装到自己的目录。
将R
packages安装到自己目录的方法：
install.packages("ggplot2",
lib="/your/own/R-packages/")
然后将指定的路径添加到R_LIBS变量即可。
gsalib的安装有点不一样。
从该网站下载gsalib：https://github.com/broadgsa/gatk
解压缩后，进入build.xml文件所在的目录，输入命令ant
gsalib。这样就可以安装gsalib
了。如果该gsalib仍然无法load，可能需要将gsllib的path告诉R。
http://gatkforums.broadinstitute.org/discussion/1244/what-is-a-gatkreport
这里有较详细的步骤。（奇怪的是我没有~/.Rprofle文件，ant完后也没有add
path，仍然work了）
此外，gaslib安装要求java的版本为1.6。
设置路径：
export
R_HOME=/your/R/install/path
export
R_LIBS=/your/R-packages/path
如何查看当前使用的R的路径？输入 which R即可。
二.流程概览
图片摘自http://www.broadinstitute.org/gatk/guide/topic?name=best-practices。话说这个Best
practices，刚开始接触GATK时每个人都推荐我看，可是我就是看不懂...许多东西要亲自做过后，再回头看才能懂。不过既然那么元老级人物都说要看看，所以如果你是新手，还是看看吧。
好消息是现在该页面多了个BroadE
Workshop的东东，里面的一些材料非常有用。http://www.broadinstitute.org/gatk/guide/events?id=2038#materials
在工作之前，最好把这些work
materials通览一遍。
seqanswers上的这个wiki也是很好的入门材料，可以让你看看每一步在做什么。
http://seqanswers.com/wiki/How-to/exome_analysis
但是呢，它比较过时了，它使用的GATK版本应该是2.0以前的，现在的流程有了一些调整，一些参数设置也有变化。所以该pipeline仅供参考，切勿照抄代码。例如它关于reference
sequence的那一部分，是自己手动generate hg19的genome
sequence。但是，由于后面要用到dbSNP和indel的各种annotation，而这些file的index很可能和自己generate的hg19的index不一样(resource
bundle上hg19除了chr1-22，X,Y，chrC，ChrM外还会多出一些没有组装的序列)，这样会导致程序出错，到时候还得重头来过。所以，务必使用GATK
resource bundle上的东西！
我后面贴出的代码也一样，切勿照抄。鬼知道过几个月GATK又会有什么变化。It keeps
changing all the
time...时刻跟上GATK网站上的update才是王道。（虽然它的网站做得...实在算不上好...好多时候我都找不到自己想要的东西）
三.流程
(1)Mapping。
因我只处理过Miseq和Hiseq的pair-end数据，以此为例。Mapping多使用BWA，因为它能允许indels。另一常用mapping软件Bowtie则不能generate
indels，多用于perfect match的场合如sRNA mapping。
首先建立reference的index。
bwa index -a
bwtsw -p hg19 hg19.fasta
这个过程耗时良久...above
5h吧。
-p是给reference起的名字。完成后会生成如下几个文件：
hg19.amb
hg19.ann
hg19.bwt
hg19.dict
hg19.pac
hg19.sa
然后是mapping喽，开始贴代码...先对pair-end数据的pair1和pair2分别进行mapping，生成相应的sai
file，然后再结合两者生成sam file。接着利用picard对sam file进行sort，并且生成sort后的bam
file。
sample=$1
data_dir=***
work_dir=***
ref_dir=~/database/hg19
sai_dir=$work_dir/saifiles
sam_dir=$work_dir/samfiles
bam_dir=$work_dir/bamfiles
log_dir=$work_dir/logs
met_dir=$work_dir/metrics
other_dir=$work_dir/otherfiles
gatk=/software/sequencing/GenomeAnalysisTK-2.3-4-g57ea19f/GenomeAnalysisTK.jar
picard=/software/sequencing/picard_latest
bwa aln $ref_dir/hg19 -t 4 $data_dir/${sample}_1.fastq.gz >
$sai_dir/${sample}_1.sai
bwa aln $ref_dir/hg19 -t 4 $data_dir/${sample}_2.fastq.gz >
$sai_dir/${sample}_2.sai
bwa sampe -r
"@RG\tID:$sample\tLB:$sample\tSM:$sample\tPL:ILLUMINA"
$ref_dir/hg19 $sai_dir/${sample}_1.sai
$sai_dir/${sample}_2.sai $data_dir/${sample}_1.fastq.gz
$data_dir/${sample}_2.fastq.gz | gzip >
$sam_dir/$sample.sam.gz
java -Xmx10g -Djava.io.tmpdir=/tmp \
-jar $picard/SortSam.jar \
INPUT=$sam_dir/$sample.sam.gz \
OUTPUT=$bam_dir/$sample.bam \
SORT_ORDER=coordinate\
VALIDATION_STRINGENCY=LENIENT \
CREATE_INDEX=true
Tips：
a.
bwa以及后面使用的一系列tools都能直接读取gzip压缩文件。为了节省硬盘空间，以及减少压缩/解压缩过程的read/write时间，最好直接用*.gz文件作为input。output也最好先在管道里使用gzip压缩一下再输出。
b.bwa
sample一步里的-r选项不可忽略。
c.Picard的SortSam需指定一个tmp目录，用于存放中间文件，中间文件会很大，above
10G.注意指定目录的空间大小。
d.对于压缩后大小约3G的fastq.gz文件，bwa aln需run
30-60min。将生成的两个sai file merge到一起生成sam file需60-90min。利用picard sort
sam需30-60min。
e.Samtools也提供了工具进行sam file的sort和bam
file的生成，并且不会生成大的中间文件。
(2) Duplicate
marking。
这一步的目的是mark那些PCR的duplicate。由于PCR有biases，有的序列在会被overpresented,
它们有着相同的序列和一致的map位置，对GATK的work会造成影响。这一步给这些序列设置一个flag以标志它们，方便GATK的识别。还可以设置 REMOVE_DUPLICATES=true
来丢弃duplicated序列，不过还未探索过只标记不丢弃和丢弃对于后续分析的影响。官方流程里只用标记就好。
java
-Xmx15g
-Djava.io.tmpdir=/tmp \
-jar $picard/MarkDuplicates.jar \
INPUT= $bam_dir/$sample.bam \
OUTPUT= $bam_dir/$sample.dedupped.bam \
METRICS_FILE= $met_dir/$sample.dedupped.metrics
\
CREATE_INDEX=true \
VALIDATION_STRINGENCY=LENIENT
(3)Local
relignment
INDEL附近的alignment通常不准，需要利用已知的indel信息进行realignment，分两步，第一步进行 RealignerTargetCreator，输出一个包含着possible
indels的文件。第一步IndelRealigner，利用此文件进行realign。
第一步需要用到一些已知的indel信息。那么该使用哪些文件呢？
参考如下的链接：http://gatkforums.broadinstitute.org/discussion/1247/what-should-i-use-as-known-variantssites-for-running-tool-x
其中这张图表尤其有用，可以看到RealignerTargetCreator和IndelRealigner步骤中推荐使用两个indel注释file，如代码中-known所示。
Tool
dbSNP
129 -
- dbSNP
>132 -
- Mills
indels -
- 1KG
indels -
- HapMap -
- Omni
RealignerTargetCreator
X
X
IndelRealigner
X
X
BaseRecalibrator
X
X
X
(UnifiedGenotyper/
HaplotypeCaller)
X
VariantRecalibrator
X
X
X
X
VariantEval
Xjava
-Xmx15g -jar $gatk \
-T RealignerTargetCreator \
-R $ref_dir/hg19.fasta \
-o $other_dir/$sample.realigner.intervals \
-I $bam_dir/$sample.dedupped.bam \
-known $ref_dir/1000G_phase1.indels.hg19.vcf \
-known $ref_dir/Mills_and_1000G_gold_standard.indels.hg19.vcf
java -Xmx15g -jar $gatk \
-T IndelRealigner \
-R $ref_dir/hg19.fasta \
-I $bam_dir/$sample.dedupped.bam \
-targetIntervals $other_dir/$sample.realigner.intervals \
-o $bam_dir/$sample.dedupped.realigned.bam
\
-known $ref_dir/1000G_phase1.indels.hg19.vcf \
-known
$ref_dir/Mills_and_1000G_gold_standard.indels.hg19.vcf
在一些老版本的pipeline中(如seqanswers上的那个)，这一步做完后还要对pair-end数据的mate
information进行fix。不过新版本已经不需要这一步了，Indel realign可以自己fix
mate。参考http://gatkforums.broadinstitute.org/discussion/1562/need-to-run-a-step-with-fixmateinformation-after-realignment-step
(4) Base
quality score recalibration
这一步简称BQSR，对base的quality
score进行校正。同样要用到一些SNP和indel的已知文件(参考上表和代码中-known选项)。第一步BaseRecalibrator生成一个用于校正的recal.grp文件和一个校正前的plot；第一步BaseRecalibrator生成一个校正后的plot；第三步PringReads利用recal.grp进行校正，输出校正后的bam
file。
java
-Xmx15g -jar $gatk \
-T BaseRecalibrator \
-R $ref_dir/hg19.fasta \
-I $bam_dir/$sample.dedupped.realigned.bam
\
-knownSites $ref_dir/dbsnp_137.hg19.vcf \
-knownSites
$ref_dir/Mills_and_1000G_gold_standard.indels.hg19.vcf \
-knownSites $ref_dir/1000G_phase1.indels.hg19.vcf
\
-o $other_dir/$sample.recal.grp \
-plots $other_dir/$sample.recal.grp.pdf
java -Xmx15g -jar $gatk \
-T BaseRecalibrator \
-R $ref_dir/hg19.fasta \
-I $bam_dir/$sample.dedupped.realigned.bam
\
-BQSR $other_dir/$sample.recal.grp \
-o $other_dir/$sample.post_recal.grp \
-plots $other_dir/$sample.post_recal.grp.pdf
\
-knownSites $ref_dir/dbsnp_137.hg19.vcf \
-knownSites
$ref_dir/Mills_and_1000G_gold_standard.indels.hg19.vcf \
-knownSites
$ref_dir/1000G_phase1.indels.hg19.vcf
java -Xmx15g -jar $gatk \
-T PrintReads \
-R $ref_dir/hg19.fasta \
-I $bam_dir/$sample.dedupped.realigned.bam
\
-BQSR $other_dir/$sample.recal.grp \
-o
$bam_dir/$sample.dedupped.realigned.recal.bam
TIPS：
这一步里如果产生不了PDF文件，很可能是Rscript的执行出了问题。检查R里面是否安装了该安装的packages。如果仍无法找出原因，可以在BaseRecalibrator那两步加上-l
DEBUG再运行，此时的log信息中会详细写明为什么Rscript执行不了。
gatk篇2:GATK best practice每个步骤耗时如何？
上次我们介绍了完整的 GATK best practice(请点击) 在我的基因组重测续数据分析流程，详细讲解了每个步骤的代码，输入输出文件，准备文件，以及耗时。 但是对同一个样本的多个lane的数据合并的问题，缺失了一个重要步骤，而且有热心的读者咨询整个流程的耗时问题，所以特出此番外篇作为补充。
再次介绍一下我的这次基因组重测续数据共5条lane，都是单独的PE150的fastq文件。
如果仅对L1这条lane数据进行处理
首先是BWA的比对耗时如下：
[main] Real time: 15870.794 sec; CPU: 77463.156 sec
picard.sam.SortSam done. Elapsed time: 45.60 minutes.
picard.sam.markduplicates.MarkDuplicates done. Elapsed time: 64.20 minutes.
picard.sam.FixMateInformation done. Elapsed time: 58.05 minutes.
然后是GATK对bam文件的一些预处理，步骤是：
RealignerTargetCreator --> IndelRealigner --> BaseRecalibrator --> PrintReads
后面我会讲到这些步骤是否是必须的。
耗时如下：
INFO 15:50:24,097 ProgressMeter - Total runtime 1165.34 secs, 19.42 min, 0.32 hours
INFO 17:21:00,917 ProgressMeter - Total runtime 4265.44 secs, 71.09 min, 1.18 hours
INFO 19:58:23,969 ProgressMeter - Total runtime 9436.69 secs, 157.28 min, 2.62 hours
INFO 23:41:00,540 ProgressMeter - Total runtime 13349.77 secs, 222.50 min, 3.71 hours
可以看到对L1这条lane数据来说，整个流程耗时才不到10个小时，还算是可接受范围内的。
接下来用HaplotypeCaller找变异
这个步骤我对realign的bam和recal的bam分别处理了，耗时如下：
INFO 20:40:49,063 ProgressMeter - Total runtime 39243.88 secs, 654.06 min, 10.90 hours
INFO 08:53:17,633 ProgressMeter - Total runtime 43939.69 secs, 732.33 min, 12.21 hours
对bam文件找变异的时间取决于reads数量的多少以及这些reads覆盖的基因组区域大小，虽然一条lane的数据量很小，但它仍然是全基因组测序，如果是全外显子测序这个耗时是不一样的。
整个L1这条lane数据(120M的reads)处理后的文件大小如下所示：
3.0M Jun 7 01:14 L1.bam.bai
8.8G Jun 7 02:33 L1_marked.bam
9.0G Jun 7 03:31 L1_marked_fixed.bam
3.3M Jun 7 03:36 L1_marked_fixed.bam.bai
2.6K Jun 7 02:33 L1.metrics
8.2M Jun 5 17:21 L1_realigned.bai
9.0G Jun 5 17:21 L1_realigned.bam
8.2M Jun 5 23:41 L1_recal.bai
27G Jun 5 23:41 L1_recal.bam
8.1M Jun 5 23:53 L1_recal.bam.bai
39M Jun 5 15:50 L1_target.intervals
320K Jun 5 19:58 L1_temp.table
因为我的这次基因组重测续数据共5条lane，这5条lane的数据其实是可以并行完成这几个步骤的，最长耗时约12小时。 每个数据处理我都分配了 5个线程， 40G的内存。
merge后不进行AddOrReplaceReadGroups处理
意味着要把5条lane的数据当做是不同的样本，这样后续处理这5个lane的bam矫正以及找变异都是独立进行的，虽然最后生成的vcf文件只有一个，但是每条lane都有独立的基因型。 realign和recal耗时如下：
INFO 17:54:59,396 ProgressMeter - Total runtime 5194.10 secs, 86.57 min, 1.44 hours
INFO 02:04:10,907 ProgressMeter - Total runtime 22558.84 secs, 375.98 min, 6.27 hours
INFO 18:48:45,762 ProgressMeter - Total runtime 60267.18 secs, 1004.45 min, 16.74 hours
INFO 21:30:54,519 ProgressMeter - Total runtime 96119.87 secs, 1602.00 min, 26.70 hours
同样还是对realign的bam和recal的bam分别用HaplotypeCaller找变异
INFO 19:36:31,960 ProgressMeter - Total runtime 201031.47 secs, 3350.52 min, 55.84 hours
INFO 22:59:21,693 ProgressMeter - Total runtime 184944.78 secs, 3082.41 min, 51.37 hours
可以看到这个时候的耗时相比仅针对一条lane已经是非常恐怖了，近100个小时。不仅仅是因为这个时候reads数量增多，而且是因为5条lane独立样本处理。
merge后需要AddOrReplaceReadGroups处理
这个才是正确的做法，因为不同的lane出来的数据都是我本人的全基因组重测续数据，后续处理应该是当做一个样本的，所有需要AddOrReplaceReadGroups处理，代码是：
### AddOrReplaceReadGroups ###
java -Djava.io.tmpdir=$TMPDIR -Xmx40g -jar $PICARD AddOrReplaceReadGroups \
INPUT=${sample}.bam OUTPUT=${sample}_tmp.bam RGID=jmzeng RGLB=lib_all RGPL=illumina RGPU=x10 RGSM=jmzeng
mv ${sample}_tmp.bam ${sample}.bam
samtools index ${sample}.bam
这个代码需要结合前面的GATK best practice一起理解。
realign和recal耗时如下：
INFO 15:52:21,739 ProgressMeter - Total runtime 3573.62 secs, 59.56 min, 0.99 hours
INFO 22:46:39,615 ProgressMeter - Total runtime 24853.46 secs, 414.22 min, 6.90 hours
INFO 16:06:14,340 ProgressMeter - Total runtime 62366.41 secs, 1039.44 min, 17.32 hours
INFO 18:18:08,004 ProgressMeter - Total runtime 94304.46 secs, 1571.74 min, 26.20 hours
这个耗时跟上面把lane当做是独立样本的没有什么区别，因为耗时取决于reads数据量。
接下来找变异，我不仅仅是对realign和recal，还加入了最原始的bam，全部耗时如下：
INFO 02:29:32,680 ProgressMeter - Total runtime 149229.64 secs, 2487.16 min, 41.45 hours
INFO 02:15:39,234 ProgressMeter - Total runtime 148379.91 secs, 2473.00 min, 41.22 hours
INFO 21:20:51,032 ProgressMeter - Total runtime 130663.06 secs, 2177.72 min, 36.30 hours
可以看到这些耗时都显著的小于把lane当做独立样本。总流程耗时仍然超过了80个小时，这也就是为什么时隔10天我才推出番外~~~
全部流程输出的文件大小如下:
59G Jun 14 14:17 jmzeng.bam
8.4M Jun 14 14:52 jmzeng.bam.bai
1.1G Jun 21 02:29 jmzeng_merge_raw.snps.indels.vcf
12M Jun 21 02:29 jmzeng_merge_raw.snps.indels.vcf.idx
8.4M Jun 14 22:46 jmzeng_realigned.bai
59G Jun 14 22:46 jmzeng_realigned.bam
1.1G Jun 21 02:15 jmzeng_realigned_raw.snps.indels.vcf
12M Jun 21 02:15 jmzeng_realigned_raw.snps.indels.vcf.idx
8.5M Jun 16 18:18 jmzeng_recal.bai
161G Jun 16 18:18 jmzeng_recal.bam
8.5M Jun 16 19:39 jmzeng_recal.bam.bai
1.1G Jun 20 21:20 jmzeng_recal_raw.snps.indels.vcf
12M Jun 20 21:20 jmzeng_recal_raw.snps.indels.vcf.idx
47M Jun 14 15:52 jmzeng_target.intervals
323K Jun 15 16:06 jmzeng_temp.table
值得注意的是，recal之后的bam文件是原始bam的3倍大小，特别耗费存储空间。
接下来我会讲解realign和recal步骤的必要性，毕竟这两个步骤也的确太耗时了，尤其是recal，不仅仅耗时还特别占硬盘存储。
猜你喜欢
基因组 | 游记 | 工作资讯
学习课程 | 好书分享
菜鸟入门
Linux | Perl | R语言 | 可视化
R包 | perl模块 | python模块
数据分析
ChIP-seq（上）| ChIP-seq（下）| RNA-seq | miRNA
WGS,WES,RNA-seq组与ChIP-seq之间的异同
编程实践
第0题 | 探索人类基因组序列 | 最后一题
直播基因组分析
我的基因组 | 解惑帖
一个标准的基因检测报告目录生信技能树
编辑：思考问题的熊
gatk篇3:GATK best practice每个步骤都是必须的吗？
GATK4这个新的版本已经发布了，我前面写的一系列教程都是基于GATK3.X的，当然，整体数据处理流程其实并没有变化，但是时间消耗，步骤选择全部需要重新评价了。
一个全基因组重测序分析实战
GATK best practice每个步骤耗时如何？
最后这个步骤的必要性必须得马上发出来了，就是因为有些步骤也的确太耗时了，尤其是recal，不仅仅耗时还特别占硬盘存储！ 那么我们就从最后找到的变异文件的比较情况来说明这些步骤的必要性吧！ 我们再回顾一下GATK best practice流程吧！
流程介绍
bwa(MEM alignment)
picard(SortSam)
picard(MarkDuplicates)
picard(FixMateInfo)
GATK(RealignerTargetCreator)
GATK(IndelRealigner)
GATK(BaseRecalibrator)
GATK(PrintReads)
GATK(HaplotypeCaller)
GATK(GenotypeGVCFs)
我上一讲说到过，我对realign和recal以及原始的bam都用了HaplotypeCaller找变异，得到的vcf文件如下：
1.1G Jun 21 02:29 jmzeng_merge_raw.snps.indels.vcf
1.1G Jun 21 02:15 jmzeng_realigned_raw.snps.indels.vcf
1.1G Jun 20 21:20 jmzeng_recal_raw.snps.indels.vcf
这是最原始的变异文件，我们需要进行过滤，拆分SNP和INDEL，这样才能更好的对它们两两比较。(这样就是比较查看realign和recal步骤对最后找变异的影响)
拆分SNP和INDEL并过滤低质量位点
代码如下：
module load java/1.8.0_91
GENOME=/home/jianmingzeng/reference/genome/human_g1k_v37/human_g1k_v37.fasta
GATK=/home/jianmingzeng/biosoft/GATK/GenomeAnalysisTK.jar
TMPDIR=/home/jianmingzeng/tmp/software
sample=$1
## for SNP
: "
"
java -Djava.io.tmpdir=$TMPDIR -Xmx40g -jar $GATK -T SelectVariants -R $GENOME \
-selectType SNP \
-V ${sample}_raw.snps.indels.vcf -o ${sample}_raw_snps.vcf
java -Djava.io.tmpdir=$TMPDIR -Xmx40g -jar $GATK -T VariantFiltration -R $GENOME \
--filterExpression "QD60.0 || MQ \
--filterName "my_snp_filter" \
-V ${sample}_raw_snps.vcf -o ${sample}_filtered_snps.vcf
java -Djava.io.tmpdir=$TMPDIR -Xmx40g -jar $GATK -T SelectVariants -R $GENOME \
--excludeFiltered \
-V ${sample}_filtered_snps.vcf -o ${sample}_filtered_PASS.snps.vcf
java -Djava.io.tmpdir=$TMPDIR -Xmx40g -jar $GATK -T VariantEval -R $GENOME \
-eval ${sample}_filtered_PASS.snps.vcf -o ${sample}_filtered_PASS.snps.vcf.eval
## for INDEL
java -Djava.io.tmpdir=$TMPDIR -Xmx40g -jar $GATK -T SelectVariants -R $GENOME \
-selectType INDEL \
-V ${sample}_raw.snps.indels.vcf -o ${sample}_raw_indels.vcf
java -Djava.io.tmpdir=$TMPDIR -Xmx40g -jar $GATK -T VariantFiltration -R $GENOME \
--filterExpression "QD200.0 || ReadPosRankSum \
--filterName "my_indel_filter" \
-V ${sample}_raw_indels.vcf -o ${sample}_filtered_indels.vcf
java -Djava.io.tmpdir=$TMPDIR -Xmx40g -jar $GATK -T SelectVariants -R $GENOME \
--excludeFiltered \
-V ${sample}_filtered_indels.vcf -o ${sample}_filtered_PASS.indels.vcf
java -Djava.io.tmpdir=$TMPDIR -Xmx40g -jar $GATK -T VariantEval -R $GENOME \
-eval ${sample}_filtered_PASS.indels.vcf -o ${sample}_filtered_PASS.indels.vcf.eval
要深度理解这个代码请点击自行阅读文档。 其实我本人不大喜欢用这个工具的各种参数，我比较喜欢自行写脚本来过滤vcf文件。
这样得到的文件如下：
801242 jmzeng_merge_filtered_indels.vcf
760890 jmzeng_merge_filtered_PASS.indels.vcf
3812094 jmzeng_merge_filtered_PASS.snps.vcf
4034168 jmzeng_merge_filtered_snps.vcf
801240 jmzeng_merge_raw_indels.vcf
4837892 jmzeng_merge_raw.snps.indels.vcf
4034166 jmzeng_merge_raw_snps.vcf
801963 jmzeng_realigned_filtered_indels.vcf
761510 jmzeng_realigned_filtered_PASS.indels.vcf
3812442 jmzeng_realigned_filtered_PASS.snps.vcf
4034797 jmzeng_realigned_filtered_snps.vcf
801961 jmzeng_realigned_raw_indels.vcf
4839256 jmzeng_realigned_raw.snps.indels.vcf
4034795 jmzeng_realigned_raw_snps.vcf
793611 jmzeng_recal_filtered_indels.vcf
754755 jmzeng_recal_filtered_PASS.indels.vcf
3784343 jmzeng_recal_filtered_PASS.snps.vcf
4010670 jmzeng_recal_filtered_snps.vcf
793609 jmzeng_recal_raw_indels.vcf
4806696 jmzeng_recal_raw.snps.indels.vcf
4010668 jmzeng_recal_raw_snps.vcf
很容易理解：
对recal的bam得到的变异vcf文件来说，总共是480万位点，拆分成401万的SNP和79万的INDEL，然后经过过滤后剩下378万的SNP和75万的INDEL。
对原始的bam得到的变异vcf文件来说， 总共是483万位点，拆分成403万的SNP和80万的INDEL，然后经过过滤后剩下381万的SNP和76万的INDEL。
对重排的bam得到的变异vcf文件来说， 总共是483万位点，拆分成403万的SNP和80万的INDEL，然后经过过滤后剩下381万的SNP和76万的INDEL。
从位点数量级来说，原始的bam和重排的bam得到的变异vcf文件没区别，所以需要用专业的工具来具体比较它们里面的每一个位点。 我这里还是选择SnpEff套件里面的SnpSift工具。
首先比较SNP位点
代码如下：
java -jar ~/biosoft/SnpEff/snpEff/SnpSift.jar concordance -v ../jmzeng_merge_filtered_PASS.snps.vcf ../jmzeng_realigned_filtered_PASS.snps.vcf 1>concordance.txt 2>SnpSift_Concordance.log
java -jar ~/biosoft/SnpEff/snpEff/SnpSift.jar concordance -v ../jmzeng_recal_filtered_PASS.snps.vcf ../jmzeng_realigned_filtered_PASS.snps.vcf 1>concordance.txt 2>SnpSift_Concordance.log
java -jar ~/biosoft/SnpEff/snpEff/SnpSift.jar concordance -v ../jmzeng_recal_filtered_PASS.snps.vcf ../jmzeng_merge_filtered_PASS.snps.vcf 1>concordance.txt 2>SnpSift_Concordance.log
比较的结果如下：
Number of samples:
1 File ../jmzeng_merge_filtered_PASS.snps.vcf
1 File ../jmzeng_realigned_filtered_PASS.snps.vcf
1 Both files
# Errors:
ALT field does not match 31
Number of samples:
1 File ../jmzeng_recal_filtered_PASS.snps.vcf
1 File ../jmzeng_realigned_filtered_PASS.snps.vcf
1 Both files
# Errors:
ALT field does not match 204
Number of samples:
1 File ../jmzeng_recal_filtered_PASS.snps.vcf
1 File ../jmzeng_merge_filtered_PASS.snps.vcf
1 Both files
# Errors:
ALT field does not match 208
可以看到对高质量的SNP位点来说，3种bam文件得到SNP信息差别很微弱，在可接受范围内。但是我们不能忽视原始的bam和重排的bam得到的变异vcf文件要比recal后的少了近两万这个事实！！！
接着比较INDEL位点
代码如下：
java -jar ~/biosoft/SnpEff/snpEff/SnpSift.jar concordance -v ../jmzeng_merge_filtered_PASS.indels.vcf ../jmzeng_realigned_filtered_PASS.indels.vcf 1>concordance.txt 2>SnpSift_Concordance.log
java -jar ~/biosoft/SnpEff/snpEff/SnpSift.jar concordance -v ../jmzeng_recal_filtered_PASS.indels.vcf ../jmzeng_realigned_filtered_PASS.indels.vcf 1>concordance.txt 2>SnpSift_Concordance.log
java -jar ~/biosoft/SnpEff/snpEff/SnpSift.jar concordance -v ../jmzeng_recal_filtered_PASS.indels.vcf ../jmzeng_merge_filtered_PASS.indels.vcf 1>concordance.txt 2>SnpSift_Concordance.log
比较的结果如下：
Number of samples:
1 File ../jmzeng_merge_filtered_PASS.indels.vcf
1 File ../jmzeng_realigned_filtered_PASS.indels.vcf
1 Both files
# Errors:
REF fields does not match 28
ALT field does not match 45
Number of samples:
1 File ../jmzeng_recal_filtered_PASS.indels.vcf
1 File ../jmzeng_merge_filtered_PASS.indels.vcf
1 Both files
# Errors:
REF fields does not match 1295
ALT field does not match 964
Number of samples:
1 File ../jmzeng_recal_filtered_PASS.indels.vcf
1 File ../jmzeng_realigned_filtered_PASS.indels.vcf
1 Both files
# Errors:
REF fields does not match 1276
ALT field does not match 947
INDEL本身对参数非常敏感，这时候不仅仅是数量的差异，而且本身找到的位点突变情况的差异也不少。
所以我的结论是，GATK的BEST PRACTISE的每个步骤都是必须的！虽然他们非常耗时，但是对结果准确性的影响的确非常显著。 如果要想把测序数据在临床上面应用，那么每一个步骤都是非常有意义的。
比如，如果我们来分析realign之后的高质量SNP为什么会有31个的ALT改变了呢？
21 10716592
21 10701512
21 10700605
20 26317761
19 15787221
18 18515822
17 25266293
16 35230302
16 33975649
16 33972478
10 42400348
10 42393437
9 66455306
4 49111623
1 142537167
Y 58977742
Y 13478115
X 61909282
X 61730552
X 61721098
简单看了几眼， 发现都是在染色体的中心粒的地方！
虽然GATK best practice中的BSQR步骤的确非常耗时(WGS数据约40小时)，但是对最后的结果影响还是蛮大的，所以是否需要省略这个步骤取决于你自己对找变异的精度要求。毕竟这个步骤说明书上面可说的是用了机器学习的方法来矫正测序误差呀！
最后，上两张GATK4的PPT吧~
又要开始学新的软件，写新的教程啦~~~
GATK4性能提升如此显著，不得不换！
尤其是在生物信息学这个日新月异的领域。
猜你喜欢
基因组 | 游记 | 工作资讯
学习课程 | 好书分享
菜鸟入门
Linux | Perl | R语言 | 可视化
R包 | perl模块 | python模块
数据分析
ChIP-seq（上）| ChIP-seq（下）| RNA-seq | miRNA
WGS,WES,RNA-seq组与ChIP-seq之间的异同
编程实践
第0题 | 探索人类基因组序列 | 最后一题
直播基因组分析
我的基因组 | 解惑帖
一个标准的基因检测报告目录生信技能树
编辑：思考问题的熊
