新一代高通量RNA测序数据的处理与分析.doc

《新一代高通量RNA测序数据的处理与分析.doc》由会员分享，可在线阅读，更多相关《新一代高通量RNA测序数据的处理与分析.doc（16页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、Reviews and Monographs 综述与专论 生物化学与生物物理进展 Progress in Biochemistry and Biophysics 2010, 37(8): 834846 新一代高通量 RNA 测序数据的处理与分析 * 1) 王 曦 1) 汪小我 1) 王立坤 1, 2) 冯智星 1) 张学工 1)* ( 生物信息学教育部重点实验室，清华信息科学与技术国家实验室 (筹 )生物信息学研究部，清华大学自动化系，北京 100084； 2) 吉林大学计算机科学与技术学院，长春 130012) 摘要 随 着新一代高 通量 DNA 测序 技术的快速 发展， RNA 测序 (

2、RNA-seq)已 成为基 因表达和 转录组 分析新 的重要 手段 RNA-seq 技术产生的 海量数据为生物信息 学带来了新的机遇和 挑战 有效地对 测序数据进行针对 性的生物信息学 处理和分 析，成为 RNA-seq 技术能否在科学探索中发挥重大作用的关键 以新一 代 Illumina/Solexa 测序平台所产生的数据为例，在扼 要介绍高通量 RNA-seq 测序流程的基础上，对 RNA-seq 数据处理和分析的方法和现有软件做一个较为全面的 综述，并对其 中有待进一步研究的问题进行展望 关键词 高通量 RNA 测序，转录组，基因表达，数据处理与分析，生物信息学 学科分类号 Q5， Q6

3、， Q7 DOI: 10.3724/SP.J.1206.2010.00151 近年来，新一代高通量测序技术得到了突飞猛 进 的 发 展 ， 在 此 基 础 上 ， 高 通 量 RNA 测 序 即 RNA-seq 也迅 速发 展 与 基因 芯片 技术 相比 ， RNA-seq 无需设计探针，能在全基因组范围内以单 碱基分辨率检测和量化转录片段，并能应用于基因 组图谱尚未完成的物种 ，具有信噪比高 、 分辨率 高 、 应用范围广等优势，正成为研究基因表达和转 录组的重要实验手段 RNA-seq 为基因组学的研究带来了高分辨率的 海量数据，如何有效处理和分析这些海量数据成为 这一新技术能否带来新的科

4、学发现的关键，一些生 物信息学方法与软件也应运而生 本文 针对当前 RNA-seq 应用的现实情况，尝试以 Illumina/Solexa 测序平台产生的 mRNA-seq 数据为例，对 RNA 测 序数据的产生过程及数据处理和分析的基本流程 、 关键方法和现有软件进行较全面的介绍，并 讨论 RNA-seq 数据分 析中存在的挑战 1 高通量测序技术简介 诞生于 20 世纪 70 年代的 Sanger 法 是最早被 广泛应用的 DNA 测序技术 7，也是完成人类基因 组计划的基础 但是，由于它测序通量低，费时费 力，科学家们一直在寻求通量更高 、 速度更快 、 价 格更便宜 、 自动化程度更高

5、的测序技术 自 2005 年以来，以 Roche 公司的 454 技术 、 Illumina 公司 的 Solexa 技术和 ABI 公司的 SOLiD 技术为标志的 新一代测序技术相继诞生 8 新一代测序技术又称 作深 度测序技术，主要特点是测序通量高 、 测序时 间和成本显著下降 把这种高通量测序技术应用到由 mRNA 逆转 录 生成 的 cDNA 上 ，从 而 获 得来 自 不 同基 因 的 mRNA 片段在特定样本中的含量，这就是 mRNA 测序或 mRNA-seq 同样原理，各种类型的转录本 都可以用深度测序技术进行高通量定量检测，统称 作 RNA-seq 或 RNA 测序 目前，在

6、已经推出的几 种新一代测序平台中， Illumina/Solexa 测序平台上 的 RNA-seq 应 用 最 广 ， 我 们 以 此 为 例 来 综 述 RNA-seq 数据处理和分析的生物信息学问题和方法 . * 国 家 自 然 科 学 基 金 资 助 项 目 (60702002, 60721003, 30873464, 60905013)和东南大学生 物电子学国家 重点实验 室开放研 究基金资 助项目 . * 通讯联系人 . Tel: 010-62794919, E-mail: 收稿日期： 2010-03-25，接受日期： 2010-04-30 15 6 2010; 37 (8) 王曦

