picard:處理HaplotypeCaller斷點的問題--gatk提取指定染色體vcf文件

使用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.
輸入文件的順序,一定是按照染色體先后順序輸入。

最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末锐秦,一起剝皮案震驚了整個濱河市咪奖,隨后出現(xiàn)的幾起案子,更是在濱河造成了極大的恐慌酱床,老刑警劉巖羊赵,帶你破解...
    沈念sama閱讀 207,248評論 6 481
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件,死亡現(xiàn)場離奇詭異扇谣,居然都是意外死亡昧捷,警方通過查閱死者的電腦和手機,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 88,681評論 2 381
  • 文/潘曉璐 我一進(jìn)店門揍堕,熙熙樓的掌柜王于貴愁眉苦臉地迎上來,“玉大人汤纸,你說我怎么就攤上這事衩茸。” “怎么了贮泞?”我有些...
    開封第一講書人閱讀 153,443評論 0 344
  • 文/不壞的土叔 我叫張陵楞慈,是天一觀的道長。 經(jīng)常有香客問我啃擦,道長囊蓝,這世上最難降的妖魔是什么? 我笑而不...
    開封第一講書人閱讀 55,475評論 1 279
  • 正文 為了忘掉前任令蛉,我火速辦了婚禮聚霜,結(jié)果婚禮上,老公的妹妹穿的比我還像新娘珠叔。我一直安慰自己蝎宇,他們只是感情好,可當(dāng)我...
    茶點故事閱讀 64,458評論 5 374
  • 文/花漫 我一把揭開白布祷安。 她就那樣靜靜地躺著姥芥,像睡著了一般。 火紅的嫁衣襯著肌膚如雪汇鞭。 梳的紋絲不亂的頭發(fā)上凉唐,一...
    開封第一講書人閱讀 49,185評論 1 284
  • 那天,我揣著相機與錄音霍骄,去河邊找鬼台囱。 笑死,一個胖子當(dāng)著我的面吹牛读整,可吹牛的內(nèi)容都是我干的玄坦。 我是一名探鬼主播,決...
    沈念sama閱讀 38,451評論 3 401
  • 文/蒼蘭香墨 我猛地睜開眼,長吁一口氣:“原來是場噩夢啊……” “哼煎楣!你這毒婦竟也來了豺总?” 一聲冷哼從身側(cè)響起,我...
    開封第一講書人閱讀 37,112評論 0 261
  • 序言:老撾萬榮一對情侶失蹤择懂,失蹤者是張志新(化名)和其女友劉穎喻喳,沒想到半個月后,有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體困曙,經(jīng)...
    沈念sama閱讀 43,609評論 1 300
  • 正文 獨居荒郊野嶺守林人離奇死亡表伦,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 36,083評論 2 325
  • 正文 我和宋清朗相戀三年,在試婚紗的時候發(fā)現(xiàn)自己被綠了慷丽。 大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片蹦哼。...
    茶點故事閱讀 38,163評論 1 334
  • 序言:一個原本活蹦亂跳的男人離奇死亡,死狀恐怖要糊,靈堂內(nèi)的尸體忽然破棺而出纲熏,到底是詐尸還是另有隱情,我是刑警寧澤锄俄,帶...
    沈念sama閱讀 33,803評論 4 323
  • 正文 年R本政府宣布局劲,位于F島的核電站,受9級特大地震影響奶赠,放射性物質(zhì)發(fā)生泄漏鱼填。R本人自食惡果不足惜,卻給世界環(huán)境...
    茶點故事閱讀 39,357評論 3 307
  • 文/蒙蒙 一毅戈、第九天 我趴在偏房一處隱蔽的房頂上張望苹丸。 院中可真熱鬧,春花似錦苇经、人聲如沸谈跛。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,357評論 0 19
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽感憾。三九已至,卻和暖如春令花,著一層夾襖步出監(jiān)牢的瞬間阻桅,已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 31,590評論 1 261
  • 我被黑心中介騙來泰國打工兼都, 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留嫂沉,地道東北人。 一個月前我還...
    沈念sama閱讀 45,636評論 2 355
  • 正文 我出身青樓扮碧,卻偏偏與公主長得像趟章,于是被迫代替她去往敵國和親杏糙。 傳聞我的和親對象是個殘疾皇子,可洞房花燭夜當(dāng)晚...
    茶點故事閱讀 42,925評論 2 344