原核生物基因组肽编码sORFs分布及功能特征

2016年，德克萨斯大学Eric教授团队在《科学》(Science)杂志报道了一种由位于lncRNA上的小开放阅读框(small open reading frame，sORF)编码的小肽DWORF[1]，长度仅为34个氨基酸，在心肌收缩中发挥重要作用.DWORF的发现引起了人们对sORFs及其编码多肽(sORFs encoded peptides，SEPs)这种长期被忽略的基因组元件的空前关注，激发了人们对非编码RNA的激烈争论，成为当前研究的重点领域[2].一方面，能否编码蛋白质一直是区分非编码RNA与mRNA的金标准，而DWORF的发现引起人们对非编码RNA是否真的不编码蛋白质的讨论[3-5]；另一方面，SEPs的发现打破了人们长期以来认为生物活性肽多是通过蛋白质前体修饰剪切而来的推论[6].几十年来，sORFs被认为不可能具备蛋白质编码能力，往往将它们归于错误预测为蛋白质编码基因的随机序列[7]，在之前数据库及蛋白质编码基因预测算法中通常将长度低于100个氨基酸的sORFs排除在外，进而降低基因预测结果的假阳性，但这样反而使得许多真正编码蛋白的sORFs被遗漏，进一步增加了预测结果的假阴性[6].因此，肽编码sORFs的发现具有重要的学术意义，尤其在以“精准”为标签的组学测序快速发展的现在，肽编码sORFs的发现似乎表明人们在某种程度上对基因组的认识转了一圈又回到了原点.目前，已有大量研究证实了肽编码sORFs在生命活动中发挥重要生物学功能[8-9]，但由于其表达水平及丰度低、序列短、实验技术缺乏等诸多因素，对sORFs的研究还处于初级阶段[6, 10-11]，能够有效识别sORFs的生物信息算法及相关数据库资源也很缺乏[12-13]，其相关序列结构等生物特征认识亟待深入[11, 14]，多数研究集中在人、鼠、拟南芥等几种真核模式生物[6, 15]，对原核生物研究还较少[16-17].因此，本文借助已有数据资源，首次针对原核生物基因组中肽编码sORF序列组成及功能分布特征开展了系统研究，并针对当前sORFs研究存在的问题开展了深入探讨，为今后原核生物sORFs研究提供了重要理论依据.

目前仅有少数几个针对模式真核生物或特定区域的sORFs数据库发表[12-13, 18]，专门针对原核生物基因组的sORFs相关数据库资源鲜有报道.尽管如此，基于转录组、蛋白质组等组学测序与生物信息技术，已有一些原核生物中的肽编码sORFs被发现并注释出来[16-17]，为本工作提供了良好的参考资源.为了保障数据可靠性，本文所需肽编码sORFs信息(即CDS序列)从RefSeq数据库[19]获取.

基于RefSeq数据库，随机选取了80余种具有不同基因组G+C含量的原核生物作为研究对象.为了进一步研究肽编码sORFs的生物学特征，我们根据RefSeq中的注释信息，将sORFs分为3个长度区间，即50~100个氨基酸、30~49个氨基酸、29个氨基酸及以下.

为了展现sORFs编码多肽序列中氨基酸使用偏好特征，我们利用序列复杂度[20]概念来定量描述：

式中，i=1，2，3，…，20，表示氨基酸种类，fi表示序列中第i种氨基酸的频率.根据统计原理，当20种氨基酸平均使用时，K值最大，为4.32；而当序列中只有1种氨基酸时，K值最小，等于0，此时说明序列中氨基酸使用偏好最强.因此，K值能够很好地反映序列中氨基酸使用信息.本文中，将利用复杂度K值来进一步分析各原核基因组中肽编码sORFs序列组成特征.

氨基酸组成是蛋白质及多肽的生物学功能基础，直接计算各种氨基酸百分含量fi是研究蛋白质序列的重要方法.为了有效描述蛋白质序列中各种氨基酸的使用特征，可直接通过计算Ci=fi-0.05来完成.由于20种氨基酸的随机使用概率为1/20=0.05，因此本文将各种氨基酸的使用偏好用上式来简单描述，若Ci大于零，则表示相应氨基酸偏好使用，反之不偏好.

为了描述氨基酸在蛋白质序列中的分布特征，本文利用我们提出的一种蛋白质序列分析圆柱体模型[21]，通过计算di=ni /N来描述，这里ni表示任意蛋白质序列中第i种氨基酸的位置，N表示序列长度，因此di∈[0, 1]，直接计算

Ni表示序列中第i种氨基酸个数，因此Di=[D1，D2，D3，…，D20]定量描述了20种氨基酸在蛋白质序列中的分布中心.

