提出问题
给出以下的已有结果:
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