picard:处理HaplotypeCaller断点的问题--g
使用GATK4 call snp时总会出错。鉴于GATK4新版的不稳定性,即使有坑也得跳。这里是填坑的地方。
重要信息:所有的vcf一定不要手工修改,否则gatk会报错。尽管你知道那个地方的正确值,但是不用用除了gatk和picard以外的工具修改,除非后续步骤不需要使用gatk.否则gatk会报错。而且gatk的版本和picard的版本必须从一而终,不能中途换版本。否则,会出现各种错误。
对于不同版本的基因组注释文件,生成的vcf,gatk有工具可以从其中一版对应到另一版本。
问题:
1. call vcf到一半,然后程序因为内存不足挂掉了。从头开始call,又得N多天。能不能从断点处继续call?
不能,但是可以查看日志,看call 到哪条染色体断掉的。比如在:chr4:1352054处断掉,输出是文件1 XXX.g.vcf.gz。那就可以直接从chr4开始,继续call 后续文件,输出XXX.2.g.vcf.gz。再合并2次的文件到一个vcf即可。
断点之后,后续执行命令时,增加参数-L XXXX.intervals
gatk --java-options "-Xmx8g -Xms8g" HaplotypeCaller -R /public/home/$genome.fa -I $sample.bam --dbsnp /public/home/$db_snp.vcf --native-pair-hmm-threads 8 -stand-call-conf 30 -L XXXX.intervals -O XXX.2.g.vcf.gz -ERC GVCF
XXXX.intervals的内容格式,根据你的gtf注释文件决定,如果是chr1格式,则文件里面应该是chrn
,如果和我的一样,直接是1这种方式表示染色体号,则如下所示即可。一定要从断掉的那条染色体从新开始,不然后续拼接不上。
XXXX.intervals文件内容
4
5
6
n
1.2 提取前面的失败的文件中可用的染色体的vcf
gatk可以提取指定染色体的vcf文件。(无论是vcf或者是vcf.gz都可以,但是文件名一定要和格式对应。是压缩的,就是gz,这样gatk会自动解压缩)
gatk操作vcf之前需要建立索引文件,picard不一定需要。
- 建立失败vcf的索引 成功的vcf会自动构建索引
#失败在第7条染色体,所以这样命名
gatk IndexFeatureFile -F Ran.7.g.vcf
- 提取前6条正确的vcf,
#$genome是基因组文件
gatk SelectVariants -R $genome -V Ran.7.g.vcf -L ran.intervals -O Ran.1_6.g.vcf
ran.intervals
是数字1到6代表前六条染色体。格式和上面的intervals格式一致。
提取vcf工具:
2. 合并vcf的方法(所有的工具输入文件必须具有相同的样本和重叠群列表。默认情况下,将创建一个索引文件并需要一个序列字典)
- GatherVcfs 此种方法,必须要求按照基因组顺序排列输入文件,而且不能有重叠。
- SortVcf 排序多个vcf,
- MergeVcfs (拯救我们的断点文件拼接工具)
Merges multiple VCF or BCF files into one VCF file. Input files must be sorted by their contigs and, within contigs, by start position. The input files must have the same sample and contig lists. An index file is created and a sequence dictionary is required by default.
拼接的命令
java -jar $picard.jar MergeVcfs I=XXX.2.g.vcf.gz I=XXX.4.g.vcf.gz O=XXX.g.vcf.gz
此处输入文件可以是多个,I
是input,·O·是output.
输入文件的顺序,一定是按照染色体先后顺序输入。