7、等：新一代高通量 RNA 测序数据的处理与分析 835 Illumina/Solexa 测序技术的基本原理是边合成 边测序 (sequencing by synthesis， SBS) ，即测序 过程是以 DNA 单链为模板，在生成互补链时，利 用带荧光标记的 dNTP 发出不同颜色的荧光来确定 不同的碱基 新加入 dNTP 的末端被可逆的保护基 团封闭，既保证单次反应只能加入一个碱基，又能 在该碱基读取完毕后，将保护基团除去，使得下一 个反应可继续进行 为了增加荧光强度， 使之更 易被成像系统所采集，该技术在测序之前还需要对 待测片段做桥式 扩增 (bridge amplification)

8、 (http:/ 初 期 的 Illumina/Solexa 测 序 技术只能在较短的测序读长上 (20 30 碱基 )保证较 高的正确率 随着技术的改进，目前的读长已经增 加到 100 碱基以上 同时 ，随着双端测序 (paired- end， PE)技术的成熟，测序长度更可达到单端测序 的 2 倍，测序通量也随之增加 这种测序 技术是 Solexa 公司发展起 来的， 2007 年被 Illumina 公 司 收购，因此现在通常被称为 Illumina/Solexa 测序技 术 近 两年 来， Illumina/Solexa 测 序平 台不 断 升 级 ，相 继推 出了 GA (Geno

9、me Analyzer)、 GA IIx、 HiSeq 2000 等测序仪 更多关于高通量测序平台的 介绍，可以查阅相关文献 9, 14 16 2 RNA-seq 测序文库制备和测序平台数据 输出 本小节针对 Illumina/Solexa 测序平台，对 RNA 测序文库制备标准和平台底层数据产生做一个简单 的介绍 2 1 RNA-seq 测序文库制备 对 于 mRNA-seq 实 验 ， 从总 RNA 到 最 终 的 cDNA 文库制备完成主要包括以下步骤 首先，用 Poly(T)寡聚核苷酸从总 RNA 中抽取全部带 Poly(A) 尾的 RNA，其中 的主要部分就是编码基因所转录 的 mR

10、NA 将所得 RNA 随机打断成片段，再用随 机引物和逆转录酶从 RNA 片段合成 cDNA 片段 然后，对 cDNA 片段进行末端修复并连接 测序接 头 (adapter)，得到将用于测序的 cDNA 在以上过 程，将 RNA 随机片段化和采用随机引物进行反转 录，都是为了使所得 cDNA 片段较均匀地 取自各 个转录本 为提高测序效率，一般还需要用电泳切 胶法获取长度范围在 200 bp(25 bp)的 cDNA 片段， 再通过 RCR 扩增，得到最终的 cDNA 文库 在上述文库制备过 程中，如果不是只抽 取带 Poly (A) 尾 的 RNA， 而 是 使 用 全 部 的 RNA， 则

11、 RNA-seq 测得的就是细胞中的全部转录本，如果把 带 Poly(A)尾的 RNA 过滤掉，也可以得到非编码的 RNA 转 录本 ， 如 果从 总 RNA 中 只 提 取 长 度为 21 23 个碱基左右的 RNA，则得到全部的 miRNA (microRNA)转录本 , 相应的方法也称作 miRNA-seq. 样品制备最终得到 的是双链 cDNA 文库 在 后续测序中，测 得的每个读段 (read)随机地来自双 链 cDNA 的某一条链，从读段序列本身 无法得知 它是与 RNA 方向相同还是倒转互补，在后续的读 段定位时需要两个方向都考虑 在新基因识别等应 用中，转录本的方向对基因注释尤

12、为重要，需要在 文库制备和测序中保留 RNA 的方向信息 最近有 文献报道了保留方向信息的 RNA-seq 样品制备方 法 17 20 2 2 测序平台数据输出 将 RNA-seq 测序文库加入流动槽 (flow cell)中 的各通道 (lane)，在桥式 PCR 扩增后，就可以进行 测序了 测序过程中，计算机软件同步地对荧光图 像数据进行处理，通过分析荧光信号来确定被测碱 基，并给出质量评分 按照图像上的位置坐标，计 算机程序将同一位置测得的碱基根据测序顺序连成 读段 (read) 由于荧光图像文件所占有的磁盘空间 很大，通常 GA IIx 平台一次实验就能产生上太字 节 (TB)的图像文

13、件，所以一般情况下不予保留原始 的荧光图像数据，而是只保留程序读出的读段数据 及对应的质量分值，这就是多数实验室委托测序中 心进行 RNA-seq 测序后得到的最原始的数据 为了便于测序数据的发布和共享，高通 量测序 数据以 FASTQ 格式来记录所测的碱基读段和质量 分数 如图 1 所示， FASTQ 格式以测序读段为单 位存储，每条读段占 4 行，其中第 1 行和第 3 行由 文件识别标志和读段名 (ID)组成 (第 1 行以 “”开 头而第 3 行以 “+”开头；第 3 行中 ID 可以省略， 但 “+”不能省略 )，第 2 行为碱基序列，第 4 行为 对应的测序质量分数 关于 FAST

14、Q 格式更多地介 绍可参考文献 21 为方便保存和共享各实验室产 生的高通量测序数据， NCBI、 EBI、 DDBJ 等数据 中 心建 立了 大容 量的 数 据库 SRA (Sequence Read Archive, http:/www.ncbi.nlm.nih.gov/Traces/sra) 来 存 放共享的测序数据 22 23 1012 13 836 生物化学与生物物理进展 Prog. Biochem. Biophys. 2010; 37 (8) 每 4 行标识为 一个测序读段 3 RNA-seq 数据的基本处理 RNA-seq 的基本应用是测量一个样本中的基 图 1 Fig. 1 读

15、段 FASTQ 数据格式示例 FASTQ format examples 读段识别码 碱基序列 +读段识别码 测序质量分数 因表达或转录组 有实验表明，新一代高通量测序 技术重复数据之间的相关度较高 (R 0.96) ，因 此，如果对同一样本在多个通道上 进行了 RNA 测 序的技术重复，我们建议可以把几个通道的数据进 行合并，这样等效地增加了测序深度 本节讨论单 个样本 RNA 测 序数据的基 本处理流 程，如图 2a 所示 (a) 全基因组 序列及基 因注释 基因表达水平 估计 (及新基 因的识别 ) RNA-seq 测序 数据 读段定位 剪接异构体的 表达水平推断 后续处理 可视化及 读

16、段注释 选择性剪接 事件识别 (b) RNA-seq 测序数据 (样本 A) RNA-seq 测序数据 (样本 B) RNA-seq 数据 处理基本流程 RNA-seq 数据 处理基本流程 基因表达 差异分析 剪接异构体 表达差异或 选择性剪接 差异分析 分类分析 其他高层 数据分析 图 2 RNA-seq 数据处理和分析流程图 Fig. 2 The flowchart for RNA-seq data processing and analysis (a)RNA-seq 数据的基本处理 , 其方法介绍见正文第 3 节 . (b)两类样本 RNA-seq 数据比较分析的框架 , 对应于正文的第

17、 4 节 . (b)中虚线框内为(a) 所示的流程 , 虚线箭头表示可选输入 . 3 1 读段定位 获得 RNA-seq 的原始数据后，首先需 要将所 些基本的预处理 例如，过滤掉测序质量较差的读 段 、 对 miRNA 测序读段数据去除接头序列等 有测序读段通过序列映射 (mapping)定位 到参考基 因组上，这是所有后续处理和分析的基础 在读段 定位之前，有时还需要根据测序数据情况对其做某 高通量测序的海量数据对计算机算法的运行时 间提出了很高的要求 针对诸如 Illumina/Solexa 等 测序平台得到的读段一般较短 、 且插入删除错误较 2 12 2010; 37 (8) 王曦等

18、：新一代高通量 RNA 测序数据的处理与分 析 837 少等特点，人们开发了一些短序列定位算法 这 位 种 子 片 段 索 引 法 的 代 表 是 Maq 26 ， 而 采 用 些 算 法 主 要 采 用 空 位 种 子 索 引 法 (spaced-seed Burrows-Wheeler 转换的 代表是 Bowtie 总 的 来 indexing)或 Burrows-Wheeler 转换(Burrows-Wheeler Transform， BWT)技术来实现 24 空位种子索引法 首先将读段切分，并选取其中一段或几段作为种子 建立搜索索引，再通过查找索引 、 延展匹配来实现 读 段 定 位

19、 ， 通 过 轮 换 种 子 考 虑 允 许 出 现 错 配 (mismatch)的各种可能的位置组合 BWT 方法通过 B-W 转换 将基因组序列按一定规则 压缩并建立索 引，再通过查找和回溯来定位读段，在查找时可通 过碱基替代来实现允许的错配 表 1 列出了目前可 免费下载使用的部分短序列定位软件 其中采用空 说，采用 BWT 的定位算法在时间效率上要优于空 位种子片段索引法 24,28 随着读长的增加，允许读 段序列中存在插入删除 (indel)的定位变得可行而重 要 由于以上两类方法对序列中插入删除的处理较 为困难，近来人 们开发了一些基于 改进的 Smith- Waterman 动

20、态 规划 算 法 29 的 序列 比 对 工具 ， 如 BFAST 、 SHRiMP 、Mosaik(http:/bioinformatics. bc.edu/marthlab/Mosaik)等 ，但 算法 速度较 慢， 大 多需采用计算机并行编程技术来解决运行时间的 问题 表 1 Table 1 适用于 Illumina/Solexa 测序平台的读段定位软件 Mappers/aligners for Illumina/Solexa sequencing data 名称 SAM1) 质量 2) 主要采用技术 网址 MAQ 26 27 否 是 空位种子 http:/ 1) Bowtie BWA3

