準(zhǔn)備訓(xùn)練集和測試集
根據(jù)Augutus的官方教程皮官，可靠的基因結(jié)構(gòu)序列的要求如下：
提供基因的編碼部分叭爱，包含上游幾KB撮躁。通常而言，基因越多买雾，效果越好馒胆，至少準(zhǔn)備200個基因以上。還得保證這些基因中要有足夠多的外顯子凝果，這樣子才能訓(xùn)練內(nèi)含子。
這些基因的基因結(jié)構(gòu)一定要足夠的準(zhǔn)確睦尽。不過器净，也不需要百分百的正確，甚至注釋都不需要特別的完整当凡，只要保證起始密碼子和終止密碼子的準(zhǔn)確是準(zhǔn)確的即可山害。
需要保證這些基因沒有冗余，也就是說不同序列如果有幾乎相同的注釋后氨基酸序列沿量，那么僅僅取其中一個（AUGUSTUS教程的建議是：保證任意兩個基因在氨基酸水平上低于70%的相似度）浪慌，這一步既可以避免過度擬合現(xiàn)象，也能用于檢驗預(yù)測的準(zhǔn)確性
一條序列允許有多個基因朴则，基因可以在正鏈也可以在負(fù)鏈权纤，但是這些基因間不能有重疊，每個基因只要其中一個轉(zhuǎn)錄本乌妒，存放格式是GenBank
之后隨機將注釋數(shù)據(jù)集分成訓(xùn)練集和測試集汹想，為了保證測試集有統(tǒng)計學(xué)意義，因此測試集要足夠多的基因（100~200個）撤蚊，并且要足夠的隨機古掏。
基因結(jié)構(gòu)集的可能來源有:
- Genbank
- EST/mRNA-seq的可變剪切聯(lián)配, 如PASA
- 臨近物種蛋白的可變剪切聯(lián)配，如GeneWise
- 相關(guān)物種的數(shù)據(jù)
- 預(yù)測基因的迭代訓(xùn)練
由于目前的轉(zhuǎn)錄組數(shù)據(jù)比較多侦啸，我先用Trinity對轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行從頭組裝槽唾，然后用PASA將Trinity組裝的轉(zhuǎn)錄本回貼到參考基因組上
Launch_PASA_pipeline.pl \
-c alignAssembly.config \
-C \ # 創(chuàng)建數(shù)據(jù)庫
-R \ # 運行alignment/assembly 流程
-g reference.fasta \
-t Trinity.fasta \
--ALIGNERS blat,gmap \ # 可以只用一個
--CPU 12 &> &> pasa_$(date +%Y-%m-%d-%H-%M).log &
alignAssembly.config可以復(fù)制PASApipeline-v2.3.3/pasa_conf的pasa.alignAssembly.Template.txt, 修改如下內(nèi)容
DATABASE=<__DATABASE__> # MySQL中的數(shù)據(jù)庫名, 也是最后文件名的前綴
# 如下調(diào)整的是PASApipeline-v2.3.3/scripts/validate_alignments_in_db.dbi參數(shù)
script validate_alignments_in_db.dbi
validate_alignments_in_db.dbi:--MIN_PERCENT_ALIGNED=<__MIN_PERCENT_ALIGNED__>
validate_alignments_in_db.dbi:--MIN_AVG_PER_ID=<__MIN_AVG_PER_ID__>
這一步得到的<prefix>.assemblies.fasta和<prefix>.pasa_assemblies.gff3
從PASA組裝中提取ORF(開放閱讀框)
PASApipeline-v2.3.3/scripts/pasa_asmbls_to_training_set.dbi \
--pasa_transcripts_fasta <prefix>.assemblies.fasta.transdecoder.genome.gff3\
--pasa_transcripts_gff3 <prefix>.pasa_assemblies.gff3
最后會生成一系列文件,如下
-
<prefix>.assemblies.fasta.transdecoder.cds/pep/gff3/bed: 雖然不再基因組上，但是根據(jù)轉(zhuǎn)錄本信息光涂，有可能是編碼區(qū)的結(jié)果 -
<prefix>.assemblies.fasta.transdecoder.genome.bed/gff3: 對應(yīng)基因組序列的基因模型
我們需要的是后者庞萍。PASA組裝的GFF3文件的屬性一欄分為：
- complete:
- 5primer_partial: 缺少起始密碼子，翻譯到5'端為止
- 3primer_partial: 缺少終止密碼子忘闻，翻譯到3'端為止
- internal: 缺少起始密碼子和終止密碼子
格式轉(zhuǎn)換
我們基于PASA的結(jié)果挂绰，將GFF3格式轉(zhuǎn)換成Genbank格式
augustus/scripts/gff2gbSmallDNA.pl \
<prefix>.assemblies.fasta.transdecoder.genome.bed/gff3 \ #gff3 from pasa_asmbls_to_training_set.dbi
reference.fa \ # reference genome sequence
1000 \ # flank
<prefix>.gene.raw.gb # output
這一步的gff2gbSmallDNA會忽略掉UTR序列。
過濾可能錯誤的基因結(jié)構(gòu)
轉(zhuǎn)錄組數(shù)據(jù)不像以前的EST序列，最后組裝的基因數(shù)不多葵蒂。如果同時使用多個組織的轉(zhuǎn)錄組交播，就會組裝出一萬條以上的基因，所以我們更加關(guān)注于質(zhì)量践付。
先創(chuàng)建初始化的物種HMM文件
augustus/scripts/new_species.pl \
--species=for_bad_genes_removing \
--AUGUSTUS_CONFIG_PATH=/path/to/augustus/config
然后嘗試訓(xùn)練秦士，捕捉錯誤
augustus/bin/etraining \
--species=for_bad_genes_removing \
--stopCodonExcludedFromCDS=false \
<prefix>.gene.raw.gb 2> train.err
提取錯誤，并進(jìn)行過濾
egrep -o 'ctg[[:digit:]]+_[[:digit:]]+-[[:digit:]]+' train.err > badgenes_list.txt
augustus/scripts/filterGenes.pl badgenes_list.txt <prefix>.gene.raw.gb > <prefix>.genes.gb
初步訓(xùn)練
將上一步過濾后的文件隨機分成兩份永高，測試集和訓(xùn)練集隧土。其中訓(xùn)練集的數(shù)目根據(jù)gb的LOCUS數(shù)目決定，至少要有200.
augustus/scripts/randomSplit.pl <prefix>.gene.gb 200 # 200為測試集的基因數(shù)
創(chuàng)建初始化HMM參數(shù)文件命爬，并進(jìn)行訓(xùn)練
augustus/scripts/new_species.pl --species aha_train1
augustus/bin/etraining --species=<prefix> <prefix>.genes.gb.train | tee train.out
用測試數(shù)據(jù)集檢驗預(yù)測效果曹傀。這里可以比較我們訓(xùn)練的結(jié)果，和近緣已訓(xùn)練物種的訓(xùn)練效果
augustus --species=<prefix> <prefix>.genes.gb.test | tee train_test1.out
augustus --species=arabidopsis <prefix>.genes.gb.test | tee ara_test.out
可以用tail -n 40 xx_test.out查看預(yù)測結(jié)果的統(tǒng)計饲宛。要看著最后的結(jié)果皆愉，需要解釋幾個統(tǒng)計學(xué)概念:
- TP(True Positive): 預(yù)測為真，事實為真
- FP(False Positive): 預(yù)測為真艇抠，事實為假
- FN(False Negative): 預(yù)測為假幕庐，事實為真
- TN(True Negative): 預(yù)測為假，事實為假
基于上述家淤，引出下面兩個概念异剥。"sensitivity"等于TP/(TP+FP)， 是預(yù)測為真且實際為真的占你所有認(rèn)為是真的比例."specificity"等于TN/(TN+FN), 是預(yù)測為假且實際為假的占你所有認(rèn)為是假的比例絮重。我們希望在預(yù)測中冤寿，盡可能地不要發(fā)生誤判，也就是沒有基因的地方不要找出基因青伤，有基因的地方不要漏掉基因疚沐。
優(yōu)化訓(xùn)練參數(shù)
很有可能的一種情況是，我們第一次的訓(xùn)練結(jié)果沒有已有訓(xùn)練的效果好潮模，所以我們需要進(jìn)行循環(huán)訓(xùn)練找到最優(yōu)參數(shù)
augustus/scripts/optimize_augustus.pl --species=<prefix> \
--cpus=12 \ #受限于--kfold=k
--rounds=12 \
<prefix>.genes.gb.train
評估優(yōu)化訓(xùn)練結(jié)果
使用optimize_augustus.pl訓(xùn)練的參數(shù)未必會比前一步的訓(xùn)練結(jié)果好亮蛔，因此需要
augustus/bin/etraining --species=<prefix> <prefix>.genes.gb.train
augustus/bin/augustus --species=<prefix> <prefix>.genes.gb.test | tee train_test2.out
比較前后的"sensitivity"和"specificity"參數(shù)，選擇其中表現(xiàn)比較好的一個擎厢，作為最終用于預(yù)測的HMM文件究流。
