结合GATK和SAMtools从头挖掘SNPs和INDELs

提出问题

给出以下的已有结果:
1. 物种species的基因组fasta文件: genome.fasta
2. 对该物种的一个名为sample的基因组DNA样,进行了Illumina hiseq2000 paired-end测序(300bp插入片段文库,~20x基因组碱基覆盖度),结果文件: reads1.fq, reads2.fq
问怎么进行SNPs和INDELs的calling?

文章有很多地方引自Nowind的博文:GATK使用简介 Part2/2
流程参照:http://www.broadinstitute.org/gatk/guide/topic?name=best-practices

以下给出具体的答案和步骤, 以上述已有结果的3个文件所在的目录为当前工作目录,各种软件的主目录以$software表示。多线程的程序以8个线程运行。

1. Raw reads的处理

使用NGSQC toolkit使用默认设置对raw reads进行过滤。

$ perl $NGSQCToolkitHome/QC/IlluQC_PRLL.pl \
  -pe reads1.fq reads2.fq 2 5 -c 8 \
  -o QC

2. 将过滤后的reads比对到基因组上

$ bowtie2-build genome.fasta genome
$ bowtie2 --rg-id sample --rg "PL:ILLUMINA" --rg "SM:sample" \
  -x geome -1 ./QC/reads1.fq_filtered -2 ./QC/reads2.fq.filtered \
  -p 8 -S sample.sam

3. 将sam文件转换为Bam文件并标记出PCR duplicates

$ samtools view -bS sample.sam > sample.bam
$ java -Xmx10g -jar $picardHome/SortSam.jar \
  INPUT=sample.bam OUTPUT=sample.sort.bam \
  SORT_ORDER=coordinate
$ java -Xmx10g -jar $picardHome/MarDuplicates.jar \
  INPUT=sample.sort.bam OUTPUT=sample.dd.bam \
  METRICS_FILE=sample.dd.metrics \
  MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000

4. 构建索引文件

$ samtools faidx genome.fasta
$ java -Xmx10g -jar $picardHome/CreateSequenceDictionary.jar \
  R=genome.fasta O=genome.dict
$ samtools index sample.dd.bam

5. 对INDEL周围进行重新比对(realignment)

$ java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \
  -R genome.fasta -T RealignerTargetCreator \
  -I sample.dd.bam -o sample.realn.intervals
$ java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \
  -R genome.fasta -T IndelRealigner \
  -targetIntervals sample.realn.intervals \
  -I sample.dd.bam -o sample.realn.bam

6. 第1遍Variant calling

6.1 使用GATK和samtools进行SNP和INDEL calling

$ java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \
  -R genome.fasta -T UnifiedGenotyper \
  -I sample.realn.bam -o sample.gatk.raw1.vcf \
  --read_filter BadCigar -glm BOTH \
  -stand_call_conf 30.0 -stand_emit_conf 0
$ samtools mpileup -DSugf genome.fasta sample.realn.bam | \
  bcftools view -Ncvg - > sample.samtools.raw1.vcf

6.2 对Variant结果进行综合和过滤

$ java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \
  -R genome.fasta -T SelectVariants \
  --variant sample.gatk.raw1.vcf \
  --concordance sample.samtools.raw1.vcf \
  -o sample.concordance.raw1.vcf
$ MEANQUAL=`awk '{ if ($1 !~ /#/) { total += $6; count++ } } \
  END { print total/count }' sample.concordance.raw1.vcf`
$ java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \
  -R genome.fasta -T VariantFiltration \
  --filterExpression "QD < 20.0 || ReadPosRankSum < -8.0 || \
  FS > 10.0 || QUAL < $MEANQUAL" \
  --filterName LowQualFilter --variant sample.concordance.raw1.vcf \
  --missingValuesInExpressionsShouldEvaluateAsFailing
  --logging_level ERROR -o sample.concordance.flt1.vcf
$ grep -v "Filter" sample.concordance.flt1.vcf > \
  sample.concordance.filter1.vcf

7. 对第5步获得的realign的bam文件的进行校正

