java queue GATK_gatk4使用总结
昨天看了gatk的官网,从2018年发布正式版的4.0.0开始,到现在已经更新到4.1.8,在速度和准确度上都有了大幅的提升。gatk4除了整合picard软件之外,在使用上与gatk3基本相同,只不过是在命令运行、功能划分及运行速度上进行了调整。gatk4软件的安装及使用,网上有很多帖子,可以自行查阅。今天主要给大家介绍gatk4与gatk3明显不同的地方,这些使得gatk4的应用更加方便。
1. 输入文件的变化
变异数据是从比对结果中获得的,gatk之前就已经支持cram的输入,但是数据去重部分的picard软件并不支持,运行时会报错。至于cram文件,也是一种比对结果文件,只不过文件体积更小,即对于存储受限的分析人员里说是一个值得考虑的方向,具体的信息可查看公众号文章cram输出及使用。
picard原装软件使用cram运行时报错信息。
Exception in thread "main" java.lang.IllegalStateException: A valid CRAM reference was not supplied and one cannot be acquired via the property settings reference_fasta or use_cram_ref_download
at htsjdk.samtools.cram.ref.ReferenceSource.getDefaultCRAMReferenceSource(ReferenceSource.java:105)
at htsjdk.samtools.SamReaderFactory$SamReaderFactoryImpl.open(SamReaderFactory.java:406)
at picard.sam.markduplicates.util.AbstractMarkDuplicatesCommandLineProgram.openInputs(AbstractMarkDuplicatesCommandLineProgram.java:220)
at picard.sam.markduplicates.MarkDuplicates.buildSortedReadEndLists(MarkDuplicates.java:535)
at picard.sam.markduplicates.MarkDuplicates.doWork(MarkDuplicates.java:264)
at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:295)
at picard.cmdline.PicardCommandLine.instanceMain(PicardCommandLine.java:103)
at picard.cmdline.PicardCommandLine.main(PicardCommandLine.java:113)
不能给定基因组,就意味着不能使用cram文件。而gatk4在整合了picard之后对其进行了优化(猜测,具体是不是整合前picrad进行了更新不清楚,有兴趣的可以去测测看),现在是支持cram的输入和去重的。
需要注意的是,软件运行结果是bam文件,而不是cram文件。
2. gatk4变异检测模块调整
使用gatk4进行变异检测时,主要包括以下5个模块
Name
Summary
CombineGVCFs
Merges one or more HaplotypeCaller GVCF files into a single GVCF with appropriate annotations
GenomicsDBImport
Import VCFs to GenomicsDB
GenotypeGVCFs
Perform joint genotyping on one or more samples pre-called with HaplotypeCaller
HaplotypeCaller
Call germline SNPs and indels via local re-assembly of haplotypes
Mutect2
Call somatic SNVs and indels via local assembly of haplotypes
一般的变异使用HaplotypeCaller进行检测,输出格式可以是vcf,也可以是gvcf
# GVCF输出
gatk --java-options "-Xmx4g" HaplotypeCaller \
-R Homo_sapiens_assembly38.fasta \
-I input.bam \
-O output.g.vcf.gz \
-ERC GVCF
# VCF输出
gatk --java-options "-Xmx4g" HaplotypeCaller \
-R Homo_sapiens_assembly38.fasta \
-I input.bam \
-O output.vcf.gz \
-bamout bamout.bam
vcf适合单样本变异检测,GVCF主要应用在群call变异检测上。对于gvcf文件的处理,可以使用GenotypeGVCFs模块输出变异位点的信息。
gatk --java-options "-Xmx4g" GenotypeGVCFs \
-R Homo_sapiens_assembly38.fasta \
-V input.g.vcf.gz \
-O output.vcf.gz
该方法在使用上有限制,只能输入一个gvcf文件,对于多个样本群call的方式,需要使用另外两个软件对gvcf文件进行处理。
另外还有一点,gatk4在命令上有很大的改进,弃用了直接使用java包的运行模式,命令行相对简单。当然,对于非常熟悉gatk3运行并且不喜欢gatk4的分析人员来说,也可以直接使用安装包下面的jar包来运行程序,使用方法与GATK3基本相同。
3. 群call策略
群call的最终目的是获得所有样本的分型结果,最终会使用GenotypeGVCFs模块转化为vcf,但是由于输入文件个数的限制,需要先使用GenomicsDBImport和GenotypeGVCFs2个模块处理。
CombineGVCFs模块是用来合并多个gvcf文件。
gatk CombineGVCFs \
-R reference.fasta \
--variant sample1.g.vcf.gz \
--variant sample2.g.vcf.gz \
-O cohort.g.vcf.gz
另外一种可以合并gvcf方式是gatk4中新出现的,GenomicsDBImport,该方法先对gvcf文件进行处理,输出一个类似于数据库的目录,后面可以直接使用该数据目录作为输入文件进行变异检测。
gatk --java-options "-Xmx4g -Xms4g" GenomicsDBImport \
-V data/gvcfs/mother.g.vcf.gz \
-V data/gvcfs/father.g.vcf.gz \
-V data/gvcfs/son.g.vcf.gz \
--genomicsdb-workspace-path my_database \
--tmp-dir=/path/to/large/tmp \
-L 20
--genomicsdb-workspace-path参数指定数据库的输出路径,该路径在运行之前不能被创建,否则会报错。合并完成之后可以进行位点变异的检测。
gatk GenotypeGVCFs \
-R data/ref/ref.fasta \
-V gendb://my_database \
-G StandardAnnotation -newQual \
-O test_output.vcf
群call最大的优势在于我们可以添加样本后重新分析,获取群体的变异位点。如果使用GenomicsDBImport进行分析,若要添加新样本的变异数据,只要将新样本的gvcf信息添加到已有的数据库中即可。
gatk GenomicsDBImport \
-V data/gvcfs/mother.g.vcf \
-V data/gvcfs/father.g.vcf \
-V data/gvcfs/son.g.vcf \
--genomicsdb-update-workspace-path existing_database
数据库中含有的文件是二进制的文件,若想要多的数据库中存在的染色体信息,或者想要提取gvcf文件,可以使用如下命令
# 提取gvcf文件
gatk SelectVariants \
-R data/ref/ref.fasta \
-V gendb://my_database \
-O combined.g.vcf
# 提取interval信息
gatk GenomicsDBImport \
--genomicsdb-update-workspace-path existing_database
--output-interval-list-to-file /path/to/output/file
总体来说,gatk4更新之后还是很好用的,相比于gatk3,只要选择合适的分析流程和策略,时间上也会大幅缩短。不过相比于samtools这些软件,gatk分析时间还是最大的问题,不过如果考虑到群call和准确度的优势,特别是大样本量的差异分析,gatk的优势还是很明显的。
4. 参考资料
#如有侵权,请告知删除#
#如有错误,欢迎指正#
java queue GATK_gatk4使用总结相关推荐
- java queue使用_使用Java使用Amazon Simple Queue Service
java queue使用 Amazon Simple Queue Service或SQS是Amazon Webservice堆栈提供的高度可扩展的托管消息队列. Amazon SQS可用于完全解耦系统 ...
- 学习Java: Queue
15 08, 2007 学习Java: Queue Java - 作者 zybing @ 15:17 Java提供了Quere,相当好用,在1.5版本中又有增强. Queue: 基本上,一个队列就是一 ...
- java Queue中 add/offer,element/peek,remove/poll区别
java Queue中 add/offer,element/peek,remove/poll中的三个方法均为重复的方法,在选择使用时不免有所疑惑,这里简单区别一下: 1.add()和offer()区别 ...
- java queue size_Java中的PriorityQueue size() 方法 - Break易站
Java.util.PriorityQueue.size()方法用于获取PriorityQueue的大小或PriorityQueue中存在的元素数. 句法: Priority_Queue.size() ...
- java queue源码_java源码解读--queue
queue接口特点:可以模拟队列行为,即"先进先出". 接口结构 queue接口继承了Collection接口,并增加了一些新方法 1 2 3 4 5 6 7 8 9 10 11 ...
- java queue 实现类 区别_Java集合11 (Queue)
java.util.Queue接口是java.util.Collection子接口. 它代表一个有序的对象列表,就像List一样,但是它的使用有略微的区别. Queue被设计成从末端插入并且从头部删除 ...
- java queue通信_Java -- 使用阻塞队列(BlockingQueue)控制线程通信
BlockingQueeu接口是Queue的子接口,但是它的主要作用并不是作为容器,而是作为线程同步的工具. 特征: 当生产者线程试图向BlockingQueue中放入元素时,如果该队列已满,则该线程 ...
- Java: Queue
Queue: 基本上,一个队列就是一个先入先出(FIFO)的数据结构 offer,add区别: 一些队列有大小限制,因此如果想在一个满的队列中加入一个新项,多出的项就会被拒绝. 这时新的 offer ...
- java Queue
队列是一个典型的先进先出(FIFO)的容器,即从容器的一端放入事物,从另一端取出,并且事物放入容器的顺序与取出的顺序是相同的,队列常常被当作一种可靠的对象从程序的某个区域传输到另一个区域,队列在并发编 ...
- java queue put take_阻塞队列的take、offer、put、add的一些比较
最近在学习<>,有很多java.util.concurrent包下的新类.LinkedBlockingQueue就是其中之一,顾名思义这是一个阻塞的线程安全的队列,底层应该采用链表实现. ...
最新文章
- s-sar命令(System Activity Reporter系统活动情况报告)
- 你和人工智能的对话,正在被人工收听
- ajax 设置Access-Control-Allow-Origin实现跨域访问
- Aroma's Search(暴力)
- [react] react非父子组件如何通信?
- ggplot2实现分半小提琴图绘制基因表达谱和免疫得分
- matplotlib基本设置
- 配置vscode作为STM32代码的编辑器(替代keil5)。实现:代码自动补全, 编译,下载。nRF52也可以编译。
- mov格式怎么在线转换成mp4格式
- md5是什么,md5的这些作用你都知道吗
- 女友的生日礼物能随便嘛?Python小妙招:制作一款出圈九宫格抽奖小程序。
- idea创建三种应用程序的方法:springboot,控制台程序,windows服务程序
- 高德Android 定位SDK示例
- 解决 “计算机中丢失gdiplus.dll”
- 管理模型 - RACI模型
- Linux系统之Xinetd服务
- Android核心功能
- 特征选择对于机器学习重要性
- 基于sqlserver数据库的学生成绩管理系统
- Python自动化与Word