转录组测序数据分析流程

转录组测序（RNA-seq）是指利用高通量测序技术对cDNA进行测序从而获取不同mRNA片段在特定样本里的含量。通过统计相关读段(reads)数计算出不同RNA的表达量，发现新的转录本；如果有基因组参考序列，可以把转录本映射回基因组，确定转录本位置、剪切情况等更为全面的遗传信息。

RNA-seq广泛应用于生物学研究、医学研究、临床研究和药物研发等。在功能基因组学方面的应用：基因（转录本）精细结构，各种器官表达模式，不同发育阶段的表达模式；在癌症和疾病方面：SNPs，选择性剪切，可选启动子，等位基因差异表达，RNA编辑，融合基因；农业和微生物方面：动植物疾病，微生物致病机理，分子育种，作物抗性研究；等等。

服务内容

RNA测序

客户提供

足够的组织材料或者RNA样品。

我们提供

测序结果、报告。

转录组分析流程

转录组学是研究生物体在特定时期或环境中所有mRNA的总体，通过对转录组的分析可以揭示基因的表达模式、功能及调控机制，对于疾病的诊断、治疗及基因工程等领域具有重要意义。转录组分析流程是对转录组数据进行处理和解读的过程，下面将介绍转录组分析的一般流程。

1. 数据获取。

转录组数据通常来源于高通量测序技术，包括RNA-seq和单细胞RNA-seq等。首先需要进行样本采集和RNA提取，然后进行测序实验，获取原始的测序数据。

2. 数据预处理。

原始测序数据通常包含低质量序列、接头序列和重复序列等，需要进行质量控制和过滤，包括去除接头序列、低质量序列和重复序列，得到高质量的测序数据。

3. 序列比对。

将高质量的测序数据与参考基因组进行比对，确定每个序列在基因组中的位置，得到比对结果，为后续的基因表达定量和差异表达分析奠定基础。

4. 基因表达定量。

通过对比对结果进行基因表达定量，包括计算基因的读数、FPKM（每百万读数下的每个基因的预期碱基数）或TPM（每百万读数下的每个基因的转录本数），得到每个基因在不同条件下的表达水平。

5. 差异表达分析。

对不同条件下的基因表达水平进行比较，筛选出差异表达的基因，通过统计学方法进行显著性分析，找出在不同条件下表达水平显著变化的基因。

6. 功能注释。

对差异表达的基因进行功能注释，包括GO（Gene Ontology）富集分析、KEGG（Kyoto Encyclopedia of Genes and Genomes）通路分析等，确定差异表达基因的功能和通路富集情况。

7. 数据可视化。

通过绘制热图、火山图、基因表达曲线等图表，直观展示转录组数据的分析结果，帮助研究人员更好地理解和解释数据。

8. 结果解释。

最后根据分析结果，解释不同条件下的基因表达模式和功能富集情况，为后续的生物学实验和功能研究提供重要参考。

