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

1、本文档下载自HYPERLINK /文档下载网，内容可能不完整，您可以点击以下网址继续阅读或下载：HYPERLINK /doc/be6b0724c324a6024ee96d8b/doc/be6b0724c324a6024ee96d8b新一代高通量 RNA 测序数据的处理与分析 _新一代高通量 RNA 测序数据的处理与分析 _生物化学与生物物理进展,37(8):834846此处图片未下载成功此处图片未下载成功新一代高通量RNA测序数据的处理与分析*王曦1)汪小我1)王立坤1,2)冯智星1)张学工1)*)(1)生物信息学教育部重点实验室，清华信息科学与技术国家实验室(筹)生物信息学研究部，清华大学自

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

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

4、a6024ee96d8b生物信息学方法与软件也应运而生本文针对当前RNA-seq应用的现实情况，尝试以Illumina/Solexa测序平台产生的mRNA-seq数据为例，对RNA测序数据的产生过程及数据处理和分析的基本流程、关键方法和现有软件进行较全面的介绍，并讨论RNA-seq数据分析中存在的挑战格更便宜、自动化程度更高的测序技术自2005年以来，以Roche公司的454技术、Illumina公司的Solexa技术和ABI公司的SOLiD技术为标志的新一代测序技术相继诞生8新一代测序技术又称作深度测序技术，主要特点是测序通量高、测序时间和成本显著下降9把这种高通量测序技术应用到由mRNA逆

5、转录生成的cDNA上，从而获得来自不同基因的mRNA片段在特定样本中的含量，这就是mRNA测序或mRNA-seq同样原理，各种类型的转录本都可以用深度测序技术进行高通量定量检测，统称作RNA-seq或RNA测序目前，在已经推出的几种新一代测序平台中，Illumina/Solexa测序平台上的RNA-seq应用最广，我们以此为例来综述RNA-seq数据处理和分析的生物信息学问题和方法.高通量测序技术简介*国家自然科学基金资助项目(60702002,60721003,30873464,60905013)和东南大学生物电子学国家重点实验室开放研究基金资助项目.*通讯联系人./p>

6、,E-mail:zhangxg收稿日期：2010-03-25，接受日期：2010-04-30诞生于20世纪70年代的Sanger法是最早被广泛应用的DNA测序技术7，也是完成人类基因组计划的基础但是，由于它测序通量低，费时费力，科学家们一直在寻求通量更高、速度更快、价;37(8)王曦等：新一代高通量RNA测序数据的处理与分析835/Solexa测序技术的基本原理是边合成边测序(sequencingb/doc/be6b0724c324a6024ee96d8bysynthesis，SBS)10-12，即测序过程是以DNA单链为模板，在生成互补链时，利用带荧光标记的dNTP发出不同颜色的荧光来确定不

7、同的碱基新加入dNTP的末端被可逆的保护基团封闭，既保证单次反应只能加入一个碱基，又能在该碱基读取完毕后，将保护基团除去，使得下一个反应可继续进行为了增加荧光强度，使之更易被成像系统所采集，该技术在测序之前还需要对待测片段做桥式扩增(bridgeamplification)13(/)初期的Illumina/Solexa测序技术只能在较短的测序读长上(2030碱基)保证较高的正确率随着技术的改进，目前的读长已经增加到100碱基以上同时，随着双端测序(paired-end，PE)技术的成熟，测序长度更可达到单端测序的2倍，测序通量也随之增加这种测序技术是Solexa公司发展起来的，2007年被Il

8、lumina公司收购，因此现在通常被称为Illumina/Solexa测序技术近两年来，Illumina/Solexa测序平台不断升级，相继推出了GA(GenomeAnalyzer)、GAIIx、HiSeq2000等测序仪更多关于高通量测序平台的介绍，可以查阅相关文献9,14-16RNA鄄seq测序文库制备和测序平台数据输出本小节针对Illumina/Solexa测序平台，对RNA测序文库制备标准和平台底层数据产生做一个简单的介绍援1RNA鄄seq测序文库制备对于mRNA-seq实验，从总RNA到最终的cDNA文库制备完成主要包括以下步骤首先，用Poly(T)寡聚核苷酸从总RNA中抽取全部带P

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

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

