目前的從頭預(yù)測軟件大多是基于HMM(隱馬爾科夫鏈)和貝葉斯理論蜈亩，通過已有物種的注釋信息對軟件進(jìn)行訓(xùn)練括享，從訓(xùn)練結(jié)果中去推斷一段基因序列中可能的結(jié)構(gòu)愉粤，在這方面做的最好的工具是AUGUSTUS它可以僅使用序列信息進(jìn)行預(yù)測杀迹，也可以整合EST, cDNA, RNA-seq數(shù)據(jù)作為先驗?zāi)Ｐ瓦M(jìn)行預(yù)測劣针。
安裝
安裝較為復(fù)雜校镐，可選用conda進(jìn)行安裝
conda install augustus
使用
1.若有被訓(xùn)練的物種(augustus --species=help查看)
#使用一下代碼進(jìn)行預(yù)測基因，以擬南芥為例
augustus --species=arabidopsis test.fa > test.gff
2.若不存在被訓(xùn)練過的物種捺典，需要進(jìn)行訓(xùn)練
-
準(zhǔn)備訓(xùn)練集和測試集
根據(jù)Augutus的官方教程鸟廓，可靠的基因結(jié)構(gòu)序列的要求如下：
a. 提供基因的編碼部分，包含上游幾KB襟己。通常而言引谜，基因越多，效果越好擎浴，至少準(zhǔn)備200個基因以上员咽。還得保證這些基因中要有足夠多的外顯子，這樣子才能訓(xùn)練內(nèi)含子
b. 這些基因的基因結(jié)構(gòu)一定要足夠的準(zhǔn)確贮预。不過贝室，也不需要百分百的正確契讲，甚至注釋都不需要特別的完整，只要保證起始密碼子和終止密碼子的準(zhǔn)確是準(zhǔn)確的即可档玻。
c. 需要保證這些基因沒有冗余怀泊，也就是說不同序列如果有幾乎相同的注釋后氨基酸序列，那么僅僅取其中一個（AUGUSTUS教程的建議是：保證任意兩個基因在氨基酸水平上低于70%的相似度）误趴，這一步既可以避免過度擬合現(xiàn)象霹琼，也能用于檢驗預(yù)測的準(zhǔn)確性
d. 一條序列允許有多個基因，基因可以在正鏈也可以在負(fù)鏈凉当，但是這些基因間不能有重疊枣申，每個基因只要其中一個轉(zhuǎn)錄本，存放格式是GenBank
之后隨機(jī)將注釋數(shù)據(jù)集分成訓(xùn)練集和測試集看杭，為了保證測試集有統(tǒng)計學(xué)意義忠藤，因此測試集要足夠多的基因（100~200個），并且要足夠的隨機(jī)楼雹。
基因結(jié)構(gòu)集的可能來源有:
Genbank
EST/mRNA-seq的可變剪切聯(lián)配, 如PASA
臨近物種蛋白的可變剪切聯(lián)配模孩，如GeneWise
相關(guān)物種的數(shù)據(jù)
預(yù)測基因的迭代訓(xùn)練
-
訓(xùn)練流程
（1）格式轉(zhuǎn)換；基于選取物種的GFF3以及ref.fa 文件將其轉(zhuǎn)換為Genbank格式
# 以菠菜為例
perl gff2gbSmallDNA.pl spinach_gene_v1.gff3 spinach_genome_v1.fa 1000 genes.raw.gb
(2) 嘗試訓(xùn)練贮缅，捕捉錯誤
etraining --species=generic --stopCodonExcludedFromCDS=false genes.raw.gb 2> train.err
(3) 過濾掉可能錯誤掉基因結(jié)構(gòu)
cat train.err | perl -pe 's/.*in sequence (\S+): .*/$1/' >badgenes.lst
filterGenes.pl badgenes.lst genes.raw.gb > genes.gb
(4) 提取上一步顧慮后的genes.db中的蛋白 (其中第4-6步驟榨咐，也有人忽視)
grep '/gene' genes.gb |sort |uniq |sed 's/\/gene=//g' |sed 's/\"http://g' |awk '{print $1}' >geneSet.lst
seqkit grep -f geneSet.lst spinach_pep_v1.fa >geneSet.lst.fa
(5) 將得到的蛋白序列進(jìn)行建庫，自身blastp比對谴供。根據(jù)比對結(jié)果块茁，如果基因間identity >= 70%，則只保留其中之一桂肌，再次得到一個過濾后的gff文件数焊，gene_filter.gff3
makeblastdb -in geneSet.lst.fa -dbtype prot -parse_seqids -out geneSet.lst.fa
blastp -db geneSet.lst.fa -query geneSet.lst.fa -out geneSet.lst.fa.blastp -evalue 1e-5 -outfmt 6 -num_threads 8
python delete_high_identity_gene.py geneSet.lst.fa.blastp Spinach_genome/spinach_gene_v1.gff3
其中，delete_high_identity_gene.py 可自行編寫
(6) 將得到的gene_filter.gff3 轉(zhuǎn)換為genbank 格式文件
perl gff2gbSmallDNA.pl gene_filter.gff3 spinach_genome_v1.fa 1000 genes.gb.filter
(7) 將上一步過濾后的文件隨機(jī)分成兩份崎场，測試集和訓(xùn)練集佩耳。其中訓(xùn)練集的數(shù)目根據(jù)gb的LOCUS數(shù)目決定，至少要有200
## 100 為測試集的基因數(shù)目谭跨，其余為訓(xùn)練集
randomSplit.pl genes.gb.filter 100
(8) 初始化HMM參數(shù)設(shè)置（在相應(yīng)～/minicode/config/species/relative name中形成參數(shù),若之前已經(jīng)存在該物種名字蚕愤，則需要刪除），并進(jìn)行訓(xùn)練
new_species.pl --species=spinach
etraining --species=spinach genes.gb.filter.train
(9) 用測試數(shù)據(jù)集檢驗預(yù)測效果饺蚊，這里可以比較我們訓(xùn)練的結(jié)果萍诱，和近緣已訓(xùn)練物種的訓(xùn)練效果
augustus --species=spinach genes.gb.filter.test | tee firsttest.out
augustus --species=arabidopsis genes.gb.filter.test | tee firsttest_ara.out
在 firsttest.out 的尾部可以查看預(yù)測結(jié)果的統(tǒng)計，首先需要解釋幾個統(tǒng)計學(xué)概念
- TP(True Positive): 預(yù)測為真污呼，事實為真
- FP(False Positive): 預(yù)測為真裕坊，事實為假
- FN(False Negative): 預(yù)測為假，事實為真
- TN(True Negative): 預(yù)測為假燕酷，事實為假
基于上述籍凝，引出下面兩個概念周瞎。"sensitivity"等于TP/(TP+FP)（預(yù)測到的百分率）， 是預(yù)測為真且實際為真的占你所有認(rèn)為是真的比例."specificity"等于TN/(TN+FN)（其中正確的百分率）, 是預(yù)測為假且實際為假的占你所有認(rèn)為是假的比例饵蒂。我們希望在預(yù)測中声诸，盡可能地不要發(fā)生誤判，也就是沒有基因的地方不要找出基因退盯，有基因的地方不要漏掉基因彼乌。
(10) 很有可能的一種情況是，我們第一次的訓(xùn)練結(jié)果沒有已有訓(xùn)練的效果好渊迁，所以我們需要進(jìn)行循環(huán)訓(xùn)練找到最優(yōu)參數(shù)慰照；
optimize_augustus.pl --species=spinach --rounds=5 --cpus=8 genes.gb.filter.train
--rounds: 進(jìn)行r輪優(yōu)化，默認(rèn)為5
--cpus：使用n個cups進(jìn)行計算琉朽，默認(rèn)為1毒租；該直與kfold相等即可，kfold默認(rèn)為8
(11) 再次進(jìn)行訓(xùn)練箱叁，并檢驗墅垮，進(jìn)行前后比較
etraining --species=spinach genes.gb.filter.train
augustus --species=spinach genes.gb.filter.test | tee secondtest.out
- 如果此時你的gene level的sensitivity還是低于20%說明Trainning set不夠大，請?zhí)砑訑?shù)據(jù)耕漱；
(12) 進(jìn)行CRF訓(xùn)練
CRF： conditional random field算色， 進(jìn)行進(jìn)行CRF時，將備份/species/target species/中的所有參數(shù)
etraining --species=spinach --CRF=1 genes.gb.filter.train
augustus --species=spinach genes.gb.filter.test | tee secondtest.out.withCRF
比較CRF前后預(yù)測的精確性孤个，若升高則使用，若降低沛简，則用上一步結(jié)果
- 如果你獲得了滿意的Trainning結(jié)果齐鲤，請開始prediction吧
下面命令可用于從 firsttest.out 中提取氨基酸序列
sed -n '/^#/p' firsttest.out | sed -n '/start/,/\]/p' | sed 's/# start gene />/g;s/protein sequence \= \[//g;s/#//g;s/\]//g;s/^\s//g' >seq.fa
參考
