注册 登录  
 加关注
   显示下一条  |  关闭
温馨提示!由于新浪微博认证机制调整,您的新浪微博帐号绑定已过期,请重新绑定!立即重新绑定新浪微博》  |  关闭

云之南

风声,雨声,读书声,声声入耳;家事,国事,天下事,事事关心

 
 
 

日志

 
 
关于我

专业背景:计算机科学 研究方向与兴趣: JavaEE-Web软件开发, 生物信息学, 数据挖掘与机器学习, 智能信息系统 目前工作: 基因组, 转录组, NGS高通量数据分析, 生物数据挖掘, 植物系统发育和比较进化基因组学

网易考拉推荐

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

2014-04-10 13:52:37|  分类: 生信分析软件 |  标签: |举报 |字号 订阅

  下载LOFTER 我的照片书  |
http://blog.sciencenet.cn/blog-252888-731246.html

结合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进行过滤。

perlNGSQCToolkitHome/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?jarpicardHome/SortSam.jar \  INPUT=sample.bam OUTPUT=sample.sort.bam \  SORT_ORDER=coordinatejava?Xmx10g?jarpicardHome/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?jarpicardHome/CreateSequenceDictionary.jar \  R=genome.fasta O=genome.dict$ samtools index sample.dd.bam
 
5. 对INDEL周围进行重新比对(realignment)
 
java?Xmx10g?jarGATKHome/GenomeAnalysisTK.jar \  -R genome.fasta -T RealignerTargetCreator \  -I sample.dd.bam -o sample.realn.intervalsjava?Xmx10g?jarGATKHome/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?jarGATKHome/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?jarGATKHome/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?jarGATKHome/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?jarGATKHome/GenomeAnalysisTK.jar \  -R genome.fasta -T BaseRecalibrator \  -I sample.realn.bam -o sample.recal_data.grp \  -knownSites sample.concordance.filter1.vcf
 
java?Xmx10g?jarGATKHome/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?jarGATKHome/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?jarGATKHome/GenomeAnalysisTK.jar \  -R genome.fasta -T SelectVariants \  --variant sample.gatk.raw2.vcf \  --concordance sample.samtools.raw2.vcf \  -o sample.concordance.raw2.vcf
 
java?Xmx10g?jarGATKHome/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

  评论这张
 
阅读(942)| 评论(0)
推荐 转载

历史上的今天

在LOFTER的更多文章

评论

<#--最新日志,群博日志--> <#--推荐日志--> <#--引用记录--> <#--博主推荐--> <#--随机阅读--> <#--首页推荐--> <#--历史上的今天--> <#--被推荐日志--> <#--上一篇,下一篇--> <#-- 热度 --> <#-- 网易新闻广告 --> <#--右边模块结构--> <#--评论模块结构--> <#--引用模块结构--> <#--博主发起的投票-->
 
 
 
 
 
 
 
 
 
 
 
 
 
 

页脚

网易公司版权所有 ©1997-2016