使用GATK4 call snp時總會出錯。鑒于GATK4新版的不穩(wěn)定性柜候,即使有坑也得跳搞动。這里是填坑的地方。
重要信息:所有的vcf一定不要手工修改渣刷,否則gatk會報錯鹦肿。盡管你知道那個地方的正確值,但是不用用除了gatk和picard以外的工具修改辅柴,除非后續(xù)步驟不需要使用gatk.否則gatk會報錯狮惜。而且gatk的版本和picard的版本必須從一而終,不能中途換版本碌识。否則碾篡,會出現(xiàn)各種錯誤。
對于不同版本的基因組注釋文件筏餐,生成的vcf开泽,gatk有工具可以從其中一版對應(yīng)到另一版本。
問題:
1. call vcf到一半魁瞪,然后程序因為內(nèi)存不足掛掉了穆律。從頭開始call,又得N多天惠呼。能不能從斷點處繼續(xù)call?
不能峦耘,但是可以查看日志剔蹋,看call 到哪條染色體斷掉的。比如在:chr4:1352054處斷掉辅髓,輸出是文件1 XXX.g.vcf.gz泣崩。那就可以直接從chr4開始,繼續(xù)call 后續(xù)文件,輸出XXX.2.g.vcf.gz洛口。再合并2次的文件到一個vcf即可矫付。
斷點之后,后續(xù)執(zhí)行命令時第焰,增加參數(shù)-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的內(nèi)容格式买优,根據(jù)你的gtf注釋文件決定,如果是chr1格式挺举,則文件里面應(yīng)該是chrn
,如果和我的一樣杀赢,直接是1這種方式表示染色體號,則如下所示即可湘纵。一定要從斷掉的那條染色體從新開始葵陵,不然后續(xù)拼接不上。
XXXX.intervals文件內(nèi)容
4
5
6
n
1.2 提取前面的失敗的文件中可用的染色體的vcf
gatk可以提取指定染色體的vcf文件瞻佛。(無論是vcf或者是vcf.gz都可以脱篙,但是文件名一定要和格式對應(yīng)。是壓縮的伤柄,就是gz绊困,這樣gatk會自動解壓縮)
gatk操作vcf之前需要建立索引文件,picard不一定需要适刀。
- 建立失敗vcf的索引 成功的vcf會自動構(gòu)建索引
#失敗在第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
是數(shù)字1到6代表前六條染色體。格式和上面的intervals格式一致笔喉。
提取vcf工具:
2. 合并vcf的方法(所有的工具輸入文件必須具有相同的樣本和重疊群列表取视。默認(rèn)情況下,將創(chuàng)建一個索引文件并需要一個序列字典)
- 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.
輸入文件的順序,一定是按照染色體先后順序輸入。