$ java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \
  -R genome.fasta -T BaseRecalibrator \
  -I sample.realn.bam -o sample.recal_data.grp \
  -knownSites sample.concordance.filter1.vcf
$ java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \
  -R genome.fasta -T PrintReads \
  -I sample.realn.bam -o sample.recal.bam \
  -BQSR sample.recal_data.grp

8. 第2遍Variant calling

$ java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \
  -R genome.fasta -T UnifiedGenotyper \
  -I sample.recal.bam -o sample.gatk.raw2.vcf \
  --read_filter BadCigar -glm BOTH \
  -stand_call_conf 30.0 -stand_emit_conf 0
$ samtools mpileup -DSugf genome.fasta sample.recal.bam | \
  bcftools view -Ncvg - > sample.samtools.raw2.vcf
$ java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \
  -R genome.fasta -T SelectVariants \
  --variant sample.gatk.raw2.vcf \
  --concordance sample.samtools.raw2.vcf \
  -o sample.concordance.raw2.vcf
$ java -Xmx10g -jar $GATKHome/GenomeAnalysisTK.jar \
  -R genome.fasta -T VariantFiltration \
  --filterExpression "QD < 10.0 || ReadPosRankSum < -8.0 || \
  FS > 10.0 || QUAL < 30" \
  --filterName LowQualFilter --variant sample.concordance.raw2.vcf \
  --missingValuesInExpressionsShouldEvaluateAsFailing
  --logging_level ERROR -o sample.concordance.flt2.vcf
$ grep -v "Filter" sample.concordance.flt2.vcf > \
  sample.concordance.filter2.vcf
$ cp sample.concordance.filter2.vcf sample.final.vcf

最后的结果文件为sample.final.vcf。关于具体的vcf的格式详解见:http://www.hzaumycology.com/chenlianfu_blog/?p=1514

15 thoughts on “结合GATK和SAMtools从头挖掘SNPs和INDELs

  1. 您好,我们测序的样本很多,需要进行SNP calling。在比对这一步时,是将每个样本的pair-end数据分别比对到基因组吗?后期什么时候在将数据合并?

    • 貌似并不一定需要合并。比如 samtools mpileup -DSugf genome.fasta smp1.bam smp2.bam smp3.bam。 GATK可能也能这样输入。 需要自己摸索下。

      • 还有一个问题,就是测序的样本中有亲本P1,P2,它们也和个体那样比对到参考基因组上吗?对P1,P2有没有什么特殊的要求?谢谢!

        • 这个就看你的目的了。软件是死的,你首先要能明白软件能干什么事情。然后你想解决什么问题,就怎么去使用软件,怎么设置参数等~~

  2. 您好,按照您上面的流程进行SNP calling,到“将sam文件转换为Bam文件并标记出PCR duplicates”这一步时,输入命令$ java -Xmx10g -jar $picardHome/MarDuplicates.jar \
    INPUT=sample.sort.bam OUTPUT=sample.dd.bam \
    METRICS_FILE=sample.dd.metrics \
    MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000
    运行完后ls命令没有看到sample.dd.bam这个文件,出现了什么问题?

  3. 楼主您好!
    感谢您的无私分享,我从您的博客中学习到了很多,成功运行了程序,但是在过程中遇到了一个小问题需要请教您一下。
    我使用的是samtool v1.3.1对拟南芥RNA数据进行SNP、INDEL分析。
    命令:
    samtools mpileup -ugf TAIR10.fa SRR671946_sorted.bam > SRR671946.bcf
    bcftools view -cgv SRR671946.bcf > SRR671946.vcf

    在查看.vcf格式中,第8列 INFO里面的显示内容为:
    DP=1;I16=1,0,0,0,22,484,0,0,60,3600,0,0,0,0,0,0;QS=1,0;MQ0F=0
    并没有DP4,导致之后的过滤没有办法进行。
    我想请教一下,这个DP4要怎么才能弄出来呢?是和我的参数设置有关吗?

    谢谢!

发表评论

电子邮件地址不会被公开。 必填项已用*标注

此站点使用Akismet来减少垃圾评论。了解我们如何处理您的评论数据