gatk过滤_快速入门GATK | Public Library of Bioinformatics
GATK,全称是Genome Anlysis Toolkit,顾名思义,是一套用于分析基因组的工具箱。主要功能是寻找变异位点和基因分型,但是实际上功能超多,导致初学者都不知道从何学习GATK。
最近因为mapping-by-sequencing要寻找variant,所以接触了GATK。我相信很多人第一眼看到GATK是茫然的,因为它的功能实在是太多了,都不知道从何开始。这里就说下我是如何在一脸茫然的情况下学习GATK。
GATK的功能虽然超级多,但主要可以归为以下几个方面:诊断和质量控制工具(Diagnostics and Quality Control Tools)
序列数据处理工具(Sequence Data Processing Tools)
变异位点探索工具(Variant Discovery Tools)
变异位点评估工具(Variant Evaluation Tools)
变异位点操作工具(Variant Manipulation Tools)
注释模块
读段(reads)过滤
资源文件解码工具
参考序列实用工具
如何快速建立GATK的心理表征
这里面的每一项点开都有好多内容,我第一次点开的时候,也是一脸茫然,不知道从何入手。
但是根据《认知学习法》,最好的学习方式就是“不要怂,直接上”,找到一个已有流程,先把代码敲上去,然后慢慢理解每一行代码的作用,建立一个模糊的心理表征,然后循序渐进,慢慢学习其他工具,最后就能熟练使用GATK了。
所以记下来主要的任务,就是带大家建立关于GATK的模糊概念。
mapping-by-sequencing其中一个重要环节就是“SNP calling”,我最初用的是samtools和bcftools,结果的variant特别多(估计很多是假阳性).虽然最后还是找到了causual mutation, 但是为了保证今后causual mutation的准确性,我发现了有文章使用了GATK。他给的代码如下:1. Add read groups (Picard tools)
AddOrReplaceReadGroups.jar I=sorted.bam_file O=s1.rg.bam RGLB=genome RGPL=ILLUMINA
RGPU=GATKv4 RGSM=sample_name VALIDATION_STRINGENCY=LENIENT
2. Mark duplicates (Picard tools)
MarkDuplicates.jar INPUT=s1.rg.bam OUTPUT=s2.dedup.bam ASSUME_SORTED=TRUE
VALIDATION_STRINGENCY=LENIENT METRICS_FILE=s2.dedup.metrics
3. Index (samtools)
samtools index s2.dedup.bam
4. Realign reads (create intervals first, then do IndelRealigner) (GATK)
GenomeAnalysisTK.jar -I s2.dedup.bam -R ref_file -T RealignerTargetCreator -o s3.intervals
GenomeAnalysisTK.jar -T IndelRealigner -I s2.dedup.bam -R ref_file -targetIntervals s3.intervals -o
s4.realn.bam
5. Unified genotyper (GATK)
GenomeAnalysisTK.jar -T UnifiedGenotyper -R ref_file -I s4.realn.bam -glm BOTH -o s5.UG1.vcf -mbq 30 -nt4
6. Base score recalibrator (GATK)
GenomeAnalysisTK.jar -T BaseRecalibrator -I s4.realn.bam -R ref_file -knownSites s5.UG1.vcf -o s6.recal
7. Print Reads (GATK)
GenomeAnalysisTK.jar -T PrintReads -R ref_file -I s4.realn.bam -BQSR s6.recal -o s7.recal.bam
8. Unified Genotype (GATK)
GenomeAnalysisTK.jar -T UnifiedGenotyper -R ref_file -I s7.recal.bam -glm BOTH -o s8.UG2.vcf -mbq 30
9. Base score recalibrator (GATK)
GenomeAnalysisTK.jar -T BaseRecalibrator -I s4.realn.bam -R ref_file -knownSites s8.UG1.vcf -o s9.recal
10. Print Reads (GATK)
GenomeAnalysisTK.jar -T PrintReads -R ref_file -I s4.realn.bam -BQSR s9.recal -o s10.recal.bam
11. Unified Genotyper (GATK)
GenomeAnalysisTK.jar -T UnifiedGenotyper -R ref_file -I s10.recal.bam -glm BOTH
第一步: 不管三七二十一直接把代码自己敲一遍。
第二步: 读每行代码的解释,理解他的意图用了picard tools的AddOrReplaceReadGroups.jar加了read group,根据以前百度GATK的经验,了解GATK要求bam文件的header必须包含@RG,所以这一步应该是前面比对时候,没有在参数中增加相应部分。所以如果我在比对的时候增加了这个参数,这一步就可以免了。
用picard tools的MarkDuplicates.jar去重。这是因为二代测序有一部桥式PCR扩增底物的过程,在100x以上出现reads重复,很大可能是PCR扩增重复了,所以可以直接去掉。
samtools的index,建立索引,提高检索效率。
先建立indel区间文件,然后对该区域进行reads重比对。这一步我没有经验,然后我在GATK的论坛搜索了这个工具,发现原因是indel会导致附近的错配,所以需要借由这一步降低indel附近的假阳性。但是最新的HaplotypeCaller or MuTect2已经不需要了,这些工具的haplotype组装步奏效果类似。这里的UnifiedGenotyper or the original MuTect.还是要的。
使用UnifiedGenotyper寻找variant。论坛搜索发现,现在推荐用HaplotypeCaller,效果更好。
根据上一步得到vcf文件对第四步得到BAM文件进行碱基重校准,得到校准所需的文件。
使用PrintReads根据第6步得到的文件,对第四步得到的BAM文件进行校准
根据重校准的BAM文件重新找variant。
根据第8步的文件进行校准
输出第二次校准后BAM文件
第三次寻找variant。
第三步:根据上述过程对拟南芥中使用GATK寻找variant建立大致的认识:常规步骤: 先比对,比对后把SAM转成BAM然后排序,建立索引,去除PCR重复
对indel附近进行重新比对
由于拟南芥没有已知的VCF文件让UnifiedGenotyper学习,所以需要通两轮碱基校准和variant calling降低假阳性
第四步: 再看看别的流程,不断修改最初的认识。我通过翻文献又找到了一个perl脚本# 第一行,UnifiedGenotyper
# 输出read-$i.snv-gatk-snp.vcf
&execsys("java -Xmx2G -jar $cfg{gatk_dir}/GenomeAnalysisTK.jar -T UnifiedGenotyper -R $cfg{genome} -I read-$i.recal.bam -D known_snp.conv -A AlleleBalance -stand_call_conf 50.0 -stand_emit_conf 10.0 -dcov 200 -glm SNP -out_mode EMIT_VARIANTS_ONLY -log log.UnifiedGenotyper-snp.$i -o read-$i.snv-gatk-snp.vcf");
# 第二行,VariantFiltration,
# 过滤read-$i.snv-gatk-snp.vcf得read-$i.snv-gatk-snp.fil.vcf
&execsys("java -Xmx2G -jar $cfg{gatk_dir}/GenomeAnalysisTK.jar -T VariantFiltration -R $cfg{genome} -V read-$i.snv-gatk-snp.vcf --clusterWindowSize 10 --filterExpression \"QUAL < 30.0 || QD < 5.0\" --filterName \"HARD_TO_VALIDATE\" -log log.VariantFiltration-snp.$i -o read-$i.snv-gatk-snp.fil.vcf");
# 第三行 UnifiedGenotyper # 输出read-$i.snv-gatk-indel.vcf&execsys("java -Xmx2G -jar $cfg{gatk_dir}/GenomeAnalysisTK.jar -T UnifiedGenotyper -R $cfg{genome} -I read-$i.recal.bam -D known_snp.conv -A AlleleBalance -stand_call_conf 50.0 -stand_emit_conf 10.0 -dcov 200 -glm INDEL -out_mode EMIT_VARIANTS_ONLY -log log.UnifiedGenotyper-indel.$i -o read-$i.snv-gatk-indel.vcf");
# 第四行 VariantFiltration
# 过滤read-$i.snv-gatk-indel.vcf得read-$i.snv-gatk-indel.fil.vcf
&execsys("java -Xmx2G -jar $cfg{gatk_dir}/GenomeAnalysisTK.jar -T VariantFiltration -R $cfg{genome} -V read-$i.snv-gatk-indel.vcf --clusterWindowSize 10 --filterExpression \"QUAL < 10.0\" --filterName \"HARD_TO_VALIDATE_1\" --filterExpression \"MQ0 >= 4 && ((MQ0 / (1.0 * DP)) > 0.1)\" --filterName \"HARD_TO_VALIDATE_2\" -log log.VariantFiltration-indel.$i -o read-$i.snv-gatk-indel.fil.vcf");
#第五行 CombineVariants
# 结合read-$i.snv-gatk-snp.fil.vcf和read-$i.snv-gatk-indel.fil.vcf
&execsys("java -Xmx2G -jar $cfg{gatk_dir}/GenomeAnalysisTK.jar -T CombineVariants -R $cfg{genome} --variant read-$i.snv-gatk-snp.fil.vcf --variant read-$i.snv-gatk-indel.fil.vcf -o read-$i.snv-gatk.vcf.tmp -log log.CombineVariants.$i");
好像和之前建立的模型不太一样呀,这里变成先找variant然后过滤了。
但是,这些流程所用的工具基本都来自于序列数据处理工具(Sequence Data Processing Tools)
变异位点探索工具(Variant Discovery Tools)
变异位点评估工具(Variant Evaluation Tools)
其实都是先找到variant,然后通过各种方法降低假阳性。
小结
总结一下:
GATK常用套路就是先找variant,然后降低假阳性(例如BQSR, VQSR或直接过滤等)
有了这样一个大致印象后,我再去学习GATK Best Practices,就稍微有点思路了。
这些就是我这段日子学习GATK的感悟,不涉及到具体的操作。希望大家能够提供更多GATK相关的pipeline,让我能够不断更新自己的心理表征。
gatk过滤_快速入门GATK | Public Library of Bioinformatics相关推荐
- 曼哈顿算法公式_距离计算方法总结 | Public Library of Bioinformatics
在做分类时常常需要估算不同样本之间的相似性度量(Similarity Measurement),这时通常采用的方法就是计算样本间的"距离"(Distance).采用什么样的方法计算 ...
- bowtie 加mn标签_Bowtie2使用方法与参数详细介绍 - Public Library of Bioinformatics
Bowtie2 使用方法与参数详细介绍 - Public Library of Bioinformatics 懒人必看 Bowtie2 -q --phred33 --sensitive --end-t ...
- java azure blob 查询_快速入门:适用于 Java 的 Azure Blob 存储客户端库 v8 | Microsoft Docs...
您现在访问的是微软AZURE全球版技术文档网站,若需要访问由世纪互联运营的MICROSOFT AZURE中国区技术文档网站,请访问 https://docs.azure.cn. 快速入门:使用 Jav ...
- java redis快速入门_快速入门Redis系列(3)——Redis的JavaAPI操作(附带练习)
作为快速入门Redis系列的第三篇博客,本篇为大家带来的是Redis的JavaAPI操作. 码字不易,先赞后看! Redis的JavaAPI操作 看完了上一篇博客,相信大家对于Redis的数据类型有了 ...
- 运动控制器编程_快速入门 | 篇二十一:运动控制器ZHMI组态编程简介一
点击上方"正运动小助手",随时关注新动态! 运动控制器ZHMI组态编程简介一 今天我们来学习一下,运动控制器的ZHMI组态编程简介.本文主要从产品概述.控制器连接触摸屏使用.HM ...
- [XML-Jsoup]Jsoup_解析_快速入门
xml常见的解析器: 1. JAXP:sun公司提供的解析器,支持dom和sax两种思想2. DOM4J:一款非常优秀的解析器3. Jsoup:jsoup 是一款Java 的HTML解析器,可直接解析 ...
- 前端_快速入门Vue.js框架
文章目录 快速入门Vue.js框架 0.前言 1.Vue.js框架 1.1.Vue简介 1.2.第一个Vue程序 1.3.el:挂载点 2.Vue指令 2.2.v-html 2.3.v-on 2.4. ...
- java akka 教程_快速入门 Akka Java 指南
快速入门 Akka Java 指南 Akka 是一个用于在 JVM 上构建高并发.分布式和容错的事件驱动应用程序的运行时工具包.Akka 既可以用于 Java,也可以用于 Scala.本指南通过描述 ...
- python构建知识库_快速入门:创建知识库 - REST、Python - QnA Maker - Azure Cognitive Services | Microsoft Docs...
您现在访问的是微软AZURE全球版技术文档网站,若需要访问由世纪互联运营的MICROSOFT AZURE中国区技术文档网站,请访问 https://docs.azure.cn. 快速入门:通过 Pyt ...
- gatk过滤_重测序2--看了不后悔的gatk-变异检测
稍微对重测序分析有点了解的都知道gatk是干嘛的,二代变异检测的金标准,就像一个百宝箱一样,里面藏掖着众多的工具包,看文档都够看一两个月的,但是没有必要.最近事情特别多,不搞废话,变异检测流程可以参考 ...
最新文章
- java 调用存储过程structdescriptor,Java调用oracle存储过程,集合入参的正确姿势
- 解决fragment replace 重叠现象
- buu [BJDCTF 2nd]rsa0
- 印度市场,圆不了二线国产手机的美梦
- boost::mpi模块实现测试mpi版本
- linux编码 form表单,Linux curl 模拟form表单提交信息和文件
- hitchhiker部署_《 Hitchhiker的React Router v4指南》:路由配置的隐藏值
- linux7 共享盘创建,使用CentOS7建立samba文件共享服务器
- resultset需要关闭吗_Java程序员都需要懂的「反射」
- PostgreSQL11.7逻辑复制压测
- SQL Server创建视图
- Atitit.跨平台预定义函数 魔术方法 魔术函数 钩子函数 api兼容性草案 v2 q216 java c# php js.docx
- html5绘制警告牌,2.10 创建自定义图形:绘制扑克牌花色 - HTML5 Canvas 实战
- 如何用运营思维,搭建会员运营体系
- 【Excel文件合并工具】
- 51视频编辑压缩官网
- 杭州电子科技大学 计算机专硕复试分数线,2020杭州电子科技大学考研复试分数线已公布...
- 实验吧 因缺思汀的绕过 By Assassin(with rollup统计)
- 计算机网络 自顶向下方法 (一) 笔记 总结 第一章 概述
- pycharm 将本地文件添加到library root