21、2 ZOOM33 ELAND SOAP234 RazerS35 Novoalign SHRiMP31 BFAST30 Mosaik 是 是 否 否 否 否 是 否 是 是 是 是 否 否 否 否 是 是 是 是 2) BWT BWT 空位种子 空位种子 BWT q-grams 过滤 Needleman-Wunsch 算法 空位种子 q-grams 过滤 Smith-Waterman 算法 Smith-Waterman 算法 并行编程 Smith-Waterman 算法 并行编程 http:/bowtie- http:/bio- http:/ http:/ http:/ http:/www.se

22、qan.de/projects/razers.html http:/ http:/compbio.cs.toronto.edu/shrimp https:/secure.genome.ucla.edu/index.php/BFAST http:/bioinformatics.bc.edu/marthlab/Mosaik SAM: 是否能以 SAM 格式输出 ; 质量 : 是否提供读段定位质量信息 ; BWT: Burrows-Wheeler 转换 . 在 RNA 测序数据的基因组定位中，一个特殊 的问 题是 跨越两 个外显 子接合 区的读 段 (junction reads)定位 在真核生物中

23、，成熟的 mRNA 是经过 由 mRNA 前体中的外显子经过剪接形成的 如果 一个读段跨越了两个外显子，那么就无法将这个读 段完整地定位到基因组序列上 而同时，这种跨两 个外显子的读段在分析转录本的剪接形式和研究选 择性剪接中有重要的作用 为了解决这一问题，人 们采取两种典型的策略来进行接合区读段的定位： 一是根据已知的基因外显子注释，构建所有可能的 外显子接合区序列，与基因组序列一并作为定位的 参考基因组；二是不依赖基因注释，而是先利用能 完整定位到基因组的读段得到粗略的外显子区域， 并结合剪接位点序列构建出可能的剪接位点，然后 将不能完整定位的读段分段定位到两个外显子可能 的结合区域 Il

