訓(xùn)練 ab initio 基因預(yù)測(cè)工具(以SNAP為例)
對(duì)于一個(gè)新的物種而言，你大概率是沒(méi)有一個(gè)高質(zhì)量的基因模型去進(jìn)行基因預(yù)測(cè)。但是我們可以利用EST序列（少部分物種估計(jì)有）庶喜、二代測(cè)序數(shù)據(jù)、同源物種蛋白序列，先直接用Maker做基因注釋?zhuān)M管得到的模型可能不是特別的完美，但可以作為輸入反復(fù)迭代運(yùn)行Maker涛漂，從而提高最終的表現(xiàn)。
這次使用的是下載的練習(xí)數(shù)據(jù)集（見(jiàn)附錄）
cd ~/maker_tutorial/example_02_abinitio
同樣检诗，讓我們先構(gòu)建配置文件匈仗，并修改如下配置
maker -CTL
vim maker_opts.ctl
# modify the following line
genome=pyu_contig.fasta
est=pyu_est.fasta
protein=sp_protein.fasta
est2genome=1
protein2genome=1
這里的"est2genome"和"protein2genome"表示直接從EST序列和同源但序列中推測(cè)基因結(jié)構(gòu)，當(dāng)然這肯定不靠譜逢慌。不過(guò)沒(méi)有關(guān)系悠轩，我們的目標(biāo)是將其作為輸入用于訓(xùn)練而已。
運(yùn)行預(yù)測(cè)程序攻泼，大約需要20分鐘
~/opt/biosoft/maker/bin/maker &> maker.log &
那么下一步就是收集所有的GFF文件火架，整理成SNAP所需的ZFF格式
mkdir snap
cd snap
~/opt/biosoft/maker/bin/gff3_merge -d ../pyu_contig1.maker.output/pyu_contig1_master_datastore_index.log
~/opt/biosoft/maker/bin/maker2zff pyu_contig1.all.gff
于是我們就會(huì)在snap文件下得到"genome.ann"和"genome.dna". 在這兩個(gè)文件的基礎(chǔ)上，我們就可以參考SNAP的文檔開(kāi)始訓(xùn)練
可以先用fathom genome.ann genome.dna -gene-stats了解基因的一些信息忙菠，比如說(shuō)這里的測(cè)試數(shù)據(jù)集就有153個(gè)基因何鸡，幾乎平均的分布在正負(fù)鏈上。
1 sequences
0.525725 avg GC fraction (min=0.525725 max=0.525725)
153 genes (plus=79 minus=74)
5 (0.032680) single-exon
148 (0.967320) multi-exon
130.782104 mean exon (min=3 max=704)
87.851593 mean intron (min=61 max=384)
此外還可以用fathom genome.ann genome.dna -validate檢查下是否有明顯的錯(cuò)誤只搁，這里的153個(gè)基因有106個(gè)warning音比，警告類(lèi)型粗略看了一眼基本都是CDS不完整俭尖。
后續(xù)就可以開(kāi)始參數(shù)預(yù)測(cè)氢惋。步驟是，先用fathom genome.ann genome.dna -categorize 1000將序列分類(lèi)稽犁，這里的1000表示基因兩側(cè)會(huì)額外有1000bp的序列焰望。該參數(shù)推薦使用基因一半的長(zhǎng)度，如果基因比較稠密則要調(diào)低已亥。這一步會(huì)生成如下文件:
- alt.ann, alt.dna (genes with alternative splicing)
- err.ann, err.dna (genes that have errors)
- olp.ann, olp.dna (genes that overlap other genes)
- wrn.ann, wrn.dna (genes with warnings)
- uni.ann, uni.dna (single gene per sequence)
這里只用最后一類(lèi)基因熊赖，也就是每個(gè)序列上只有一個(gè)基因。用fathom uni.ann uni.dna -export 1000 -plus只輸出unigene中正鏈基因虑椎，這一步同樣會(huì)生成四個(gè)文件
- export.aa 每個(gè)基因的蛋白序列
- export.ann 正鏈的基因結(jié)構(gòu)
- export.dna 正鏈的DNA
- export.tx 每個(gè)基因的轉(zhuǎn)錄本
接著讓forge負(fù)責(zé)預(yù)測(cè)參數(shù)震鹉， 由于輸出會(huì)很多，所以建議創(chuàng)建個(gè)文件夾
mkdir params
cd params
forge ../export.ann ../export.dna
cd ..
最后是hmm-assembler.pl構(gòu)建HMM捆姜，即基因模型文件, hmm-assembler pyu params > pyu1.hmm
完成SNAP的模型構(gòu)建后传趾，修改"maker_opts.ctl"用以增加該文件，并不再用est和protein直接推測(cè)基因結(jié)構(gòu)泥技。
snaphmm=pyu1.hmm
est2genome=0
protein2genome=0
再一次運(yùn)行maker
~/opt/biosoft/maker/bin/maker &> maker.log &
這次結(jié)果會(huì)比上一次有很明顯的提升浆兰，你可以重復(fù)上面的代碼從而進(jìn)一步提高SNAP的模型。
