(宏)基因组编码基因预测

2022-05-05 13:36:10 浏览数 (1)

基因预测是指通过对组装的基因组序列进行分析,根据已知生物的基因结构知识或数据库序列来识别其所包含的基因等功能区域。编码基因预测,就是识别基因组序列上所包含的蛋白质编码区域(Coding sequence,CDS),通过在基因组序列上寻找开放阅读框(Open Reading Frame,ORF)实现。

ORF是指从序列5'端的一个起始密码子(ATG)到3'端的一个终止密码子(TTA、TAG、TGA)之间的片段,可以理解为理论上的编码区(不一定所有识别的ORF均是完整的CDS),一般通过计算机进行序列分析得到。当面对一条陌生的DNA序列(尤其是不完整的contigs),由于对其遗传信息完全不清楚,可以有6种方法来尝试解读序列,分别是序列第1、2、3个碱基开始以及反向互补序列的第1、2、3个碱基开始,因此每一个基因有6种框架阅读模式,通常情况下选择中间没有被终止密码子隔开的最大ORF作为基因预测的正确结果。

目前,基因预测的基本方法有2种,基于序列相似性的搜索和基于模式序列特征的从头预测。基于序列相似性的搜索方法思路是将待预测的基因组序列在6种模式的阅读框中进行翻译并与蛋白质数据库中的序列进行比对,如blastx,或者对EST数据库中同一生物的cDNA序列进行比对分析,如blastn,然后确定基因的数目和对应的CDS序列。该方法很依赖于数据库里的数据,对于数据库中没有收录相关基因序列的新物种或者无法确定同源关系的物种,不适用该方法。而且如何界定基因序列的起始、终止位置,尤其真核生物基因的外显子和内含子边界以及筛选比对结果也很重要。

基因的从头预测方法依据人们对已知基因结构特征的认识,如启动子区的TATA box、密码子偏好性等,采用统计学方法,如隐马尔可夫模型、决策树方法、神经网络分析法等,对基因组作基因预测。原核生物基因的各种信号位点(如启动子和终止子信号位点)特异性较强且容易识别,因此相应的基因预测方法已经基本成熟。而真核生物的基因预测工作的难度则大为增加。首先,真核生物中的启动子和终止子等信号位点更为复杂,难以识别。其次,真核生物中广泛存在可变剪切现象,使外显子和内含子的定位更为困难。因此,预测真核生物的基因结构需要运用更为复杂的算法,常用的有隐马尔科夫模型等。

原核生物基因预测软件常用的有GeneMark、Glimmer、Prodigal等,真核生物基因预测软件有GENSCAN、Augustus、GlimmerHMM、PASA等。在这里我们主要介绍原核生物基因预测。

Prodigal