24、lumina/Solexa 平台提供的 RNA- seq 软件分析包 GApipeline 采用了第一种策略 采用 第二种策 略的软 件有 Tophat 36和 G-Mo.R-Se37 等， 最新的 Tophat 软件增加了利用已知外显子边界注 释信息的选项 27 25 30 31 838 生物化学与生物物理进展 Prog. Biochem. Biophys. 2010; 37 (8) 不论是哪种测序平台，测序中都不可避免地存 在一定的错误，基因组中又存在单核苷酸多态性等 引起的序列变化，所以在读段定位时通常允许一定 数量的错配，可以根据不同应用调节允许错配的程 度 另一方面，由于基因组中重复

25、序列和高相似度 序列的影响，某些读段会出现定位到基因组多个位 置的情况 这些因素影响了各个读段到基因组的定 位质量，在一些新的读段定位算法中，同时给出每 个读段与基因组匹配质量 通常在后续处理前，人 们将多定位的 读段都过滤掉，也有人尝试用适当的 策略把多定位读段 “分配 ”到其中某些位置上 . 读段定位到基因组后通常采用 SAM(Sequence Alignment/Map)格式或其二进制版本 BAM 格式 来存储 二进制版本可大大节省存储空间，但不能 直接用普通文本编辑工具显示 关于 SAM 格式的 详 细 介 绍 ， 可 查 阅 (http:/ SAM1.pdf) 3 2 基因表达水平估