1 / 17 RNA-seq(转录组学)的分析流程和原理 在开始详细讲解RNA测序之前，我们先来了解一下它的基本步骤： 1. 建库：提取RNA，富集mRNA或消除rRNA，合成cDNA和构建测序文库。 2. 测序：然后在高通量平台（通常是Illumina）上进行测序（每个样本测序reads在DNA测序中，读数是对应于单个DNA片段的全部或部分的碱基对（或碱基对概率）的推断序列。深度为10-30 Million reads。） 3. 分析：先比对/拼装测序片段到转录本，通过计数、定量，样本间过滤和标准化，以进行样本组间基因/转录本统计差异分析。 大致了解这个过程之后，我们就先从建库开始了解 建库的难点在于提纯出mRNA, 一般在我们抽离出的RNA中rRNA占比很大，其他还会有tRNA、microRNA等。我们需要从抽离出的RNA中提取出mRNA，并建立cDNA文库。 这里以应用最广泛的Illumina公司的Truseq RNA的建库方法为例来进行介绍。 首先，利用高等生物的mRNA通常有poly(A)尾的（使mRNA更稳定，翻译不容易出错）特点，用带有poly(T)探针的磁珠与总RNA进行杂交，这样磁珠就和带poly(A)尾巴的mRNA结合在一起了。 接下来，就回收磁珠，把这些带poly(A)的mRNA从磁珠上洗脱下来。再用镁离子溶液（或者超声波）进行处理，把mRNA打成小段。 然后，利用这些被打断的mRNA片段，以随机引物进行逆转录，得到第一链cDNA。 再根据第一链cDNA合成出ds-cDNA。 对cDNA在平末端进行3’端加A碱基（腺苷酸）（adapter接头上带了T碱基头，为了和adapter配对） 2 / 17 在双链cDNA的两端加分别上Y型接头 再经PCR扩增经筛选的目的基因，就得到可以上机的测序文库了。 这个建库方法对RNA的完整度有较高的要求。也就是说，只有在mRNA大部分是完整的状态下，才能得到比较好的效果。 因为带Poly(T)的磁珠，它所吸附的是带有Poly(A)的那些序列。那么如果mRNA发生了降解，也就是mRNA断掉了，那么磁珠所吸附下来的片段，都是那些靠近3'端的那些断片，而那些5'端的断片，是吸附不下来的。会在富集过程中被洗脱掉。这样进行数据分析的时候，就会发生一定的数据偏差。 同时，为了保证能够测到尽可能完整的mRNA序列，Illumina公司建议：先对总RNA进行一次质量检测。 一般是用Agilent公司出品的Bioanalyzer 2100毛细管电泳仪，对总RNA样本进行一次电泳质检。根据18S和28S这两个核糖体RNA的电泳峰是否高、是否尖，来判断RNA的质量。这两个峰越高、越尖，也就说明RNA的降解就越少，完整度就越高。 同时系统会自动打分。这个分值，叫“RIN”值。也就是RNA的完整度评分值。是“RNA Integrity Number”的英文首字母缩写。这两个峰越高、越尖，那么打分，也会越高。（之前说过，越高越尖说明RNA完整度越好）RIN值最高是10分，最低是0分。 Illumina公司推荐用RIN值在8.0以上的RNA进行建库和测序。 Illumina测序即二代测序技术通常是RNA-seq测序中最常用的。接下来，我们简单了解一下它的原理。先说flow cell ，它是含有管道的一片玻璃，是测序反应发生的位置。 然后在接上接头之后，我们首先要进行桥式PCR。 3 / 17 首先是右图第一步接好接头，然后将DNA样品调整到合适的浓度加入到flowcell中，再加入特异的化学试剂。就可以使得序列的一端与flowcell上面已经存在的短序列通过化学键十分强健地相连。如左图。图中不同的颜色表示的是两种不同的接头，分别对应序列之前加入的两种接头。 连接以后就正式开始桥式PCR。首先进行第一轮扩增，将正向链序列补成双链。加入NaOH强碱性溶液破坏DNA的双链，并洗脱。留下从固定的接头上延续出来的反向链，正向链被洗脱。得到右图的第二步。 加入缓冲溶液，这时候序列自由端的部分就会和旁边的接头进行配对，形成桥状。得到右图的第三步。 进行一轮PCR，在PCR的过程中，序列是弯成桥状，所以叫桥式PCR，一轮桥式PCR可以使得序列扩增1倍。也就是右图的第四步。 然后用碱液进行解链，由于合成是各有一端是延续在接头上的，不会被洗脱。从而得到两条固定好的序列。如右图的第五步 我们一直重复第3步到第5步，就可以得到就会得到一个具有完全相同序列的簇，一般叫cluster。就是第6步 利用多个序列进行桥式PCR就会得到多个簇。 形成簇之后，加入带不同荧光标记的含有叠氮基团的dNTP和聚合酶，由于叠氮基团的存在（又称为可逆阻断终止技术），一个循环只能延长一个碱基，即一个碱基结合到测序链上为一个循环，当结合上一个碱基后，把其他dNTP和酶用水冲掉，然后进行激光扫描，根据激光扫描颜色判断是哪个碱基，根据互补原则反推出结果。再切掉叠氮基团，进入新的循环，加入新的dNTP和酶，再延长碱基，冲掉dNTP和酶，再激光扫描。如此往复，就可以测出序列的内容。如图，我们每轮读一个碱基，将每一轮的碱基按顺序组合，就可以得到待测序列了。 我们回到文中的图片。这个流动池图，就是刚才我们所说的原理，不同颜色代表不同碱基。我们根据颜色就可以读出这一轮循环所接上去的碱基。 4 / 17 说完了测序，我们回到数据分析。第一步，一般需要进行Data Cleaning。就是把从原始数据（Raw Data）到干净数据（Clean Data）的过程。Illumina测序仪下机的数据通常为Bcl格式，是将同一个测序通道（Lane）所有样品的数据混杂在一起的，所以公司一般不会提供Bcl文件。测序公司使用Illumina官方出品的Bcl2FastQ软件，根据Index序列分割转换成每个样品的FastQ文件，打开长这样：每一条序列（read）包含四行，第一行是read的ID，第二行是序列，第四行是序列中每个碱基的测序质量。 原始数据没法直接分析，是因为部分reads测序质量较低，可能会误导后续结果，因此需要对低质量碱基太多或N（未能识别的碱基）太多的reads进行去除； 此外，部分测序文库的插入片段太短，导致测到的是两侧的接头序列，这些序列接头也需要从reads中去除。最后，我们也会对清洗前后的Raw Data和Clean Data进行评估，评估内容包括碱基质量、序列长度、碱基比例、GC含量、重复序列等。一般测序公司都会提供clean data。就不详细赘述了。 然后呢，需要把测到的RNA片段，先mapping(比对)到基因组上。 这张图我们可以看到，在外显子表达的时候，由于中间有内含子存在，使得两端外显子表达之后再连接。这个带有junction site 的部分，我们再测序后称为junction reads。了解到的几种mapping方法，都是先放置non-junction reads 部分再根据对 junction reads进行处理（根据junction Library或切断为seed）。也就是我们最后看到的这样，把序列比对上基因组。 可是对于用于分析没有参考基因组和基因注释的物种，我们是没法进行比对的。所以我们需要进行序列组装。以Cufflinks软件为例。 首先，先不要看下面的三个转录本，假设只知道上面的已经比对好的reads。Cufflinks会识别这些不会出现在同一个转录本的reads即不相容。例如比较清晰的蓝色和黄色，如果在一个转录本，那么，黄色应该在相同的位置断开，而不是跨过去。如图中所示，经识别不同颜色片段两两互不相容的，同一颜色的片段， 5 / 17 彼此是相容的。 那么，通过将每个相容的片段作为节点并和它最近且相容的片段相连。就可以得到所谓的重叠图。 在此基础上，我们基于精简原则。Cufflinks在能够覆盖所有reads的路径中选择出互不相连且最少的一组路径，作为最优路径。我们怎么理解这个最优路径呢，互不相连且最少。那不符合这个原则，我们可以走出什么样的路径呢？比如这样的，这样的。路径就会变得很复杂，没有这个必要。选择最优途径，即覆盖了所以reads，又减少了不必要的重复工作。 最终呢，我们就得到了这样的三个转录本集合。 我们可将序列组装想像成从大量片段的文字中拼凑出一整篇文章的过程：被测序的分子就是那篇文章，而测序片段就是那段文章中，随机切取出来的句子。其中一种重建出这段文章的方式，就是找到句子中重叠的部分，因为一旦找到够多重叠的部分，我们就有机会将每个句子连接到一起，进而得到原始的文章。 好了，说完组装之后呢，我们回到这个数据分析。 那么在比对或者是组装完了之后，可以先看一下，有多少的RNA片段，是在靠近基因的5'端的位置，又有多少片段在是靠近基因的3'端的位置。 那么这张图上，就是把所有的基因，都按其外显子的长度，拉直。然后来看，比对上的片段。 这样一个比对的结果，就可以让我们看见前面Poly(T)磁珠在抓mRNA的时侯。捕获下来的这些mRNA是不是完整的，如果捕获下来的这些mRNA大部分是完整的话，那么这个图形靠近5'端的曲线就会显得比较饱满。它的高度会和3'端的高度差不多。 6 / 17 反之，如果这根曲线的3'端是很高的，而5'端是比较低的，我们就可以初步判断，这个RNA有一定程度的降解。 因此，我们可以推断在捕获过程当中，有相当一部分（mRNA），它的5'片段因为与3'片段的Poly(A)片段的尾巴断开了，所以，没有被捕获下来。所以，这个RNA，是有一定程度降解的。 在知道了测序的质量之后，接下大家来要关注的就是不同样本之间、各个基因的mRNA的表达量的差异。 在做这些比较的过程当中，目前最常用的，对基因表达量进行相对定量的一个指标（也叫归一化方法），就是RPKM值。那么RPKM，是Reads Per Kilobase of exon model perMillion mapped reads的英文的首字母缩写。 大概意思就是每一百万条可以比对到基因组上的，Read当中，有几条是可以比对到某个特定基因的，这个数值再除以该基因的外显子的长度，得到的这样一个最终的比值。好，再来一遍，每一百万条比对到全基因组的read中，有几条可以比对到特定基因上面去。几条这个数值再除以该基因的外显子的长度，得到一个最终的比值。 基本构成我们是知道了，那我们接下来，在详细的分步地对这个公式进行一下解释。 首先，就是比对到某个基因的外显子上的Read数，去除以这次所测到的、全部可以比对到基因组上的Read数。这个比较容易理解就是：这个基因所表达出来的mRNA它所被测到的片段，来和所有被测到的、可以Mapping(比对）到基因组上的片段来进行比较。这个的意义是什么呢？因为在是随机抽样的情况下，可以认为这个mapped到特定转录本的read数正比于该基因（也就是特定转录本）的表达量，那么我们就可以得到该特定基因的表达量在总表达量中的占比。 这点是比较容易理解的。 7 / 17 那么比较费解的是，为什么还要除以第二项，就是“除以这个外显子的长度”。这是因为建库过程当中，RNA被打断（并逆转录）成若干个180-200BP左右的小片段，如果一个基因的长显子越长，那么它所产生的mRNA就越长，那么mRNA越长，被打出来的小片段就越多，那么它的read数就越多。 我们来假设，一个A基因，它的mRNA的长度，假设它是500bp，那么它的mRNA可能被打成2 个，200Bp左右的小片段；还有一个B基因，如果这个B基因的mRNA是200bp长，那么它刚好就是自成一个200bp的小片段。 我们来看，A基因是2个小片段，而B基因是1个小片段，所以，A基因在测序过程当中，它被测到的概率就会比B基因整整大出去一倍。那当我们检测到时，就可能会认为A基因的表达量是B基因的2倍。那么结合我们之前的公式来看，我们的mapped read读数也就是分子的确是A是B的两倍。 但其实我们知道这是不对的，我们可以看出，他们的表达量是一样的，只是A基因比较长。这就是我们为什么要把刚才第一项比出来的比值，然后再除以这个外显子的长度。通过上面的解释，我们就可以理解：除以这个外显子的长度，它的目的，是修正这个mRNA长度所引起的mRNA的Read数的偏差。 通过这种修正，就能够还原出一个比较真实的、原始的表达拷贝数状态。 这个，就是“RPKM”指标的原理。 那么作为一种针对全转录组的分析，我们希望是一次看到一个整体的样本（表达）差异的情况。而不仅仅是看少数几个基因的表达差异。 所以，科学家做了一种叫“火山图”的一个图形，来比较形象地来说明样本之间的表达差异。 那么我们来看这张图。 看最上面，我们知道对照样本为ALD，检测样本为ABF。然后，我们先明确一下。在差异表达火山图中的每一个点，都表示一个基因。一个什么样的基因呢？就是

博学笃行 自强不息

1

转录组数据分析

引言：

转录组数据分析是研究生物体在特定条件下细胞内mRNA的表达情况的一种方法。转录组数据分析的目的是识别差异表达的基因，揭示其在生物学过程中的功能和调控机制。随着高通量测序技术的发展，转录组数据的获取变得更加容易和可行，因此转录组数据分析也成为生命科学研究中的重要工具之一。本文将介绍转录组数据分析的基本流程和常用的分析方法。

一、转录组数据的获取

转录组数据的获取通常使用高通量测序技术，例如RNA-Seq和microarray。RNA-Seq是一种基于测序的转录组分析技术，可以直接测定细胞或组织中的mRNA的序列，无需参考基因组序列，因此能够检测到新的转录本和未注释的基因。Microarray是一种基于杂交的转录组分析技术，通过固定的探针阵列检测RNA样本中的mRNA的水平。两种技术各有优劣，研究人员可以根据自己的研究目的选择适合的方法。

二、转录组数据分析的基本流程

转录组数据分析的基本流程包括数据预处理、差异表达分析和功能注释。数据预处理是对原始转录组数据进行筛选、过滤和归一化处理，以得到可靠的表达量信息。差异表达分析通过统计学方法找出博学笃行 自强不息

2

在不同条件下差异表达的基因。功能注释则针对差异表达基因进行功能分析，识别其参与的生物学过程和通路。

三、数据预处理

数据预处理是转录组数据分析的第一步，其目的是筛选和过滤掉噪声数据，并对数据进行归一化处理。常见的数据预处理步骤包括：读取原始数据、去除低质量读段、去除接头序列和低质量碱基、去除rRNA序列、对碱基进行质量修剪、进行序列比对和表达量计算、基因表达量归一化等。数据预处理的目的是为后续的差异表达分析做好准备。

四、差异表达分析

差异表达分析是转录组数据分析的核心步骤，它通过统计学方法找出在不同条件下差异表达的基因。常见的差异表达分析方法有DESeq2、edgeR和limma等。这些方法利用统计学模型和假设检验来计算基因表达的显著性差异，并生成差异表达基因列表。差异表达基因的筛选标准通常是根据p值和fold change的阈值进行选择。