使用HISAT StringTie和Ballgown进行的转录本表达水平分析（中译）(二)

实验设计

用HISAT进行reads比对

RNA-seq分析开始时，将reads与参考基因组进行映射，以确定其基因组位置。这种映射信息使我们能够收集每个基因对应的reads子集，然后进行组装，并对这些reads所代表的转录本进行量化。已经开发了几种将reads与基因组对齐的快速算法，其中Bowtie和Burrows-Wheeler对齐法(Burrows–Wheeler aligner, BWA)是最广泛使用的算法之一。这两种方法都使用一种叫做Burrows-Wheeler变换（Burrows–Wheeler transform, BWT）的数据结构。它以高度压缩的形式存储参考基因组。通过使用一个特殊的，称为费拉吉娜-曼奇尼（FM）索引的方式，这些程序可以极其迅速地搜索基因组，其比对速度以每小时数百万个reads计算。 RNA-seq映射器（mapper）需要解决一个额外的问题，它不存在于在纯DNA比对中：许多RNA-seq 产生的reads将横跨内含子。RNA-seq流程捕获并测序成熟的mRNA分子进行，其中内含子已被剪接去除。因此，一个短的reads可能与相隔10,000 bp或更多的两个位置对齐（人类内含子的平均长度为 >6,000 bp ，一些内含子的长度为>1 Mbp）。对于一个典型的使用100bp reads的人类RNA-seq数据集，>35%的reads将跨越多个外显子。对齐这些跨外显子的reads比对齐一个外显子内的reads要困难得多。 HISAT使用两种类型的索引（index）进行比对：一个全局的全基因组索引和数以万计的小型局部索引。这两种类型的索引都是使用与Bowtie2相同的BWT/FM索引构建的。而HISAT系统甚至使用了一些Bowtie2的代码。由于使用了这些高效的数据结构和算法，HISAT生成拼接排列的速度比Bowtie和BWA快了数倍，而使用的内存却只有Bowtie的两倍。HISAT是被设计为作为TopHat 和TopHat2（后者是由本协议的一些作者开发的）的继承者而设计的，具有兼容的输出，但它的运行速度是TopHat和TopHat2的50倍。快50倍。对于与人类基因组的比对，该索引需要 8千兆字节（GB）的内存，这使得它可以在传统的台式机上运行；一台电脑处理20个人类RNA-seq样本，每个样本有1亿个读数，应该需要不到一天时间。对于具有较小基因组的物种，内存和处理量的要求也会相应减少。在本协议中，我们使用了HISAT2，它包含了算法上的改进，产生了比原来的HISAT系统略高的准确性。 RNA-seq数据的比对可以用来检测新的剪接点、转录起始点和转录终止点。与TopHat一样，如果这些事件（event）存在于数据集（data）中，HISAT将检测所有这些事件。用户还可以提供基因注释文件，指定已知基因的位置，包括其所有的外显子和内含子边界。使用高质量的注释文件可以改善HISAT在已知基因附近的比对，但这种输入是可选的。

用StringTie进行转录本组装和量化

RNA-seq实验的分析需要准确重建每一个基因所表达的异构体，以及对这些异构体的相对丰度进行估计。因为注释仍然不完整，即使是对研究最深入的物种，如人类，我们也不能为我们的程序提供所有基因和所有异构体的完整列表。尽管如此，准确的量化得益于准确了解哪些reads来自于哪些异构体，由于reads比转录本短得多，所以无法完美计算。StringTie将RNA-seq reads中的转录本组装起来，这些reads已经与基因组进行了比对。首先将reads分组为不同的基因位点，然后将每个基因位点组装成不同的基因组。然后将每个基因座组装成所需数量的异构体来解释这些数据。StringTie使用一种网络流算法，从表达量最高的转录本开始，它同时对其进行组装和定量。然后删除与该转录本相关的reads，并重复这一过程，组装更多的异构体，直到所有的reads都被使用，或者直到剩余reads的数量低于设定的转录噪音（transcriptional noise）水平（用户可自行调整）。