基于RefSeq中的注释信息，附件表S1给出了77种不同G+C含量的原核生物基因组中肽编码sORFs分布信息，其中百分比表示各基因组中肽编码sORFs占该基因组所有蛋白质编码基因的比例.可以看到，各个基因组中均有小于100个氨基酸的sORFs被不同程度注释出来.从长度分布来看，已经有注释信息的sORFs还多集中在50个氨基酸以上，其次是30~49个氨基酸之间的sORFs，而且有些模式基因组中甚至已注释出了长度在29个氨基酸以下的sORFs.例如，在大肠杆菌E. coli str. K-12(NC_000913)中，有35个长度小于29个氨基酸的sORFs被注释以来，其中最短的sORFs仅有14个氨基酸.可见，随着基因组信息研究的深入，sORFs不再是之前认为的不具备蛋白编码能力的随机序列，sORFs也是普遍存在的基因组元件.由于目前对于原核生物基因组中肽编码sORFs的相关研究很少，为了揭示原核生物中sORFs的相关特征规律，我们进一步以6种研究相对较为广泛、深入的大肠杆菌基因组为例开展了序列分析.从表 1给出的6种大肠杆菌菌株sORFs的注释情况来看，每个大肠杆菌菌株中都有长度小于29个氨基酸的sORFs.根据表 1，计算了各种长度分布sORFs对应G+C含量与其编码的多肽序列复杂度之间的散点图，结果见图 1.可以发现，不同长度肽编码sORFs相应的G+C含量分布差别不大，但与基因组G+C含量(~50%)相比，分布比较离散.而对序列复杂度而言，不同长度区间sORFs分布区域具有明显差异，50~100个氨基酸sORFs序列复杂度最大，其次是30~49个氨基酸sORFs，0~29个氨基酸sORFs序列复杂度整体最低，3个长度区间序列复杂度平均值及标准偏差见表 2，各样本间方差分析P值均 < < 0.05，表明相应序列中氨基酸使用具有不同程度偏好特征.

从信息熵的数理统计原理很容易理解图 1中不同长度分布肽编码sORFs展现出的序列复杂度差异，而从生物学角度，则需要进一步揭示其背后的序列特征.为此，接下来本文又深入分析了各长度分布肽编码sORFs中氨基酸组成及分布特征，结果见图 2.为了便于观察和对比分析，图 2横坐标中对各氨基酸按照其理化特征依次排列，即非极性氨基酸(A、V、L、I、F、W、M、P)、极中性氨基酸(G、S、T、C、Y、N、Q)、极正电氨基酸(H、K、R)和极负电氨基酸(D、E).由图 2a可以看出，不同长度sORFs在氨基酸使用偏好上具有一定程度的差别，就本文选择的6种大肠杆菌而言，3种长度分布sORFs在非极性氨基酸使用上偏好性较为明显，且长度小于29个氨基酸的肽编码sORFs氨基酸使用频率浮动范围较广，展现的更明显的氨基酸使用偏好特征.而由图 2b，3种长度分布sORFs中氨基酸分布差别相对较大.

通常情况下，序列长度的降低增加了分子生物学实验中PCR引物设计的难度，实验研究sORFs难度也相应增加，因此对其功能的研究是目前生命科学的一项重要课题.在附件表S2中，我们整理了不同长度分布肽编码sORFs的具体功能信息，从目前注释结果来看，3种长度区间sORFs分别有35%(6337/18118)、25%(380/1498)和53%(58/110)具有明确功能注释.从分布上来看，长度越大，涉及到的功能类型也越多.这里，我们依然以大肠杆菌中长度在29个氨基酸以下的sORFs为例，来进一步研究其功能特征.表 3中给出了来自表 1中所列6种大肠杆菌基因组的79个肽编码sORFs(< 30氨基酸)功能分布信息，除掉没有功能注释的30个sORFs，其余49个sORFs的生物功能可以大致分为6大类，其中近50%为前导肽(leader peptide).

表 3中的sORFs呈现出相对集中的生物学功能特征，在对真核生物学研究中发现sORFs编码肽具有一定程度的序列保守性[9]，接下来我们分别以具有leu operon leader peptide、thr operon leader peptide、tryptophanase leader peptide及toxic membrane protein功能的4组sORFs编码肽为例进一步分析其序列保守特征.借助多序列比对程序Clustal，可以得到图 3中的序列分析结果.可以看到，尽管具有前导肽功能的序列来自不同大肠杆菌菌株，但展现出高度的序列保守性，其中leu operon leader peptide和tryptophanase leader peptide对应的sORFs在不同菌株中完全一致.同样，在具有toxic membrane protein功能的sORFs序列中也展现出一定程度的序列保守特征.在图 3e中，我们对上述4种不同功能的sORFs编码肽进行了序列比对，其中每条序列前标记的数字1~4分别表示图 3a~d的功能类别，可以看出在这些序列中缺少保守序列.因此，图 3的结果表明尽管sORFs序列较短，在依然具备了同蛋白质编码基因类似的保守性序列特征.

截止到2017年2月，GenBank收录已完成基因组计划70 036个[22]，其中超过90%为细菌(62 720个)和古细菌(713个).作为与日常生活及工农业生产最为密切的原核生物，一直都是科学研究及生物技术应用的重点领域.本文从全基因组角度证明了原核生物基因组同样普遍存在具有重要功能的肽编码sORFs，而且数量分布也非常广泛，这可为今后基于多肽的抑制剂、药物设计等提供丰富的资源.同时，肽编码sORFs的普遍存在也进一步说明之前基因预测研究中将长度小于100个氨基酸的sORFs排除掉缺少理论依据，反而加剧了原核生物蛋白质编码基因错误注释的进一步积累.因此，原核生物肽编码sORFs还有待于系统研究、开发，这也是当前测序技术快速发展背景下提升基因组注释质量急需解决的重要课题.