11、的RNA-seq样品制备方法17-20援2测序平台数据输出将RNA-seq测序文库加入流动槽(flowcell)中的各通道(lane)，在桥式PCR扩增后，就可以进行测序了测序过程中，计算机软件同步地对荧光图像数据进行处理，通过分析荧光信号来确定被测碱基，并给出质量评分按照图像上的位置坐标，计算机程序将同一位置测得的碱基根据测序顺序连成读段(read)由于荧光图像文件所占有的磁盘空间很大，通常GAIIx平台一次实验就能产生上太字节(TB)的图像文件，所以一般情况下不予保留原始的荧光图像数据，而是只保留程序读出的读段数据及对应的质量分值，这就是多数实验室委托测序中心进行RNA-seq测序后得到的

12、最原始的数据为了便于测序数据的发布和共享，高通量测序数据以FASTQ格式来记录所测的碱基读段和质量分数如图1所示，FASTQ格式以测序读段为单位存储，每条读段占4行，其中第1行和第3行由文件识别标志和读段名(ID)组成(第1行以“”开头而第3行以“ ”开头；第3行中ID可以省略，但“ ”不能省略)，第2行为碱基序列，第4行为对应的测序质量分数关于FASTQ格式更多地介绍可参考文献21为方便保存和共享各实验室产生/doc/be6b0724c324a6024ee96d8b的高通量测序数据，NCBI、EBI、DDBJ等数据中心建立了大容量的数据库SRA(SequenceReadArchive,/Tr

13、aces/sra)来存放共享的测序数据22-23836生物化学与生物物理进展Prog.Biochem.Biophys.2010;37(8)每4行标识为一个测序读段读段识别码碱基序列此处图片未下载成功此处图片未下载成功此处图片未下载成功 读段识别码测序质量分数鄄seq数据的基本处理seq的基本应用是测量一个样本中的基因表达或转录组有实验表明，新一代高通量测序技术重复数据之间的相关度较高(R2抑0.96)1-2，因此，如果对同一样本在多个通道上进行了RNA测序的技术重复，我们建议可以把几个通道的数据进行合并，这样等效地增加了测序深度本节讨论单个样本RNA测序数据的基本处理流程，如图2a所示图1读段

14、FASTQ数据格式示例Fig.1FASTQformatexamples(a)(b)图2所示的流程,虚线箭头表示可选输入.鄄seq数据处理和分析流程图鄄seqdataprocessingandanalysis(a)RNA-seq数据的基本处理,其方法介绍见正文第3节.(b)两类样本RNA-seq数据比较分析的框架,对应于正文的第4节.(b)中虚线框内为(a)援1读段定位获得RNA-seq的原始数据后，首先需要将所有测序读段通过序列映射(mapping)定位到参考基因组上，这是所有后续处理和分析的基础在读段定位之前，有时还需要根据测序数据情况对其做某些基本的预处理例如，过滤掉测序质量较差的读/do

15、c/be6b0724c324a6024ee96d8b段、对miRNA测序读段数据去除接头序列等高通量测序的海量数据对计算机算法的运行时间提出了很高的要求针对诸如Illumina/Solexa等测序平台得到的读段一般较短、且插入删除错误较;37(8)王曦等：新一代高通量RNA测序数据的处理与分析837少等特点，人们开发了一些短序列定位算法这些算法主要采用空位种子索引法(spaced-seedindexing)或Burrows-Wheeler转换(Burrows-WheelerTransform，BWT)技术来实现24空位种子索引法首先将读段切分，并选取其中一段或几段作为种子建立搜索索引，再通过查

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

17、发了一些基于改进的Smith-Waterman动态规划算法29的序列比对工具，如BFAST30、SHRiMP31、Mosaik(/marthlab/Mosaik)等，但算法速度较慢，大多需采用计算机并行编程技术来解决运行时间的问题表1Table1名称MAQ26Bowtie27BWA32ZOOM33ELANDSOAP234RazerS35NovoalignSHRiMP31BFAST30Mosaik)否是是否否/doc/be6b0724c324a6024ee96d8b否否是否质量2)是是是否否否否是是适用于Illumina/Solexa测序平台的读段定位软件Mappers/alignersforI

18、llumina/Solexasequencingdata主要采用技术空位种子BWTBWT空位种子空位种子BWTq-grams过滤Needleman-Wunsch算法空位种子q-grams过滤Smith-Waterman算法:/shrimp:/index.php/BFAST/marthlab/Mosaik网址:/:/index.shtml/bwa.shtml:/products/zoom:/software/genome_analyzer_software.ilmn:/www.seqan.de/projects/razers.html是是是是Waterman算法并行编程Smith-Waterma

