谁用过Augustus基因预测软件

Posted 2023-05-11

参考技术A 　　Augustus的安装和使用参数

　　AUGUSTUS is a program that predicts genes in eukaryotic genomic sequences.

　　1. Augustus的安装

　　Augustus下载：bioinfuni-greifswaldde/augustus/binaries/

　　$ wget bioinfuni-greifswaldde/augustus/binaries/augustus.2.7.tar.gz

　　$ tar zxf augustus.2.7.tar.gz

　　$ cd augustus.2.7

　　$ cd src

　　$ make -j 8

　　$ export AUGUSTUS_CONFIG_PATH=$PWD/../config/ (可以加入到.bashrc中）

　　2. Augustus使用方法

　　2.1 基因预测例子

　　$ augustus --strand=both --genemode=partial --singlestrand=false --hintsfile=hints.gff --extrinsicCfgFile=extrinsic.cfg --protein=on --introns=on --start=on --stop=on --cds=on --codingseq=on --alternatives-from-evidence=true --gff3=on --UTR=on ----outfile=out.gff --species=human genome.fa

　　$ augustus --noprediction=true --species=SPECIES sequences.gb

　　2.2 Augustus使用参数

　　Usage:

　　augustus [parameters] --sepcies=SPECIES queryfilename

　　重要参数：

　　--strand=both, --strand=forward or --strand=backward report predicted genes on both strands, just the forward or just the backward strand.default is 'both'

　　--genemodel=partial, --genemodel=intronless, --genemodel=complete,

　　--genemodel=atleastone or --genemodel=exactlyone partial : allow prediction of incomplete genes at the sequence boundaries (default) intronless : only predict single-exon genes成都网络整合运营推广http://www.yingtaow.com/wltg?like in prokaryotes and some eukaryotes

　　complete : only predict complete genes atleastone : predict at least one complete gene exactlyone : predict exactly one complete gene

　　--singlestrand=true predict genes independently on each strand, allow overlapping genes on opposite strands. This option is turned off by default.

　　--hintsfile=hintsfilename When this option is used the prediction considering hints (ex trinsic information) is turned on. hintsfilename contains the hints in gff format.

　　--extrinsicCfgFile=cfgfilename Optional. This file contains the list of used sources for the hints and their boni and mali. If not specified the file "extrin sic.cfg" in the config directory $AUGUSTUS_CONFIG_PATH is used.

　　--maxDNAPieceSize=n This value specifies the maximal length of the pieces that the sequence is cut into for the core algorithm (Viterbi) to be run. Default is --maxDNAPieceSize=200000.

　　AUGUSTUS tries to place the boundaries of these pieces in the intergenic region, which is inferred by a preliminary prediction. GC-content dependent parameters are chosen for each piece of DNA

　　if /Constant/decomp_num_steps > 1 for that species. This is why this value should not be set very large, even if you have plenty of memory.

　　--protein=on/off

　　--introns=on/off

　　--start=on/off

　　--stop=on/off

　　--cds=on/off

　　--codingseq=on/off Output options. Output predicted protein sequence, introns, start codons, stop codons. Or use 'cds' in addition to 'initial', 'internal', 'terminal' and 'single' exon.

　　The CDS excludes the stop codon (unless stopCodonExcludedFromCDS=false) whereas the terminal and single exon include the stop codon.

　　--AUGUSTUS_CONFIG_PATH=path path to config directory (if not specified as environment var iable)

　　--alternatives-from-evidence=true/false report alternative transcripts when they are suggested by hints

　　--alternatives-from-sampling=true/false report alternative transcripts generated through probabilistic sampling

　　--sample=n --minexonintronprob=p --minmeanexonintronprob=p --maxtracks=n --proteinprofile=filename Read a protein profile from file filename. See section 7 below.

　　--predictionStart=A, --predictionEnd=B A and B define the range of the sequence for which predictions should be found. Quicker if you need predictions only for a small part.

　　--gff3=on/off output in gff3 format.

　　--UTR=on/off predict the untranslated regions in addition to the coding sequence. This currently works only for human, galdieria, toxopl asma and caenorhabditis.

　　--outfile=filename print output to filename instead to standard output. This is useful for computing environments, e.g. parasol jobs, which do not allow shell redirection.

　　--noInFrameStop=true/false Don't report transcripts with in-frame stop codons. Otherwise, intron-spanning stop codons could occur. Default: false

　　--noprediction=true/false

　　If true and input is in genbank format, no prediction is made.

　　Useful for getting the annotated protein sequences. Augustus也可以以 genebank格式文件为输入文件，进行基因预测，并将预测结果和genebank的结果进行比较后 得出一个精确性的统计结果。

　　当然，由于genebank格式文件中有些sequences没有cds的注释结果，因此可以使用该 参数进行检测，从而得到没有cds的序列号，在人为去去除这些没有cds注释的序列，再去进行 预测准确性的评估。

　　--contentmodels=on/off If 'off' the content models are disabled (all emissions unif ormly 1/4). The content models are; coding region Markov chain (emiprobs),

　　initial k-mers in coding region (Pls), intron and int ergenic regin Markov chain. This option is intended for special applications that require judging gene structures from the signal models only,

　　e.g. for predicting the effect of SNPs or mutations on splicing. For all typical gene predictions, this should be true. Default: on

　　--paramlist For a complete list of parameters, type "augustus --paramlist"

Augustus 进行基因注释

 

目前的从头预测软件大多是基于HMM(隐马尔科夫链)和贝叶斯理论，通过已有物种的注释信息对软件进行训练，从训练结果中去推断一段基因序列中可能的结构，在这方面做的最好的工具是AUGUSTUS它可以仅使用序列信息进行预测，也可以整合EST, cDNA, RNA-seq数据作为先验模型进行预测。

 

安装

具体看看这个http://www.chenlianfu.com/?p=1224

 

 

使用

（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. 预测基因的迭代训练

 

流程如下

这一部分有时间可以看看

http://www.chenlianfu.com/?p=2307 还有这个https://www.jianshu.com/p/679bd6bb0ea4

（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

　　

来源：https://www.cnblogs.com/zhanmaomao/p/12359964.html