长链非编码 RNA( Long non-coding RNA，Ln-cRNA) 是一类长度超过 200 nt 的长链非编码 RNA分子，是 RNA 聚合酶 II 转录的副产物，它可在多个层面上( 表观遗传调控、转录调控以及转录后调控等) 调控基因的表达［1-4］。据统计，哺乳动物蛋白编码基因占总 RNA 的 1%，而长链非编码 RNA 占总RNA 的比例可达4% ～ 9% ，这些长链非编码 RNA现已成为继 MicroRNA 后的研究热点［2，5］。目前发现 LncRNA 的序列保守性较低，相比编码蛋白基因具有较强的可塑性，进化速度快，且在基因组上的位置往往与其功能有一定的相关性［6-9］。

　　LncRNA 的作用机制比较复杂，目前发现的LncRNA的作用机制包括: 在蛋白编码基因上游启动子区发生转录，干扰下游基因的转录表达; 抑制RNA 聚合酶 II 活性、介导染色质重构或组蛋白修饰等，影响相关基因表达; 与蛋白编码基因的 mRNA结合，干扰 mRNA 的剪切成熟过程; 与蛋白编码基因的 mRNA 结合，在 Dicer 酶作用下产生内源性的siRNA，从而调控相关基因的表达; LncRNA 还可结合到特定蛋白质上调节相应蛋白的活性，通过结合到特定蛋白上，改变该蛋白的胞质定位［4，7-10］。总体而言，LncRNA 可通过表观遗传学调控、转录调控、转录后调控、蛋白活性调控等多种方式调控相关基因［4，7-10］。

　　目前，LncRNA 的相关研究大多集中在人类和哺乳动物中。截至 2013 年 7 月 1 日，NONCODE v3．0非编码数据库中已收录 73 327条 LncRNA，主要分布在人 类 ( 33 831，占 46. 14%) 和 小 鼠 ( 37 047，占50. 52% ) 中［11］。Cabili 等利用 RNA-seq 技术测定人体24 种细胞与组织中的 RNA，结合已知的 LncRNA数据，凭借序列、结构、转录、同源性等 30 个属性对LincRNA 进行大规模鉴定，得到 8 195 条 LincRNA( Long intergenic non-coding RNA)［12］。同时发现 Lin-cRNA 在不同的组织有着差异表达，LincRNA 还与其两旁的蛋白编码基因存在着共表达的关系［12］。Brunner 等完成了首个大型的癌症 LncRNA 表达谱分析，对64 个肿瘤样品高通量 RNA-seq 测序，在各种肿瘤类型之间找出差异表达的 1 065个LncRNA［13］。Li等使用 RNA-Seq 首次系统地验证了鸡骨骼肌发育过程中的 LncRNAs，证实了鸡基因组中 281 个新的基因组间 LncRNAs，发现和蛋白编码基因相比，新证实的LncRNAs 的保守性更差，但和随机非编码序列相比，这些 LncRNAs 较为保守［14］。Pauli 等对斑马鱼早期发育过程中的 8 个阶段进行时间序列的 RNA-seq 分析，重新构建了 28 912个位点上的 56 535个转录本，除了覆盖大多数已表达的 RefSeq 转录本外，还找到数千个新的亚型和表达位点，定义了一系列( 共1 133条) 在胚胎发育过程中表达的非编码多外显子转录本［15］。

　　较人类和哺乳动物中大量 LnRNA 的鉴定及其大量调控机制的功能解析相比，植物中 LncRNA 的研究则相对滞后［16］。目前，在新建立的植物 Ln-cRNA 数据库 PLncDB( Plant long noncoding RNA da-tabase) 中，利用不同技术手段 ( ESTs、Tiling array、RNA-Seq 等 ) 鉴 定 发 现 了 13 000 多 个 拟 南 芥LncRNA［17］。Swiezewski 等在拟南芥冷处理的种子中发现一个 LncRNA，命名为 COLDAIR，COLDAIR在春化过程中的对重要开花抑制因子 FLC 的沉默起着重要的作用［16］。此外，Boerner 等在玉米中利用全长 cDNA 序列( flcDNA) ，采用生物信息学方法鉴定了 1 011的 LncRNA，并对其进行了相关结构分析及基因组定位［18］。

