使用MAKER進(jìn)行注釋: 訓(xùn)練SNAP基因模型

準(zhǔn)備階段

訓(xùn)練SNAP模型蛉拙，需要準(zhǔn)備三個(gè)文件，分別是參考基因組序列画舌，組裝的轉(zhuǎn)錄本序列和同源蛋白序列。

對(duì)于參考基因組序列已慢，我們要保證N50需要超過(guò)預(yù)期基因長(zhǎng)度的中位數(shù)曲聂，否則注釋效果不好。不過(guò)目前的基因組在三代測(cè)序的加持下佑惠，基本上都不是問(wèn)題朋腋。

對(duì)于同源蛋白，建議只用UniProt/Sprot的人工檢查過(guò)的高質(zhì)量蛋白序列膜楷，而不是盲目參考文獻(xiàn)旭咽，使用近源物種的所有蛋白。除非你的近源物種都是模式物種赌厅，蛋白序列可靠性高穷绵，否則用錯(cuò)誤的輸入進(jìn)行訓(xùn)練，數(shù)據(jù)越多反而錯(cuò)的越多特愿。

我們可以在ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/taxonomic_divisions/選擇合適的uniprot_sprot數(shù)據(jù), 然后將其輸出為fasta格式仲墨。以植物為例

wget ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/taxonomic_divisions/uniprot_sprot_plants.dat.gz
zcat uniprot_sprot_plants.dat.gz |\
    awk '{if (/^ /) {gsub(/ /, ""); print} else if (/^AC/) print ">" $2}' |\
    sed 's/;$//'> protein.fa

對(duì)于轉(zhuǎn)錄本，我們通常會(huì)測(cè)一些轉(zhuǎn)錄組數(shù)據(jù)揍障，有三種策略可以得到轉(zhuǎn)錄本目养。（這里暫時(shí)不考慮三代全長(zhǎng)轉(zhuǎn)錄本）

Trinity重頭組裝轉(zhuǎn)錄本
使用STAR + Trinity 獲取轉(zhuǎn)錄本
使用STAR + StringTie + gffread 獲取轉(zhuǎn)錄本

對(duì)于這三種策略，不推薦策略一毒嫡，因?yàn)樵谟袇⒖蓟蚪M的情況下癌蚁，策略二不但計(jì)算效率高，而且能避免組裝錯(cuò)誤（多倍體等位基因之間相似度高）兜畸。對(duì)于策略二和策略三努释，我會(huì)推薦策略三。因?yàn)閷?duì)于靠的比較近的基因膳叨，Trinity很可能會(huì)把這兩個(gè)基因裝成一個(gè)洽洁。

Fig1

并且利用該轉(zhuǎn)錄本作為輸入訓(xùn)練SNAP模型痘系，之后以SNAP模型作為輸入菲嘴，將轉(zhuǎn)錄組和同源蛋白作為證據(jù)而不是直接用作模型，我們?cè)贆z查maker的結(jié)果, 也會(huì)發(fā)現(xiàn)使用StringTie進(jìn)行組裝的結(jié)果才是對(duì)的汰翠。

Fig2

因此使用策略二不但計(jì)算量大龄坪，而且有些情況下還會(huì)導(dǎo)致過(guò)近的轉(zhuǎn)錄本錯(cuò)誤融合，反而影響了最終效果复唤。盡管你可以通過(guò)調(diào)整參數(shù)健田，或者使用鏈特異性測(cè)序的方式來(lái)來(lái)優(yōu)化結(jié)果，但STAR + StringTie 的組合同樣也可以這樣子做佛纫，因此最終我還是推薦策略三妓局。

當(dāng)然总放，這是我定性通過(guò)IGV瀏覽結(jié)果得出的結(jié)論，樣本小好爬，結(jié)論未必可靠局雄，僅供參考。

訓(xùn)練階段

假設(shè)我們準(zhǔn)備的三個(gè)文件分別命名為, genome.fa, protein.fa 和不同組織來(lái)源的stringtie組裝后的轉(zhuǎn)錄本或GFF文件存炮。

注炬搭，如果是StringTie組裝的GTF文件，需要做如下的轉(zhuǎn)換(感謝上海歐易生物-鮑志貴提供的幫助）

gffread -E sample.gtf -o - | sed -e "s#transcript#match#g" -e "s#exon#match_part#g" > sample.gff

接著使用maker -CTL新建配置文件, 設(shè)置如下選項(xiàng)

genome=genome.fa
est=組織1.fa,組織2.fa,組織3.fa
est_gff=組織1.gff,組織2.gff,組織3.gff
protein=protein.fa
est2genome=1
protein2genome=1

然后使用mpiexec -n 線程數(shù) maker &> run.log運(yùn)行程序穆桂。

處理結(jié)果后宫盔，我們新建一個(gè)snap目錄訓(xùn)練模型

mkdir snap && cd snap
gff3_merge -d ../genome.maker.output/genome_master_datastore_index.log

使用makerzff構(gòu)建輸入文件

maker2zff -c 0.8 -e 0.8 -o 0.8 -x 0.2 genome.all.gff

maker2zff會(huì)根據(jù)AED(-x)和QI值進(jìn)行過(guò)濾，其中QI值一共有9項(xiàng)享完，每一項(xiàng)的含義如下

Length of the 5' UTR
Fraction of splice sites confirmed by an EST alignment (-c)
Fraction of exons that overlap an EST alignmetn(-e)
Fraction of exons that overlap EST or Protein alignments(-o)
Fraction of splice site confrimed by a ab-initio prediction(-a)
Fraction of exons that overlap a ab-initio prediction(-t)
Number of exons in the mRNA
length of the 3' UTR
Length of the protein sequence produced by the mRNA (-l)

如果QI值第二項(xiàng)為-1灼芭，表示沒(méi)有支持該剪切位點(diǎn)的EST證據(jù).

接著構(gòu)建模型

fathom -categorize 1000 genome.ann genome.dna
fathom -export 1000 -plus uni.ann uni.dna
forge export.ann export.dna
hmm-assembler.pl snap . > ../snap.hmm

修改配置，然后重新運(yùn)行驼侠。MAKER會(huì)自動(dòng)處理沖突的部分姿鸿，避免重復(fù)序列屏蔽等的一些重復(fù)計(jì)算。

genome=genome.fa
est=Trinity-GG.fasta
protein=protein.fa
snap=snap.hmm
est2genome=0
protein2genome=0

根據(jù)輸出結(jié)果再一次訓(xùn)練模型

mkdir snap2 && cd snap2
gff3_merge -d ../genome.maker.output/genome_master_datastore_index.log
fathom -categorize 1000 genome.ann genome.dna
fathom -export 1000 -plus uni.ann uni.dna
forge export.ann export.dna
hmm-assembler.pl snap . > ../snap2.hmm

通常迭代2-3次就夠了倒源，畢竟我們可能還會(huì)訓(xùn)練AUGUSTUS和GeneMark模型苛预，通過(guò)比較多個(gè)模型來(lái)得到最終結(jié)果。

關(guān)于SNAP軟件細(xì)節(jié)可以參考使用MAKER進(jìn)行基因注釋(高級(jí)篇之SNAP模型訓(xùn)練）