GATK4 多个样本GenotypeGVCFs前用 CombineGVCFs还是GenomicsDBImport

Posted 2023-04-23

参考技术A

我们知道， GATK 4 多个样本joint genotyping用模块 GenotypeGVCFs , 目前 GenotypeGVCFs 只支持以下三种形式的输入文件：
1）a single single-sample GVCF
2）a single multi-sample GVCF created by CombineGVCFs
3）a GenomicsDB workspace created by GenomicsDBImport.
即：单个样本的GVCF文件；由 CombineGVCFs 模块将多个样本的GVCF文件生成在一起的文件；由 GenomicsDB 模块将多个样本GVCF处理生成一起的工作空间。当然这里的GVCF 文件是由 HaplotypeCaller 模块的 -ERC GVCF 或者 -ERC BP_RESOLUTION 参数产生， 如果是其他工具生成的GVCF可能会因为缺少某些 GenotypeGVCFs 需要的重要信息导致出错。

因此对于多个样本joint genotyping，在用GenotypeGVCFs前需要将多个样本的gvcf文件用 CombineGVCFs 方式或者 GenomicsDBImport 方式合并成一个文件，前者是一个总的gvcf文件，后者是一个GenomicsDB工作目录。

按照GATK的官方说明 ， CombineGVCFs 效率比较低，需要许多内存，推荐使用 GenomicsDBImport 。不过需注意， GenomicsDBImport 是需要分开interval计算的，如分染色体计算，其每次只能处理一个interval。鉴于GATK极力推荐 GenomicsDBImport ，我们以染色体 chr10 为例测试 CombineGVCFs 和 GenomicsDBImport 对一个trio家系的外显子数据效果，这两个模块的命令分别如下：
CombineGVCFs ：

GenomicsDBImport :

分别运行上述两个命令，但是 问题来了... ，发现 CombineGVCFs 不到5分钟即可完成，而 GenomicsDBImport 已经超过24个小时还在运行中，而且在其genomicsdb-workspace-path产生的chr10_database.db目录下生成超过185,752个子文件，并还在继续生成中！！！这只是对一个染色体的运行，产生如此多的子文件对于有文件数目limit的存储是一个很大的限制，更重要的是运行时间如此之久。 这是为什么呢？

首先检查了上述命令是没有问题的，Google粗略搜索了下没找到原因，直接去GATK 论坛上询问，给出的 解释 如下，比较简单直接：
GenomicsDBImport is used for samples in the order of thousands. For <1000 samples it is better to use CombineGVCFs.

也就是说 GenomicsDBImport 更适用于1000个样本以上的joint genotyping！好吧，这点在GATK的官方使用文档中并没有说明。带着这个问题的疑虑，我又搜索了下发现其实先前已有很多人问过相同的问题并在GATK论坛上深入讨论过，大体总结如下：

结合GATK和samtools以及picardtools call snp

刚开始学生物信息学，老师给了个以snp为标记来画遗传图的课题，研究了一段时间，开始用bwa+samtools来call snp，师姐以前用这套做过，她建议我用另外的方法来做，于是准备学下用GATK来做snp calling。

call snp首先要有个比较准确的参考基因组，然后有样本，我的样本是杂交产生的F2，下面使自己的一些使用过程和心得体会，

这套流程相对于bwa和samtool来说有所不同，先需要对sample的fasta进行筛选，我自己找了下主要是用NGSQC toolkit下的perl脚本筛选，命令为:

$ IlluQC_PRLL.pl -pe read1.fq reads2.fq 2 5 -c 8 -o QC

这个是加入了环境变量的命令

然后要用bowtie2对参考基因组建立索引文件：

$ bowtie2-build genome.fasta genome

genome为索引文件前缀，此索引文件生成目录默认为当前目录而不是ref genome所在的目录

然后将sample mapping到ref genome上命令如下：

bowtie2 --rg-id sample --rg "PL:ILLUMINA" --rg "SM:samplename" -x genome -1 IlluQC_Filtered_files/sampleR1.fastq.gz_filtered -2 IlluQC_Filtered_files/sampleR2.fastq.gz_filtered -p 8 -S sample.sam

--rg-id 是设定read group ID到sample上 -rg是加上“@RG”的头文件，这个类似于bwa mem中的 -R参数

下一步将得到的sam文件转换成二进制的bam文件，这不用到samtools，命令如下：

$ samtools  view -bS sample.sam  > sample.bam

samtools view是专门处理bam文件和sam文件的，-bS命令是输入文件格式为sam 输出文件格式为bam。

接下来是标记出PCR扩增产生的大量重复片段：

这步用到了picardHome里面的SortSam.jar

命令如下：

$ java -Xmx10g -jar picardHome/SortSam.jar INPUT=sample.bam OUTPUT=sample.sorted.bam SORT_ORDER=coordinate

排序

$ java -Xmx10g -jar picardHome/MarDuplicates.jar INPUT=sample.sorted.bam OUTPUT=sample.dd.bam METRICS_FILE=sample.dd.metrics MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000

再次对genome.fasta索引

$ samtools faidx genome.fasta

$ java -Xmx10g -jar picardHome/CreateSequenceDictionary.jar R=genome.fasta O=genome.dict

$ samtools index sample.dd.bam

此处有一个需要注意的地方就是genome.fasta的前缀genome必须和genome.dict的前缀genome一样，否者下一步会报错

接下来对INDEL周围进行重新比对（realignment）：

开始用刀GATK这个软件了

java -Xmx10g -jar 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

接下来是第一遍call snp：

$ 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 -q 20 -Q 25 -ug -t AD -t DP -f genome.fasta sample.realn.bam | bcftool call -c -V indels -v - > sample.samtools.raw1.vcf