开放阅读框本是分子生物学中一个非常基础概念，其定义是从起始密码子到终止密码子结束的一段碱基序列.虽然定义很简单，但从已有参考书中却很难找到对ORF的深入解释，在一些专业论坛及学术资料中，许多研究人员甚至认为ORF就是CDS序列.而实际上，只有很少数ORF才具备蛋白质编码能力，这也是利用ORFfinder等预测程序得到的ORF数量要远多于基因预测程序得到的蛋白质编码基因数量的主要原因.因此，在非编码RNA中存在ORF序列也已是人们习以为常的现象，但很少有人想到过这些ORF能够编码蛋白质，因而非编码序列中肽编码sORFs的发现引起人们很大兴趣[2, 23].在首个公开发表的sORFs数据库sORF.org中[12]，收集了人、鼠和果蝇等几种真核生物中共计348 844条由核糖体谱识别得到的肽编码sORFs.其中，仅人类基因组中收集的新sORFs就达到了190 195个，而一直以来，人们对人类基因组蛋白质编码基因数目估计仅为20 000~25 000个.尽管原核基因组基因预测工作已持续了30余年，近年来大量工作表明，目前原核生物基因组蛋白质编码基因普遍存在错误注释[24]，已有的基因预测程序还有很大发展空间，而对具有肽编码能力的sORFs的有效预测算法还有很多难点.因此，刚刚处于起步阶段的肽编码sORFs研究带来的挑战及面临的系列问题还需要进一步系统、深入研究.

继基因组测序之后，转录组、蛋白质组(质谱分析)等测序手段已成为当前生物医学领域研究的标配技术[25-26].基因组测序为人们打开了基因组这本天书，而转录组、蛋白质组等测序技术为解读这本天书提供了有力工具.长期以来，转录组与蛋白质组在mRNA与非编码RNA研究中发挥了极为重要作用[16]，因而人们对借助这些组学测序技术来研究sORFs一度寄予厚望[26].然而，最终测序结果表明距离预期存在很大差距[14].这其中的主要原因就在于sORFs整体表达水平低、丰度低，有些仅在某种条件下选择性表达[27].借助转录组测序，人们无法进一步判断转录产物是停留在RNA水平还是继续被翻译，而以质谱分析为代表的蛋白质组测序虽然能够检测到mRNA翻译产物，但对于低丰度、低表达、不表达的sORFs没有效果，因而实际应用有限[16].近几年发展的核糖体谱(ribosomeprofiling)技术被认为是继质谱分析之后检测RNA翻译与否的更有效测序技术[28-29]，并被广泛应用于肽编码sORFs识别[12, 30-31]，其原理是通过分析结合在核糖体的分子来判断mRNA的翻译情况，但近期发现许多lncRNA等非编码分子也会与核糖体结合[4, 27]，因而依然存在许多问题.因此，肽编码sORFs的发现为当前各种组学测序技术带来了更大的挑战，将测序技术与生物信息技术有效结合成为今后sORFs研究的必然趋势[2, 10, 23, 32].

缺少有效的高通量研究方法是目前肽编码sORFs遇到的首要问题[33].通过转录组测序可以获得大量的RNA转录产物，但无法判断这些转录产物是否能够翻译.尽管某些分子生物学和生物化学方法(如实时荧光定量PCR等)也可以不同程度用于sORFs研究，但也仅能停留在转录层次[27]，也无法实现高通量实验研究.因此，计算方法研究sORFs同样成为当前的主要研究手段[26-27].目前对sORFs计算研究包括两个层次，一个是通过分析sORFs序列组成[23, 34]，包括碱基组成、密码子和二联密码子(六聚体碱基)[35]等；另一种是借助BLAST等相似性分析程序来获取sORFs特征[11, 14, 36].然而，由于sORFs序列相对较短，有些传统的序列分析方法可能无法直接应用，但一些应用于其他类型短肽或RNA的计算方法可以为sORFs研究提供有效参考[37].另一方面，通过进一步发展序列结构分析方法来深入揭示肽编码sORFs的固有特征[38]，与实验研究sORFs有效互补[10]，有望为肽编码sORFs研究提供新思路.

总之，作为一种长期被忽略的基因组“新元件”，从其在基因组中的识别定位到功能研究，肽编码sORFs承载了极为重要的学术意义，也为多肽应用提供了丰富的生物资源[39].面对肽编码sORFs带来的机遇与挑战，充分融合生物信息技术、各种组学测序技术及生物实验技术将是今后切实可行的研究趋势.因此，希望本文能够为今后原核生物基因组生物活性多肽资源挖掘及应用提供可靠的理论支持.

附件  表S1和表S2见本文网络版附录(http:// www.pibb.ac.cn).