準(zhǔn)備訓(xùn)練集和測試集

根據(jù)Augutus的官方教程，可靠的基因結(jié)構(gòu)序列的要求如下：

• 提供基因的編碼部分，包含上游幾KB。通常而言，基因越多，效果越好，至少準(zhǔn)備200個(gè)基因以上。還得保證這些基因中要有足夠多的外顯子，這樣子才能訓(xùn)練內(nèi)含子。

• 這些基因的基因結(jié)構(gòu)一定要足夠的準(zhǔn)確。不過，也不需要百分百的正確，甚至注釋都不需要特別的完整，只要保證起始密碼子和終止密碼子的準(zhǔn)確是準(zhǔn)確的即可。

• 需要保證這些基因沒有冗余，也就是說不同序列如果有幾乎相同的注釋后氨基酸序列，那么僅僅取其中一個(gè)（AUGUSTUS教程的建議是：保證任意兩個(gè)基因在氨基酸水平上低于70%的相似度），這一步既可以避免過度擬合現(xiàn)象，也能用于檢驗(yàn)預(yù)測的準(zhǔn)確性

• 一條序列允許有多個(gè)基因，基因可以在正鏈也可以在負(fù)鏈，但是這些基因間不能有重疊，每個(gè)基因只要其中一個(gè)轉(zhuǎn)錄本，存放格式是GenBank

之后隨機(jī)將注釋數(shù)據(jù)集分成訓(xùn)練集和測試集，為了保證測試集有統(tǒng)計(jì)學(xué)意義，因此測試集要足夠多的基因（100~200個(gè)），并且要足夠的隨機(jī)。

基因結(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 \ # 運(yùn)行alignment/assembly 流程
      -g reference.fasta \
      -t Trinity.fasta \
      --ALIGNERS blat,gmap \ # 可以只用一個(gè)
      --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

最后會(huì)生成一系列文件,如下

• <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會(huì)忽略掉UTR序列。

過濾可能錯(cuò)誤的基因結(jié)構(gòu)

轉(zhuǎn)錄組數(shù)據(jù)不像以前的EST序列，最后組裝的基因數(shù)不多。如果同時(shí)使用多個(gè)組織的轉(zhuǎn)錄組，就會(huì)組裝出一萬條以上的基因，所以我們更加關(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)練，捕捉錯(cuò)誤

augustus/bin/etraining \
    --species=for_bad_genes_removing \
    --stopCodonExcludedFromCDS=false \
    <prefix>.gene.raw.gb 2> train.err

提取錯(cuò)誤，并進(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)練

將上一步過濾后的文件隨機(jī)分成兩份，測試集和訓(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àn)預(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)計(jì)。要看著最后的結(jié)果，需要解釋幾個(gè)統(tǒng)計(jì)學(xué)概念:

• TP(True Positive): 預(yù)測為真，事實(shí)為真
• FP(False Positive): 預(yù)測為真，事實(shí)為假
• FN(False Negative): 預(yù)測為假，事實(shí)為真
• TN(True Negative): 預(yù)測為假，事實(shí)為假

基于上述，引出下面兩個(gè)概念。"sensitivity"等于TP/(TP+FP)， 是預(yù)測為真且實(shí)際為真的占你所有認(rèn)為是真的比例."specificity"等于TN/(TN+FN), 是預(yù)測為假且實(shí)際為假的占你所有認(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ù)未必會(huì)比前一步的訓(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)比較好的一個(gè)，作為最終用于預(yù)測的HMM文件。