26、计 在深度测序技术出现之前，高通量测量不同基 因表达水平的主要手段是基 因芯片，在此基础上可 以对不同组织或者不同发育阶段的基因表达差异和 模式进行分析 mRNA-seq 数据最基本的应用也是 检 测 基 因的 表 达 水 平 ，与 基 因 芯 片 数 据 相 比 ， RNA 测序得到的是数字化的表达信号，具有灵敏 度高 、 分辨率高 、 无饱和区等优势 RNA 测序数据是对提取出的 RNA 转录本中随 机进行的短片段测序，如果一个转录本的丰度高， 则测序后定位到其对应的基因组区域的读段 也就 多，可以通过对定位到基因外显子区的读段计数来 估计基因表达水平 很显然，读段计数除 了与基因 真实表

27、达水平成正比，还与基因长度成正比，同时 也与测序深度即测序实验中得到的总读段数 正相 关 为了保持对不同基因和不同实验间估计的基因 表达值的可比性，人们提出了 RPM 和 RPKM 的概 念 RPM(reads per million reads) 即每百万读段中 来自于某基因的读段数，考虑了测序深度对读段计 数 的 影 响 RPKM(reads per kilo bases per million reads)是每百万读段中来自于某基因每千碱基长度 的读段数，公式表示为： 基因长度 测序深度 RPKM 不仅对测序深度作了归一化，而且对基 因长度也作了归一化，使得不同长度的基因在不同 测序深度

28、下得到的基因表达水平估计值具有了可比 性 ，是 目前 最常 用 的基 因 表达 估 计方 法 软 件 rSeq43、 DEGseq 软件包 44和 Cufflinks45等都提供了 用上述方法进行基因表达水平计算的功能 根据 RNA-seq 文库制备标准，在 不考虑基因 结构的理想情况下，读段会均匀地分布在基因上 而实际上，通过对实际数据的可 视化分析很容易发 现，读段在基因上的分布有着自身的一些模式，呈 现出不均匀性 (图 3) 这一问题已经引起很多学者 的关注 造成读段分布出现偏好的原因可能有 多个方面：在制备 cDNA文库时，反转录所采用的 随机 引物对 RNA 序列具 有一定 的偏好

29、性，使 得 cDNA 片段不能够完全 均匀地取自各 转录本； 在 PCR 扩增中，扩增效率与序列的 GC 含量等特征相 关，可导致 GC 含量高的 cDNA 片段在文库中拷贝 数增加超过其他片段；舍弃多定位的读段也可能导 致读段的非均匀分布；等等 如果 能对读段分布的 不均匀性进行建模并加以校正，可以提高 RNA-seq 推断基因表达量的准确度 但根据对实际数据的观 察，对于较长转录本，读段非均匀分布带来的误差 很大程度上可相互抵消，用 RPKM 来估计基因的 表达水平可以得到比较满意的结果 3 3 选择性剪接事件识别和剪接异构体表达水平 推断 在真核生物中，选择性剪接现象普遍存在 基 因转录

30、形成的 mRNA 前体 (pre-mRNA)在剪接过程 中因去掉不同的内含子区域或保留不同的外显子区 域，可形成不同的剪接异构体 根据 RNA-seq 原 理，只要 测序深度足够深，就能检测到所有转录本 的全部序列，包括来自剪接接合区的序列 利用考 虑到接合区的读段定位方法，就有可能系统地研究 某一组织或某一条件下的基因选择性剪接事件 前面已经提到， Tophat 等软件定位剪接 接合 区读段的策略能标定出剪接事件中的两 个剪接位 点：供体位点和受体位点 通过比较供体位点和受 体位点的组合，就能识别选择性剪接事件 图 3 中包含 了选择性剪接识 别的一个例子 进 一步， 通过对供体和受体位点的

31、读段计数，结合外显子其 他区域的读段数据，还能定量地计算选择性剪接 事 件之间的比例 50 51 对于每一个剪接异构体， RNA-seq 数据能在一 2, 38 39 4042 2 基因区读段计数 RPKM= 10 2010; 37 (8) 王曦等：新一代高通量 RNA 测序数据的处理与分析 839 定程度上推断其表达水平 比如，可以根据已知外 显子组成和各外显子长度对剪接异构体建立数学模 3 5 读段的可视化及注释 对于复杂的组学数据，能尽可能方便地直接 观 型，在测序读段转录本上均匀分布的假设下，利用 各外显子上的读段数和接合区读段数求解异构体的 表 达 值 Jiang 等 43 的 方

