-
安裝
安裝較為復雜，可選用conda進行安裝
-
使用
（1）若存在已經被訓練的物種（augustus --species=help查看），則直接使用一下代碼進行預測基因，以擬南芥為例：
1 augustus --speices=arabidopsis test.fa > test.gff
（2）若不存在被訓練過的物種，則需要進行訓練
-
准備訓練集和測試集
根據Augutus的官方教程，可靠的基因結構序列的要求如下：
a. 提供基因的編碼部分，包含上游幾KB。通常而言，基因越多，效果越好，至少准備200個基因以上。還得保證這些基因中要有足夠多的外顯子，這樣子才能訓練內含子
b. 這些基因的基因結構一定要足夠的准確。不過，也不需要百分百的正確，甚至注釋都不需要特別的完整，只要保證起始密碼子和終止密碼子的准確是准確的即可。
c. 需要保證這些基因沒有冗余，也就是說不同序列如果有幾乎相同的注釋后氨基酸序列，那么僅僅取其中一個（AUGUSTUS教程的建議是：保證任意兩個基因在氨基酸水平上低於70%的相似度），這一步既可以避免過度擬合現象，也能用於檢驗預測的准確性
d. 一條序列允許有多個基因，基因可以在正鏈也可以在負鏈，但是這些基因間不能有重疊，每個基因只要其中一個轉錄本，存放格式是GenBank
之后隨機將注釋數據集分成訓練集和測試集，為了保證測試集有統計學意義，因此測試集要足夠多的基因（100~200個），並且要足夠的隨機。
基因結構集的可能來源有:
a. Genbank
b. EST/mRNA-seq的可變剪切聯配, 如PASA
c. 臨近物種蛋白的可變剪切聯配，如GeneWise
d. 相關物種的數據
e. 預測基因的迭代訓練
-
流程如下
（1）格式轉換；基於選取物種的GFF3以及ref.fa 文件將其轉換為Genbank格式
1 perl ~/miniconda2/bin/gff2gbSmallDNA.pl ./Spinach_genome/spinach_gene_v1.gff3 ./Spinach_genome/spinach_genome_v1.fa 1000 genes.raw.gb
（2）嘗試訓練，捕捉錯誤；
1 etraining --species=generic --stopCodonExcludedFromCDS=false genes.raw.gb 2> train.err
（3）過濾掉可能錯誤掉基因結構
1 cat train.err | perl -pe 's/.*in sequence (\S+): .*/$1/' >badgenes.lst 2 filterGenes.pl badgenes.lst genes.raw.gb > genes.gb
（4）提取上一步顧慮后的genes.db中的蛋白 (其中第4-6步驟，也有人忽視)
1 grep '/gene' genes.gb |sort |uniq |sed 's/\/gene=//g' |sed 's/\"//g' |awk '{print $1}' >geneSet.lst 2 python extract_pep.py geneSet.lst Spinach_genome/spinach_pep_v1.fa
（5）將得到的蛋白序列進行建庫，自身blastp比對。根據比對結果，如果基因間identity >= 70%，則只保留其中之一，再次得到一個過濾后的gff文件，gene_filter.gff3
1 makeblastdb -in geneSet.lst.fa -dbtype prot -parse_seqids -out geneSet.lst.fa 2 blastp -db geneSet.lst.fa -query geneSet.lst.fa -out geneSet.lst.fa.blastp -evalue 1e-5 -outfmt 6 -num_threads 8 3 python delete_high_identity_gene.py geneSet.lst.fa.blastp Spinach_genome/spinach_gene_v1.gff3
（6）將得到的gene_filter.gff3 轉換為genbank 格式文件
1 perl ~/miniconda2/bin/gff2gbSmallDNA.pl gene_filter.gff3 ./Spinach_genome/spinach_genome_v1.fa 1000 genes.gb.filter
（7）將上一步過濾后的文件隨機分成兩份，測試集和訓練集。其中訓練集的數目根據gb的LOCUS數目決定，至少要有200
1 ## 100 為測試集的基因數目，其余為訓練集 2 randomSplit.pl genes.gb.filter 100
（8）初始化HMM參數設置（在相應～/minicode/config/species/relative name中形成參數,若之前已經存在該物種名字，則需要刪除），並進行訓練
1 new_species.pl --species=spinach 2 etraining --species=spinach genes.gb.filter.train
（9）用測試數據集檢驗預測效果，這里可以比較我們訓練的結果，和近緣已訓練物種的訓練效果
1 augustus --species=spinach genes.gb.filter.test | tee firsttest.out 2 augustus --species=arabidopsis genes.gb.filter.test | tee firsttest_ara.out
在 firsttest.out 的尾部可以查看預測結果的統計，首先需要解釋幾個統計學概念
- TP(True Positive): 預測為真，事實為真
- FP(False Positive): 預測為真，事實為假
- FN(False Negative): 預測為假，事實為真
- TN(True Negative): 預測為假，事實為假
基於上述，引出下面兩個概念。"sensitivity"等於TP/(TP+FP)（預測到的百分率）， 是預測為真且實際為真的占你所有認為是真的比例."specificity"等於TN/(TN+FN)（其中正確的百分率）, 是預測為假且實際為假的占你所有認為是假的比例。我們希望在預測中，盡可能地不要發生誤判，也就是沒有基因的地方不要找出基因，有基因的地方不要漏掉基因。
（10）很有可能的一種情況是，我們第一次的訓練結果沒有已有訓練的效果好，所以我們需要進行循環訓練找到最優參數；（運行會非常費時間，而且最終的效果一般只能提高准確度幾個百分點，慎重使用）
1 optimize_augustus.pl --species=spinach genes.gb.filter.train
（11）再次進行訓練，並檢驗，進行前后比較
1 etraining --species=spinach genes.gb.filter.train 2 augustus --species=spinach genes.gb.filter.test | tee secondtest.out
- 如果此時你的gene level的sensitivity還是低於20%說明Trainning set不夠大，請添加數據；
- 如果你獲得了滿意的Trainning結果，請開始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
參考
------END-------
關注下方公眾號可獲得更多精彩
