目前的从头预测软件大多是基于HMM(隐马尔科夫链)和贝叶斯理论，通过已有物种的注释信息对软件进行训练，从训练结果中去推断一段基因序列中可能的结构，在这方面做的最好的工具是AUGUSTUS它可以仅使用序列信息进行预测，也可以整合EST, cDNA, RNA-seq数据作为先验模型进行预测。
安装
安装较为复杂，可选用conda进行安装
conda install augustus
使用
1.若有被训练的物种(augustus --species=help查看)
#使用一下代码进行预测基因，以拟南芥为例
augustus --species=arabidopsis test.fa > test.gff
2.若不存在被训练过的物种，需要进行训练
-
准备训练集和测试集
根据Augutus的官方教程，可靠的基因结构序列的要求如下：
a. 提供基因的编码部分，包含上游几KB。通常而言，基因越多，效果越好，至少准备200个基因以上。还得保证这些基因中要有足够多的外显子，这样子才能训练内含子
b. 这些基因的基因结构一定要足够的准确。不过，也不需要百分百的正确，甚至注释都不需要特别的完整，只要保证起始密码子和终止密码子的准确是准确的即可。
c. 需要保证这些基因没有冗余，也就是说不同序列如果有几乎相同的注释后氨基酸序列，那么仅仅取其中一个（AUGUSTUS教程的建议是：保证任意两个基因在氨基酸水平上低于70%的相似度），这一步既可以避免过度拟合现象，也能用于检验预测的准确性
d. 一条序列允许有多个基因，基因可以在正链也可以在负链，但是这些基因间不能有重叠，每个基因只要其中一个转录本，存放格式是GenBank
之后随机将注释数据集分成训练集和测试集，为了保证测试集有统计学意义，因此测试集要足够多的基因（100~200个），并且要足够的随机。
基因结构集的可能来源有:
Genbank
EST/mRNA-seq的可变剪切联配, 如PASA
临近物种蛋白的可变剪切联配，如GeneWise
相关物种的数据
预测基因的迭代训练
-
训练流程
（1）格式转换；基于选取物种的GFF3以及ref.fa 文件将其转换为Genbank格式
# 以菠菜为例
perl gff2gbSmallDNA.pl spinach_gene_v1.gff3 spinach_genome_v1.fa 1000 genes.raw.gb
(2) 尝试训练，捕捉错误
etraining --species=generic --stopCodonExcludedFromCDS=false genes.raw.gb 2> train.err
(3) 过滤掉可能错误掉基因结构
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/\"//g' |awk '{print $1}' >geneSet.lst
seqkit grep -f geneSet.lst spinach_pep_v1.fa >geneSet.lst.fa
(5) 将得到的蛋白序列进行建库，自身blastp比对。根据比对结果，如果基因间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 转换为genbank 格式文件
perl gff2gbSmallDNA.pl gene_filter.gff3 spinach_genome_v1.fa 1000 genes.gb.filter
(7) 将上一步过滤后的文件随机分成两份，测试集和训练集。其中训练集的数目根据gb的LOCUS数目决定，至少要有200
## 100 为测试集的基因数目，其余为训练集
randomSplit.pl genes.gb.filter 100
(8) 初始化HMM参数设置（在相应～/minicode/config/species/relative name中形成参数,若之前已经存在该物种名字，则需要删除），并进行训练
new_species.pl --species=spinach
etraining --species=spinach genes.gb.filter.train
(9) 用测试数据集检验预测效果，这里可以比较我们训练的结果，和近缘已训练物种的训练效果
augustus --species=spinach genes.gb.filter.test | tee firsttest.out
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) 很有可能的一种情况是，我们第一次的训练结果没有已有训练的效果好，所以我们需要进行循环训练找到最优参数；
optimize_augustus.pl --species=spinach --rounds=5 --cpus=8 genes.gb.filter.train
--rounds: 进行r轮优化，默认为5
--cpus：使用n个cups进行计算，默认为1；该直与kfold相等即可，kfold默认为8
(11) 再次进行训练，并检验，进行前后比较
etraining --species=spinach genes.gb.filter.train
augustus --species=spinach genes.gb.filter.test | tee secondtest.out
- 如果此时你的gene level的sensitivity还是低于20%说明Trainning set不够大，请添加数据；
(12) 进行CRF训练
CRF： conditional random field， 进行进行CRF时，将备份/species/target species/中的所有参数
etraining --species=spinach --CRF=1 genes.gb.filter.train
augustus --species=spinach genes.gb.filter.test | tee secondtest.out.withCRF
比较CRF前后预测的精确性，若升高则使用，若降低，则用上一步结果
- 如果你获得了满意的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
参考
