结合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

匿名进行回复 取消回复

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

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