　　随着高通量测序和生物信息学的快速发展，RNA 测序( RNA-seq) 成为研究基因表达、调控网络和转录组分析的重要技术手段［19］。目前，利用RNA-Seq 数据鉴定 LncRNAs 已逐渐成为代替传统微阵列、cDNA 测序的重要技术手段。然而，相较于人类和哺乳动物，植物因其功能基因组学发展的落后性，造成现有的基于哺乳动物发展的 LncRNA 鉴定方法在植物中产生较高的假阳性。本研究结合植物基因组特点以及高通量 RNA-Seq 测序的优势，拟设计优化一套基于 RNA-Seq 数据的植物 LncRNA鉴定方法。并进一步对其进行全基因组定位，为后续功能机制研究奠定基础。

　　1、材料方法

　　1． 1、玉米自交系 B73 顶端分生组织 RNA-Seq 数据

　　从 NCBI SRA 数据库(检索下载玉米 B73 顶端分生组织的 RNA-Seq 数据，编号分别为 SRR424649 和SRR424650。该数据是通过 Illumina 测序平台，长度为 76 bp 的单端测序获得，碱基数约计 2. 35 Gb。

　　1． 2、序列预清理

　　B73 顶端分生组织的 RNA-Seq 测序数据，首先通过 NCBI SRA 数据库中的 SRA Toolkit 软件转化产生碱基及质量值。随后，利用SolexaQA 软件包对原始测序数据的质量 ( Q20，Phred-Score≥20 即 1% 的错误率) 和测序长度( L20，长度≥20 bp) 进行过滤［20］。

　　1． 3、转录组重建

　　采用基因组指导的参考性组装，从而获得转录本集合。清理后的 B73 顶端分生组织 RNA-Seq 数据通过 Tophat v2． 0． 8［21-22］比对到参考基因组。随后，对那些没比对上的序列在重新比对之前清理到50 nt 长。利用 Cufflinks v2． 1． 1 程序［22］程序组装这些比对上的序列。在转录组组装时，用参考基因组序列信息与参考的基因注释集相结合的方法，减少因序列覆盖不足而造成的转录本不完整的情况，进而获得可靠的 B73 顶端分生组织转录本。

　　1． 4、Coding / non-coding 编码预测

　　对组装获得的 B73 顶端分生组织转录本进行编码能力区分及 LncRNA 鉴定与分类。主要是利用NCBI 非冗余蛋白数据库( nr) 和 Pfam 蛋白家族数据库，通过 CPC 程序［23］综合考虑 ORF 长度以及与已知蛋白的同源性来鉴定转录本是否编码，随后进一步过滤低表达转录本，从而获得候选的 B73 顶端分生组织 LncRNA 集合。

　　1． 5、基因组定位

　　获得候选 B73 顶端分生组织 LncRNA 数据集后，进一步通过自行编写的 Perl 脚本程序从基因组比对结果中释放基因组定位信息，从而获得 Ln-cRNA 在全基因的分布情况，进一步解析其类型和相关的基因信息。

　　1． 6、同源比对分析

　　为了 进 一 步 明 晰 基 于 RNA-Seq 数 据 预 测LncRNA的可行性及可靠性，将基于玉米 B73 顶端分生组织全长cDNA序列( flcDNA) 和 RNA-Seq 不同数据预测的 LncRNA 利用 Blastn［24］程序做了同源性分析。同源性分析采用参数为 E 值≤1e-15，最小匹配碱基≥100 bp，相似度≥80%。

　　2、结 果

　　2． 1、基于 RNA-Seq 数据预测植物 LncRNA 的流程设计

　　结合 RNA-Seq 数据分析方法以及植物基因组、转录组特点，优化了植物 LncRNA 的分析策略，主要分为 3 个部分: ①RNA-Seq 序列的预处理。主要利用 SolexaQA 对低质量碱基和长度进行过滤; ②转录组重建。目前的转录组重建方法主要有两类，一类是基因组指导的，另一类是不依赖于基因组的从头组装( de novo assembly) 。De novo 组装对测序深度、长度均有较大的依赖性，因此这里我们采用基因组指导的转录组重建( Genome-guided transcriptome re-construction) 方法，主要通过 Tophat2 程序利用剪切方式比对，采用外显子优先的比对算法，利用已知的基因集信息作为参考，将测序 Reads 映射到基因组上，进一步构建转录本; ③coding/non-coding 编码鉴定。通过 nr 和 Pfam 蛋白数据库过滤具有编码能力的转录本后，利用 CPC 程序综合考虑 ORF 长度以及与已知蛋白的同源性来鉴定转录本编码能力，过滤低表达转录本，获得 LncRNA 集合。

图 1 基于RNA-Seq数据预测植物LncRNA流程

　　2． 2、玉米 B73 顶端分生组织的 RNA-Seq 数据预测LncRNA

　　为了考察上述流程的可行性和可靠性，我们利用玉米自交系 B73 顶端分生组织的 RNA-Seq 数据来预测 LncRNA。选取的约 2. 35 Gb 测序碱基的B73 illumina RNA-Seq 数据( 登录号: SRR424649 和SRR424650) 清理后共 28 130 060条序列，碱基数约2. 1 Gb，均长为 72 bp。利用 Tophat 2 转录组比对程序将测序序列定位到玉米 B73 参考基因组( V2) 上。

　　28 130 060 条 测 序 Reads 序 列 中，18 297 964( 65. 11%) 条序列被特异地定位在 B73 参考基因组上，1 693 776( 6. 05%) 条序列定位在多个位点，而8 138 320( 28. 93% ) 条序列则未能定位到参考基因组。随后，利用 Cufflinks 转录组组装程序，采用基因组指导性组装( Genome-guided assembly) 策略，结合B73 的基因注释信息，组装共获得157 058 个转录本。通过 FPKM 表达量和长度过滤获得候选的43 865个转录本，用于后续的编码鉴定。

　　进一步通过 CPC 程序预测，综合考虑了转录本ORF 大小以及与 nr、Pfam 蛋白数据库的序列同源性，来鉴定 coding/non-coding 编码能力。结果发现，6 122( 14. 00% ) 转录本不具备编码能力，为非编码RNA( non-coding RNA) ，而 37 743( 86. 00% ) 转录本与已知的功能基因具有较高的同源性而被定义为编码 RNA( coding RNA) 。因此，基于 B73 RNA-Seq 数据共获得 6 122个长非编码转录本( LncRNA) ，平均长度为 619 bp( 表 1) 。同时，利用 Cufflinks 程序计算 6 122 LncRNA 的表达量( FPKM) ，结果显示 Ln-cRNA 总体表达量普遍较低( 图 2) 。

图 2 玉米B73顶端分生组织中候选 LncRNA FPKM表达量分布

表1 玉米B73 LncRNA预测结果汇总

　　2． 3、玉米 B73 LncRNA 的全基因组定位

　　利用 B73 参考基因组和 RNA-Seq 比对信息，对上述预测的 LncRNA 进行了全基因组定位，揭示其LncRNA 的在玉米染色体上的分布及邻近可能调控的基因信息。总体而言，6 122个 LncRNA 中，6 103( 99. 69%) 定 位 在 已 知 的 染 色 体 区 域，19 个( 0. 31%) 定位未知区域。其中，2 431个 LncRNA 定位在已注释的基因区域，而3 691个 LncRNA 则定位在已注释的基因区域外。

　　2． 4、基于 RNA-Seq 与 flcDNA 预测玉米 B73 Ln-cRNA 集合的同源分析

　　为了进一步考察 RNA-Seq 技术预测植物 Ln-cRNA 的可行性，将我们利用基于 RNA-Seq 预测的玉米 B73 的 LncRNA 集合与 Boerner 等利用 flcD-NA 预测的 LncRNA［18］做了比较分析( 图 3) 。两种途径发现的 LncRNA 通过同源比对分析发现，基于 flcDNA 预测的 B73 1011 个 LncRNA 中，101 个与基于 RNA-Seq 预测的 B73 LncRNA 序列存在较高的相似性。未能比对上的 910 条序列，应用本研究设计的方法进行编码能力鉴定，结果 Boerner等报道的 688 条序列［18］与目前报道的功能基因具有较高的同源性，因此应被鉴定为编码 RNA( cod-ing RNA) 。这一结果证实，基于 RNA-Seq 数据预测 LncRNA 相较传统低通量测序数据具有较高的准确性及可靠性。

图3 基于RNA-Seq和flcDNA不同数据预测的LncRNA集合比较

　　3、讨 论

　　利用生物信息学方法来鉴定筛选 LncRNA 是目前研究最为经济有效，同时也是主流的方法。在鉴定筛选 LncRNA 中，转录组重建和 Coding/non-cod-ing 识别是其中最为关键的两个部分。首先，转录组重建往往受序列读长、测序深度、文库大小等多个因素的影响。传统 cDNA Sanger 测序读长较长有助于组装，但是其通量低、测序错误率高以及组装常常产生嵌合体，从而会造成后续分析的误差［25］。RNA-Seq 读长短，但其通量高、测序质量较高、基因覆盖度高，有助于获得全面的转录组信息，揭示转录本UTR 长度、外显子边界定位、选择性剪切机制以及动态的表达调控机制。此外，双末端 Paired-end 测序方式，同时结合参考基因组信息，将极大有助于转录组的准确重建［25-26］。本研究采用的参考基因组指导性组装，利用已知基因组注释信息，采用外显子优先算法有效保证了转录组的准确重建。另外，coding / noncoding 编码能力鉴定则是关系到 LncRNA能否发现的重要步骤。本研究采用综合性方法，利用 CPC 程序综合考虑 ORF 长度以及与已知蛋白的同源性，进而鉴定其编码能力，这一方法已广泛应用在人类、哺乳动物以及一些模式植物如拟南芥、玉米的 LncRNA 鉴定中［18，27-28］。但该鉴定方法不可避免地受蛋白数据库大小的影响，蛋白数据库越大，收录的蛋白序列也越多，鉴定的 LncRNA 数目则越少。

　　目前，Swissport 数据库和 NCBI RefSeq 数据库是通过文献支持以及专家审阅的高质量非冗余蛋白数据库，收录的蛋白序列虽然少，但较为可靠。而 NCBI的 nr 蛋白数据库则是目前收录蛋白序列最多的数据库，但其中包含大量机器翻译推测的蛋白序列，其准确性相对较低。因此，选用蛋白数据库需要根据试验自身的需求来调整。

　　此外，植物基因组因其自身结构特点以及基因组学发展的滞后性，造成了不能简单套用人类和哺乳动物 LncRNA 的鉴定方法。人类和哺乳动物 Ln-cRNA 鉴定的流程大致可分 3 个步骤: ①转录组重建; ②过滤已知基因集，选取非基因集或者与已知基因集存在变异的转录本; ③编码能力鉴定。然而，植物基因组学发展相对滞后，其注释基因集含有大量机器 de novo 从头预测的基因。因此，植物中的参考基因集并不能作为准确可靠的功能基因去过滤，而是需要重新评价或者仅作为参考基因定位数据。鉴于此，本研究充分结合了植物基因组特点以及高通量 RNA-Seq 测序数据的优势，通过生物信息学手段发展、优化了一套基于 RNA-Seq 数据的植物 Ln-cRNA 筛选及鉴定方法，为植物中鉴定筛选 LncRNA提供重要的思路和方法。

TAG标签： RNA-Seq     LncRNA     植物