StringTie的网络流算法可以比以前的一些方法更准确地重建转录本，这是因为它同时计算丰度和外显子结构。该算法还具有计算运行效率，使它的运行速度比其前身Cufflinks快很多倍，而且使用的内存也少很多。大多数RNA-seq的运行可以由StringTie在只有8GB随机存取内存（RAM）的传统台式电脑上处理，甚至大型数据集也可以在一小时内运行。用户可以选择向StringTie提供一个包含参考基因模型的文件；这个注释文件包含 “已知 “基因的外显子结构说明，以及这些基因的名称。如果用户提供了这个注释文件，StringTie将使用它作为组装的指南，它也将用该文件中的名称来标记它组装的基因。注释对重建低丰度基因很有帮助，对于这些基因来说reads数量太少，无法准确组装。请注意，注释文件是可选的，StringTie会根据需要组装额外的、没有注释的基因来解释数据。在用StringTie组装转录本后，用户可以使用gffcompare工具（Box1；表1）来确定有多少组装的转录本与注释的基因完全或部分匹配，并计算出有多少基因是全新的。另外，用户可以跳过对新基因和转录本的组装，只用StringTie对注释文件中提供的所有转录本进行量化（图1）。

• Box1 | 将一组转录物与已知基因列表进行比较

  比较基因和转录物的清单首先要求清单的格式相同。许多生物信息学程序利用 GFF（gene finding format, 基因搜索格式）来表示基因和转录物。GFF是一种简单的以制表符分隔的文本格式，描述基因、转录本、外显子和内含子在基因组上的属性。StringTie使用GTF，它是GFF的一个特殊版本，包含一些额外的数据类型。详见此链接 我们已经开发了gffcompare工具（可从JHU或GitHub下载）来比较一个或多个由StringTie产生的GTF文件和一个GFF或GTF的参考注释文件。假设StringTie的输出是transcripts.gtf，参考注释文件是chrX.gtf，可以用以下命令运行gffcompare软件。

  1	$ gffcompare –G –r chrX.gtf transcripts.gtf

  在-r选项之后是要用作参考的注释文件，而-G选项告诉gffcompare比较输入的transcripts.gtf文件中的所有转录本，尽管那些可能是多余的。 gffcompare程序是基于CuffCompare工具的，它是Cufflinks/Tuxedo套件的一部分，CuffCompare手册中记载的许多使用选项和输出结果也适用于gffcompare程序。所有由gffcompare生成的文件的前缀都是gffcmp，除非用户在 -o 选项中选择不同的前缀。当如上所示使用时，gffcompare会产生一个输出文件，称为gffcmp.annotated.gtf，它在每个转录本上添加一个 “类代码”（在表1中描述）和参考注释文件中的转录本名称。这使用户能够快速检查预测的转录本与注释文件的关系。这里显示的gffcompare命令还将在gffcmp.stats输出文件中计算不同基因特征（如外显子、内含子、转录本和基因）的灵敏度和精确度统计。敏感性被定义为正确重建的来自注释的基因的比例，而精确性（positive predictive value, 也被称为阳性预测值）捕捉到与注释重叠的输出的比例。因此，这两个衡量标准都是相对于输入注释而言的，可能不能完全反映出组装的基因和转录物的真正准确性。gffcmp.stats文件还包含其他几个衡量标准，如StringTie输出中包含的新外显子、内含子或基因的总数。

• Table1 | 类别代码——用于描述转录本集合与参考注释的比较情况

  Class Code	Description
  =	预测的转录本与参考转录本有完全相同的内含子
  c	预测的转录本包含在参考转录本中
  j	预测的转录本是一种潜在的新异构体，它与参考转录本至少有一个剪接结点
  e	预测的单外显子转录本与一个参考外显子和一个参考内含子的至少10bp重叠，表明可能是前体mRNA（possiblepre-mRNA）片段
  i	预测的转录本完全落在参考内含子中
  o	预测转录本的外显子与参考转录本重合
  p	预测的转录本位于参考转录本的2kb范围内（可能是聚合酶活跃的片段）
  r	预测的转录本有>50%的碱基与软屏蔽的soft-masked（重复的）参考序列重叠
  u	与已知参考转录本相比，预测的转录本是基因间区（intergenic）的
  x	预测转录本的外显子与参考文献重叠，但位于相反的链上
  s	预测转录本的内含子与相反链上的参考内含子重叠

在包括此协议在内的大多数实验设计中，实验包括多个RNA-seq样本.存在于一个样本中的基因和异构体很少与存在于所有其他样本中的基因和异构体相同，但它们需要以一致的方式进行组装，以便可以进行比较。我们通过使用StringTie的合并功能将所有的装配合并在一起来解决这个问题，该功能合并了在全部样本中发现的所有基因。因此，在一个样本中由于缺乏覆盖率而缺少外显子的转录本可能会被恢复到它的全长；关于合并（merging）如何改善组装(assemblies)的说明，见图2。

