RNAseq-表达图谱 清华大学本科论文.pdf

清 华 大 学 综 合 论 文 训 练 综 合 论 文 训 练 题目 用 RNA Seq 数据估计剪接异构体 表达的方法调研与程序实现 系 别 自动化系 专 业 自动化 姓 名 马鑫云 指导教师 张学工 教授 2011 年 6 月 10 日 于学位论文使用授权的说明 本人完全了解清华大学有关保留 使用学位论文的规定 即 学校有权保留 学位论文的复印件 允许该论文被查阅和借阅 学校可以公布该论文的全部或部 分内容 可以采用影印 缩印或其他复制手段保存该论文 涉密的学位论文在解密后应遵守此规定涉密的学位论文在解密后应遵守此规定 签名 导师签名 日期 I 中文摘要 RNA Seq 技术是利用深度测序技术来进行转录组水平研究的新兴技术 是当 前生物信息领域的一个研究热点 RNA Seq 技术凭借其高精度 高通量 高检测 范围等一系列优点冲击着传统的转录组研究的方法 RNA Seq 技术在转录组研究方面有很多应用 其中很重要的一个应用是估计 基因和剪接异构体的表达值 本文研究了几种利用 RNA Seq 技术进行表达值估计 的方法 对各方法的特点进行了一定的分析 并且编程实现了其中的一种基于 读 段非均匀分布假设 的算法 关键词 RNA Seq 表达值估计 非均匀分布 II ABSTRACT RNA Seq is a recently developed technology to study transcriptome which is based on deep sequencing RNA Seq has shaken the status of traditional transcriptome researching methods because of its advantages such as unprecedented precision throughput and high dynamic range over them RNA Seq can be applied to many transcriptome researching fields and one of the most important applications is to estimate the expression level of genes and isoforms This paper conducts a study about some methods based on RNA Seq to estimate the expression level of transcriptome analyze the characteristics of each method and realize one of these methods which is based on the hypothesis of reads start position follow non uniform distribution by programming Keywords RNA Seq estimation of expression level non uniform distribution III 目录 第 1 章 引言 1 1 1 研究背景 1 1 2 RNA Seq 原理 1 1 3 RNA Seq 技术研究现状 2 1 4 本文的主要工作 3 第 2 章 基因表达值的估计 4 2 1 基因表达与基因表达值 4 2 2 利用基因芯片和 RNA Seq 进行基因表达值的估计 4 2 3 RPKM 简介 5 2 4 几种利用 RNA Seq 估计表达值的方法 6 2 4 1 基于读段均匀分布假设的基因表达值计算 6 2 4 2 基于读段均匀分布假设的剪接异构体表达值计算 7 2 4 3 基于读段非均匀分布假设的剪接异构体表达值计算 7 2 4 4 几种方法的联系 8 第 3 章 基于读段非均匀分布假设模型 N URD 简介 9 3 1 符号说明 9 3 2 基于读段均匀分布假设模型 URD 简介 9 3 3 基于读段非均匀分布假设模型 N URD 简介 10 3 3 1 在 URD 模型的基础上引入了表征非均匀分布的参数 12 3 3 2 GBC Global Bias Curve 12 3 3 3 LBC Local Bias Curve 13 3 3 4 GBC 矩阵和 LBC 矩阵 14 3 3 5 b 矩阵 15 3 4 似然函数 16 第 4 章 N URD 模型的算法实现 17 4 1 程序概况 17 4 2 程序框架 18 IV 4 3 获取数据 18 4 4 数据预处理 19 4 5 核心算法的实现 20 4 5 1 获取带权值的基因结构矩阵 20 4 5 2 求解最大似然估计 21 4 6 结果输出 23 4 7 程序总流程图 25 4 8 仿真实验 26 4 9 实验结果 28 4 10 程序有待改进之处 28 第 5 章 结论与展望 29 5 1 结论 29 5 2 展望 29 插图索引 31 表格索引 32 参考文献 33 致谢 34 声明 35 附录 A 外文资料的书面翻译 36 附录 B CisGenome browser 介绍 51 1 第1章 引言 1 1 研究背景 随着测序技术的不断发展 测序平台产生的数据量不断增大 表 1 1 如何 从这些海量数据中发掘其背后隐藏的生物意义成为了当今生物信息学的一个挑 战 表1 1 测序技术的发展 1 2 测序技术 通量 轮 第一代测序技术 Sanger ABI3730 DNA Analyzer 56Kbp 第二代测序技术 Roche GS FLX Titanium 500Mbp Illumina HiSeq 2000 HiScan 200Gbp ABI SOLiD 5500 xl 100Gbp 第三代测序技术 Helicos Heliscope 21 37Gbp 高通量 RNA 测序技术 RNA Seq 是利用深度测序技术进行生物信息研究的 一个非常成功的典范 RNA Seq 利用深度测序技术对转录组 transcriptome 进 行测序 进而对所关心的细胞或组织进行转录组水平的研究 RNA Seq 在估计基 因表达 发现单核苷酸多态性 SNP 新基因的检测等诸多重要领域都有非常重 要的作用 以估计基因表达为例 与传统的基因芯片 microarray 技术相比 RNA Seq 技术不需要事先设计探针 对于基因组图谱还未完成的物种同样有效 能够在全基因组范围内进行检测 具有单核苷酸分辨率的高检测精度 而且信噪 比高 应用范围广 已经成为估计基因表达和研究转录组的重要手段 3 本文的研究重点是对当下利用 RNA Seq 数据进行基因和剪接异构体表达值 估计的几种模型进行研究 并对其中一种基于 读段非均匀分布假设 的算法进 行程序实现 1 2 RNA Seq 原理 RNA Seq 是利用深度测序技术对转录组进行测序的技术 RNA Seq 的整个流 程如图 1 1 所示 fr 库 RNA 1 3 下生 相关 发展 在 RNA S fragment 利用深度 A Seq 的过 3 RNA Se RNA Seq 生物信息学 关的文献数 年份 2007 2008 2009 2010 2011 从表1 2可 展之迅速以 Seq 中 被 在这些片段 度测序技术对 过程 eq 技术研 技术凭借其 学的一个研究 数量 表1 2 至 2011 年 可以看到 近 以及在当前生 图1 1 测序的单链 段的单端或 对这些 cDN 研究现状 其突出的优 究热点 表 PubMed 近几 6 月 6 日 近几年与R 生物信息研 2 RNA Seq 流 链 RNA 首先 或者双端拼装 NA 文库进行 优点 迅速成 表 1 2 列出了 几年收录的 R 所有文章 0 6 33 152 139 RNA Seq相 研究领域的影 流程 4 先被打断并 装上接头 行测序得到 成为基因组 了 PubMed 近 RNA Seq 相关 篇 综 0 0 7 17 5 关的文献呈 影响力 并反转录成 c adaptor 制 到读段序列的 组学研究的重 近几年收录 关文献 综述性文章 7 呈指数增长 cDNA 片段 制备成 cDN 的过程就是 重要手段 录的与 RNA 篇 长 可见RNA 段 NA 文 是 是当 A Seq A Seq 3 在利用 RNA Seq 进行表达值估计方面 很多生物信息学实验室都进行了广泛 而深入的研究 2008 年 Mortazavi 等人提出了 读段在基因上均匀分布 的假设 5 在此假 设的基础上提出了利用 RNA Seq 技术来估计基因表达的方法 并提出了一个能够 定量表征基因表达值的概念 RPKM 值 5 2009 年 Hui Jiang 等人基于 读段在基因上均匀分布 的假设 结合基因的 结构信息 提出了利用 RNA Seq 技术进行剪接异构体水平的表达值估计的算法 6 2010 年 武征鹏等人基于 读段在基因上非均匀分布 的假设 提出了一个 新的模型 该模型能够利用读段在基因上的实际分布的信息 进行剪接异构体水 平的表达值估计 7 2010 年 Bo Li 等人利用生成模型 从测序的机理出发建立统计模型 并且 利用统计推断的方法进行剪接异构体水平的表达值估计 8 该方法由于在建模时 利用统计的方法处理读段定位的不确定性 因此可以处理读段多定位的问题 1 4 本文的主要工作 本文的主要工作是研究近几年在生物信息领域提出的关于基因表达值估计和 剪接异构体表达值估计的几种方法 在此基础上 本文还对武征鹏的基于 读段 非均匀分布假设 的算法进行了程序实现 4 第2章 基因表达值的估计 2 1 基因表达与基因表达值 基因表达 gene expression 是指生物细胞在生命周期中 依照中心法则把储 存在 DNA 中的遗传信息通过一定的生物过程转移到具有生物功能的基因产物的 过程 基因的表达值是指由基因生成的基因产物的数量 在生物信息学领域 通 常以一个基因所转录出来的转录本 transcript 的数量来表征该基因的表达值 由于转录本的数量对于细胞当前和将来的行为有着很大的影响 5 因此 对基因 表达值的正确估计是十分重要的 在基因转录成信使 RNA mRNA 的过程中可能发生选择性剪切 同一个基 因可能转录出多个 mRNA 这些 mRNA 被称为 剪接异构体 isoform 对某 个基因中的各个剪接异构体丰度的估计称为剪接异构体表达值估计 由于剪接异 构体的表达值的变化与许多复杂疾病的发展有着密切联系 9 因此剪接异构体的 表达值估计也是非常重要的 2 2 利用基因芯片和 RNA Seq 进行基因表达值的估计 估计基因表达的传统方法是用基因芯片 microarray 技术 基因芯片基于杂 交 hybridization 原理 利用人工合成的特定核苷酸序列作为探针 probe 来 检测与之互补配对的核苷酸序列的丰度 其流程是在一个很小的芯片上种上大量 特定的 DNA 探针 将带有荧光标记的待测样品加入到该芯片中 令其与探针进 行杂交 通过对基因芯片进行荧光成像 可以获取与每个探针杂交的核苷酸序列 的浓度 进而便可获得相应的基因表达值 基因芯片曾经是基因表达值估计的主流方法 但是由于该技术只能用于检测 已知基因的表达值 而且需要事先设计探针 检测范围较小 现在在基因表达值 估计领域中基因芯片的地位正逐渐被新兴的 RNA Seq 技术所动摇 RNA Seq技术相对于基因芯片而言是一个新兴的技术 它以深度测序为基础 可以进行很多生物信息学的分析 估计基因表达是 RNA Seq 技术应用的一个方面 利用 RNA Seq 技术估计基因表达值的做法是先利用深度测序技术获取转录组的 序列信息 该信息以一定长度的读段 read 序列形式来表示 得到读段序列信 息之 即可 测序 在基 位到 2 3 具 于 2 越多 之后 利用 可得到读段 序的结果 因 基因上的读 到某个基因 3 RPKM 简 RPKM 是 RPKM 是 2008 年第一 显然 测 多 因此 用读段定位软 段在基因上的 因此转录本 读段数 在利 因上的读段数 图2 1 利 简介 是利用 RNA 是 Reads Per 一次提出 5 RPK 序深度越深 简单地利用 软件将各读 的分布信息 本的丰度在 R 利用 RNA S 数进行某种 利用 RNA Seq Seq 技术用 r Kilobase p RPKM 的 KM 被定 总的 深 基因的 用被定位到 5 读段定位 m 息 由于 RNA RNA Seq 中 Seq 估计基 种形式的归一 q 进行基因表 用来定量估 per Million m 的计算公式 定位到某基因 的测序深度 长度越长 到某个基因内 mapping 到 A Seq 的测 中就体现为 因表达时 一化来作为 表达值估计的 估计基因表达 mapped rea 式如下 因的读段数 基因长度 被定位到 内部的读段 到参考基因 测序读段是对 为能够被定位 通常采用的 为其表达值的 的流程 4 达值的一个 ads 的缩写 9 10 数 到该基因内部 段数来表示其 因的相应位置 对转录本打 位到该转录 的方法就是 的定量表示 个非常有效的 由 Mortaz 部的读段数 其表达值是 置 打断后 录本所 是将定 示 的工 zavi 数也会 是不合 6 理的 RPKM 通过将定位到某个基因内的读段数对总测序深度以及该基因的长度 进行归一化 来表示该基因的表达值 由于 RPKM 对基因的长度和测序深度进行了归一化 因此 它具有了在不同 样本和同一个样本中的不同基因 剪接异构体之间的可比性 6 2 4 几种利用 RNA Seq 估计表达值的方法 本节介绍近几年发表的几种利用 RNA Seq 技术进行基因和剪接异构体表达 值估计的模型和算法 2 4 1 基于读段均匀分布假设的基因表达值计算 Mortazavi 等人于 2008 年提出了一个 读段在基因上均匀分布 的假设 5 并且基于这个均匀分布的假设提出了一个简单而有效的估计基因表达值的方法 利用基因的 RPKM 值作为其表达水平的一个估计 5 在 读段均匀分布 的假设下 被定位到某个基因内的读段数与该基因的长 度和测序深度呈正比 因此 为了估计某基因的表达值 应当将所有被定位到该 基因的读段数对该基因的长度和测序深度进行归一化 这便是 RPKM 的概念 当所有的读段都被唯一地定位到某个基因时 可以用简单地用 RPKM 来估计 基因表达值 但是当读段存在多定位的时候 尤其是能定位到多个基因的时候 由于该读段到底来自于哪个基因并不能从读段序列和读段定位信息来重现 因此 涉及到这些读段的基因的表达值就不能简单地套用 RPKM 的概念 在 Mortazavi 等人的工作中对这个问题进行了一定的分析 并且提出了一种解决办法 首先利 用唯一定位的读段进行基因表达值的估计 之后对于那些定位到多个基因的读段 以当前基因表达值的比例将其分配到所涉及到的基因 再进行基因表达值的估计 此外 Mortazavi 等人在 2008 年的论文中还对 读段均匀分布 这一假设的 合理性进行了一定的分析 Mortazavi 发现在实验中得到的真实数据并没有呈现很 好的均匀分布特性 读段分布有一定的偏向 Mortazavi 认为其原因可能与 RNA 的二级结构会使引物结合位点的选择产生偏好有关 Mortazavi 还发现先在反转录 之前先进行 RNA 打断会一定程度上减轻读段分布的非均匀性 不过并不能完全 消除 7 2 4 2 基于读段均匀分布假设的剪接异构体表达值计算 由于在转录的过程中有选择性剪接的存在 一个基因很有可能转录出多个剪 接异构体 研究各个剪接异构体的表达值也是一个十分重要的课题 Trapnell 等 人的实验 2010 年 证明 基因的各个剪接异构体表达值的相对高低对老鼠的肌 肉细胞生成有着决定性的作用 10 Humbert 等人的实验 2007 年 和 Beyer 等人 的实验 2008 年 都显示出剪接异构体表达值的变化与许多复杂疾病的发展有着 密切的联系 9 11 Hui Jiang 等人在 2009 年提出了一个基于 读段均匀分布假设 的模型估计 剪接异构体表达值的算法 6 Hui Jiang 证明 当读段在基因上均匀分布时 落入 该基因各外显子内的读段数服从泊松分布 基于基因各外显子之间相互独立的假 设可以得到读段在整个基因上的分布模型 我们称该模型为 URD Uniform Read Distribution 模型 7 得到读段的在基因上的分布模型之后 便可以利用最大似 然估计的方法估计各剪接异构体的表达值 2 4 3 基于读段非均匀分布假设的剪接异构体表达值计算 Dohn 于 2008 年发现 高通量测序产生的数据分布有明显的偏差 12 Mortazavi 于同年发现读段分布的非均匀性 5 现在认为造成读段在基因上分布非均匀的原因不止一个 一方面可能是由于 RNA 的二级结构造成反转录时引物对 RNA 的结合位点具有偏好性 另一方面 GC 含量可能会影响扩增效率 进而导致读段分布的非均匀性 武征鹏等人在 2010 年提出了一个基于 读段非均匀分布假设 的模型 N URD 7 该模型在 Hui Jiang 的 URD 模型的基础上引入了表征非均匀分布的特 征 使得该模型能够更好地刻画真实数据的分布特性 其估计的剪接异构体表达 值更加接近真实值 7 同年 Bo Li 等人基于深度测序的原理建立了一个生成模型 8 该模型的建立 过程模拟了深度测序的几个主要步骤 并且考虑了每个步骤中存在的不确定因素 该模型在所有剪接异构体的表达值和利用 RNA Seq 得到的所有读段序列之间建 立了联系 并利用统计推断中的 EM 算法给出了估计各个剪接异构体表达值的方 法 由于该模型很好地刻画了测序机理 因此能够天然地解决读段分布非均匀的 问题 此外 由于该模型使用统计模型来刻画测序中的不确定性 因此还能够解 决读段多定位的问题 8 2 4 4 几种方法的联系 上述几种方法之间虽然具有很大的差别 研究问题的侧重点也有所区别 但 是还是有着很明显的联系 Mortazavi 等人的方法较为朴素 并且只能估计基因的表达值 但是他们提出 的 利用归一化的读段数来定量表示表达值 的思想为之后几种表达值估计方法 奠定了基础 Hui Jiang 等人的方法将表达值估计的问题从基因水平上升到剪接异 构体水平 并且在表达值估计问题中引入了统计推断的思想 武征鹏等人的方法 对 Hui Jiang 等人的模型进行了改进 使其能够更好地反映真实数据的信息 从而 改善表达值估计的结果 Bo Li 等人的工作是对 均匀分布假设 比较彻底的改 进 不需要依赖于先前的分布假设 能够较好地反映表达值和测序读段序列之间 的关系 同样是处理 读段非均匀分布 的问题 武征鹏等人的方法虽然并没有建立 一个完全基于 读段非均匀分布 假设的模型 理论上不够完备 但是由于该方 法部分借鉴了 URD 模型的很好的性质 如似然函数的凸性 因此能够在不明显 增加问题复杂度的前提下对表达值的估计效果进行较为明显的改进 Bo Li 等人 的工作虽然理论上比较完备 不依赖于分布假设 但是在实际操作中还是需要很 多经验性的知识 如读段起始位置的经验分布 来帮助模型的求解 因此会在一 定程度上影响估计的精度 武征鹏等人利用仿真数据证明 在绝大多数情况下 N URD 模型的估计精度要略高于 Bo Li 的生成模型 7 武征鹏等人的工作更加偏 重经验性 Bo Li 等人的工作更加看重模型的普适性 长远来说 Bo Li 等人的模 型因其比较强的适应能力可能会愈加重要 但是在经验性还很重要的现在 武征 鹏等人的工作也是十分重要的 9 第3章 基于读段非均匀分布假设模型 N URD 简介 3 1 符号说明 为了便于对问题进行描述 本节先介绍一些本文中对于模型描述需要用到的 符号 研究基因 G 该基因的外显子数 n 该基因的剪接异构体数 m 该基因各个外显子的长度 bp 12 n l ll 落在该基因的各个外显子中的读段数 12 n x xx 落在整个基因上的所有读段数 w 该基因各个剪接异构体的表达值 12 m 为了反映该基因各剪接异构体的结构信息 本文引入两个结构矩阵a和b 其 中 a是一个 0 1 指示矩阵 元素1 ij a 表示该基因的第 i 个剪接异构体中含有第 j 个外显子 0 ij a 表示该基因的第 i 个剪接异构体中不含第 j 个外显子 b是一 个浮点类型的指示矩阵 能够比a矩阵反映更加精细的基因结构信息 对于剪接 异构体 i ij b的值反映了定位到该剪接异构体的第 j 个外显子内的读段数的期望 ij b的值越大 表明被定为到该外显子内的读段数的期望也越大 3 2 基于读段均匀分布假设模型 URD 简介 Hui Jiang 等人在 2009 年提出了一个基于 读段在基因上均匀分布 的假设 来计算剪接异构体表达值的模型 URD Uniform Read Distribution 6 Hui Jiang 证明 如果读段被均匀地定位到基因上 则被定位到该基因的一个外显子内的读 段数是一个近似服从泊松分布的随机变量 6 以该基因的第 j 个外显子为例 该 随机变量的参数为 1 m jjiji i l wa 3 1 其中各参数的含义请见第 3 1 节 由此很容易得到估计剪接异构体表达值的最大似然函数为 10 jj x j j j e Lx x 3 2 进一步可以假设该基因的各个外显子上的读段数之间是相互独立的 则可以 得到联合的对数似然函数为 12 11111 log loglog nmnmn njijijjijij jijij Lx xxwl axl wax 3 3 解此优化问题 找到一组使得对数似然函数最大的参数 12 m 则此 12 m 就是最优的剪接异构体表达值 可以证明 联合对数似然函数 式 3 2 是凸的 6 具有良好的优化性质 保 证局部最优解等价于全局最优解 3 3 基于读段非均匀分布假设模型 N URD 简介 虽然基于均匀分布假设的模型在数学上有很好的性质 而且已经有前人基于 该假设建立了有效的模型 但是实际的实验数据却显示 读段在基因上的分布有 着很明显的非均匀性 a EI24 11 b KNG1 c SLCO2B1 图3 1 读段在基因上分布的真实情况 7 图 3 1 是利用 Hui Jiang 等人开发的 CisGenome browser 软件得到的基因读段 分布的可视化表示 该图给出了人类基因组中的三个基因上的读段分布情况 每 幅图中最上方显示的是基因的坐标 中间显示的是该基因各个剪接异构体中的外 显子分布情况 最下方显示了被定位到该基因的读段的分布情况 可以看到 读段在基因上面的分布呈现出明显的非均匀性 3 端 右 的读 段密度明显高于 5 端 左 可以想象 在估计基因表达值的时候将读段分布的非 均匀性也考虑进来理应会提升估计的精度 武征鹏等人于2010年发表了一个 基于读段非均匀分布假设模型 的算法 7 该算法在估计基因表达值的时候引入了能够表征读段非均匀分布的参数 使得模 型更加接近真实情形 从而提升了基因表达值估计的性能 7 12 3 3 1 在 URD 模型的基础上引入了表征非均匀分布的参数 武征鹏的模型是建立在 Hui Jiang 的均匀分布模型的基础之上的一个模型 通 过在均匀分布模型中引入两个能够表征非均匀分布的参数来提升估计性能 这两 个参数分别是 全局偏差曲线 GBC Global Bias Curve 和局部偏差曲线 LBC Local Bias Curve 这两个参数分别表征了所有基因中读段分布的总体特性和 某个特定基因中的读段分布特性 3 3 2 GBC Global Bias Curve GBC 是一条从基因的 5 端到 3 端的曲线 来自于对整个测序样本读段的统计 综合考虑了所有基因上面的读段分布情况 反映了整个测序样本中读段的分布情 况 本文采用归一化的统计直方图作为对 GBC 的估计 本文估计 GBC 曲线的具体流程如下 1 选择只有一个剪接异构体的基因 由于那些有多个剪接异构体的基因其读段的分布还受到该基因结构的影 响 因此本文在计算GBC曲线时只考虑那些只有一个剪接异构体的基因 2 筛去低表达基因 由于那些低表达的基因有较大的不确定性 因此我们滤去那些读段数少 于一定阈值 本文的实验选取的阈值为 100 的基因 3 统计读段数 通过前两步筛选之后得到符合估计 GBC 曲线要求的基因 利用读段定位 信息 获取这些基因上的读段分布情况 在本文的实验中 用统计直方 图来表示读段的分布情况 为了减少读段数的波动性 也为了便于综合 不同基因之间的读段分布信息 在统计直方图的时候需要先将个基因分 成等长的几个区间 区间数不可过多 在本文的实验中 各基因被等分 为 10 个区间 得到各基因的区间之后统计每个基因落在各个区间内的 读段数 通过统计所有的基因落在各个区间内的读段数的总和 可以绘 制成总的直方图如图 3 2 a 所示 4 归一化 为了便于后续的计算 本文进一步将 3 2 a 中的曲线进行归一化 使 得各区间的平均值为 1 如图 3 2 b 所示 图 3 2 b 即是本文的 GBC 曲线 13 a 读段分布直方图 b 偏差曲线 图3 2 偏差曲线示意图 7 3 3 3 LBC Local Bias Curve LBC 和 GBC 的不同在于它是针对某个特定的基因 所刻画的是某个基因自 身的读段分布情况 本文中的 LBC 曲线为一个阶梯函数 该阶梯函数在该基因的 每个外显子的范围内维持恒定值 该恒定值为 1 1 2 m jjij i xlajn 3 4 式 3 4 相当于将落在每个外显子中的读段数对该外显子在 URD 模型中的 泊松分布参数 见 3 2 节对 URD 模型的说明 进行归一化的结果 仅差一个常 数w w为所有的读段数 在某次测序中为常数 因此对读段分布情况没有影响 而对于服从泊松分布的随机变量而言 表示该变量的期望值 因此式 3 4 反 映了各外显子中实际的读段数与预期读段数之间的比值 又由于泊松分布是基于 读段均匀分布 假设的 因此式 3 4 也反映了实际读段的分布情况和均匀分 布假设的差距 注意到 该式中的 值在估计 LBC 曲线时是未知的 在第一次计算 LBC 时 可以令 12 1 1 1 m 来求解 LBC 求得各个 i 之后再代入式 3 4 中 更新 LBC 曲线的值 进而求得新的 i 如此往复 这是一个迭代的过程 14 3 3 4 GBC 矩阵和 LBC 矩阵 在 3 3 2 和 3 3 3 小节中得到的 GBC 曲线和 LBC 曲线均是一维的曲线 并不 能反映基因各剪接异构体的结构 而在实际的算法中 需要的是一个能够表示基 因结构的带权值的指示矩阵 图 3 3 给出了利用偏差曲线获取带权指示矩阵的过程 a 将区间分成三部分 b 对每个部分取平均值 图3 3 利用偏差曲线得到带权值的指示矩阵 7 以 GBC 矩阵为例 用第 3 3 2 节中的方法得到 GBC 曲线之后 图 3 3 a 实 线 为了得到某个基因的 GBC 矩阵 首先 选定该基因的某个剪接异构体作为 研究对象 将 GBC 曲线分割成若干段 其段数恰好为该剪接异构体的外显子数 且每段的长度之比恰好为各个外显子的长度之比 之后计算每段区间中 GBC 曲 线下方所围成的面积 并将该面积对该段区间的长度进行归一化 得到的结果即 为该剪接异构体的外显子在 GBC 矩阵中所对应的元素的值 对于那些不在该剪 接异构体中的外显子 只需令 GBC 矩阵对应的元素为 0 即可 对该基因的每个 剪接异构体都进行如上的计算之后即可得到最终的 GBC 矩阵 从 LBC 曲线获得 LBC 矩阵L有两种算法 第一种算法是由武征鹏提出的 其计算方法与本小节所述的 GBC 矩阵的计算过程类似 只需将 GBC 曲线替换成 LBC 曲线即可 第二种算法是由本文提出的一种新的算法 新的算法并没有将区 间分段然后求取各区间内 LBC 曲线所围面积的步骤 而是简单地对剪接异构体中 所包含的外显子的 LBC 值进行按比例放大 这个比例是整条 LBC 曲线下方所围 面积和该剪接异构体中的外显子的 LBC 曲线下方所谓面积的比值 新算法和武征鹏算法的示意图如下 15 a 武征鹏的算法 b 本文的新算法 图3 4 新算法和武征鹏提出的算法比较 图 3 4 中所示的基因一共有 4 个外显子 各外显子的长度为 25 20 25 30 其 LBC 曲线为 4 段 阶梯函数 所考察的剪接异构体有两个外显子 即该基因的前 两个外显子 2 段 阶梯函数的函数值对应于该剪接异构体在 LBC 矩阵中的非零 元素值 由于武征鹏的算法直接将剪接异构体长度延拓到整个基因的长度进行求解 而不是简单地依赖于所考察的剪接异构体所含的外显子读段分布情况 这会导致 计算某个外显子在 LBC 矩阵中的权值时 很大程度依赖于其他外显子内的读段密 度 而事实上这些外显子并没有为该外显子的读段密度提供任何贡献 如图 3 4 所示 LBC 曲线中第一个外显子的读段密度明显高于第二个外显子 但是利用武 征鹏的方法求解各外显子在 LBC 矩阵中的权值时 第二个外显子的权重与其在 LBC 曲线上的值没有直接的关系 而是直接地依赖于第三 四个外显子的读段密 度 这导致了在 LBC 矩阵中第二个外显子的权值反而大于第一个外显子的反常现 象 而新的方法由于是对各外显子读段密度按相同比例放大 因此能够很好地保 持外显子密度之间的大小关系 如果某基因的 LBC 曲线的各段函数值确实存在逆 序并且某剪接异构体的外显子分布恰好又满足类似于图 3 4 中的分布 那么用本 文的新算法得到的估计值理应会比武征鹏的算法更加合理 3 3 5 b 矩阵 在 3 1 节中 URD 模型用一个由 0 1 元素构成的指示矩阵a来代表基因的结 构 表示各个剪接异构体包含外显子的情况 由于a矩阵只能表示某个剪接异构 16 体是否含有某个外显子 而不能提供更多的信息 诸如读段的分布特性等 因此 URD 模型对基因表达值和剪接异构体表达值的估计受到了限制 在 N URD 模型中引入一个带权值的指示矩阵 b矩阵 b矩阵的计算需要用到 3 3 4 小节中得到的 GBC 矩阵G和 LBC 矩阵L 在本 文中 b矩阵是这两个矩阵的线性加权和 1 ijijij bGL 3 5 其中权值 的取值不定 是一个可以人为调节的量 当0 时 b矩阵特殊 化为L矩阵 当1 时 b矩阵特殊化为G矩阵 当0 5 时 b矩阵为L矩阵 和G矩阵的平均 b矩阵综合考虑了读段分布的全局偏差和局部偏差 能够较好的反应读段在 各基因上的真实分布情况 用其替代 URD 模型中的a矩阵 即可得到 N URD 模 型 3 4 似然函数 由 3 1 节可以知道 URD 模型的对数似然函数为 12 11111 log loglog nmnmn njijijjijij jijij Lx xxwl axl wax 其中的a为一个由 0 1 元素构成的指示矩阵 将其换成 N URD 模型中的b矩 阵 即可得到 N URD 模型的对数似然函数 12 11111 log loglog nmnmn njijijjijij jijij Lx xxwl bxl wbx 3 6 可以证明 将似然函数中的a矩阵替换为b矩阵不会影响似然函数的凸性 7 因此 N URD 模型中的优化问题与 URD 模型的优化问题拥有同样好的性质 17 第4章 N URD 模型的算法实现 4 1 程序概况 本程序的功能是利用基因注释和读段在基因上的定位信息来估计各基因的表 达值 表 4 1 是本程序的一些信息概要 表4 1 程序概况 版本 1 0 运行环境 Linux 使用语言 Perl C 文件形式 Perl 脚本 C 源码 输入 基因注释 支持 refFlat 格式 读段定 位文件 支持 eland2 格式和 SAM 格式 输出 各基因的表达值 程序的输入数据是基因注释和读段定位文件 本程序支持两种格式的读段定 位文件 eland2 和 SAM 由于本程序计算中所需的读段都是唯一定位读段 因此 若采用 SAM 作为读段定位文件 需要在利用读段定位工具生成 SAM 文件时指明 只选择唯一定位的读段 例如 利用bowtie做读段定位时 可以选择参数 m 1 由于 eland2 文件本身就包含了读段唯一定位的信息 因此在做读段定位时不需要 专门处理 程序的输出是基因表达值文件 该文件记录了基因的表达值和其他相 关信息 本程序主要由两部分构成 数据预处理部分和表达值估计部分 数据预处理 部分由 perl 语言实现 实现从输入数据中提取并整合信息的功能 表达值估计部 分由 C 语言实现 实现对单基因或者全基因组的表达值估计的功能 本程序由 perl 脚本文件和 C 源码构成 程序的入口是一个名为 NURD pl 的 perl 脚本文件 其语法是 perl NURD pl SAM eland2 map file gene anno file g gene name 其中各参数的含义是 18 1 SAM eland2 选择读段定位文件的格式 本程序支持 SAM 格式和 eland2 格式 2 map file 读段定位文件 3 gene anno file 基因注释文件 4 g gene name 可选项 如果用户在命令中使用 g 选项 则表示只对特定 的某个基因估计表达值 该基因名为紧跟 g 的参数 gene name 如果用 户未选该选项 默认对基因注释文件中的所有基因进行表达值估计 该命令的功能是 完成对读段定位信息和基因注释信息的数据预处理 C 文件的编译 对指定基因或全基因组的表达值估计 程序运行结束之后 将计算结果输出到文件 如果选择对特定基因求表达值 则将该基因的表达值信息输出到名为 singleGeneExpression expr 的表达值文件 如果对全基因组求表达值 则将所有基因的表达值信息输出到名为 allGeneExpression expr 的表达值文件 4 2 程序框架 程序的主要流程如图 4 1 所示 图4 1 程序主要流程 4 3 获取数据 本程序的输入数据是两个数据文件 一个是读段的定位 mapping 文件 另 一个是基因注释信息 本程序支持两种格式的定位文件 一种是 SAM 另一种是 eland2 由于本程序只处理唯一定位的读段 因此只支持读段唯一定位 SAM 格式 的文件 19 4 4 数据预处理 数据预处理的目的是将已经获得的数据进行一定的整合 提取出对后续计算 有用的信息 输入是之前获取的读段定位信息和基因注释信息 输出是经过整合 之后的一些简单统计信息以及各基因的基本信息 本程序的数据预处理部分主要 分为如下几个步骤 1 提取基因注释中所有基因的基本信息 包括 各个基因的名字 各基因 剪接异构体的数目和各剪接异构体的名字 各基因所有外显子的起始位 置 结束位置和长度 各基因的结构 即其各个剪接异构体中包含哪些 外显子等信息 2 从所有基因中提取出那些只有一个剪接异构体的基因 利用读段定位信 息获取定位到该基因上的读段的数量和起始位置 将各基因分成等长的 若干个区间 本文中为 10 个区间 统计被定位到各个区间内的读段数 得到各基因在各区间内的读段数之后 筛去那些读段数过少 本文中的 阈值为 100 个读段 的基因 以减少随机性 将这些基因在各个区间内 的读段数进行 归一化 使其呈均值为 1 的一个分布 并对各个基因分 布情况取平均 即可得到本次测序的全局偏差曲线 GBC 3 统计所有定位到基因上的读段总数 并且统计被定位到各个基因的各个 外显子内的读段数 4 将 GBC 曲线 所有被定位的读段总数 各基因的详细信息以及被定位到 各基因的各外显子内的读段数等信息保存到一个文件中 以便后续计算 20 图4 2 数据预处理部分流程 4 5 核心算法的实现 核心算法是本文工作的核心内容 是利用数据预处理得到的所有信息进行基 因和剪接异构体表达值估计的过程 本程序实现了 6 种算法来对表达值进行估计 核心方法都是利用斐波那契法求解优化问题 4 5 1 获取带权值的基因结构矩阵 求解 N URD 模型的先决条件是获取读段在基因上的分布信息 这一信息在 本文中用b矩阵来表示 由 3 3 5 小节可知 b矩阵由 GBC 矩阵G和 LBC 矩阵L 以及权重 决定 而利用 3 3 节中的方法 则可以直接计算得到各基因的 GBC 矩 阵G和 LBC 矩阵L 因此 在 GBC 矩阵和 LBC 矩阵确定之后 b矩阵只与权值 有关 不同的权值 会产生不同的b矩阵 进而对表达值的估计产生不同的影响 本文实现了 3 种对应于不同 值的算法 当0 时 b矩阵特殊化为L矩阵 对 应的算法为 LN URD 算法 当1 时 b矩阵特殊化为G矩阵 对应的算法为 GN URD 算法 当0 5 时 b矩阵为L矩阵和G矩阵的平均 对应的算法为 MN URD 算法 21 3 3 3 小节曾经提到 对 LBC 曲线的估计是一个迭代的过程 因此 不同的 迭代次数会导致 LBC 曲线的不同 从而造成 LBC 矩阵和b矩阵的不同 并影响 基因表达值的估计 本文实现了 2 种对应于不同的迭代次数的算法 1 M 和 5 M 分别表示在0 5 时进行 1 次迭代和 5 次迭代 此外 本程序还实现了基于读段均匀分布假设的 URD 算法 图4 3 获取带权值的基因结构矩阵 4 5 2 求解最大似然估计 求解 N URD 模型的核心是求解最大似然估计 本文中求解最大似然估计所 采用的方法是斐波那契法 22 斐波那契 Fibonacci 法 近似等价于 0 618 法 优选法 是求解优化 问题的一个十分有效的方法 能够用尽可能少的试探次数尽快找到最优解 其主 要思想是利用直线搜索中的区间压缩原理快速定位优化问题的最优解 由于本文的 N URD 模型的对数似然函数是凸函数 因此局部最优解就是全 局最优解 这一良好的性质可以保证用斐波那契法求解最优值的正确性 本文采用斐波那契法求解最大似然估计的具体流程如下 图4 4 用斐波那契法解最大似然估计的流程 23 其中涉及到两重循环 外循环 循环 1 用于控制算法收敛到最优的表达值 其循环停止条件是循环前后的似然函数值增加的幅度小于某个预先设定值或者 迭代次数超过某个预先设定的值 内循环 循环 2 用于控制斐波那契数列的区 间压缩过程 该循环的停止条件是区间长度小于某个预先设定的值 4 6 结果输出 利用 4 5 节中的算法求解最大似然估计 可以估计各基因的各个剪接异构体 表达值的最优解 求得的最优解被输出到本地文件进行保存 输出文件的数据格 式如下 以基因 EI24 为例 EI24 2 1568 NM 001007277 NM 004879 UR 166 529 218 165 43 2888 56 7112 GN 84 5251 296 07 22 2086 77 7914 MN 91 2242 289 722 23 9468 76 0532 LN 97 0417 284 209 25 4535 74 5465 1M 111 237 270 7 29 1244 70 8756 5M 117 429 264 814 30 7211 69 2789 第一行 基因名 该基因的剪接异构体数 被定位到该基因的总读段数 由 于有些基因表达值非常低 被定为到该基因内的读段数很少 此时对该基因的表 达值估计并不一定可靠 用户可以根据读段数来决定是否采纳本程序对该基因表 达值的估计结果 第二行 该基因的各个剪接异构体的名字 第三行至第八行 用 6 种算法得到的各个剪接异构体的表达值 每一行表示 一种算法的结果 以 EI214 的第 5 行为例 MN 91 2242 289 722 23 9468 76 0532 该行表示用 MN 算法得到的 EI 基因各剪接异构体的表达值 其中第一列是算 法的缩写 第2 3列分别表示EI24基因的两个剪接异构体的表达值 单位为RPKM 顺序对应于第二行的顺序 第 4 5 列表示两个剪接异构体表达值占基因表达总 量的百分比 24 6 种算法的缩写及其对应的含义如下所示 表4 2 算法缩写表 算法缩写 含义 UR URD 均匀分布模型 GN N URD 模型中 1 无迭代 MN N URD 模型中 0 5 无迭代 LN N URD 模型中 0 无迭代 1M N URD 模型中 0 5 迭代一次 5M N URD 模型中 0 5 迭代五次 25 4 7 程序总流程图 图4 5 程序总流程图 26 4 8 仿真实验 为了验证 N URD 模型对 URD 模型的改进 本文进行了仿真实验 仿真实验是指利用计算机生成符合模型要求 如读段分布的非均匀性 的伪 数据 并利用这些伪数据进行实验的过程 由于生成仿真数据的时候 实验者可 以事先得到真实的参数 对后续的实验结果形成比照 因此可以很容易地评估实 验结果的优劣 本文所进行的仿真实验的数据是 伪基因 的详细信息 即 人为生成一个 基因 并生成其相关信息 如基因结构 外显子长度 定位到各外显子内的读段 数等 这些信息是用随机采样和生成特定分布的随机数的方式得到的 此外 基 因中各剪接异构体的表达值是事先知道的 为了不同方法估计基因表达值的性能优劣 本文引入 MIRR 和 DS 这两个概 念来对各方法的性能进行定量的比较 MIRR 是 major isoform recovery rate 的缩写 指的是某种剪接异构体表达值估 计算法是否能够正确判断各基因中主要剪接异构体 major isoform 的能力 所 谓的主要剪接异构体 是指在该基因的各个剪接异构体中表达值最高的剪接异构 体 MIRR 指标越高 说明该算法的估计结果越接近真实情况 DS 是 difference score 的缩写 是对估计的表达值和真实表达值之间距离的一 个度量 当真实表达值和估计表达值分别为 12 m 和 12 m 时 其定 义式如下 1 11 100 m ii mm i jj jj DS 4 1 本实验中 利用所有基因的平均 DS 值做为评价一个算法性能的定量指标 DS 越小 说明算法的精度越高 本文的仿真实验的数据是 1000 个 伪基因 的数据 每个基因都有相同的剪 接异构体数和外显子数 不同的是各基因的基因结构 外显子长度和读段定位情 况 27 仿真实验的结果如下 表4 3 各方法仿真结果的 MIRR 值 单位 基因结构 URD GN URD MN URD LN URD 1 M 5 M 2 6 73 6 78 8 79 7 78 7 78 2 78 1 73 6 78 8 79 2 75 5 79 2 79 2 2 8 76 6 80 2 80 7 81 0 80 0 80 0 76 6 80 2 79 9 78 1 79 9 79 9 3 6 63 4 70 6 71 5 67 7 69 9 70 1 63 4