32、法 及 软 件 IsoInfer 52 和 cufflinks 都采用了这种思路来实现剪接异构体的 表达推断 需要指出的是，某些形式的剪接异构体 表达水平在这种方法框架下不可辨识 53 3 4 新基因的检测 在对 RNA-seq 数据的分析中，人们发 现 ，往 往不是所有读段都能定位到已有注释的基因区，说 明除了转录噪声或测序错误等的影响外，可能还存 在尚未被注释的基因 这里，我们把这种尚未注释 的基因称为新基因，包括新的蛋白质编码基因和非 编码 RNA 基因 能检测新基因，尤其是低表达基 因是 RNA-seq 技术优于基因芯片的特点之一，因 为它不需要利用已知基因注释来设计检测探针 RNA-

33、seq 技术灵敏度高，但样品污染 、 测序错 误等仍可能带来背景噪声 从基因组未注释区域的 RNA 测序读段信号中检测新基因是典型的信号检 测问题 如何控制新基因识别的误 发现率 (FDR)是 检测方法的关键 Useq 软件包 将 ChIP-seq 数据 分析的方法移植到 RNA-seq 数据上，用滑窗的方 法来识别测序读段定位富集的区域，给出反映滑窗 所在区域读段富集显著程度的 P 值 (P-value)及新基 因误发现率，通过设定 P 值或误发现率的阈值，可 筛选出读段富集的区域，再将相邻区域合并或根据 剪接接合区读段将相应区域连接，完成新基因的 检测 察数据对于数据的分析和解释都非常重要

34、，对新一 代测序数据的可视化和交互展示是一个非常重要但 容易被人忽 视的问题 不深入考查数据的细节，而 是满足于对数据的统计分析，是高通量数据应用中 经常容易陷入的误区，方便有效的可视化工具能够 帮 助避 免这 样的 误 区 表 2 列 出 了部 分 适用 于 RNA-seq 数据的全基因组浏览器，其中比较具有代 表 性 的 有 UCSC Genome Browser、 CisGenome Browser 和 IGV(Integrative Genomics Viewer)等 这 些浏览器具有如下特点： a 能在不同尺度下显现 单个或多个读段在基因组上的位置，包括来源于 剪 接接合区的读段；

35、b 能在不同尺度下显示不同区 域的读段丰度，以反映不同区域的转录水平或测序 效率； c 能显示基因及其剪接异构体的注释信息； d 能显示其他注释信息，例如物种间基因组序列 保守性 、 序列 GC 含量等； e 能直接或间接支持 SAM/BAM 读段定位数据存储格式 . UCSC Genome Browser 属于基 于网络模 式的全基 因组浏览 器， 所有数据都需要上载到远程服务器，经过处理后将 图形返回客户端显示 图 3 中的例子就是从 UCSC Genome Browser 的显示截取的 . CisGenome Browser56 是典型的本地版基因组浏览器，所有读段数据 、 注 释信息都

36、存于本地文件，因此不需要网络连接，方 便内部考查数据用 IGV(http:/www.broadinstitute. org/igv)可以说是以上两种模式的融合 ，既可以从 远程服务器端下载各种注释信息，又可以从本地加 载注释信息 表 2 适用于 mRNA-seq 数据的全基因组浏览器 名称 Table 2 The genome browsers/viewers applicable to mRNA-seq data viewing 支持的数据格式 网址 IGV UCSC Genome Browser 55 GFF/GFF3, BED, SAM/BAM, WIG, BED, bigBed, BE

37、DGRAPH, GFF, GTF, WIG, bigWig, BAM, http:/www.broadinstitute.org/igv http:/genome.ucsc.edu CisGenome Browser MochiView SeqMonk Gambit GenomeGraphs57 56 BAR, BED,refFlat, FA, Wig, BED, GFF, FASTA, ELAND, GFF, BED, MAQ, SAM BAM, BED, GFF2, GFF3, FASTA, VCF Expression data, Annotation tracks http:/biog

