Augustus进⾏基因注释
⽬前的从头预测软件⼤多是基于HMM(隐马尔科夫链)和贝叶斯理论，通过已有物种的注释信息对软件进⾏训练，从训练结果中去推断⼀段基因序列中可能的结构，在这⽅⾯做的最好的⼯具是AUGUSTUS它可以仅使⽤序列信息进⾏预测，也可以整合EST, cDNA, RNA-q数据作为先验模型进⾏预测。
安装
安装较为复杂，可选⽤conda进⾏安装
使⽤
（1）若存在已经被训练的物种（augustus --species=help查看），则直接使⽤⼀下代码进⾏预测基因，以拟南芥为例：
1 augustus --speices=arabidopsis test.fa > test.gff
（2）若不存在被训练过的物种，则需要进⾏训练
准备训练集和测试集
根据Augutus的官⽅教程，可靠的基因结构序列的要求如下：
moutain top
a. 提供基因的编码部分，包含上游⼏KB。通常⽽⾔，基因越多，效果越好，⾄少准备200个基因以上。还得保证这些基因中要有⾜够多
的外显⼦，这样⼦才能训练内含⼦
b. 这些基因的基因结构⼀定要⾜够的准确。不过，也不需要百分百的正确，甚⾄注释都不需要特别的完整，只要保证起始密码⼦和终⽌
密码⼦的准确是准确的即可。
c. 需要保证这些基因没有冗余，也就是说不同序列如果有⼏乎相同的注释后氨基酸序列，那么仅仅取其中⼀个（AUGUSTUS教程的建
议是：保证任意两个基因在氨基酸⽔平上低于70%的相似度），这⼀步既可以避免过度拟合现象，也能⽤于检验预测的准确性
d. ⼀条序列允许有多个基因，基因可以在正链也可以在负链，但是这些基因间不能有重叠，每个基因只要其中⼀个转录本，存放格式是
插叙的作用GenBank
之后随机将注释数据集分成训练集和测试集，为了保证测试集有统计学意义，因此测试集要⾜够多的基因（100~200个），并且要⾜够的随机。
基因结构集的可能来源有:
迈克英文
legacya. Genbank
b. EST/mRNA-q的可变剪切联配, 如PASA
c. 临近物种蛋⽩的可变剪切联配，如GeneWi
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=fal genes.raw.gb 2>
（3）过滤掉可能错误掉基因结构
1 | perl -pe 's/.*in quence (\S+): .*/$1/' >badgenes.lst
stress是什么意思
2 filterGenes.pl badgenes.lst genes.raw.gb > genes.gb
（4）提取上⼀步顾虑后的genes.db中的蛋⽩ (其中第4-6步骤，也有⼈忽视)
1 grep '/gene' genes.gb |sort |uniq |d 's/\/gene=//g' |d '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 -par_qids -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 格式⽂件vela
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=spinachkapo
2 etraining --species=spinach genes.ain
（9）⽤测试数据集检验预测效果，这⾥可以⽐较我们训练的结果，和近缘已训练物种的训练效果
1 augustus --species=spinach genes.st | tee firsttest.out
2 augustus --species=arabidopsis genes.st | tee firsttest_ara.out
在 firsttest.out 的尾部可以查看预测结果的统计，⾸先需要解释⼏个统计学概念今年考研时间
TP(True Positive): 预测为真，事实为真
gag
FP(Fal Positive): 预测为真，事实为假
FN(Fal Negative): 预测为假，事实为真
TN(True Negative): 预测为假，事实为假
基于上述，引出下⾯两个概念。"nsitivity"等于TP/(TP+FP)（预测到的百分率），是预测为真且实际为真的占你所有认为是真的⽐例."specificity"等于TN/(TN+FN)（其中正确的百分率）, 是预测为假且实际为假的占你所有认为是假的⽐例。我们希望在预测中，尽可能地不要发⽣误判，也就是没有基因的地⽅不要找出基因，有基因的地⽅不要漏掉基因。
（10）很有可能的⼀种情况是，我们第⼀次的训练结果没有已有训练的效果好，所以我们需要进⾏循环训练找到最优参数；（运⾏会⾮常费时间，⽽且最终的效果⼀般只能提⾼准确度⼏个百分点，慎重使⽤）
1 optimize_augustus.pl --species=spinach genes.ain
（11）再次进⾏训练，并检验，进⾏前后⽐较
1 etraining --species=spinach genes.ain
2 augustus --species=spinach genes.st | tee condtest.out
如果此时你的gene level的nsitivity还是低于20%说明Trainning t不够⼤，请添加数据；
gre词汇下载如果你获得了满意的Trainning结果，请开始吧
下⾯命令可⽤于从 firsttest.out 中提取氨基酸序列
d -n '/^#/p' firsttest.out | d -n '/start/,/\]/p' | d 's/# start gene />/g;s/protein quence \= \[//g;s/#//g;s/\]//g;s/^\s//g' >q.fa
参考
------END-------
关注下⽅公众号可获得更多精彩