Prodigal也即Prokaryotic Dynamic ProgrammingGenefinding Algorithm(原核生物动态规划基因查找算法,https://github.com/hyattpd/Prodigal/),主要应用于细菌和古生菌基因预测,如果要对meta样品做基因预测,prodigal还专门提供了meta的版本。相对与glimmer基因预测工具,prodigal步骤简单更加好用,而且软件可以直接输出基因的核酸序列并翻译出的相应的氨基酸序列。

基本参数如下所示:

代码语言:javascript复制
-a  输出预测蛋白质的序列文件名
-c  不允许基因一边断开,也就是要求完整的ORF,有起始和终止结构
-d  输出预测基因的序列文件名
-f  选择输出文件格式,有gbk、gff、和sco格式可供选择
-g  指定密码子翻译表(不同生物密码子稍有不同),普通细菌/古菌选11,支原体、螺原体选4,如果为auto则先尝试11然后4,也可以其他(支持第1-25套密码子)
-i  输入文件,即需要预测的基因组序列文件
-m  屏蔽基因组中的N碱基(对于有gap的scaffolds)
-o  预测结果输出文件名,默认为屏幕输出
-p  选择项目性质,是单菌'single'还是宏基因组'meta'
-q  不输出错误信息到屏幕
-t  指定训练集,不指定则使用自身数据创建训练集
-s  输出所有潜在基因及其分值到一个文件中

使用Prodigal对组装的基因组序列进行基因预测:

代码语言:javascript复制
prodigal -a scaffolds.protein.fa -d scaffolds.gene.fa -f gff -g 11 -i new.scaffolds.fasta -o scaffolds.gff -s scaffolds.stat -p single -q -m

运行结束后产生的结果文件中scaffolds.gene.fa和scaffolds.protein.fa分别为预测的基因序列和对应的蛋白质序列,如下所示:

scaffolds.gff为预测结果文件,也即不同scaffold序列预测到的基因及其基本信息,如下所示:

使用Prodigal预测宏基因组编码基因,方法如下所示:

代码语言:javascript复制
prodigal -a prod.meta.protein.fa -d prod.meta.gene.fa -f gff -g 11 -i new.spades.contig.fasta -o prod.meta.gff -s prod.meta.stat -p meta -q

GeneMark

GeneMark(http://topaz.gatech.edu/GeneMark/)是一款主流的基因预测软件,可以使用蛋白质编码序列和非编码序列的Markov模型(及启发式算法Heuristic Model),以及起始位置核苷酸频率矩阵来提高基因预测的准确性,广泛适用于细菌、古菌、宏基因组、宏转录组的基因预测。GeneMark Suite一般有以下几个主程序:

①GeneMark:gm和gm.pl

②GeneMark.hmm:gmhmmp和gmhmmp.pl

③GeneMark.hmm(Heuristic Model):gmhmmp_heuristic.pl

④GeneMarkS:gmsn.pl (GeneMarkS pipeline之一)

除此之外,还有gc(计算GC含量)、mkmat(产生核苷酸频率矩阵模型)、viewmat(查看矩阵的完整内容)、matinfo(展示GeneMark所使用的转移矩阵中的文本信息)等辅助程序。

GeneMark程序基于编码区和非编码去的马尔科夫模型,并采用滑动窗口的方法,预测一条DNA序列中潜在的蛋白质编码区。该方法对编码可能性之间的局部变化非常敏感,但能生成一幅展示编码可能性分布的细节图。GeneMark与GeneMark.hmm程序都需要利用序列中核酸使用的频率矩阵作为基础,来预测序列中潜在的编码区域,这些矩阵都是物种特异的。如果没有合适的矩阵模型,需要使用该物种或近缘物种的编码序列与非编码序列利用软件包里的mkmat命令创建一个新矩阵,要么使用一个近缘物种的矩阵。事实上对于de novo测序可以直接使用GeneMarkS。不同于GeneMark、GeneMark.hmm则将序列作为一个整体,GeneMarkS基于一个可以反映基因组织规则的隐状态网络,采用隐马尔可夫模型,预测基因和基因间的区域。GeneMarkS程序使用方法如下所示:

代码语言:javascript复制
gmsn.pl --gcode 11 --fnn --faa --prok --format GFF3 --output output.gff3 input.fasta
--format 结果文件的格式(default: LST; supported: LST, GFF and GFF3)
--output 输出的结果文件名
--clean 运行结束后删除所有临时文件
--fnn 生成预测基因的核酸序列文件
--faa 生成预测基因的蛋白序列文件
--gcode 密码子表(default: 11; supported: 11, 4, 25 and 1)
--strand 在序列正向还是反向互补预测基因(default: 'both'; supported: direct, reverse and both )
--prok 原核生物基因预测(相当于--combine --clean --gm)
--euk 真核生物内含子较少的生物基因预测(也即低等真核生物)
--virus 真核生物病毒基因预测
--est 使用EST格式序列进行预测
--phage 噬菌体(也即原核生物病毒)基因预测(相当于--combine --clean --gm)

使用GeneMarkS对组装的基因组序列进行基因预测:

代码语言:javascript复制
cp …/genemark/gmsuite/gm_key ~/.gm_key  #拷贝key文件到自己主文件夹
#注意,此license下载后有400天有效期,到期后需要续期或者重新下载安装软件。
gmsn.pl --gcode 11 --fnn --faa --combine --clean --gm --format GFF3 --output gm.scaffolds.gff3 new.scaffolds.fasta

运行结束后,会生成基因序列与对应的蛋白序列。

基因组分析中使用了GeneMarkS预测编码基因,在宏基因组则使用MetaGeneMark。MetaGeneMark利用GeneMark.hmm主程序(gmhmmp)基于自带的核苷酸频率矩阵模型MetaGeneMark_v1.mod进行基因预测,其范围是范围是细菌和古菌。其使用方法如下所示:

代码语言:javascript复制
gmhmmp [options] [sequence file]
其中sequence file为输入的contig序列,fasta格式。
必需参数:
-m  物种矩阵模型,可以是宏基因组,一般为自带的MetaGeneMark_v1.mod
选项参数:
-o  输出结果文件的文件名全名,默认为fasta文件名加.lst
-f  输出结果格式,可选L(LST)和G(GFF),默认为L
-a  输出预测基因的蛋白质序列(默认输出到总结果文件)
-d  输出预测基因的核酸序列(默认输出到总结果文件)
-A  预测基因蛋白质序列单独输出到文件的文件名
-D  预测基因核酸序列单独输出到文件的文件名
-k  输出RBS(核糖体结合位点Ribosomal Binding Site)得分和位置
-K  基因RBS得分和位置信息文件
-b  输出每条contigs序列的一个最佳预测
-r  确定是否使用核糖体绑定位点(RBS)模型来预测基因的开始。如果使用的话,必须给出一个包含RBS模型的文件
-s  预测基因的链,d为正向,r为反向互补链,默认为'.'也即正反向均预测
-p  允许基因之间有重叠,1为允许,0为禁止,默认为1,0用于预测没有内含子的真核生物基因组
-g  密码子表编号(详见2.1.1.2)
-M  宏基因组的模型文件
-R  RBS模型的文件文件名
-i  非编码状态开始和结束的可能性,默认为0.5
-n  不预测gap处不完整的基因

使用MetaGeneMark预测宏基因组基因,如下所示:

代码语言:javascript复制
cp …/mgm/gm_key ~/.gm_key  #拷贝key文件到自己主文件夹
nohup gmhmmp -m MetaGeneMark_v1.mod -A mgm.meta.protein.fa -D mgm.meta.gene.fa -o mgm.meta.gff -g 11 -f G new.spades.contig.fasta &

运行结束后,结果如下所示:

在基因组、宏基因组项目中,一般序列组装完成后的第一个步骤就是编码基因预测,这也是后续功能注释分析的基础。需要注意的是,真核生物基因结构与原核生物完全不同,其基因预测原理也不相同,通常我们使用原核生物基因预测工具预测宏基因组序列,获得的均为原核生物基因。

END

0 人点赞