19、n算法并行编程):是否能以SAM格式输出;2)质量:是否提供读段定位质量信息;BWT:Burrows-Wheeler转换.在RNA测序数据的基因组定位中，一个特殊的问题是跨越两个外显子接合区的读段(junctionreads)定位在真/doc/be6b0724c324a6024ee96d8b核生物中，成熟的mRNA是经过由mRNA前体中的外显子经过剪接形成的如果一个读段跨越了两个外显子，那么就无法将这个读段完整地定位到基因组序列上而同时，这种跨两个外显子的读段在分析转录本的剪接形式和研究选择性剪接中有重要的作用为了解决这一问题，人们采取两种典型的策略来进行接合区读段的定位：一是根据已知的基因外

20、显子注释，构建所有可能的外显子接合区序列，与基因组序列一并作为定位的参考基因组；二是不依赖基因注释，而是先利用能完整定位到基因组的读段得到粗略的外显子区域，并结合剪接位点序列构建出可能的剪接位点，然后将不能完整定位的读段分段定位到两个外显子可能的结合区域Illumina/Solexa平台提供的RNA-seq软件分析包GApipeline采用了第一种策略采用第二种策略的软件有Tophat36和G-Mo.R-Se37等，最新的Tophat软件增加了利用已知外显子边界注释信息的选项838生物化学与生物物理进展不论是哪种测序平台，测序中都不可避免地存在一定的错误，基因组中又存在单核苷酸多态性等引起的序