38、ibbs.stanford.edu/ jiangh/browser http:/johnsonlab.ucsf.edu/sj/mochiview-start http:/www.bioinformatics.bbsrc.ac.uk/projects/seqmonk http:/ http:/bioconductor.org/packages/release/bioc/html/GenomeGraphs.html 以上列出的全基因组浏览器均可在 Windows、 Linux 和苹果公司的 Mac OS 等计算机平台下运行 . 45 54 840 生物化学与生物物理进展 Prog. Biochem

39、. Biophys. 2010; 37 (8) 剪接接合区 样本 A 测序标签分布 剪接接合区 chr22: 50 0 37860000 可变剪接 10kb 37865000 37870000 37875000 测序标签的非均匀分布 样本 B 测序标签分布 基因注释 50 CBX7 该基因在 A、 B 样本中差异表达 内含子区的测序标签 图 3 mRNA-seq 数据可视化示例 Fig. 3 An example for mRNA-seq data visualization 图示区域为人类基因 CBX7. 图中红色表示样本 A 的数据，蓝色表示样本 B. 各轨道 (track)依次为 : 基

40、因组坐标 、 样本 A 的剪接接合区 、 样本 A 的读段分布 、 样本 B 的剪接接合区 、 样 本 B 的读段分布 、 UCSC 基因注释 . 图中还标识了： 因受体位点 不同而形成的选择性剪接 ; 基 因的 5端出现读段的非均匀分布 ; 在两个样本中 , 差异表达基因的读段信号强度不同 ; 在基因标注的内含子 (intron)区域存在少量不连续的 读段 . 除对读段的可视化外，用描述统计学方法对实 验数据进行分类别统计也十分重要 例如，统计读 段在各个染色体上的分布情况和在注释的外显子 、 内含子 、 剪接接合区 、 基因间区的分布情况等 目 前，已经有一些用于测序数据注释的生物信息学软

41、 件，比如 SAMtools 、 BEDtools 等，但由于测序 技术发展迅速，用户需求因人而异，用户经常还需 要根据需求编写一定的程序或脚本完成或完善注释 分析的任务 对于熟悉图 形用户界面的研究人员， 还可以 利用 UCSC Table Browse56 和 Galaxy59 60 来 配合完成注释分 析 由于 UCSC Table Browser 集 成了大量基因组尺度上的注释信息，而 Galaxy 又 为用户提供了书写简单 、 接口明晰和直观的数据处 理流程，这种方法十分方便有效，也为很多学者在 展示研究成果时所采用 以上描述了对于基因组已知的物种进行 RNA- seq 数据处理基本

42、流程，图 2a 给出了其主要步骤 若研究对象尚未完成基因组测序，则需要采用读段 的从 头 拼 装 (de novo assembly) 来代替 读段 定 位，后续流程也须做相应的调整 若 RNA-seq 实 验在文库制备时保留了 RNA 的方向信息，则应分 别研究来自正链和反链的转录产物，并通过与基 因 注释 比较 来检 测反 义 转录 本 最近 ， RNA- MATE 软件在 其分析流 程中加 入如此 的处理 策 略 此外，通过分析定位到外显子接合区的读段， 还可以获取转录本结构，这为研究基因的剪接调控 机 理提供了重要信息 5 而利用 RNA-seq 数据提供 的序列信息，通过与 DNA

43、序列的细致比较可分析 转录组的序列差异 (如 SNP 等 ) ，从而研究等位基 因的表达模式 及 RNA 编辑 等 最后需 要指 出， 由于 miRNA 在序列 和结构 上具有 一定的 特 点， miRNA-seq 数据的基本处理流程也与本节所述 有 所 不 同 ， 感 兴 趣 的 读 者 可 参 考 软 件 工 具 miRDeep69提供的处理策略 4 多类样本 mRNA-seq 数据间的比较分析 很多 RNA-seq 实验的目的是为了 比较两种或 多种样本中基因表达或整个转录组的差异，如比较 癌症组织和正常组织的转录组差异等 这些差异既 包括通常意义下的差异表达基因，也主要包括选择 性剪接模式的差异 、 剪接异构体表达的差异 、 非编 码转录本的差异等 这些差异一般可以用一些统计 假设检验方法检测，但这种检验有时会受到测序深 度 、 基因长度等因素的影响 70 71，需要对结果进行 仔细分析，消除可能的混杂因素，必要时可以用读 段 的绝 对 表达 值倍 数 变化 (fold-change) 来 作为 补 充