论文-小RNA高通量测序数据分析方法.pdf

- 1 - 小小 rnarna 高通量测序数据分析方法高通量测序数据分析方法 彭 骅 , 王胜华 , 陈 放 四川大学生命科学学院，成都（610064） e-mail： 摘要摘要 本文应用 perl 语言和 mysql 数据库构建了小 rna 高通量测序数据分析平台， 以 5 个水稻数 据集为分析对象，详细介绍了小 rna 高通量测序数据的处理方法和流程. 我们以 msu 6. 1 水稻基 因组为参考，构建了该版本的全基因组结构及已知 ncrnas 位点信息数据库，结合 perl 脚本可以实 现小 rna 在基因组上的详细定位与统计，同时我们从数据库中提取已知 pre-mirnas 表达特征，设 计了一个新的 mirnas 挖掘方法，该方法可以筛选出大量的新 mirnas，其中已知 mirnas 命中率 可以达到 98%. 针对水稻小 rna 种类的多样性， 我们对 mirnas 和 endo-sirnas 的鉴别也给予了探 讨和说明. 本文设计的高通量测序数据分析平台，方法简单高效，以数据库作为存储和查询媒介， 能够实现多位点 reads 的分析，可以得到灵活多样的数据统计结果. 依照本文的方法同样可以构建其 他模式物种的小 rna 数据分析平台， 在高通量测序逐渐普及的将来， 本文的方法对中小实验室建立 自己的数据分析平台具有实践指导意义. 关键词关键词 小 rna(small rnas)， 非编码 rna(ncrnas)， 微小 rna(micrornas，mirnas) 内源小干扰 rna(endo-sirnas)， 小 rna 高通量测序(small rna-seq) 中图分类号：中图分类号： q531 q751 tp391 1 引言引言 小 rna(small rnas)主要指长度在 18-30nt 的一类非编码 rna(ncrnas)，在真核生物中，具有 基因表达调控功能的小rna主要有微小rna(micrornas， mirnas)、 内源小干扰rna(endo-sirnas) 和 piwi 干扰 rna(pirnas). pirna 长度集中在 26-31nt， 目前只在动物的生殖系细胞及干细胞中被发 现，其主要功能是参与转座子的沉默1,2. mirnas 和 endo-sirnas 长度主要集中在 20-24nt. mirnas 在动植物和微生物中都普遍存在， 目 前在 mirbase 14 数据库中已包含 115 个物种的 12627 条记录3. 在细胞质中，mirnas 与 ago1 等 蛋白形成 risc 复合体(rna-induced silencing complex), risc 通过 mirna 与特定的 mrna 靶基因 互补配对，在配对区域的中间位置，ago1 通过对 mrna 的切割促使其降解或者通过翻译抑制实现 转录后调控4-6. 据估计一个物种中约1/3的基因会受到mirna的调控7， 大量的实验也表明mirnas 参与了诸多生命过程的调控，例如细胞周期，细胞分化，组织器官的发生，营养代谢，信号途径以 及对外界生物的非生物的环境的反应9-12；同时，mirnas 在生产实践与临床治疗上也具有很大的应 用前景13,14. 小分子干扰 rna(sirnas)最初在植物转录后基因沉默现象中被报道,长度在 20-25nt,来 自外源的双链核糖核酸(dsrna)切割产生15. 随着小 rna 研究的深入发展，大量的 endo-sirnas 被 发现，在植物体内，endo-sirnas 目前可以分为三种类型：(1)trans-acting sirnas(ta-sirnas)，21nt， 功能与 mirnas 类似，与 ago1 或 ago7 组成 risc 复合体参与转录后调控. (2) natural-antisense sirnas(nat-sirnas)，21nt 和 24nt, 21nt 的 nat-sirnas 功能与 mirnas 类似，参与转录后调控. (3)repeat-associated sirnas(ra-sirnas)， 24nt， 参与染色体水平的基因沉默、 抑制转座子的转座16-18. 显然, 生命体内还存在许多具有重要功能的small rnas. 如何鉴定与发现这些 rna是值得人们思索 中国中国科技论文在线科技论文在线 - 2 - 的重要问题. 以往用于寻找 mirnas 等小 rna 的方法有实验克隆法， 计算机预测法19-23. 克隆法可以直接用 于鉴定新小 rna，是初期发掘小 rna 的常用方法，不足之处是实验周期较长，对低表达的小 rna 的发现能力十分有限. 计算机预测法多是针对某一已知的小 rna 特征设计算法， 从全基因组或 est 数据库中快速发掘大量潜在的小 rna，一定程度上弥补了克隆法的缺点，然而，预测的小 rna 最 终还需要实验证明，同时计算机预测法对新类型小 rna 的发掘能力十分有限. 随着第二代高通量测 序技术的问世，以测序为中心的功能基因组学研究开始全面展开24,25，其中的小 rna 高通量测序 (small rna-seq)技术开始逐渐取代原始的小 rna 发掘法方法，该法具有速度快、成本低、覆盖度深 等多方面的优点， 对鉴定与发现生命体内的小分子rna及其功能与机理研究起极大的推动作用26-28. 从最初应用 454 技术建立的谷类小 rna 数据库(csrdb，cereal small rnas database) 29，到应用 illumina/solexs 技术建立的拟南芥小 rna 计划30 (asrp，arabidopsis small rna project)，small rna-seq 已被广泛应用于特定发育阶段全基因组水平上的小 rna 的鉴定与发现31-34. 然而， 在真核 生物内小 rna 种类繁多，高通量测序产生的数据量巨大，例如，一个水稻小 rna 样本的 illumina 测序数据可达 2g 大小，含有约 1200 万个测序片段(reads)读数，如何有效快速的处理这些数据，深 入挖掘未知小 rna，是小 rna 高通量测序数据分析的一个主要问题35. 目前小 rna 高通量测序数据的挖掘能力还十分有限，已发表文献中所挖掘出的小 rna 也只是 数据中很小的一部分，而用于数据挖掘的综合软件也不多见，多是针对某些特殊需求设计的单一功 能软件. 目前的开源软件主要有 learn36、mirdeep37、cashx34等，learn 是基于 perl 语言的 已知 ncrnas 注释程序，mirdeep 是基于 perl 语言的 mirnas 挖掘程序，cashx 是拟南芥小 rna 计划组开发的，是基于 perl 和 mysql 的存储程序，它的功能和数据库结构比较简单. 而目前少数几 种可用于高通量数据分析的商业软件价格相对昂贵，中小实验室难以负担，并且单纯依赖商业软件 的分析结果往往又不能满足客户的各种要求。本文从实际应用出发，为科研人员详细介绍了高通量 数据的分析方法. 我们采用 perl 语言和 mysql 数据库构建了小 rna 高通量测序数据分析平台. 该 平台设计简单科学， 效率和功能强大， 能够被大多数的中小实验所采用； 可以实现多位点 reads 分析， 而这部分数据在以往的分析中常常是被过滤掉的；可以实现数据的综合分析，得到灵活多样的统计 结果；同时数据库本身具有很强的扩展能力，可以为各种小 rna 及未知小 rna 的挖掘提供支持， 利用该平台我们也成功地挖掘出了一批新 mirnas. 2 2 材料和方法材料和方法 2.12.1 实验平台及软件实验平台及软件 本实验平台配置：ibm p615 小型服务器，2 核 power 1. 45g 处理器，8g 内存，fedora 10 ppc 64 位操作系统，mysql 5. 1 数据库，用于数据库存储与查询；hp dx2390 台式机，2 核 inter 2. 8g 处 理器，4g 内存，fedora 10 x86 64 位操作系统，perl 1. 6，主要用于数据处理. 使用的软件有：perl 脚本，用于多个步骤中的数据处理；perl-dbi 模块主要用于 mysql 数据表 之间查询与更新；phpmyadmin 用于 mysql 数据库管理和可视化查询；soap38软件用于映射 tags 到参考基因组；rfam/infernal 用于已知 ncrna 收集；dchip39用于芯片数据分析；rnahybrid40 用于 mirnas 靶基因筛选；rnalfold41及 mipred41 用于 mirnas 二级结构分析. 2.22.2 小小 rnarna 高通量数据高通量数据 .1 实验数据集实验数据集 本实验使用的数据集(表 1)均为 illumina 平台的水稻小 rna 测序数据，sample1-3 来源于 ncbi 的 geo 数据库， sample4 个来源我们实验室的测序结果， 数据集在选择上并没有实验设计上的考虑， 中国中国科技论文在线科技论文在线 - 3 - 原因是能够获得到的原始数据集的数量非常有限，截止到目前，geo 数据库中仅有 13 个水稻高通 量测序数据集. table 1 small rna-seq datasets table 1 small rna-seq datasets sample gsm subspecific source raw format data size submission date sample1 gsm278571 japonica group grains/1-5daf* scarf 516m 2008-8-21 sample2 gsm278572 japonica group grains/6-10daf scarf 517m 2008-8-21 sample3 gsm361264 indica group flag leafs/120days sanger-fastq 1. 5g 2009-7-15 sample4 - indica group pollen/trinuclear illumina-fastq 2. 1g *daf,days after fertilization. 2.2.2 原始数据格式原始数据格式 实验样品涉及的三种原始数据格式 （表 1） ， 其均可以由 illumina 测序平台生成， 不同之处如下： scarf 格式(图 1a)： 每个 read 用单行表示， 序列前面的部分表示实验编号和 read 所处的阵列位置； 序列 53排列，测序结果直接对应待测定的样本 rna 或 dna 序列，低质量的核苷酸用 n 表示； 序列后边数字(-540)表示对应的序列质量，数值越小表示质量越不可靠，由于 sbs 测序法技术上 的限制，随着测序长度的增加测序结果就越不可靠，所以 illumina 平台的 read 长度一般控制在 40nt 以下. scarf 格式是 fastq 的压缩格式，数据较小便于传输，序列质量便与观察. fastq 格式：每个 read 都用四行表示(图 1b)，第 1 行和第 3 行为序列标识，含义同上；第 2 行为 reads 序列；第四行表示序列质量，sanger-fastq 序列质量采用 33 位至 126 位的 ascii 码编码， illumina-fastq 序列质量采用 59 位至 126 位的 ascii 码编码. 2.2.3 原始数据特点原始数据特点 （1） 在一个实验的数据集中所有 reads 长度均相同，长度一般在 40nt 以下. （2） reads 3端包含接头序列片段，长度不定. 例如，实验样本分离 18-28nt 范围的小 rna， 测序长度为 33nt，则每个 read 的 3端有 5-15nt 长度不等的接头序列. （3） 相同 read 的数量可以直接代表该序列的表达量. （4） reads 数据量大，一般一个 small rna-seq 数据集包含的 reads 数不低于 500 万个. 2.2.4 高通量数据处理中的相关概念高通量数据处理中的相关概念 由于 reads 的数量巨大，去除接头后的，我们可以把序列相同的 reads 合并，称为一个序列标签 (tags)，tags 映射(mapped)到基因组后，位点相互重叠的 tags 我们可以聚集成簇(cluster)，在一定距离 内相互邻近的 cluster 我们可以进一步归为一个假设转录单元(unit)(图 3). read 映射到基因组的位点我们称为 hits，大多数 reads 的位点是唯一的(即 hits=1)，这些 reads 在后续分析中比较容易处理，而 hits1 的 reads 在以往的分析中往往被过滤掉，实际上很多 mirnas 在基因组中是多位点的， 其它类型 ncrnas 中多位点现象更是普遍， 目前多位点 reads 的定位和处理 也是高通量测序中急需解决的一个基本问题， 本文根据表达有效性对多位点 reads 进行过滤， 根据上 下游表达环境进行位点定位，可以实现大部分多位点 reads 的处理. (a) scarf (b) fastq fig.1 small rna-seq raw data format 中国中国科技论文在线科技论文在线 - 4 - 2.32.3 参考基因组数据库参考基因组数据库 水稻有两个亚种： japonica 和 indica， 目前的基因组注释版本有 japonica： irgsp build5(2008-6-26 更新)，msu v6. 1(2009-6-3)；indica：bgi-9311(2006-12-28). 由于 msu v6. 1 注释较为完整，在本文 作为参考基因组43 . msu v6. 1 共有 56797 个基因位点，其中 32%的基因位点有 go 注释，28. 6%的 基因位点与转座元件(tes)相关. 我们根据 msu v6. 1 的基因组注释信息，从新划分了基因组结构(图 3)，构建了用于高通量数据分析的 rice 数据库(图 4). 2.3.1 基因组结构信息表基因组结构信息表的的构建构建 如果将 dna单链分别对待， 例如以“+”链为例， 可将单链分为基因间隔区(图 3a)和基因编码区(图 3b)，在某些基因间隔区中还包含“-”链基因的反义链区(图 3c)；如果把 dna 双链为整体对待，可将 dna 分为基因位点间隔区(图 3d)和基因位点区(图 3e). 在 msu v6. 1 的注释文件中，我们根据已知的(b)、(d)位点信息构建了(a)的位点信息，分别导入 rice(23) genes genes_cross intergenic fb_intergenic feature mir14 ncrna ncrna_typelist mule organelle repeat fst gff_typelist sequence map_genes map_fbintergenic map_fst map_mir14 map_mir14_ant map_mule map_ncrna map_organelle map_repeat table 2 fb_intergenic table structure fields comment id record id, primary key seqid chromosome number start start sites end end sites strand +/- strand length intergenicss length f_gene_id adjacent front genes id, foreign key b_gene_id adjacent back genes id, foreign key fig.4 rice genome database structure - - + + (a) (b) (c) (d) (e) fig.3 illustration:genomic structure (a)intergenics on ssdna. (b) genes on ssdna. (c)gene_antisense on ssdna. (d)intergenics on dsdna. (e) genes on dsdna (a) (e) read1 read2 read3 read4 tag1 tag3 tag2 tag4 cluster1 cluster3 unit1 cluster2 fig.2 illustration: read, tag, unit 中国中国科技论文在线科技论文在线 - 5 - 到本地 rice 数据库的 genes 表、 intergenic 表和 fb_intergenic 中表. 由于某些基因存在重叠情况(图 3)， 我们创建了重叠基因信息表 genes_cross 表. 对于基因内部的结构信息(内含子、外显子、3、5端)， 将其存储在 feature 表中. 其中 genes 表和 fb_intergenic 表(表 2)是基因组数据库的核心表，与核心表 对应的 map_genes 表和 map_fbintergenic 表用于存储不同样本中 tags 所映射的 genes 及 fb_intergenic 的表达情况. 2.3.2 已知已知 mirnas、ncrnas 和和 marker 信息信息表表的的构建构建 如果基因组的注释比较全面，则有关的注释信息可以直接从注释文件中获取，如果注释不够全面 或者达不到实验目的要求，则需要自己收集相应的序列信息再通过 blastn 映射到基因组，本实验建 立了 6 个相关注释信息表，如(图 4)中所示，包括 mir14 表、ncrna 表、repeat 表、fst 表、mule 表、organelle 表，其中将 repeat、fst、mule、organelle 归类为基因标识(marker)，三个表的信息来 源及说明见(表 3). table 3 mirna、ncrna and maker tables description classes type label total rows source and description mirnas pre-mirna os01 684 414 sequences come from mirbase_14,blastn mapped mirna os02 731 451 sequences come from mirbase_14, blastn mapped ncrnas snorna os03 843 rfam_8. 1 and plant snorna db, blastn mapped snrna os04 196 rfam_8. 1 and noncode, blastn mapped srrna os05 40 rfam_8. 1and irgsp5, blastn mapped lrrna os06 15 rfam_8. 1and irgsp5, blastn mapped rrna_region os07 7 irgsp5, blastn mapped srp os08 14 srpdb, blastn mapped trna os09 851 rfam_8. 1 and irgsp5, blastn mapped trnascan(1) os10 94 msu v6. 1, perl script extracted intron os11 240 rfam_8. 1，blastn mapped other os12 16 rfam_8. 1and noncode，blastn mapped marker repeat(2) gf01 264562 all come from msu v6. 1 annotation, perl script extracted fst(3) gf02 58032 mule(4) gf03 13365 organelle(5) gf04 165 (1)trnascan, identified using the trnascan-se (2)repeat, come from msu oryza repeat database v3. 1 (3)fst, flanking sequence tags. (4)mule, mutator like elements. (5)organelle, organellar insertions repeat 表存储水稻重复序列位点信息，来自 msu oryza repeat database v3. 1；fst 表存储侧翼 序列标签(fst,flanking sequence tags)，来自 tos17， t-dna 或 ac/ds 的插入突变位点侧翼序列； mule 表存储基因突变元件(mule，mutator like elements)，来自突变基因或基因片段；organelle 表 存储细胞器插入片段(organellar insertions)，来自水稻叶绿体和线粒体在基因组上的同源序列. 我们 选择这些序列作为 marker 不但可以对已知 mirnas 的做起源分析， 还可以用于未知的 tags 的深入分 析. 与 6 个注释信息表相对应的 6 个 map_*表(图 4)，用于存储不同样本在每个注释位点上表达差异 信息. 如 map_mir14 表中， 每个已知的 mirnas 分别存在 4 个 tc 字段、 4 个 rc 字段、 4 个 h1rc 字段， 分别代表 4 个样本中的 tags 数量、reads 数量、位点唯一的 reads 的数量. 2.42.4 样品样品 tagstags 数据数据库库 本实验建立了四个样品数据库，每个数据库用于存储各自样品中的 tags 信息和 tags 在基因组上 的映射信息，数据库结构如(图 5)所示. 结合(图 3)所示的基因结构对 tags 映射的基因组类型进行分 类并存储在不同的映射表中（表 4). 中国中国科技论文在线科技论文在线 - 6 - 样品数据库的核心表为 tags_map 表、 tags_to_cluster 表和 cluste_to_unit 表. tags_map 表中核心字 段为 r_count，hits，map_label，map_ncrna，gff_label(表 5)，其字段说明如. table 5 tags_map table structure fields comment id record id, primary key tagsid_idg storage tagsid table id, foreign key seq_type storage tags sequence and mismatches type seqid chromosome number start start sites end end sites strand +/- strand hits storage tags alignment genome sites number r_count the tags contains reads number map_label storage tags sites information of genome structure，which come from sample. table (3)-(8) ncrna_label storage tags sites information of ncrna, which come from table rice. map_ncrna gff_label storage tags sites information of marker, which come from marker tables on rice database tags_to_cluster 表中核心字段为 tc，rc，h1rc，map_label，ncrna_lable，gff_label，all_tagsid_rc. 其中 tc 存储 cluster 在一起的 tags 计数；rc 是对所有 reads 计数；h1rc 是对 hits=1 的 reads 计数； map_label，ncrna_lable，gff_label 来自 tags_map 的 map_lalel，ncrna_lable，gff_label 的推导； all_tagsid_rc 存储聚集在一起的 tags 的 id 等信息，我们根据 tags 聚集情况对多位点 tags 进行打分， 评价其真实的转录位点. 除了四个样本数据库之外，本文还构建了一个 hts 数据库，用于存储 4 个样品中 reads 数量大 于 3 的 tags 记录，该数据库可以用于 tags 在不同样本中的表达差异分析. sample1 (14) tagsid(1) tags_map(2) tags_map_fbintergenic(3) tags_map_genes(4) tags_map_genes_3cross(5) tags_map_genes_5cross(6) tags_map_genes_ant(7) tags_map_intergenic(8) tags_map_mrna_cds(9) tags_map_mrna_cds_cross(10) tags_map_mrna_intron(11) tags_to_cluster(12) tags_to_cluster_ant(13) cluster_to_unit(14) fig.5 sample database structure table 4 sample tables function description db.tables function description (1)* storage tags sequence,read count,hits count and so on. (2) storage tags mapping information,contain genome structure and ncrna/marker information. (3)/(8) tags mapped to intergenic relation (n:1). (4)/(5)/(6) tags mapped to genes relation (n:1). (7) tags mapped genes_antisece relation (n:1). (9)/(10)/(11) tags mapped genes exon/intron structure relation (n:1). (12) storage cluster site information,tags-count,read-count,so on. (13) storage cluster_antisece site information (14) storage unit site information,cluster-count,tags-count,so on. *(1)- (14) tables name are same as the fig.5 tables 中国中国科技论文在线科技论文在线 - 7 - 2.52.5 高通量数据处理流程与方法高通量数据处理流程与方法 .1 处理流程处理流程 trimmed adapter reads convert into tags perl script tags mapped to genome soapaligner loaded to myql perl script rice genes intergenic fb_intergenic tags_map_genes tags_map_genes_ant tags_map_intergenic tags_map_ sample2 sample3 sample1 raw data tags mapping data tags data tagsid and tags_map tables tags site classification mysql tags_map_mrna_cds tags_map_mrna_intron tags_to_cluster tags_to_cluster_ant tags genes site classification tags convert into cluster cluster convert into unit perl script rice mir14 ncrna repeat rice map_genes map_mir14 map_ncrna mysql ncrnas express feature hyp-mirna structure-filtered hyp-pre-mirna target-filtered hyp-mirna new-mirna sample1、2 tags_map tags_to_cluster cluster_to_unit (a) (c) (b) (c) (d) (b) (a) fig.6 data flow (a) basic processing. (b) new-mirnas data mining. 1694 hpy-mir 142 osa-mir 662 hpy-mir 103 osa-mir 135 new-pre-mir 102 osa-pre-mir tags mapped to mirna/ncrna/marker 中国中国科技论文在线科技论文在线 - 8 - .2 处理方法处理方法 1)数据基本处理. 如(图 6a)所示，包括步骤-，概述如下： 步骤：reads 接头去除. 由于部分 reads 3端序列质量较低，错误的碱基较多，为去除接头带来 一定难度，本实验采用的动态接头去除法，对末尾低质量的碱基进行递归以确保最大程度的去除接 头. 与样品 1、2 的原作者处理结果相比较32，我们的处理结果要优于前者. 步骤：tags 映射到基因组. 目前针对高通量测序开发的映射软件主要有 soap(soapaligner)、 eland、maq、bowtie、genomemapper 等，我们使用了华大基因组开发的 soap 2. 0，其优点是可 以输出多位点 tags 的所有映射位点信息，参数设置允许 2 个错配无 gaps. 步骤：tags 映射信息导入数据库. 由于多位点的 tags 的位点信息量过大，我们根据 tags 的表达 的有效性进行筛选，hits 1-6 的 tags 映射信息全部导入，hits 7-20 的 tags，将所含 reads 数量大于 3 的位点信息导入数据库，hits 21-50 的 tags，将所含 reads 数量大于 10 的位点信息导入数据库. 已知 trna(met_cat)在基因组上的位点数最多(hits=50)， 所以我们将 hits 值大于 50 的 tags 不再导入数据 库. 经过对多位点 tags 的筛选，映射到基因组的 reads 中约 90%的数据被导入到数据库(表 6). 步骤：tags 位点结构归类. 将 tags_map 表中的位点信息与 rice 数据库中的 genes、intergenic 及 fb_intergenic 的位点信息进行匹配查询，匹配结果分别导入到不同的结构表中，同时更新 tags_map 表的 map_label 字段，该字段用不同的标识代表 tags 所在位点的基因结构属性. 我们按 map_lable 字 段分组汇总即可获得 tags 在基因组上的分布情况. 步骤：tags 位点注释. 将各样本的 tags_map 中的位点信息与 rice 数据库的 mir14 和表 ncrna 分别进行匹配查询，查询结果导入表 map_mir14 和表 map_ncrna，同时更新 tags_map 表的 ncrna_label 字段；对 4 个 marker 表匹配查询结果分别导入 map_repeat、map_fst、map_mule、 map_organelle 4个表， 同时将各类ncrna的标识写入tags_map的gff_label字段. 我们对ncrna_lable 字段的汇总即可获得 tags 在已知 ncrna 中的表达情况(图 10)， 对 gff_label 字段的汇总即可得到 tags 在已知 marker 上的分布情况. 步骤： tags 聚簇. 我们将表 tags_map 的位点信息聚集成簇，结果导入 tags_to_cluster 表；在 将 cluster 信息聚集成单元，结果导入 cluster_to_unit 表. 生成 cluster 之后, 我们根据 cluster 包含 的 tags 聚集特点对 hits1 的 tags 进行定位处理. 2)新 mirnas 挖掘. 如(图 6b)所示，包括步骤(a)-(d)，概述如下： 步骤(a)：初步筛选. 通过对已知 pre-mirnas 的表达特征的分析，设计算法，对 tags_to_cluster 表中没有映射到已知 ncrnas 的 cluster 进行筛选，再从筛选的 cluster 中对 tags 进行筛选，筛选出来 的 tags 称为假设的 mirnas(hyp-mirnas). 对 cluster 筛选的主要条件：18=5 or rc/tc=50)； 对 tags 筛选的主要条件：1810 ； 步骤(b)：靶基因过滤. 如对应的小 rna 样本阶段存在相应的 mrna 表达谱数据(如转录组芯片 数据，数字标签测序数据)，表达谱中下调的 mrna 可以假设为 mirnas 的靶基因(target)，以这些 靶基因为过滤材料，通过靶基因预测我们可以过滤出一大批符合 mirnas/target 配对条件的 hyp-mirnas. 步骤(c)：二级结构过滤. 提取 hyp-mirnas 所在 unit 的 rna 序列，如果 unit 长度小于 80nt，则 两边各取 150nt. 我们采用 rnalfold 程序对提取的 rna 进行了筛选，该程序可以在指定的长度(-l) 内递归搜索二级结构最稳定的子序列， 我们从最大长度开始以 20nt 为单位递减， 但最短不小于 80nt， 该值作为 -l 的参数，筛选满足无内环的发夹结构并含有 tags 位点的最长子序列，该子序列作为假 设的 mirnas 前体(hyp-pre-mirnas). 中国中国科技论文在线科技论文在线 - 9 - 步骤(d)：pre-mirnas 鉴定. 我们采用 mipred 程序对 hyp-pre-mirnas 序列进行打分，一般得分 在 52 以上的序列，mipred 就认为其是真实的 pre-mirnas，而真实的 pre-mirnas 上处于发夹结构 茎区的表达量最高的 tags 我们就认为其是新的 new-mirna. 3 3 结果结果 按照图 6 所示的流程和方法，分别对四个样本数据集进行了处理，应用 mysql 查询语句可以从 不同角度对每个数据表进行统计，这里主要从以下这几个方面展示处理结果。 3.13.1 基本数据处理结果基本数据处理结果 3.1.1 样本样本 reads/tags 数量数量统计统计 由于处理过程多处优化，导入数据库中用于分析的数据比例比原始文献都要高32. 如果进行样本 间的表达差异比较，则 reads 计数需要做标准化处理，由于 4 个样本分离的小 rna 长度有所不同， 我们选取长度在 17nt28nt， tags 包含的 reads 数大于 3 的所有 reads 计数做平均值标准化处理(表 6). table. 6 sample reads / tags count sample raw date mapped data load to mysql data nf tags reads tags reads tags reads cluster unit per* per* 1 1593541 3428861 1239641 2963375 1069018 2727537 589267 69511 92. 0% 76. 6% 1. 05 2 1547469 3412945 1224883 2857403 1041850 2610522 617266 62068 91. 4% 68. 0% 1. 00 3 2522046 7956301 1712339 6487265 1362016 5791197 1054495 199711 89. 3% 47. 1% 0. 43 5 3795643 13665389 2628717 10656649 2032599 9633619 1643696 243367 90. 4% 32. 0% 0.56 per*,percentage of load_db_reads/mapped_reads. per*, hits=1 reads percentage same as per*. nf, normalization factor. 3.1.2 样本样本 reads/tags 长度分布长度分布 样品 1， 种子发育 1-5 天， 分离 18-25nt 范围的小 rna. 将去除接头后的数据和导入数据库的数据 按 reads/tags 数量对长度做图(图 7)，图中 raw_rc 曲线比较规则，说明 rna 样本完整性较好；峰值 出现在 21、24，说明此长度的小 rna 大量表达. 样品 4，花粉发育 3 核阶段，分离 12-40nt 范围的小 rna. 将不同处理阶段的数据做图(图 8)， 图中 raw_rc 曲线相对不是很规则，说明样本完整性一般；将已知 ncrnas 过滤后得到 free_rc 曲线。 比较 db_rc 曲线，可以看出峰值 16、19 及 30-39 范围的 reads 基本消失，说明这些 reads 大多是已知 的 ncrna 降解片段；而 free_rc 曲线中仍然存在 21-23nt 及 23-25nt 两个峰，尤其 24nt 峰值很高，说 明还有大量的未知 small rnas 待挖掘. fig.7 sample1 reads/tags length distribution. raw_rc, reads count in raw data. db_tc, tags count of load mysql 中国中国科技论文在线科技论文在线 - 10 - 3.1.3 reads 位点分布图位点分布图 根据样本数据库 tagsid 表中的 hits 字段分组，对 r_count 字段汇总，即可得到 reads 在基因组上 的位点分布情况(图 9)，由图可知，不同组织的样本小 rna 表达存在显着的差异. 如样本 1、2 为同 一组织， 差异并不明显， 而叶片与成熟花粉多位点小 rna 比例很高， 这说明该组织中大量位于 repeat 区域的小 rna 表达，或者某些多位点小 rna 大量表达，也可能为转录后扩增现象所致。 3.1.4 reads 在基因组及在基因组及已知已知 ncrna 分布情况分布情况 根据tags_map表中的ncrna_label字段分组， 对r_count字段汇总， 可以得到reads在已知ncrnas 中的表达情况(图 10). 由图可知， 不同组织中的 ncrnas 的表达存在显著差异， 如叶片样本 3 中已知 的 mirnas 大量表达， 占总 reads 的 47. 6%； 在花粉样本 4 中， 已知的 trna 大量表达， 约占总 reads 的 51. 8%；图中 intergenics、genes 和 gene_antisence 部分是没有映射到已知 ncrnas 的 reads 部分， 其中 intergenics 和 gene_antisence 的比例可以代表样本中未知的有待挖掘的小 rna 部分. 对样本 4 分析发现，大量的 trna 降级片段存在与 30-39nt 范围内(图 8)，因其他 3 个样本分离 的小 rna 长度都不超过 30nt，所以我们对样本 4 中小于 30nt 的 reads 重新做了一个统计(图 11)，由 图可见已知的 trna 仍然占了总量的 34%，这说明花粉中 trna 确实存在大量表达，或者 trna 在 样本 4 中发生了大量降解. fig.9 sample1-4 reads hits distribution. fig.8 sample4 reads/tags length distribution. free_rc, reads count of no mapped known ncrnas. fig.11 sample4 reads distribution on ncrnas. na, no mapped known ncrnas reads. 中国中国科技论文在线科技论文在线 - 11 - 3.23.2 已知已知 mirnasmirnas 表达特征分析表达特征分析 3.2.1 已知已知 ncrnas 的的表达特征表达特征 如果一些 tags 片段是来自一个较长前体的随机降解， 那么这些 tags 头尾重叠的概率较大， cluster 处理后就很有可能聚集在一个 cluster 中， 这样的一个 cluster 所包含的 reads/tags 值应该较小. 如果这 些短的 tags 片段是来自前体的规则剪切，那么 tags 间重叠的概率就较小，剪切点两侧的序列在一个 cluster 中的概率就很