21、列变化，所以在读段定位时通常允许一定数量的错配，可以根据不同应用调节允许错配的程度另一方面，由于基因组中重复序列和高相似度序列的影响，某些读段会出现定位到基因组多个位置的情况这些因素影响了各个读段到基因组的定位质量，在一些新的读段定位算法中，同时给出每个读段与基因组匹配质量通常在后续处理前，人们将多定位的读段都过滤掉，也有人尝试用适当的策略把多定位读段“分配”到其中某些位置上2,38.读段定位到基因组后通常采用SAM(SequenceAlignment/Map)格式或其二进制版本BAM格式39来存储二进制版本可大大节省存储空间，但不能直接用普通文本编辑工具显示关于SAM格式的详细介绍，可查阅(

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

23、与基因真实表达水平成正比，还与基因长度成正比，同时也与测序深度即测序实验中得到的总读段数正相关为了保持对不同基因和不同实验间估计的基因表达值的可比性，人们提出了RPM和RPKM的概念2RPM(readspermillionreads)即每百万读段中来自于某基因的读段数，考虑了测序深度对读段计数的影响RPKM(readsperkilobasespermillionreads)是每百万读段中来自于某基因每千碱基长度的读段数，公式表示为：=基因区读段计数伊伊109不仅对测序深度作了归一化，而且对基.Biochem.Biophys.2010;37(8)因长度也作了归一化，使得不同长度的基因在不同测序深

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

25、率与序列的GC含量等特征相关，可导致GC含量高的cDNA片段在文库中拷贝数增加超过其他片段；舍弃多定位的读段也可能导致读段的非均匀分布；等等如果能对读段分/doc/be6b0724c324a6024ee96d8b布的不均匀性进行建模并加以校正，可以提高RNA-seq推断基因表达量的准确度但根据对实际数据的观察，对于较长转录本，读段非均匀分布带来的误差很大程度上可相互抵消，用RPKM来估计基因的表达水平可以得到比较满意的结果援3选择性剪接事件识别和剪接异构体表达水平推断在真核生物中，选择性剪接现象普遍存在基因转录形成的mRNA前体(pre-mRNA)在剪接过程中因去掉不同的内含子区域或保留不同的

26、外显子区域，可形成不同的剪接异构体根据RNA-seq原理，只要测序深度足够深，就能检测到所有转录本的全部序列，包括来自剪接接合区的序列利用考虑到接合区的读段定位方法，就有可能系统地研究某一组织或某一条件下的基因选择性剪接事件前面已经提到，Tophat等软件定位剪接接合区读段的策略能标定出剪接事件中的两个剪接位点：供体位点和受体位点通过比较供体位点和受体位点的组合，就能识别选择性剪接事件4,49图3中包含了选择性剪接识别的一个例子进一步，通过对供体和受体位点的读段计数，结合外显子其他区域的读段数据，还能定量地计算选择性剪接事件之间的比例50-51对于每一个剪接异构体，RNA-seq数据能在一;3

27、7(8)王曦等：新一代高通量RNA测序数据的处理与分析839定程度上推断其表达水平比如，可以根据已知外显子组成和各外显子长度对剪接异构体建立数学模型，在测序读段转录本上均匀分布的假设下，利用各外显子上的读段数和接合区读段数求解异构体的表达值Jiang等43的方法及软件IsoInfer52和cufflinks45都采用了这种思路来实现剪接异构体的表达推断需要指出的是，某些形式的剪接异构体表达水平在这种方法框架下不可辨识533援4新基因的检测在对RNA-seq数据的分析中，人们发现，往往不是所有读段都能定位到已有注释的基因区，说明除了转录噪声或测序错误等的影响外，可能还存在尚未被注释的基因这里，我

28、们把这种尚未注释的基因称为新基因，包括新的蛋白质编码基因和非编码/doc/be6b0724c324a6024ee96d8bRNA基因能检测新基因，尤其是低表达基因是RNA-seq技术优于基因芯片的特点之一，因为它不需要利用已知基因注释来设计检测探针seq技术灵敏度高，但样品污染、测序错误等仍可能带来背景噪声从基因组未注释区域的RNA测序读段信号中检测新基因是典型的信号检测问题如何控制新基因识别的误发现率(FDR)是检测方法的关键Useq软件包54将ChIP-seq数据分析的方法移植到RNA-seq数据上，用滑窗的方法来识别测序读段定位富集的区域，给出反映滑窗所在区域读段富集显著程度的P值(P-

29、value)及新基因误发现率，通过设定P值或误发现率的阈值，可筛选出读段富集的区域，再将相邻区域合并或根据剪接接合区读段将相应区域连接，完成新基因的检测援5读段的可视化及注释对于复杂的组学数据，能尽可能方便地直接观察数据对于数据的分析和解释都非常重要，对新一代测序数据的可视化和交互展示是一个非常重要但容易被人忽视的问题不深入考查数据的细节，而是满足于对数据的统计分析，是高通量数据应用中经常容易陷入的误区，方便有效的可视化工具能够帮助避免这样的误区表2列出了部分适用于RNA-seq数据的全基因组浏览器，其中比较具有代表性的有UCSCGenomeBrowser、CisGenomeBrowser和I

30、GV(IntegrativeGenomicsViewer)等这些浏览器具有如下特点：a能在不同尺度下显现单个或多个读段在基因组上的位置，包括来源于剪接接合区的读段；b能在不同尺度下显示不同区域的读段丰度，以反映不同区域的转录水平或测序效率；c能显示基因及其剪接异构体的注释信息；d能显示其他注释信息，例如物种间基因组序列保守性、序列GC含量等；e能直接或间接支持SAM/BAM读段定位数据存储格式.UCSCGenomeBrowser55属于基于网络模式的全基因组浏览器，所有数据都需要上载到远程服务器，经过处理后将图形返回客户端显示图3中的例子就是从UCSCGenomeBrowser的显示截取的.C

31、isGenomeBrowser56是典型的本地版基因组浏览器，所有读段数据、注释信息都存于本地文件，因此不/doc/be6b0724c324a6024ee96d8b需要网络连接，方便内部考查数据用IGV(/igv)可以说是以上两种模式的融合，既可以从远程服务器端下载各种注释信息，又可以从本地加载注释信息表2名称IGV适用于mRNA鄄seq数据的全基因组浏览器/viewersapplicabletomRNA鄄seqdataviewing支持的数据格式网址:/igv/GFF3,BED,SAM/BAM,WIG,55BED,bigBed,BEDGRAPH,GFF,GTF,bigWig,BAM,5657

32、,BED,refFlat,FA,Wig,BED,GFF,FASTA,ELAND,GFF,BED,MAQ,SAMBAM,BED,GFF2,GFF3,FASTA,VCFExpressiondata,Annotationtracks:/jiangh/browser/sj/mochiview-start:/www.bioinformatics.bbsrc.ac.uk/projects/seqmonk/p/gambit-viewer:/packages/release/bioc/html/GenomeGraphs.html以上列出的全基因组浏览器均可/doc/be6b0724c324a6024ee96d

33、8b在Windows、Linux和苹果公司的MacOS等计算机平台下运行.840生物化学与生物物理进展Prog.Biochem.Biophys.2010;37(8)剪接接合区样本A测序标签分布该基因在A、B样本中差异表达内含子区的测序标签接测序标布剪接接合区样本B测序标签分布基因注释图3鄄seq数据可视化示例鄄seqdatavisualization图示区域为人类基因CBX7.图中红色表示样本A的数据，蓝色表示样本B.各轨道(track)依次为:基因组坐标、样本A的剪接接合区、样本A的读段分布、样本B的剪接接合区、样本B的读段分布、UCSC基因注释.图中还标识了：因受体位点不同而形成的选择性剪

34、接;基因的5忆端出现读段的非均匀分布;在两个样本中,差异表达基因的读段信号强度不同;在基因标注的内含子(intron)区域存在少量不连续的读段.除对读段的可视化外，用描述统计学方法对实验数据进行分类别统计也十分重要例如，统计读段在各个染色体上的分布情况和在注释的外显子、内含子、剪接接合区、基因间区的分布情况等目前，已经有一些用于测序数据注释的生物信息学软件，比如SAMtools39、BEDtools58等，但由于测序技术发展迅速，用户需求因人而异，用户经常还需要根据需求编写一定的程序或脚本完成或完善注释分析的任务对于熟悉图形用户界面的研究人员，还可以利用UCSCTableBrowse56和Ga

35、laxy59-60来配合完成注释分析由于UCSCTableBrowser集成了大量基因组尺度上的注释信息，而Galaxy又为用户提供了书写简单、接口明晰和直观的数据处理流程，这种方法十分方便有效，也为很多学者在展示研究成/doc/be6b0724c324a6024ee96d8b果时所采用以上描述了对于基因组已知的物种进行RNA-seq数据处理基本流程，图2a给出了其主要步骤若研究对象尚未完成基因组测序，则需要采用读段的从头拼装(denovoassembly)6,61-62来代替读段定位，后续流程也须做相应的调整若RNA-seq实验在文库制备时保留了RNA的方向信息，则应分别研究来自正链和反链的

36、转录产物，并通过与基因注释比较来检测反义转录本63最近，RNA-MATE64软件在其分析流程中加入如此的处理策略此外，通过分析定位到外显子接合区的读段，还可以获取转录本结构，这为研究基因的剪接调控机理提供了重要信息5而利用RNA-seq数据提供的序列信息，通过与DNA序列的细致比较可分析转录组的序列差异(如SNP等)65，从而研究等位基因的表达模式66-67及RNA编辑68等最后需要指出，由于miRNA在序列和结构上具有一定的特点，miRNA-seq数据的基本处理流程也与本节所述有所不同，感兴趣的读者可参考软件工具miRDeep69提供的处理策略多类样本mRNA鄄seq数据间的比较分析很多RN

37、A-seq实验的目的是为了比较两种或多种样本中基因表达或整个转录组的差异，如比较癌症组织和正常组织的转录组差异等这些差异既包括通常意义下的差异表达基因，也主要包括选择性剪接模式的差异、剪接异构体表达的差异、非编码转录本的差异等这些差异一般可以用一些统计假设检验方法检测，但这种检验有时会受到测序深度、基因长度等因素的影响70-71，需要对结果进行仔细分析，消除可能的混杂因素，必要时可以用读段的绝对表达值倍数变化(fold-change)来作为补充图2b给出了两类样本数据分析的框架;37(8)王曦等：新一代高通量RNA测序数据的处理与分析841援1差异表达基因的识别虽然新一代测序相对第一代测序的单

38、位成本大大降低，但是，利用RNA测序进行基因表达研究的成本仍很高，因此，很多实验室没有条件进行样本重复如果两类样本均没有生物重复，例如只/doc/be6b0724c324a6024ee96d8b对两个细胞系各进行一次mRNA样本测序，则可以用随机采样模型通过假设检验来分析差异表达对于某个基因，如果一个读段来自于这个基因，我们称事件A发生对于一次RNA-seq实验，事件A发生的概率可以用这个基因上的读段数n除以所有基因上的读段总数N来估计，即RPM事件A发生的概率反应了这个基因的表达水平如果要判断(a)某个基因在两个样本中的表达水平是否一致，就可以通过检验事件A在两种条件下发生的概率是否一致来实

39、现，采用似然比检验1、Fisher精确检验72以及基于MA图的统计检验方法44等同样，也可用RPKM作为统计量来进行假设检验分析，由于是比较同一个基因在两个样本间的差异，基因长度的影响被抵消，用RPKM和用RPM得到的结果相似对无生物重复的RNA-seq数据进行差异表达基因分析，已经有几个公开发表的软件，包括DEGseq44、Useq54、Cufflinks45中的Cuffdiff模块等图4展示了我们开发的DEGseq软件提供的多种差异表达基因识别方法的应用例子(b)-A(d)A-0A8log2(readcountsforeachgene)inB(c)-AvsBAvsB图4用DEGseq软件包

40、识别差异表达基因的结果(a)各基因在样本A和样本B中表达水平的散点图.(b),(c),(d)图中红点表示分别用FET、LRT和MARS方法得到的差异表达基因.FET:FishersExactTest,Fisher精确检验.LRT:LikelihoodRatioTest,似然比检验.MARS:MA-plot-basedmethodwit/doc/be6b0724c324a6024ee96d8bhRandomSamplingmodel,基于MA图的随机采样模型.如果每一类样本都包含了若干生物重复，如病人和正常人对照研究，则可以沿用基因芯片数据分析中的很多方法比如，可以用t检验结合倍数变化的方法来分

41、析差异表达如果两类样本具有配对的信息，也可以通过整合每对样本分析结果来实现其步骤为，先在每对样本中识别出差异表达的基因，再寻找这若干组差异表达基因之间的相同者，或用投票的方法来为基因的差异程度打分针对某些RNA-seq数据生物样本量小，R软件包DEGseq44和edgeR73等还专门提供了基于改进模型的统计方法此外，一类将分类器与特征选择包裹在一起的方法也同样适用于此类问题(见4.3)842生物化学与生物物理进展4援2差异表达剪接异构体的识别差异表达剪接异构体的识别方法与差异表达基因的识别相似如果把剪接异构体看成是独立的基因，那么前面讨论的用于识别差异表达基因的方法对剪接异构体完全适用但是，注

42、意到来自于同一个基因的剪接异构体并不独立，某些假设检验的基本条件并不满足，得到的结果就不一定正确此外，由于现在剪接异构体表达推断的方法还不够成熟，加之在基因结构不可辨识的剪接异构体上作表达推断会出现病态结果，差异表达剪接异构体的识别问题还处于探索的阶段目前，在剪接异构体表达水平可辨识且读段覆盖度较高的基因上，BASIS74方法通过贝叶斯模型来推断差异表达的剪接异构体换一个角度，剪接异构体由选择性剪接造成，如果剪接异构体的表达有差异，那么导致这些异构体的选择性剪接事件及异构体特异的外显子的表达也会有差异因此，对差异表达剪接异构体的识别可以转变为分析选择性剪接事件和外显子表达的差异75外显子表达差

43、异的分析可以完全利用基因表达差异的分析方法而剪接接合区也可以看成是一个较短的“外显子”(长度一般与测序长度相当)不过，由于外显子长度较基因的长度短，对应的读段数量较少，差异识别的敏感度会有所下降Solas方法75就是根据类似的原理，采用统计学假设检验的方法来识别差异表达的剪接异构体4援3对样本的分类分析:/doc/be6b0724c324a6024ee96d8b通过统计方法识别出来的差异表达基因及剪接异构体能否有效地区别两类样本，可以通过分类分析进一步证实如果把每个基因(或剪接异构体)的表达值作为特征，则差异表达基因(或剪接异构体)的选取也就是特征筛选的过程把前面用统计方法等检测出来的差异表达

44、基因(或剪接异构体)用于分类分析，常被称为过滤法另一类基于分类器的包裹法，例如R-SVM76、SVM-RFE77等，可以根据每个特征在分类器中所占的权重来筛选特征，因此也可以用于差异表达基因(或剪接异构体)的识别分类的性能可以用交叉验证(cross-validation，CV)方法来评估需要特别注意的是，交叉验证应该包括对特征选择步骤的交叉验证，防止发生信息泄露而导致评估结果过于乐观具体做法是：将样本按一定的策略分成两份，一份(通常是样本数多的一份)用于特征选取和分类器训练，而用余下的样本进行分类器性能的估计；重复以上步骤多次，就得.Biochem.Biophys.2010;37(8)到交叉验

45、证错误率必要时还可以用随机置换检验(permutationtest)来推断所得错误率的统计显著性76当样本数较小时，可以采用留一法交叉验证(leave-one-outcross-validation，LOOCV)4援4其他高层分析方法检测差异表达的基因或差异表达异构体是人们认识所研究的生物问题机理的第一步，接下来需要从功能上研究这些差异转录现象的分子机理这与在基因芯片应用中所面临的是同样的生物学问题，对芯片数据分析结果的后续处理方法，都可以借鉴到测序数据上来如何进一步地从机理来解释结果，还需结合已知生物学知识进行后续分析人们对基因芯片得到的基因表达数据进行分析的很多方法都可以用到RNA-seq

46、数据上来，比如利用机器学习方法进行分类和特征选择，对差异表达的基因进行GO(geneontology)78类别富集分析、信号通路富集分析等，一些常用的分析工具包括GoMiner79、DAVID80和VisANT81等需要说明的是，在各种以差异表达基因为基础的分析中，由于基因表达水平都是通过读段计数来估计的，表达水平较高或/doc/be6b0724c324a6024ee96d8b转录本较长的基因拥有更多的读段，更容易被多数统计方法识别为差异表达基因70这种偏好可能对后续分析带来影响以GO类别富集分析为例，这种偏好将导致长基因占主导的功能类别更有可能被识别为富集的功能这将对生物机理的研究带来误导最

47、近，Young等71发展了一种GOseq方法，针对这一偏好对GO类别富集分析做了改进RNA鄄seq数据处理中的生物信息学挑战高通量测序技术的发展十分迅速，这要求相应的数据处理与分析方法快速跟进正是这些方法，架起了高通量实验数据与科学问题之间的桥梁这种桥梁作用正日趋重要，也为生物信息学带来了挑战7,71这里，我们重点讨论两方面的挑战：a如何实现剪接接合区读段的准确定位？b在数据处理各阶段中，如何对RNA-seq数据的系统误差和固有偏好建模或补偿，以消除它们可能带来的错误推断及结论？援1剪接接合区读段的定位测序技术的一个发展趋势是测序长度不断增加随着读长的增加，RNA-seq中来自剪接接合区的读段

48、会越来越多我们粗略估算，按照人类基因组refSeq基因注释，一般情况下，如果测序读长;37(8)王曦等：新一代高通量RNA测序数据的处理与分析843为50个碱基，则约有10%的读段来自剪接接合区而当测序长度达到100个碱基时，这个比例将达到25%左右对这些剪接接合区读段的分析，将使我们能够更准确地检测剪接事件和推断剪接异构体的表达水平，大大推进人们对选择性剪接的研究在RNA-seq出现的早期，人们没有意识到剪接接合区读段的重要性因为当时的读长只有2030个碱基，来自剪接接合区的读段所占比例甚小当时读段定位的通常做法是，先将读段与全基因组序列做映射定位，再考虑不能定位的读段是否来自于剪接接合区2

49、这种做法虽然在一定程度上保证了读段定位的比率，但由于基因组中重复序列和相似序列的存在，部分接合区读段有可能在容许错配的情况下被定位到基因组上其他位置，从而失去了定位到正确的剪接接合区的机会在读段定位时，如果要同时/doc/be6b0724c324a6024ee96d8b考虑基因组序列和剪接接合区序列，就要利用已知的剪接事件注释，这是目前软件通用的方法然而，包括人类在内的各物种的基因注释信息都还有待完善，也没有较完整的剪接组(splicome)数据库，能够不依赖注释信息和对剪接机理的现有认识，高效、准确地定位所有已知和未知的接合区读段，仍然是对读段映射定位算法的一个挑战援2系统噪声和偏好的分析虽

50、然深度测序技术的准确性较以前的技术有了很大提高，但仍然存在错误和噪声比如从图3中可以看到，内含子区内有一些不连续的读段，很可能由系统噪声造成，如样品污染、测序错误和不恰当的读段定位策略等从图3还能看出，外显子区域内的读段信号分布也很不均匀有文献报道，序列组成尤其是GC含量46、RNA二级结构2等也有可能是导致读段不均匀分布的原因这些噪声和分布偏好将影响新基因的识别和对剪接异构体形式和表达水平的推断合理地建模RNA-seq数据中的系统噪声和偏好是解决上述问题最有效的办法基本的思路可以是：首先根据实验原理寻找可能产生系统噪声或偏差的因素，并尽可能将这些因素转化成可量化的特征，如序列特征、二级结构等

51、；然后，将用实验数据对这些特征做统计分析，构造和训练模型，用模型来对数据进行校正需要注意的是，某些偏好是由当前的测序技术和分析方法共同造成的，难以完全消除71在这种情况下，后续处理和解释时需要充分意识到这种偏好可能对生物学结论带来的影响，必要时通过补充其他实验来验证和修正通过高通量测序得到的生物学结论总结与展望本文以Illumina/Solexa测序平台为例，尝试对新一代测序技术的RNA-seq数据处理和分析方法做了较为全面的梳理，并对各个环节上可用的软件进行了汇总高通量测序是正在飞速发展的技术，相应的生物信息学方法也在快速发展，这里讨论的是RNA-seq中一些代表性的方法和问题，希望能对正在

52、或即将采用RNA-seq实验进行科学研究的学者和进行RNA测序数据处理的同行提供参考.测序和基因芯片有很多共同的应用领域，尽管相对还不是很成熟，RNA-seq技术在很多方面已经表现/doc/be6b0724c324a6024ee96d8b出了优势，有人甚至预言基因芯片时代即将结束36但也有报道认为，RNA-seq数据在基因表达水平的估计上和基因芯片相比没有明显的优势60，加上测序的成本目前还远高于芯片实验的成本，所以更多人认为测序和基因芯片将长期共存，以各自不同的特点在现代组学研究中发挥作用新一代高通量测序技术的应用面非常广82，RNA-seq只是其中一个方面，除此之外，基因组的从头测序和重测

53、序83-84、染色质免疫沉淀测序(ChIP-seq)85-86、甲基化测序(Methyl-seq)87-88等技术都同样有着广泛的应用尤其是，用ChIP-seq研究蛋白质与DNA的相互作用，能够得到高分辨率的转录因子结合数据和组蛋白修饰等表观遗传学数据发展有效的生物信息学方法，将ChIP-seq数据与RNA-seq得到的转录组数据进行综合分析，将大大推进人们对复杂的基因转录调控系统的认识致谢感谢本实验室刘霖曦、谢芃、孟璐等同学对本工作有意义的讨论，感谢斯坦福大学WingHWong教授、HuiJiang博士和JunLi同学等的讨论和帮助参考文献1MarioniJC,MasonCE,ManeSM,

54、etal.RNA-seq:anassessment.GenomeRes,2008,18(9):1509-15172MortazaviA,WilliamsBA,McCueK,etal.MappingandSeq.NatMethods,2008,5(7):621-628844生物化学与生物物理进展3NagalakshmiU,WangZ,WaernK,etal.Thetranscriptional:/doc/be6b0724c324a6024ee96d8bing.Science,2008,320(5881):1344-13494SultanM,SchulzMH,RichardH,etal.Aglob

55、alviewofgene.Science,2008,321(5891):956-9605WangET,SandbergR,LuoS,etal.Alternativeisoformregulation.Nature,2008,456(7221):470-4766BirzeleF,SchaubJ,RustW,etal.Intotheunknown:expression.NucleicAcidsRes,2010,doi:10.1093/nar/gkq1167SangerF,NicklenS,CoulsonAR.DNAsequencingwithchain-terminatinginhibitors.

56、ProcNatlAcadSciUSA,1977,74(12):5463-54678MarguliesM,EgholmM,AltmanWE,etal.Genomesequencingdensitypicolitrereactors.Nature,2005,437(7057):376-3809ShendureJ,JiH.Next-generationDNAsequencing.NatBiotechnol,26(10):1135-114510RuparelH,BiL,LiZ,etal.Designandsynthesisofa3忆-O-ally/doc/be6b0724c324a6024ee96d8

57、bl.ProcNatlAcadSciUSA,2005,102(17):5932-593711SeoTS,BaiX,KimDH,etal.Four-colorDNAsequencingby.ProcNatlAcadSciUSA,2005,102(17):5926-593112JuJ,KimDH,BiL,etal.Four-colorDNAsequencingbysynthesis.ProcNatlAcadSciUSA,2006,103(52):19635-1964013FedurcoM,RomieuA,WilliamsS,etal.BTA,anovelreagentforphaseamplifi

58、edDNAcolonies.NucleicAcidsRes,2006,34(3):e2214ShendureJA,PorrecaGJ,ChurchGM.OverviewofDNA/AusubelFM,BrentR,KingstonRE,A:JohnWileyandSons,Inc.,2008:Unit7.115MardisER.Next-generationDNAsequencingmethods.AnnuRev,2008,9:387-40216Full/doc/be6b0724c324a6024ee96d8berCW,MiddendorfLR,BennerSA,etal.Thechallen

59、gesof.NatBiotechnol,2009,27(11):1013-102317CroucherNJ,FookesMC,PerkinsTT,etal.Asimplemethodfor.NucleicAcidsRes,2009,37(22):e14818ParkhomchukD,BorodinaT,AmstislavskiyV,etal.TranscriptomespecificsequencingofcomplementaryDNA.NucleicAcidsRes,2009,37(18):e12319PerkinsTT,KingsleyRA,FookesMC,etal.Astrand-s

60、pecificSeqanalysisofthetranscriptomeofthetyphoidbacillusSalmonellatyphi.PLoSGenet,2009,5(7):e100056920MamanovaL,AndrewsRM,JamesKD,etal.FRT-seq:free,strand-specifictranscriptomesequencing.Nat.Biochem.Biophys.;37(8),2010,7(2):130-13221CockPJ,FieldsCJ,GotoN,etal.TheSangerFASTQfileformat,andtheSolexa/Il