图2 | 使用StringTie的合并功能合并转录本组合。在这个例子中，来自四个不同样本的四个部分集合被合并成两个转录本A和B。样本1和2都与参考注释一致，在此用于合并和扩展以创建转录本A。样本3和4彼此一致，但与注释不一致，这些被合并以创建转录本B。

在合并步骤之后，再运行一次StringTie来重新估计合并后的转录本的丰度。重新估计的步骤使用与初始步骤中使用的相同的丰度估计算法。

许多（不是全部）RNA-seq的用户都对发现有差异表达的基因感兴趣，而且许多人都专注于研究得很好(well-studied)的基因。因此，使用一种完全依赖注释的，不会去发现新的基因和新的异构体的方法（特别是对于像人类、小鼠和果蝇这样研究良好的物种），可能是很诱人的。然而，这些物种可用的标准注释文件可能遗漏了成千上万的异构体，以及许多基因。例如，最近的一项研究使用核糖体足迹分析（ribosome footprint profiling）验证了7273个以前未被注释的人类基因，其中大部分是研究人员在一组非常大的RNA-seq数据（28亿读数）上运行StringTie时发现的。StringTie自动检测新的基因和新的异构体，如果它们存在于你的实验数据中，不管它们是否出现在标准注释文件中，这里描述的协议可以发现影响这些基因的差异表达。

用Ballgown 进行差异性表达分析

分析、可视化和统计建模是收集和量化转录本后的下一个步骤。R编程语言和Bioconductor软件套件拥有一套全面的工具来处理从绘制(plotting)原始数据到归一化(ormalization)再到下游统计建模的任务。我们的流水线(pipeline)使用Ballgown软件包，它被设计为成为SAT和StringTie等上游命令行软件与Bioconductor社区的下游功能之间的桥梁。

StringTie可以产生一组链接的表格（a linked set of tables），可以使用Ballgown包中的函数直接读入R。在R中产生的数据包括三个表中的信息。(i) 表型数据（phenotype data）-关于正在收集的样本的信息；(ii) 表达数据（expression data）-每个样本中每个外显子、连接点、转录本和基因表达量的标准化和非标准化措施；(iii) 基因组信息（genomic information）-给出外显子、内含子、转录本和基因的位置的坐标，以及包括基因名称等信息的注释。

一旦信息被读入R会话（session），分析就会以交互方式进行。换句话说，R和Bioconductor中的任何包或函数都可以使用。大多数差异表达分析都是沿着这个路径进行的：(i) 数据可视化和检查，(ii) 差异表达的统计测试，(iii) 多重测试校正，(iv) 下游检查和结果总结。Ballgown为这些步骤中的每一步都提供了函数，当数据集出现特定问题（包括RNA-seq分析中通常出现的问题）时，可以与其他R命结合使用。

协议的Ballgown部分从向R中加载数据开始。这需要加载StringTie产生的丰度数据和描述样本的表型信息。这个过程中的一个关键步骤是确保基因组样本的标识符（ID）与表型数据的ID相匹配。如果样本似乎并不匹配，Ballgown将抛出一个错误。装入数据后，下一步是检查转录本的丰度估计分布。在这一点上，丰度估计（以FPKM值表示，即每百万映射读数中每千碱基的转录物片段）已经相对于文库大小进行了归一化处理，因此，样品之间整体分布的任何极端差异都应引起人们关于样品、排列或丰度估计出现问的担心。

差异表达分析采用线性模型。附在转录本上的FPKM值通常是高度倾斜（skewed）的。因此，为了稳定方差，Ballgown的内置函数会应用对数转换，然后拟合标准的可用于测试差异表达的线性模型。Ballgown允许进行时间序列和固定条件（time-course and fixed-condition）的差异表达分析。RNA-seq分析中出现的一个常见问题是没有考虑到混杂因素–诸如批次效应（batch effects）–这些因素与分析并无直接关系，但可能影响基因的表达水平。使用Ballgown stattest函数，你可以直接指定任何已知的混杂因素。Ballgown允许你在基因、转录本、外显子或连接点水平上进行分析。其结果是一个表格（table），它包含了差异性表达测试的特征信息；条件之间的折叠变化（若存在）；P值；差异性表达的q值。 这些结果可用于下游基因组分析软件的进一步挖掘，对重要的结果进行可视化检查（visual inspection），或将结果导出到表格中，以便共享并使用Excel等程序进行人工检查。