Expression Analysis

第10章 Expression Analysis

10.1 Introduction

范围：PDF page 299；对应印刷页码 279。下一小节 Step 0: Choose an Expression Analysis Technology 已排除。

---

2003 年人类基因组完成测序，使我们第一次获得了一个相对完整的人类全部基因初步目录（Lander et al. 2001；Venter et al. 2001）。尽管自第一版草图序列发布以来，基因组本身以及其中编码的基因集合又经历了大量演化，但围绕基因表达如何被调控、以及这些表达数据如何用于刻画不同表型状态并探索其性质的问题，至今仍有许多尚未解决。事实上，我们已经知道：即便在同一个个体体内，同一套基因组也会在不同细胞类型中呈现出截然不同的表现形式，而这些基因表达谱还会随着不同条件而变化，包括健康状态与疾病状态之间的变化。

科学家在基因组完成测序之前就已经意识到这些问题的重要性，并发展出多种方法来检测不同表型之间 RNA 表达的差异。早期技术一次只能检测一个或少数几个基因，但 DNA microarray 技术的出现，使研究者能够同时检测大量基因，从而把表达分析扩展到全基因组尺度（Schena et al. 1995）。此后，DNA microarray 被广泛用于探索模式生物和人类疾病中的基因表达模式（DeRisi et al. 1996；Spellman et al. 1998；Golub et al. 1999；Perou et al. 1999；Callow et al. 2000；Konstantinopoulos et al. 2011）。

不过，基于 microarray 的早期基因表达分析也带来了明显挑战，许多早期研究都受到可重复性不足问题的困扰（Ioannidis et al. 2009；Ishmael et al. 2009）。随后，计算生物学家与实验生物学家在实验流程和分析方法上投入了大量工作，逐步建立起更加稳定的实验室操作规范和分析流程，使 DNA microarray 研究结果的一致性显著提高，同时也进一步凸显出严谨实验设计与重复设置的重要性（Hegde et al. 2000；Simon et al. 2002；Irizarry et al. 2003, 2005；Bolstad et al. 2004；Larkin et al. 2005；Quackenbush 2005）。

随后，ultra-high-throughput sequencing 技术的引入，又推动了 RNA sequencing (RNA-seq) 的发展。这类实验不再像 microarray 那样强烈依赖预先设定“要测哪些分子”的假设，而是能够以更开放的方式检测转录本（Kahvejian et al. 2008；Nagalakshmi et al. 2008）。虽然 RNA-seq 在表达测定的技术路径上与 microarray 明显不同，但稳健 RNA-seq 分析方法的建立，实际上仍然继承了许多在 DNA microarray 发展过程中总结出来的基本经验教训。

本章的目标，是按步骤介绍开展基因表达评估时需要考虑的问题与常用方法：从实验设计出发，逐步进入数据标准化、样本比较与结果解释等核心环节。虽然本章属于入门性介绍，但作者希望这里提供的内容能够成为读者后续深入研究的起点，并引导读者进一步系统学习本章涉及的各种方法。

Ch10 Step 0: Choose an Expression Analysis Technology

第10章 表达分析 / 第0步：选择表达分析技术

---

引言

尽管 RNA-seq 目前是主流技术，看似从这个角度切入不太寻常，但 DNA microarray（DNA 微阵列）至今仍被广泛使用，并且在开发实验和分析方案时，它相对于 RNA-seq 的一些优势值得纳入考量。这两种技术的应用流程都始于从感兴趣样本中提取和纯化 RNA，并利用来自逆转录病毒的逆转录酶（reverse transcriptase）将这些 RNA 逆转录为互补 DNA（cDNA）。随后，通过测序或杂交的方式测定基因组中各基因的相对丰度。这些丰度数据再被用于下游分析，以理解不同生物状态之间基因表达模式的变化规律，并揭示这些变化如何帮助我们认识所研究系统的生物学特性。尽管现有技术在原理上存在共性，但它们之间也确有值得关注的差异。

---

DNA 微阵列（DNA Microarrays）

DNA 微阵列是首个实现全基因组层面基因表达分析的技术。DNA 微阵列的工作原理是检测溶液中标记的 cDNA 靶分子与固定在固相载体表面单链基因特异性 DNA 探针之间的杂交事件。尽管 DNA 微阵列在早期曾受到噪声问题的困扰，经常出现结果不可重复的情况，但技术的进步和分析方法的完善已大幅提升了其所产生数据的质量。DNA 微阵列还具有其他多项值得考虑的优势。

第一，DNA 微阵列的基因内容（gene content）是明确已知的，即每个被检测表达情况的基因或转录本都由一个或多个探针（或探针组）来代表。这种设计的一个好处是，我们事先就知道哪些基因会被检测；因此，通常可以相当有把握地判断某个特定基因是否表达以及其相对表达水平。虽然可能存在交叉杂交或其他伪影（如不同的杂交效率），但商业化阵列的质量和当今实验室操作流程的稳健性已经使检测质量和可重复性大幅提升。

第二，由于该技术已经相当成熟，几乎在微阵列分析的各个环节都有经过充分验证的分析方法可用。DNA 微阵列分析工具最核心的单一资源库是 Bioconductor，其中拥有大量、成熟的微阵列分析全流程方法。

第三，通过 Gene Expression Omnibus（GEO）和 ArrayExpress 数据库，可以获取海量的 DNA 微阵列数据资源。这些数据库提供了额外的独立数据集，可用于估算验证单个实验结果所需的样本量。

值得注意的是，尽管 Bioconductor 也包含许多 RNA-seq 数据分析方法，但在最佳实践的共识形成方面，RNA-seq 远不如微阵列分析成熟。此外，虽然 GEO 和 ArrayExpress 中也包含 RNA-seq 数据，但这些数据仅占这两个资源库总数据量的一小部分。虽然可以通过比较微阵列数据与 RNA-seq 数据来判断大致趋势，但无法直接将微阵列的杂交信号强度与 RNA-seq 的 reads 计数进行对比。

---

RNA-seq

RNA-seq 的历史可以追溯到 1990 年代的 cDNA 测序，以及 2000 年代初的基因表达系列分析（Serial Analysis of Gene Expression, SAGE）——一种允许对短 cDNA 片段进行测序的技术。然而，RNA-seq 真正发展成为独具特色的方法，是在过去十年间随着超高通量测序（ultra-high-throughput sequencing）的出现而实现的：从极少量的 RNA 出发，即可产生数千万条（甚至更多）序列 reads。

RNA-seq 之所以迅速成为基因表达谱分析的主流技术，原因有以下几点。

第一，RNA-seq 不仅能够测量"基因"的表达水平，还能提供单个等位基因（allele）和转录变体（transcript variant）表达水平的信息。尽管这种灵活性使得该技术具有独特性和开放性，但实际上很少有研究真正利用了这一能力，大多数分析仍回归到考察转录本的总表达水平。此规律的最显著例外是在癌症中利用 RNA-seq 识别融合转录本（fusion transcripts）——这是 DNA 微阵列完全无法实现的应用。

第二，RNA-seq 方法随着时间推移已显著改进，如今只需更少量的 RNA 即可作为起始材料。RNA-seq 的应用还包括来自单个细胞的转录组分析，这需要专门的分析方法来应对数据的稀疏性（见 Single-Cell Sequencing）。使用少量起始材料还意味着可以从单个生物样本中采集小型活检样本，或获取多种不同类型的数据（如 RNA-seq 和 DNA 甲基化数据），从而使整合分析变得可行。

第三，RNA-seq 不受预定义转录本集合的限制。RNA-seq 实验能够揭示尚未被描述的"新"基因的表达情况，包括非编码 RNA（non-coding RNA）的转录水平。RNA-seq 数据也可用于检测聚腺苷酸化病毒转录本序列（polyadenylated viral transcript sequences）。

---

选择权在你手中（The Choice is Yours）

尽管微阵列仍是一种可行的替代方案，但微阵列与 RNA-seq 之间的成本差距已降至使 RNA-seq 通常成为默认选择的程度。有鉴于此，本章将重点介绍 RNA-seq 分析，并在读者有兴趣了解微阵列表达分析方法时，引用本书的前版内容。然而，表达分析的许多通用原则是相通的，可以作为思考各类大规模基因组研究的通用模板，因此在下面的讨论中也会包含一些相关案例。

---

> 边界说明：PDF page 301 底部出现的 "Step 1: Design the Experiment" 及后续正文已排除，不属于本节范围；下一节将处理实验设计相关内容。

Ch10 Step 1: Design the Experiment

第10章 表达分析 / 第1步：设计实验

---

实验设计的基本逻辑

在生物学中，大多数成功的实验都建立在成熟的假设检验（hypothesis testing）思想之上。我们首先识别一个问题，并提出一个可能的机制。随后，设计一种能够检验该假设的扰动方式，对系统施加干预，并收集数据，以观察是否出现与所提出机制一致的变化。我们在系统中观察到的响应，将支持或否定原先的假设。

在这类实验中，我们通常会尽量严格控制变量，从而精确测量它们的影响，并且一次只扰动一个参数。良好的实验设计还要求有足够的重复，以便估计我们希望测量的效应大小。

全基因组基因表达技术改变了我们处理生物学问题的方式。过去我们常常只关注单个基因，而现在可以同时考察某个系统中数千个基因的响应，并寻找与表型变化相关的表达模式改变。我们既可以利用这类大规模实验来检验既有假设，也可以根据基因表达模式的变化生成新的假设，再在后续实验中进一步验证。不过，基因组尺度技术所带来的观察范围和数据规模，并不意味着我们可以忽视实验设计与数据分析的严谨性。

像所有实验一样，一个基因表达谱实验应当始于一个定义明确的问题，并且实验必须收集足以回答该问题的数据。最常见的实验设计，是比较两组实验组别（groups 或 cohorts），例如处理组与对照组，或患病人群与健康人群。设计这类研究时，一个关键要素是确保实验具有足够数量的独立生物学重复（independent biological replicates），从而使处理组和对照组在样本规模上足以支持合理比较。

---

功效、验证与样本量

在大规模转录表达谱实验中，功效与样本量（power size）计算 notoriously difficult，也就是出了名地难做，原因很大程度上在于：表达水平本身波动很大，而相对效应大小在实验开展前通常并不清楚。一个可行策略，是先做一个小规模预实验（pilot experiment），从中识别潜在信号，用来估计效应大小，然后据此进行更严格的统计功效计算，并设计完整实验。

另一种策略，则是把视野放到原始实验之外，在整体方案中预先纳入验证阶段（validation stage）。例如，可以使用一种独立技术——如 reverse transcription polymerase chain reaction（逆转录聚合酶链式反应，RT-PCR）——去验证一个较小的“显著基因集”（significant gene set）；更理想的做法，是纳入一个独立的验证人群（validation population），再对其单独进行表达谱分析，以评估原始结果是否成立。

---

混杂因素、批次效应与元数据

另一个重要考虑，是设计一种能够避免混杂（confounding）并消除批次效应（batch effects）的实验策略。这既包括样本采集的实验安排，也包括基因表达数据采集的实施策略。具体而言，应尽量确保“处理组”和“对照组”样本在同一时间、相同条件下被采集；并且在 RNA 提取、文库构建以及测序数据生成等环节中，对样本进行合理混合（mixing）。

一个重要却常被忽视的问题是：待分析样本是否具有充分的元数据（metadata）。例如，如果分析的是乳腺癌样本，就必须了解每个样本的疾病亚型（disease subtype），并在实验设计阶段考虑各亚型的分布情况。如果缺少这些信息，再考虑到 RNA-seq 实验中样本数相对于被检测基因数往往很少，就很容易出现这样的情况：观察到的表达差异，其实只是样本被分配到不同组别时某种偏倚的结果。

作者举了一个具体例子：他们曾分析过一个据称能够根据原发肿瘤表达谱预测乳腺癌肺转移（lung metastasis）的基因特征（gene signature），后来却发现，用于识别该特征的所有发生转移的样本都属于 basal subtype（而这正是最容易向肺转移的亚型）。换句话说，这个特征实际上预测的是 basal subtype，而不一定是真正预测转移本身。因此，在分析数据之前，必须考虑处理组与对照组之间是否存在人口统计学差异（demographic differences），不同亚群患者所接受的治疗是否不同，患者是否来自不同医院或国家，以及样本是否以不同方式采集和处理。

信不信由你，这些混杂因素都曾在已发表研究中真实出现过，而且几乎都本可以事先避免。与其在数据分析阶段再费力解释这些问题，不如在实验开始之前就主动识别潜在混杂因素，这样更值得投入精力。

---

从分析策略反推实验设计

作者提出了一种非常实用的方法：先从你在数据收集完成后打算采用的分析策略（analytical strategy）出发，再反向推导实验设计。这样可以确保你拥有足够数量的样本，以及恰当且完整的元数据，从而保证所获得的数据和信息真正能够回答你的实验问题。

---

> 边界说明：本节止于 PDF page 302 中 Step 2 标题前；Step 2: Collect and Manage the Data – and Metadata 及其正文已排除，不属于本节范围。

Ch10 Step 2: Collect and Manage the Data – and Metadata

第10章 表达分析 / 第2步：收集并管理数据——以及元数据

---

为什么这一步不能被忽视

一次转录表达谱实验，通常包括对对照系统施加某种扰动、采集生物学样本，然后生成最终要被分析的数据。乍看之下，收集并管理这些相关数据似乎是理所当然的事情，但现实中，这恰恰是一个经常被忽视、并且日后很容易“反噬”研究者的环节。

---

数据存在哪里不是最关键，关键是要有组织

处理这类问题的方法有很多。最简单的做法，是把数据存放在共享硬盘的某个文件夹中；更系统的做法，则是建立数据库，把最终数据统一纳入其中。无论采用哪一种策略，最重要的一点始终不变：必须保持组织有序，并清楚记录每一份数据分别对应哪个项目。

换句话说，这一步的核心并不在于你选用了多复杂的存储平台，而在于数据管理是否可追踪、可核对、可复现。如果样本、文件、项目和后续分析之间的对应关系没有被明确记录，那么即便前面的实验做得再认真，后面的分析也会很容易陷入混乱。

---

> 边界说明：本节从 PDF page 302 底部真实标题 Step 2: Collect and Manage the Data – and Metadata 开始，包含其在 PDF page 303 顶部延续的正文；PDF page 303 中部再次出现的 Step 3: Data Pre-Processing 标题及其后正文均已排除，不属于本节范围。

第 10 章 表达分析

第 3 步：数据预处理（Step 3: Data Pre-Processing）

在数据分析的初始阶段，人们可以采用多种策略来管理数据：既可以将数据简单地存放在共享磁盘的某个文件夹中，也可以创建一个数据库来统一存放所有数据。无论选择哪种策略，最重要的是保持有序的组织习惯，并清楚地记录每份数据对应哪个项目。

在数据能够用于实验组之间的比较或其他任何用途之前，必须首先将原始数据映射到具体的基因或基因转录本。尽管这一步看起来相当简单，但实际上可以采用多种不同的方法——而这些不同的方法，自然也可能导致不同的最终结果。虽然存在许多公认的进行数据预处理的方法，但研究者应当注意并仔细记录自己在从原始数据中识别基因转录本时所做的选择。

DNA 微阵列

对于 DNA 微阵列而言，原始数据的映射似乎很简单，因为人们通常认为微阵列不过是由一系列针对每个目标基因设计的固定探针组成的阵列。然而，许多微阵列——最著名的是 Affymetrix GeneChip——实际上使用成组的探针（即"探针组"）来共同定义一个基因的表达。事实上，Affymetrix 芯片的设计不仅包含利用参考基因序列设计的"完全匹配"（perfect match, PM）探针组，还包含了"错配"（mismatch, MM）探针——这些探针的序列中央仅有一个碱基与参考序列不同。PM 探针用于提供杂交信号的估计值，而 MM 探针则用于估计由非特异性杂交和背景荧光造成的背景信号。

探针到基因的映射信息通常包含在"芯片设计文件"（chip design file, CDF）中，该文件作为输入进入任何分析流程的早期阶段，用于建立荧光强度与基因表达水平之间的映射。研究界对于究竟应该使用哪些数据来进行探针映射，存在相当大的争议：有人主张仅使用 PM 探针，有人开发了非标准的 CDF 文件，而大多数研究者则使用 Affymetrix 提供的 CDF 文件。与基因表达分析的许多其他方面一样，这里没有绝对正确的答案——只需要做出一个合理的选择，并记录下这一决定，以便他人能够重现分析。然而，研究者必须始终清楚探针映射到的是哪些基因标识符（以及这些标识符的哪个版本）——无论这些标识符是官方基因名、RefSeq ID、Ensembl ID 还是其他类型，因为这些决策都会影响后续涉及将表达数据映射到生物通路或功能分类系统（如 Gene Ontology, GO）以及使用基因集富集分析等技术时的分析结果。

RNA-seq

RNA-seq 面临着类似的一系列挑战，不过其映射过程的神秘性要低一些。RNA-seq 的原始输出是一组序列 reads，这些 reads 被映射到一组基因或基因转录本上。为此，最常用的方法是"基于参考序列的组装"（assembly on reference）：先将 reads 映射到基因转录本，然后进行组装，最后量化每个基因的整体代表性。在此过程中，参考数据库的选择决定了映射结果。研究者可以选择 RefSeq、Ensembl 基因或其他任何合适的参考序列。目前已开发出大量用于 reads 映射、组装和量化的算法，包括 Burrows-Wheeler aligner（BWA）（Li and Durbin, 2009）、Bowtie/Bowtie2（Langmead et al., 2009; Langmead and Salzberg, 2012）以及 STAR aligner（Dobin et al., 2013）。

近年来，为了应对日益庞大的 RNA-seq 数据集，研究者开发了伪比对（pseudo-alignment）和拟映射（quasi-mapping）方法，包括 Salmon（Patro et al., 2017）、Sailfish（Patro et al., 2014）和 Kallisto（Bray et al., 2016）。这些方法专门设计用于削弱大数据带来的计算复杂性，包括内存限制问题。伪比对和拟映射绕过了全比对的过程，其核心思想是用 k-mer 来表示转录组，然后通过哈希表将 k-mer 映射到 de Bruijn 图表示（一种表示 k-mer 之间重叠关系和映射路径的图形化方法）或后缀数组（suffix array，即一个 k-mer 的所有扩展或后缀的有序数组）。通过合理定义哈希函数，可以忽略参考序列的大部分内容，从而将 reads 查询映射到有限的潜在目标上。

小结

与基因组数据分析的许多方面一样，对于最优方法的选择目前尚无明确共识，而且方法本身也在不断演进。重要的是从标准方法中进行选择，对希望分析的数据始终如一地应用同一方法，并以确保分析可被重现的方式记录自己的选择——包括记录所使用的软件版本和数据库版本。

Step 4: Quality Control

作为科学家进行的任何测量都包含误差。有些误差是随机的，分析所用的统计方法正是用来在自然变异中估计真实信号；有些误差是系统性的，同样可以用统计方法进行估计和处理。然而，有些误差源于失败的实验，最佳处理方式是识别并剔除这些失败实验产生的数据。在基因表达分析实验过程中，这类误差可能来自 RNA 样本中的污染物、劣质的实验试剂，或者简单的实验室操作失误。

一旦获得原始数据，最重要的问题之一就是判断这些数据是否具有足够的质量，能够进入后续的分析流程。尽管生物变异是任何数据集中都希望体现的特征，但毫无疑问，失败的实验应当从待分析数据集中剔除。实验室中的表达分析涉及许多复杂步骤，从降解的输入 RNA 到劣质试剂再到简单的操作失误，任何环节出问题都可能产生噪声水平极高的数据，使任何有意义的分析都无法进行。幸运的是，现已有一系列工具可用于分析微阵列表达分析和 RNA-seq 实验产生的数据，为微阵列和基于测序的数据提供了成熟的指标体系。与该领域的所有内容一样，用于分析这些数据的工具也在快速发展，因此建议读者关注相关文献综述或向正在进行基因表达分析的同行了解可能出现的新方法。

Quality Control Tools

Bioconductor 软件包 arrayQualityMetrics 提供了一套广泛的工具（包括许多从其他 Bioconductor 软件包整合而来），用于评估单色和双色微阵列数据的质量。作为 arrayQualityMetrics 包的输入，用户需提供一个微阵列强度矩阵，并可选地提供样本和探针信息，这些信息以 Bioconductor 对象形式组织，如 AffyBatch、ExpressionSet、NChannelSet 或 BeadLevelList 对象——这些对象将不同技术的表达数据与表型信息相关联。arrayQualityMetrics 的输出包括每个阵列的伪彩色图像和 MA 图，用于评估数据质量。在 MA 图中，M 值是两个强度的对数比，A 值是强度的对数平均值。对于双色阵列，这些图使用每个通道的强度；对于单色阵列，M 值使用每个样本的中位强度作为比值的分母。图 10.1 展示了一个在标准化前后 MA 图的示例，其中可见标准化过程消除了水平轴以下的系统性曲率。

还有其他多种诊断图可用于识别单个劣质阵列或整体较差的数据集。这些包括来自 affy 包的 RNA 降解图（Gautier et al. 2004）、来自 affyPLM 包的相对对数表达（RLE）箱线图和归一化未缩放标准误差（NUSE）箱线图（Brettschneider et al. 2008），以及来自 simpleaffy 包的 QC 统计图（Wilson and Miller 2005）。这些质量控制工具的结果被整合到一个 HTML 文档中，为理解已收集的原始数据提供了宝贵的资源。

对于 RNA-seq，FastQC 是一个广泛使用的软件包，以与微阵列 arrayQualityMetrics 包非常相似的方式，为原始高通量测序数据提供了一系列简单的质量控制检查。FastQC 包含多个分析模块，允许用户探索序列质量的各个方面，提供多个汇总图和表格，并将结果导出为基于 HTML 的报告。FastQC 中的模块可提供基本统计信息，包括 reads 数量、读长和 GC 含量。用户还可以查看箱线图，显示沿所有 reads 每个位置上序列质量分数的逐位评估，以及每个序列质量分数的分布图。这两种图都可以对测序运行的整体质量提供良好的评估。

另一个用于评估整体序列质量的有用图是 per base sequence content（逐碱基序列组成）。对于任何基因组，GC 含量应沿任意随机序列 reads 的长度保持一致，即 %A = %T 且 %G = %C。然而，文库制备方案通常会在待测 DNA 的 5' 端连接短的引物和接头序列，而这正是预期看到 GC 分布显著偏离的位置。相关的度量是 per base N content（逐碱基 N 含量），用于量化因无法以足够置信度调用碱基而将定义核苷酸替换为 N 的频率；该信息有助于识别测序反应中的失败循环。

图 10.2 展示了一组 RNA-seq 样本相对于参考序列的碱基错配（MM）率直方图示例。虽然大多数样本的 MM 率较低，但仍有一些离群值可从下游分析中移除。还有其他工具可识别异常的序列重复水平、过度代表的序列、未检出的接头和过度代表的 k-mers。FastQC 还包含用于分析 microRNA、宏基因组序列和表观遗传学检测（如甲基-seq）的工具。

提供研究整体质量概览的一个令人兴奋的新工具是 MultiQC。该工具将多个样本的质量控制报告（来自 FastQC 以及其他工具）聚合在一起，在单个易于阅读和理解的 HTML 报告中呈现，可帮助识别并随后从分析中移除低质量样本。

Screening for Misidentified Samples: PCA on Y Chromosome Expression

我们尚未触及的一个要素是与每个样本关联的注释质量。任何分析最终都依赖于将样本分配到不同组别、比较组间表达水平，以及确保没有可能扭曲分析的混杂因素。分析的质量（以及我们对结论的信心）取决于个体被分配到特定组的可靠性；这反过来又取决于我们是否能准确地将每个样本与适当的注释元数据（如个体的性别、年龄、治疗状态和其他表型数据）相关联。虽然这看起来相对简单，但样本的错误注释是一个比人们预期更常见的问题。例如，发现 GEO 中 46% 的研究存在劣质或错误的注释——这些错误在提交前通过简单的验证步骤就可以很容易地识别出来（Toker et al. 2016）。

通常很难测试样本注释的准确性，因为大多数实验的目的通常是发现组间差异，而不是利用已知差异将样本分配到组别。不过，对于几乎任何数据集都可以运行一项测试来了解样本注释的质量：男性和女性的注释是否正确。如果简单地查看 Y 染色体基因的表达并进行主成分分析（PCA；见 Principal Component Analysis），则预期会发现两个明显不同的组，因为女性不表达 Y 染色体基因（Paulson et al. 2017）。

作为一个示例，我们分析了来自 The Cancer Genome Atlas（TCGA；来自 Genomic Data Commons [GDC] Data Portal）的结直肠癌基因表达数据，以及来自 GEO 的五个结直肠癌数据集（GSE14333、GSE17538、GSE33113、GSE37892 和 GSE39582）。当我们使用 PCA（见 Principal Component Analysis）分析 Y 染色体基因的表达时，发现了两组明显不同的样本——一组表达这些基因（因此可能是男性），另一组仅有背景噪声水平的表达（因此可能是女性）。然而，我们发现 TCGA 的 456 个样本中有 11 个（2%）存在性别错标：有些样本虽标注为女性，却与男性样本聚类（因此表达 Y 染色体基因）；另一些则应为相反情况。相比之下，在 GEO 数据集中，我们发现 1376 个样本中有 85 个（6%）按性别错误聚类。虽然我们可以从下游分析中移除这些错标样本，但我们放弃了其中一个 GEO 研究，因为其错标率接近 15%，这让我们对其余样本注释的真实性也产生了怀疑。

Step 5: Normalization and Batch Effects

The Importance of Normalizing and Batch-Correcting Data

任何基因表达分析的输出，都可以表示为一个由正值填充的表达矩阵，这些数值代表每个样本中每个探针或基因的观测表达水平。对于 RNA-seq 数据，这类数据也可表示为“计数矩阵（count matrix）”C：其中每一行对应一个基因，每一列对应一个样本，而矩阵中每个位置的条目则表示在该特定样本中比对到该基因的 reads 数目。

可将其形式化写为：C = [c11 c12 … c1n; c21 c22 … c2n; … ; cm1 cm2 … cmn]。

如果世界足够完美，我们只需逐基因比较不同样本组之间的表达值或 read counts，就能直接比较它们的表达水平。但现实中有许多因素会扭曲这些数据，例如起始 RNA 量不相等、文库制备中的伪差、杂交或测序效率差异，以及其他各种问题。

标准化（normalization）是数据分析中的常规步骤，它使我们能够在不同样本或不同数据集之间进行可比的比较。比如在 northern blot 分析中，常会使用 actin 或 GAPDH 这类 housekeeping genes 来校正每个样本中其他基因的杂交信号，以补偿实验中样本量的变化。这里隐含的假设是：一个或多个基因在所有样本中的表达水平保持一致，因此可以把这种“基线”用作校正其他基因测量值的参照。总体而言，存在若干类标准化过程，可用于实现样本间或基因间的比较。

FPKM and Count Data

在 DNA microarray 中，每个基因通常由一组大体等价的探针表示。例如，在 Affymetrix 芯片上，基因由 probe sets 表示。每个 probe set 中的探针长度均为 25 bp，并包含一一对应的 PM 和 MM probes，它们只在一个碱基上不同（而且恰好位于探针正中间）。此外，每个 probe set 所含探针数相同，并定位于其目标转录本的 3' 端附近，以尽量保证从 RNA 制备杂交文库时具有近似一致的效率。

但在 RNA-seq 数据分析中，按基因进行标准化时考虑的问题会有所不同。RNA-seq 需要先制备文库并测序；如果只是简单统计每个基因获得的 reads 数，那么较长的基因天然更可能获得比短基因更多的 reads。因此，原始计数数据通常会被转换为 reads per kilobase million（RPKM）、fragments per kilobase million（FPKM）或 transcripts per million（TPM）等数值。三者的计算方式略有差别，但都以将 reads 比对到转录本并对结果进行缩放为起点。

RPKM 的计算方式是：先统计比对到某个转录本的 reads 数目，再除以“每百万比对到基因组的 reads 数”，并进一步按该基因长度进行标准化。最后这一步是为了校正这样一个事实：一个 2 kb 的基因天然会比一个 1 kb 的基因多获得大约两倍的 reads。例如，若有 400 万条 reads 比对到某个基因组，其中 5000 条比对到一个长度为 2 kb 的特定基因，则该基因的 RPKM 值为 625（[5000/2]/4）。FPKM 是 RPKM 的扩展形式，适用于 paired-end sequencing，因为在这种情形下，一个 RNA-seq 文库片段的两端都会被测序。FPKM 的基本思想与 RPKM 相同，只是额外考虑到两条 reads 可能来自同一个片段。

TPM 与上述指标相似，但它被标准化到统一的“每百万”尺度上，因此更便于样本间比较。TPM 的计算方法是：先对每个基因，用 reads（或转录本数）除以该基因的长度（单位为 kb），得到每个基因的 reads per kilobase（RPK）值；然后将基因组中所有基因的 RPK 值求和，得到累计 RPK；再将该累计值除以 100 万，得到一个“scaling value”；最后，再用每个基因的 RPK 除以这个 scaling value，即得到每个基因的 TPM。概念上，TPM 与 microarray 的测量值最为接近，因为它既考虑了每个基因的长度，也比较了标准化后的转录本计数。

Sample and Quantile Normalization

随着使用 DNA microarrays 产生的高通量基因表达数据分析越来越多，人们很快发现，所谓存在“恒定不变”的 housekeeping genes 这一假设并不成立；事实上，所有基因的表达水平都会发生变化。在缺乏可靠参照的情况下，新一代标准化方法的重点转向：考察一个样本中所有基因表达水平的分布，再将不同样本的分布调整得彼此相似。

最初，人们曾采用对单个样本的平均表达水平或中位表达水平进行校正的方法，但这类方法无法补偿由实验伪差导致的分布差异。如果我们假设一个细胞所能产生的 RNA 总量是有限的，那么当某些基因表达升高时，其他基因的表达就应相应降低，从而使相关样本的整体表达分布保持一致。

从概念上讲，一种简单做法是观察基因表达值的分布，并将其切分为更小的区间，也就是 quantiles（分位区间）。随后，可以按 quantile 逐段调整数据，使实验中的所有样本都具有相同的分布，从而能够比较任意特定基因在不同样本之间的表达变化。这个思路听起来相对简单，但仍值得更细致地看一看它的具体过程。

图 10.3 用一个包含 4 个样本和 6 个基因的例子展示了 quantile normalization 的过程。我们把这些测量值表示为一个 genes-by-samples 矩阵，并用三个简单步骤实现标准化。第一步，按照图中的蓝色箭头，对每个基因计算其在所有样本中的中位值，并将这些中位值从低到高排序。第二步，按照绿色箭头，对每个待分析样本，将每个基因的表达值替换为它在该样本内从低到高的秩（rank），从而得到一个 rank matrix。第三步，按照金色箭头，将 rank matrix 与上述排序后的中位值集合结合起来，用对应秩次的中位值替换 rank matrix 中的秩值。

在该图示例中，所有 rank 1 的基因都被赋予表达值 6，所有 rank 2 的基因都被赋予表达值 9，依此类推，于是原始数值就被这些按秩排序后的中位值所替代。这种方法背后的一个重要假设是：所选参考对于所有样本都适用；但在底层生物学差异较大的情况下（例如同时分析多种组织类型），这一假设未必成立。尽管这一过程乍看之下有些反直觉，但研究表明，这种方法能产生稳健且具有可重复性的表达值，从而可在样本间进行比较（Bolstad et al. 2003）。

Figure 10.3 展示了 quantile normalization 的总体思路：左上角方框中给出了 4 个样本、6 个基因的表达数据；首先计算每个基因的中位表达值并从低到高排序（蓝色箭头）；然后根据各样本中的表达高低对基因排序（绿色箭头）；最后，用排序后中位值中对应秩次的数值替换 rank matrix 中的秩（黄色箭头）。

Additional Methods of Sample Normalization

标准化方法的选择几乎是任何基因表达分析讨论中最容易引发争议的问题之一，因为每种方法都依赖一组特定假设，而这些假设对具体数据集而言，可能有人认为合理，也可能有人认为并不充分。下面介绍 4 种已被广泛使用的补充标准化方法（Li et al. 2015）。

Counts per Million

counts per million（CPM）缩放方法与 TPM 相似，都是把计数值标准化到统一的“每百万”尺度上，以便更容易在样本之间进行比较。CPM 及其他缩放方法的基本假设是：每个细胞产生的 RNA 总量大致相同，因此总 reads 数或总 counts 应该保持恒定。计算 CPM 时，只需先汇总每个样本的 reads 总数，再对所有样本的 read counts 做缩放，使其总量相等。CPM、TPM 以及按基因长度标准化后的数值，都是 RNA-seq 早期基因表达分析中最常见的标准化方法。

Upper Quantile Normalization

这种缩放方法假定：在低到中等表达水平范围内，各样本的计数分布彼此相似，但在高于第 75 个分位点后开始出现偏离。该方法的做法很直接：将各样本中低于第 75 分位点的 counts 规模调整为彼此一致，并令标准化后的计数 yij 按如下方式缩放：yij = cij / q75j。

Relative Log Expression

这种方法假设：各样本中的计数值与跨样本的基因表达几何均值高度一致，并且 read count 的频率会随着测序深度呈指数增长。RLE 使用每个样本的几何均值，并对样本中的 reads 进行缩放，使所有样本的几何均值相同。对于每个样本，其标准化因子被定义为：特征 read counts 与所有样本对应 read counts 几何均值之比的中位数。DESeq（Anders and Huber 2010）和 DESeq2（Love et al. 2014）都采用该方法作为默认标准化策略，后文还会继续介绍这两个框架。

Trimmed Mean of M Values

trimmed mean of M values（TMM）方法建立在“绝大多数基因并不存在差异表达”这一假设之上。TMM 选取一个样本作为参考，然后将其余每个样本与该参考进行比较，计算 log fold-change（在 microarray 分析中通常称为 “M” 值），去掉最外侧 30% 的 M 值，再计算平均 M0 值，并将其缩放到所有样本之间相等。TMM 被用于 edgeR（Robinson et al. 2010）的检验框架中，后文也会进一步讨论 edgeR。

Batch Correction

batch correction（批次校正）是任何大规模基因组学分析中的另一个关键环节；更广泛地说，凡是以分组方式采集测量值、或者在不同时间采集不同样本的科学研究，都会面临这一问题。所谓 batch effects，是指由系统性误差带来的偏差来源，例如使用了不同批次的试剂、处于不同实验室条件、处理了不同时间采集但相似的样本，或者仅仅因为不同实验人员执行同一实验时存在内在差异。由于每个样本通常都会产生大量测量值，因此在高通量实验中，batch effects 尤其明显。

抵御 batch effects 的第一道防线，是良好的实验设计。理想情况下，应尽量在同一时间采集所有样本，并在同一时间、同样条件下、使用同一批试剂完成检测。对于样本量很大的研究，这通常并不现实；这时次优但非常重要的策略，是在实验的每个环节都混合安排 cases 和 controls，避免样本分组与批次效应发生混杂。举例来说，如果周一测 controls、周二测 cases，那么任何可能出现的 batch effect 都会与 case/control 状态完全纠缠在一起，使得差异根本无法分辨。而若将 cases 与 controls 混合处理，就能在很大程度上抵消潜在的批次影响。

检测 batch effects 的一个非常简便的方法，是基于基因表达数据绘制 PCA 图（见 Principal Component Analysis），其中图中的每个点代表一个样本。首先，可以按实验条件给点着色（例如 cases 与 controls）。在理想情况下，我们希望看到不同条件之间有清晰分离；但现实中看到的结果，会取决于信号、噪声以及潜在 batch effects 的共同作用。接着，还可以根据其他相关变量重新给这些点着色，比如样本采集日期、RNA 提取日期、标记或文库构建日期、芯片批次或测序运行批次等，以观察是否出现某种模式。理想状态下，PCA 图中除条件分离外不应出现其他结构；但事实上，经常会看到同一批次的样本聚在一起，或不同批次彼此分离。批次效应校正试图控制的，正是这些差异。若希望更高通量地识别问题，则突出显示 PCA 成分与表型变量、批次变量之间的相关性，通常是一种很有用的技术。图 10.4 就给出了一个逐步展示这一过程的示例。

需要指出的是，尽管 PCA 非常有用，但它通常只能捕捉较大的 batch effects。实际上，单个基因或某些基因子集也可能受到实验条件的不利影响，从而改变任何下游分析的结论。

幸运的是，目前已经发展出一些方法，能够处理 batch artifacts：它们通过识别与批次相关的实验信号，并对这些信号进行校正，从而更好地识别真正与实验分组相关的差异表达基因。

两种广泛使用的批次校正方法是 COMBAT（Johnson et al. 2007）和 surrogate variable analysis（SVA；Leek et al. 2012）。这两种方法都会检查表达数据，寻找那些与 batch 相关的基因；而在 SVA 中，相关对象还可以是其他非生物学变量。然而，正如前文所述，只有在相关生物学变量没有与这些非生物学变量发生混杂时，才能识别并纠正这些关联。例如，如果所有 cases 都在一个批次中运行，而所有 controls 都在下一个批次中运行，那么 batch effects 就会与表型完全混杂；这样一来，某些本来由非生物学因素导致的差异，就会被误判为真实的生物学差异。要发现并消除 batch effects，就必须尽可能让每一个实验组都在每一个实验室批次中得到代表。

Figure 10.4 展示了一个 batch effects 的 PCA 示例。图中使用 BatchQC 软件及其 vignette 所模拟的数据：(a) 为前两个主成分的 PCA 散点图，点代表样本，并按 phenotype condition 着色；(b) 为模拟基因表达数据在各样本中的箱线图，同样按 phenotype condition 着色；(c, d) 则展示了相同数据，但按测序“batch”着色，此时可以看到第一主成分以及基因表达值都出现了显著的批次差异。

第六步：探索性数据分析

在大规模基因组数据的分析过程中，一个常见的问题是：某一群体中是否存在由独特基因表达模式所定义的亚群？这个问题只有在拥有足够数据来搜索可用于识别和区分群体的模式时，才能得到合理的回答。这类无偏倚搜索所使用的方法称为无监督方法（unsupervised methods），因为它们是去发现模式，而非去验证是否能区分预先定义的组别。在这一大类方法中，最常用的三种是：层次聚类、主成分分析（PCA）和非负矩阵分解（NMF）。

如前所述，表达数据的一种便捷表示方式是表达矩阵（expression matrix）——一个基因×样本矩阵，其中每一行是一个"基因向量"（gene vector），代表某一基因在所有样本中的表达水平；每一列是一个"样本向量"（sample vector），代表单个样本中所有基因的表达水平；矩阵中的每个元素代表某一基因在某一样本中的表达值。热图（heatmap）是该矩阵的一种可视化形式，其中每个单元格的着色取决于其信号的强度。层次聚类与 NMF 都根据共有的表达模式对样本和/或基因的子集进行分组，并在热图背景下可视化结果。PCA 则对样本或基因的表达向量（即矩阵的行或列）进行运算。

层次聚类

层次聚类（Hierarchical Clustering）已成为基因表达数据分析中应用最广泛的技术之一，其优势在于方法简单、结果易于可视化（Eisen et al. 1998; Michaels et al. 1998; Wen et al. 1998）。初始时，以 N 个聚类开始，其中 N 为目标数据集中的样本数（或基因数）。层次聚类是一种聚合式方法（agglomerative approach）：单个表达谱被逐步连接形成节点，这些节点又进一步连接，直到过程完成，最终形成一棵层次树。

层次聚类的核心逻辑是：找出彼此最接近的向量，然后根据它们之间的距离将样本归为一组。当然，根据表达谱来测量样本（或基因）之间的距离有多种方式，其中最常用的两种是：欧氏距离（Euclidean distance measure，适用于基因表达的绝对水平重要的情况）和皮尔逊相关系数距离（Pearson correlation distance measure，适用于关注相关性模式的情况）。

选择好距离度量方法后（图 10.5），算法以直接的方式推进。以下说明假设是对样本进行分组，但聚类过程同样适用于基因。

1. 计算所有待聚类样本之间的成对距离矩阵。
2. 在距离矩阵中搜索最相似的两个样本或两个聚类。初始时，每个聚类仅包含一个样本。如果若干对样本恰好具有完全相同的相似度得分，通常随机选择其中一对，但也可以用其他方法来决定选择哪一对。
3. 将选中的两个聚类合并，生成一个新的聚类，该聚类现在包含两个或更多对象。
4. 计算这个新聚类与所有其他聚类之间的距离。无需重新计算所有距离，因为只有涉及新聚类的距离发生了变化。
5. 重复步骤 2–4，直到所有对象归入一个聚类。

图 10.5 展示了层次聚类过程的简化示意。(a) 反复搜索彼此"接近"的样本或样本组，并将它们逐步合并形成更大的聚类。(b) 对应的层次聚类树状图（dendrogram），根据评估得到的相似性将样本连接在一起。

层次聚类方法有多种变体，反映了计算新定义聚类与其他基因或聚类之间距离的不同策略（即所谓的聚合方法，agglomeration methods）：

• 单连锁聚类（single linkage clustering）：使用一个聚类与另一聚类之间的最短距离。
• 完全连锁聚类（complete linkage clustering）：取任意两个聚类之间的最大距离。
• 平均连锁聚类（average linkage clustering）：使用两个聚类之间的平均距离。

通常，样本之间的关系使用树状图（dendrogram）来表示，树中的分支根据聚类过程中的连接关系构建。为了可视化样本之间的关系，用树状图重新排列表达矩阵热图中的行（或列，视情况而定），以显示数据集中的模式（图 10.6）。

图 10.6 热图展示了三个不同心脏组织中 100 个变异性最高基因的基因表达聚类。该热图中的表达水平经过行 Z-score 标准化，以最佳呈现表达差异。低表达用蓝色表示，高表达用红色表示。顶部的树状图使用欧氏距离进行层次聚类，表明三个组织（顶部颜色条）的样本聚为三组。左侧组织（绿色）对这组基因的表达较低，而中间组织（灰色）和右侧组织（橙色）的表达较高。基因集中有一小部分在右侧组织高表达，而在其他两个组织中低表达。

树状图的树形结构使其有助于识别在外观上可以将样本分成若干聚类的位置。然而，与其使用"肉眼检验"，不如使用客观方法来确定聚类的数量及其成员——幸运的是，有多种方式可以搜索这样的分组。

一种方法是直接利用构建聚类时计算的距离，作为每个单独聚类连通性的度量。当从单个元素沿着树状图向上移动时，聚类之间的距离增加；因此，随着距离阈值的增加，有效聚类数减少。另一种方法是使用自助法（bootstrapping）或刀切法（jack-knife）来测量树状图中关系的稳定性，并以此作为聚类数的度量。有多种自助法可以使用，但最简单的一种是：带放回地对数据集进行抽样，每次计算一棵新的层次聚类树状图，然后简单地统计树状图中每个分支被恢复的频率；以树状图上的百分比截断值来设定聚类数。在对基因聚类稳定性进行自助法估计时，应对生物样本集合进行重抽样；而在估计生物样本中的聚类数时，则应对基因表达向量进行自助法。刀切法类似，但并非重抽样，而是依次留出相应的向量，计算新的树状图，直到所有向量都被考虑过。同样，根据每个分支关系在树状图中被恢复的频率来估计每个聚类的稳定性。

许多层次聚类方法存在一个潜在问题：随着聚类的增长，用于计算距离的聚类表达向量可能不再准确地代表聚类内部的任何元素。例如，在对基因进行聚类时，每个聚类的"中心"通常是该聚类中所有基因的平均值；由此产生的基因表达向量的线性组合有时被称为元基因（metagene）。因此，随着聚类的推进，基因本身的实际表达模式变得不那么相关。此外，如果在过程早期做出了错误的聚类分配，该错误将被固定下来且无法纠正。一种可以避免这些问题的替代方案是使用分裂式聚类方法（如 k-means），将数据（基因或样本）划分为具有相似表达模式的组。

尽管聚类方法适用于任何数据集，但在实践中，对于其中许多基因在样本间并不变异的大型数据集，聚类往往效果不佳。因此，可以首先对数据应用统计过滤，选择在实验组之间变异性最大的基因。然而，这种过滤可能会给聚类带来偏差。如果选择能够区分两个实验组的基因，那么过滤最可能产生的结果将是两个聚类恰好落入预先定义的两个实验组中。更无偏的方法是：直接计算每个基因在样本群体中的方差，然后删除在数据集中变化不显著的基因——因为这些基因最不可能揭示样本收集中存在的亚类。但这种方法依赖于各实验组之间有良好的样本平衡。

主成分分析

如果我们观察一个样本×基因的表达矩阵，可以想象每个样本在每个被评估的约 25 000 个基因上都有其独特的表达水平。那么，每个样本都可以表示为该 25 000 维"基因表达状态空间"中的一个点。不难想象，样本集合可以可视化为该 25 000 维空间中的一团点云。然而，许多基因的表达水平可能是相关的，因此它们并不能提供太多可用于分离样本和区分点云中不同样本组的信息。

PCA（图 10.7）是一种降维方法（dimensionality reduction method），它搜索表达水平（即基因表达水平）变量的线性组合，以最佳解释样本之间的方差，然后将数据转换，使基因表达的特征向量（eigenvector。线性变换的特征向量是指当该变换作用于该非零向量时，所有值按相同的标量因子发生变化）按最佳分离数据的顺序排列。这样，数据被转换后，第一个"主成分"解释了数据集中最大的变异量，第二个主成分解释次大的变异量，依此类推。使用这些特征向量坐标来绘制数据，通常更容易可视化样本分离为不同组的情况。这反过来有助于理解表达数据集中的样本是否聚为具有基因表达巨大差异的特定子集。PCA 也是一个良好的质量控制工具，因为技术变异（如批次效应）可以通过直观检查 PCA 图来轻松检测。

图 10.7 对来自 GTEx（Genotype-Tissue Expression）项目 v6 版本的皮肤、全血及由这些组织衍生的细胞系的归一化基因表达矩阵进行主成分分析（PCA）得到的前两个主成分散点图。每个点代表一个样本，并根据其来源着色。第一主成分（PC1）分离了组织类型，第二主成分（PC2）将组织与细胞系分开。

PCA 建立在底层基因×样本表达矩阵的一系列简单线性代数变换之上。PCA 的工作流程示意图如图 10.8 所示。

1. 首先对矩阵（这里是矩阵的行，即基因）进行标准化，使每个基因的表达范围处于相同的尺度。
2. 计算协方差矩阵（covariance matrix），其中条目 ij 是基因 i 与基因 j 之间的协方差。两个基因之间的协方差基本上衡量的是它们在群体中所有样本的表达偏差是否相关。
3. 计算协方差矩阵的特征向量和特征值。特征向量是与协方差矩阵相乘后返回相同向量的向量，其中每个值乘以一个标量因子（即相应的特征值）。特征向量（或主成分）是矩阵的不变量，是基因的线性组合（因此有时被称为"特征基因"，eigengenes）。
4. 用特征向量来重新表述原始数据。这通过将原始表达矩阵乘以特征向量矩阵来实现。
5. 在新特征向量（彼此正交，类似于 x–y–z 轴）的基上绘制结果。

图 10.8 主成分分析（PCA）是一种降维方法，用于识别数据中捕获最大变异的变量组合，然后在主成分空间中绘制数据。此处，点代表在更高维"表达空间"中的实验，每个样本都有自己独特的表达谱（因此有独特的坐标）。PCA 识别数据变异最大的正交轴，并计算作为各基因线性组合的新坐标轴。然后将样本投影到"PC 空间"，通常只绘制前几个主成分。从数学上讲，PCA 将基因×样本矩阵分解为基因×模式矩阵（其列为各主成分）和模式×样本矩阵。在这种情况下，主成分有时被称为"元基因"（metagenes），因为它们是基因的线性组合。

在这一框架下，第一个特征向量解释了数据中最大的变异量；第二个特征向量解释次大的变异量，依此类推。在二维或三维图中，通常会检查样本的分布，按批次、受试者性别或处理组等着色，观察各种系统性因素和生物因素如何影响数据中的变异。

非负矩阵分解

NMF 是另一种降维方法。它将基因表达数据建模为两个非负矩阵的乘积，通过将基因汇总为更少数量的所谓"元基因"（meta-genes）。在 NMF 中，我们从一个 n × m（基因×样本）的表达矩阵 R 开始。我们使用一个 n × k 的特征矩阵 P，它包含每个基因在 k 个聚类中的质心值。用一个 n × k 的权重矩阵 Q 乘以 P，Q 提供了将 R 的列表示为 P 列的非负线性组合的权重。得到的乘积 R̂ 是原始矩阵 R 的近似：

R ≈ P × Q^T = R̂

上述说明中隐藏了一个关键点：NMF 需要事先了解数据中可能存在多少个实验组（k）。

在实践中，往往并不知道数据集中预期有多少个实验组，因此常见的做法是用多个 k 值运行该方法，然后选择最能解释所研究系统生物学的划分方式。在此探索性分析中，一个有用的定量指标是共表型系数（cophenetic coefficient），它衡量的是基因必须有多相似才能被归入同一聚类——本质上是组内距离与组间距离之比。一种广泛使用的方法是绘制共表型系数图，然后寻找其值的急剧下降（表明真正的聚类正在被过度分割），以此来选择最优的 k 值。这种图的示例如图 10.9 所示。

图 10.9 NMF 共表型系数图——见 PDF 第 296/297 页。

---

第六步 完

（下一节：第七步：差异表达分析——PDF 第 316 页）

第七步：差异表达分析

虽然探索性数据分析非常有用，但大多数基因表达实验的设计目的，仍然是检验这样一个假设：表型差异是否与功能相关基因的表达差异有关。检验这一假设最直接的方法，是判断是否存在一批基因在不同样本组之间表现出显著不同的表达水平。

在基因表达分析的早期应用中，人们常常先施加一个生物学过滤条件（biological filter），即寻找在不同条件之间平均变化达到 2 倍或以上的基因，这就是所谓的倍数变化过滤（fold-change filter）。这种做法对很多生物学家来说较为直观，但它忽略了基因表达的变化幅度与波动性。因此，差异表达的统计学评估很快就成为判断转录差异的标准方法。

然而，基因表达统计分析面临的最大问题在于：我们通常同时测量的是成千上万个基因，而样本数往往只有几十到几百个。这就带来了多重检验（multiple testing）问题——当测量项目远多于样本数时，即便完全出于随机波动，也很可能会出现一些看似在组间不同的基因。幸运的是，针对这一问题已经发展出了多种校正方法。为了理解如何更可靠地识别显著差异，首先需要了解常用差异检验方法及其背后的基本假设。

Student's t-Test: The Father of Them All

Student's t-test，通常简称 t 检验，是各学科中最常用的两组差异检验方法之一。它本质上可以看作一种信号与噪声的比较：先比较两组样本的平均表达水平，再利用标准差来判断两组均值差异是否足够显著。换句话说，t 检验衡量的是：均值差异相对于数据波动是否足够大，并据此估计观察到的差异由随机因素造成的概率。

根据具体问题不同，t 检验有多种形式：

• 双侧 t 检验（two-sided t-test）：检验某基因在两组之间是否存在高低任一方向的差异；
• 单侧 t 检验（one-sided t-test）：检验某基因在一组中是否显著高于另一组；
• 配对 t 检验（paired t-test）：用于配对样本设计，判断组间表达差异是否超过随机预期。

对于多于两组的比较，t 检验可以推广为 F-test 或 方差分析（analysis of variance, ANOVA）。

在表达分析中，如果实验具有配对设计（例如患者与匹配对照），配对 t 检验尤其常用。传统 t 检验默认数据近似服从正态分布，因此它不能直接修正基因表达研究中常见的均值—方差依赖关系（mean-variance dependency）。一种改进办法是使用经验 t 检验（empirical t-test）：在不同组之间反复置换样本，针对每个基因反复计算 t 统计量，再判断真实数据得到的 t 统计量是否显著高于置换数据的预期水平。

不过，t 检验的所有变体都面临同一个核心困难：多重检验问题。由于被检测的基因数远多于样本数，即便基因在真实生物学上并无差异，也可能仅因偶然而显得显著。因此，人们投入了大量工作来开发更稳健的方法，以识别实验组之间真正存在差异表达的基因。下面概览几类目前仍广泛使用的方法。

Limma

Limma（linear models for microarray and RNA-seq data）最初是为微阵列数据分析开发的，后来也被广泛用于 RNA-seq 数据分析。Limma 的做法是：为每个基因的表达水平拟合一个线性模型，然后使用moderated t-test 来识别显著差异表达基因。

这里的 moderated t-test 建立在一种经验贝叶斯方法（empirical Bayesian method）之上。与标准 t 检验不同，它会根据表达水平对方差进行缩放，并利用合并方差估计（pooled estimate of variance）和自由度（degrees of freedom）来更稳定地估计显著性。因此，即使样本量较小，Limma 也通常能够给出更稳定的统计推断。

此外，Limma 默认倾向于报告假发现率（false discovery rate, FDR）的估计，而不是只给出简单的 p 值。对于高维基因表达数据来说，这一点尤其重要，因为 FDR 比单个检验的 p 值更能反映整体显著结果中的错误发现风险。

Voom

负二项模型（negative binomial model）是处理计数型数据的重要方法之一，但在样本量较大时，直接计算负二项分布的代价可能较高，这会使一些方法（如 DESeq 和 edgeR）在某些场景下显得较为笨重。虽然对 RNA-seq 原始计数做对数转换（log transformation）有助于标准化数据，但这样做往往又会进一步扭曲方差估计。

Voom（variance modeling at the observational level）的核心思想，是为每个基因经验性建模其均值—方差关系：它把每个基因的 counts per million 的对数值（log-cpm）的标准差，建模为平均对数计数的函数。随后，Voom 将这一均值—方差趋势转化为每个观测值的精度权重（precision weight），并把这些权重纳入 Limma 的分析流程中。

这样做的好处在于：RNA-seq 数据仍然可以放进一个经验贝叶斯的线性建模框架中分析，同时又能无缝接入许多既有的 Limma 工作流。因此，Voom 因其兼容性和实用性，已经成为 RNA-seq 差异表达分析中最常用的方法之一。

Negative Binomial Models

DESeq（Anders and Huber 2010; Love et al. 2014）和 edgeR（Robinson et al. 2010）是两种广泛使用的 RNA-seq 差异表达分析方法。它们都直接处理 RNA-seq 的计数数据，而这类数据通常具有过度离散（over-dispersion）的特点，也就是其波动程度高于理想无偏计数数据的预期。

DESeq 与 edgeR 都会把一个基因表达水平的方差拆分为两部分：

• 生物学变异（biological variability）
• 技术变异（technical variability）

生物学变异是不同基因或不同条件之间真实差异的来源；技术变异则来自实验流程中的多个环节，例如样本采集、RNA 提取、文库构建、测序误差及其他实验条件，这些因素都会扭曲真实 RNA 计数。

这两类方法都假定：生物学变异的大小与过度离散程度成正比。因此，DESeq/DESeq2 与 edgeR 一般都假定数据服从负二项分布（negative binomial distribution），而不是泊松分布（Poisson distribution），并进一步拟合广义线性模型（generalized linear model），以估计并修正这种过度离散，从而更稳健地识别差异表达基因。

Fold-Change

如前所述，在 DNA 微阵列全基因组表达分析的早期阶段，倍数变化（fold-change）曾被广泛使用。统计学方法能够较为无偏地评估某个基因差异表达的证据强弱，但它们往往不符合许多生物学家的直觉：即某些基因表达水平仅有非常小的变化时，未必足以解释明显的表型差异。

当然，也有人会指出：对于转录因子（transcription factor）、激酶（kinase）等调控节点分子，即便很小的表达变化，也可能带来显著的下游效应。但对于大多数基因而言，绝对变化幅度很小时，其生物学解释通常并不容易。因此，很多研究会把统计显著性与倍数变化阈值结合起来使用：先找出统计学上显著的基因，再进一步筛选其中倍数变化大于 2 倍（或其他阈值）的基因，用于后续分析。

Correcting for Multiple Testing

如前所述，在任何 RNA-seq 实验中，被检测的基因数通常都非常大，而样本数相对较少，因此单纯由于偶然因素而找到“差异表达基因”的概率会显著升高。举例来说，假设我们分析 25 000 个基因，并按某种“组间差异显著性指标”（例如 t 统计量）对它们排序。每个基因对应一次单独检验，因此总共相当于进行了 25 000 次统计检验。

如果我们仅仅取其中“最显著”的前 5%，那么就相当于挑出了 1250 个基因，但我们并不能有把握地说这些基因中任何一个都是真正差异表达的。因此，差异表达分析必须引入多重检验校正。

Family-Wise Error Rate

家族错误率（family-wise error rate, FWER）用于估计：在进行多次统计检验时，出现一个或多个假阳性（false positives，也称 type I errors）的概率。

如果我们总共进行 c 次检验，并为每次检验设定显著性水平 \(\alpha\)，则有：

FWER ≤ 1 − (1 − α)^c

回到前面的 25 000 个基因示例，若显著性水平设为 0.05，那么几乎必然会至少出现一个假阳性，即：

FWER ≤ 1 − (3.8 × 10^-55)

处理多重检验问题的经典 FWER 方法之一是 Bonferroni 校正（Bonferroni correction）。其思路非常直接：把原本的 p 值阈值除以检验次数，即将 \(\alpha\) 替换为 \(\alpha / c\)。

不过，对于基因表达分析来说，Bonferroni 往往过于严格。例如，当检测 25 000 个基因时，原本的 p < 0.05 会被压缩成 p < 2 × 10^-6。如此苛刻的阈值常常会导致：即便比较的是生物学上差异非常明显的样本，也很难找到哪怕一个满足阈值的基因。虽然也有一些 Bonferroni 的修正版本试图缓和这种过严问题，但实际使用中，更常见的是基于 FDR 的方法。

False Discovery Rate

Benjamini and Hochberg (1995) 提出了假发现率（false discovery rate, FDR）的概念，用于处理多重检验带来的统计问题；这一思想后来又由 Benjamini and Yekutieli (2001) 进一步扩展。与 FWER 试图估计“是否至少会出现一个假阳性”不同，FDR 承认假阳性不可避免，并进一步估计：在一组被判定为显著的结果中，有多大比例可能是假阳性。

换句话说，FDR 本质上是在回答：显著结果中预计有多少比例是错误发现。因此，研究者可以根据研究目的自行设定一个可接受的错误比例阈值。这里讨论的 FDR 专指多重检验中的 FDR-controlling procedure，它与二分类任务中其他语境下的 FDR 含义相似，但并不相同（参见 Box 5.4）。

在最简单的实现中，从任意统计检验的 p 值出发计算 FDR 并不困难。设我们有 N 次检验，并为每次检验都得到了一个 p 值，则可以按如下步骤计算：

1. 将所有 p 值按从小到大排序：p1, p2, …, pk, …, pN。
2. 对于给定的 FDR 水平 q，寻找满足下式的第 k 个 p 值：

pk ≤ (i/N)(q/c(N))，其中

c(N) = Σ(i=1 到 N) (1/i)

然后把排名位于 1 到 k 的那些基因判定为在 FDR = q 水平下显著。

1. 许多计算 FDR 的方法还会为每个基因报告一个 q value，它可以写作：

qi = (piN/i)c(N)

其中 c(N) 的定义同上。

1. 但按此方式直接计算得到的 q 值并不一定是 p 值的单调函数，因此 Benjamini 和 Yekutieli 引入了调整后的 q 值，定义为：

qi = min qk for k ≥ i

理解 FDR 及其正确用法，对于基因表达数据分析极其重要，因为绝大多数差异表达分析方法默认报告的都是 FDR 或 q 值。如果某个方法默认不提供 FDR，研究者也可以使用 R 中 stats 包的 p.adjust 等函数进行计算。

---

第七步 完

（下一节：第八步：通过功能富集分析探索机制——PDF 第 320 页）

第八步：通过功能富集分析探索机制

在识别出一组“显著”差异表达基因之后，下一步就是利用这份基因列表来探索研究对象背后的生物学机制。如果你对所研究的系统已经有相当了解，那么直接查看显著基因列表，并挑出一个或几个基因来描述它们可能参与的过程，往往并不困难。然而，这种做法很容易沦为所谓的 "bio-poetry" ——也就是基于零散经验知识拼接出的解释。它虽然听起来合理，却未必能够真正捕捉数据中的整体趋势。

因此，与其逐个讨论单个基因，不如进一步追问：哪些生物学过程在不同状态之间整体上发生了表达模式改变？

幸运的是，我们可以借助多种资源为基因赋予更高层级的功能注释（见第 13 章）。其中最常用的注释体系之一是 GO（Gene Ontology）。GO 是一个建立成熟、人工审查较充分的功能注释系统，它结合生物学文献和其他信息来源，为每个基因在以下三个方面分别赋予分类（见第 7 章）：

1. cellular component (CC)：基因产物所在的细胞区域或细胞外环境；
2. molecular function (MF)：基因产物在分子层面执行的主要功能，例如运输或结合；
3. biological process (BP)：基因产物所参与的整体生物学过程，例如细胞生长、信号传导或能量代谢。

在这三类信息中，GO biological process 往往最具解释力。除了 GO 之外，还有许多 pathway database 和其他 gene set 数据库也可用于类似分类；在分析方法上，它们的核心思想基本一致。我们真正关心的问题是：在区分研究人群的显著基因中，是否有一个或多个生物学过程（或通路）出现了过度代表（over-represented）？

这类方法统称为 gene set enrichment analysis，也常被称作 functional enrichment analysis。目前已有许多不同的分析策略、R packages 和 online tools。最常用的方法大体可分为两类：

• list-based methods：基于预先定义好的差异表达基因列表；
• rank-based methods：不只看显著基因，而是使用按某种显著性指标排序后的完整基因列表，例如 p value 或 q value。

List-Based Methods

当我们得到一组差异表达基因后，一个自然问题是：映射到某个特定功能类别的基因数，是否多于随机情况下的预期？

大多数回答这一问题的方法，都基于 Fisher's exact test。它也常被称为 hypergeometric test，不过这里之所以称为“exact”，是因为它并不依赖对超几何函数的近似。

Fisher's exact test 需要两个集合：

• 一个是被选中的差异表达基因集合；
• 另一个是背景集合（background set）。

对于 microarray，背景集合通常是芯片上所有出现的基因；对于 RNA-seq，背景集合则常常取整个基因组。

该方法的核心是：针对每一个待检验的功能类别（例如每一个 GO BP term），构造一个 contingency table。具体来说，设总共有 n 个基因，其中：

• a：显著且映射到该 BP term 的基因数；
• b：不显著但映射到该 BP term 的基因数；
• c：显著但不映射到该 BP term 的基因数；
• d：既不显著也不映射到该 BP term 的基因数。

此时，就可以用超几何分布来估计观察到这种分布是由随机造成的概率。

直观地说，如果基因组中有 10% 的基因被注释到某个 BP term，那么在显著基因集合中也看到大约 10% 映射到该术语，并不奇怪；但如果显著基因中有 20% 都映射到同一 BP term，那么这就可能具有统计学意义。

这种方法的价值，在于它能帮助我们避免过度解读某些功能类别的出现频率。不过，由于通常会同时检验很多个功能类别，因此这些 p values 仍然必须进行 multiple testing 校正。

基于（修正版）Fisher's exact test 的一个经典 list-based pathway enrichment 工具是 DAVID（Database for Annotation, Visualization and Integrated Discovery；da Huang et al. 2009a,b）。DAVID 的做法是在列联表中把 a 替换为 (a − 1)，从而让检验更保守一些。DAVID 是一个较为友好的 online tool，可以针对给定的输入基因列表和背景基因列表，测试多种 pathway annotation databases，包括 GO terms、Biocarta pathways 以及 KEGG（Kyoto Encyclopedia of Genes and Genomes）pathways。它既提供常见 microarray platforms 的背景基因集，也支持 RNA-seq 分析所需的 whole genome 背景集合。

另一类执行相似分析的工具是 R package topGO（Alexa et al. 2006）。除超几何检验外，topGO 还能使用专门设计的算法，显式考虑不同 GO terms 之间的依赖关系。一个典型例子是其中的 elim algorithm：它会把已经注释到某个显著富集节点（即 GO term）的基因，从该节点所有祖先节点中移除。这样做的结果是：分析会减少对笼统、上层细胞功能术语的偏重，而更强调 GO 层级结构中位置更低、功能定义更具体的术语。

Rank-Based Methods

list-based methods 的一个主要问题，是它们对“阈值设在何处”非常敏感。实际分析中常会发现，只要把显著性的 FDR cut-off 稍微调高或调低，结果就可能明显改变，因为在阈值上下附近，也许恰好聚集着一批来自同一功能类别的基因。

rank-based approaches 正是为了解决这一问题。它们不依赖单一显著性阈值，而是先按照某种显著性度量，对实验中所有基因进行排序，例如：

• p value
• q value
• t-statistic
• signal-to-noise metric
• 或其他可反映差异程度的指标

随后，rank-based 方法会检验：某一特定功能类别对应的基因（例如某个 GO BP term 或 KEGG pathway），是否在排序列表的两端出现了过度聚集。

这一思路最早的经典实现是 Gene Set Enrichment Analysis (GSEA)（Subramanian et al. 2005）。GSEA 既可以作为 online tool 使用，也提供 stand-alone Java program，可被 R 等其他编程语言调用。GSEA 使用 weighted Kolmogorov-Smirnov test 来计算某个 gene signature 的 enrichment score，随后再通过置换（permutations）来判断该富集分数是否显著（见第 13 章）。

另一个相关工具是 R package GSEAlm（Oron et al. 2008）。它通过线性模型先计算 p values，然后检验：映射到某个特定注释类别的基因，其 p value 分布是否不同于背景基因集合的 p value 分布。

---

第八步 完

（下一节：第九步：构建分类器——PDF 第 322 页）

Step 9：构建分类器

许多分析的最终目标是识别和探索驱动表型差异的生物学过程。然而，基因表达谱分析还有另一个常见应用，特别是在临床或转化医学场景中——利用数据构建分类模型，将新样本分配到研究中的某一表型组。要构建这样的分类器，首先要选择一组特征（在我们的例子中即基因），以区分不同的生物学类别，然后拟合模型的参数，使分类器能够基于所选特征集的基因表达来准确分类样本。正如本章前述的基因表达分析其他步骤一样，在特征选择、分类器训练与测试方面都有大量可选方案，其中涉及多种可用于分类的统计学和机器学习方法，但该领域尚未就哪种方法最优达成明确共识。也就是说，生物标志物（包括特征集和分类方法）需要使用独立数据集进行仔细验证，这一点是明确的。

成功且可重复的分类器开发的一个关键要素是从良好的实验设计开始。迄今为止已发表了成千上万的基因表达生物标志物，但大多数在首次发表的研究之外再未被使用过。尽管这些问题在前述实验设计讨论中已有部分涉及，但仍有一些额外的标准需要考虑。

首先，所分类的不同组之间应该有良好的平衡。如果研究的是一种罕见疾病，人群中只有 10% 患有该病，那么只需简单地得出"没有人患有该病"的结论，就能构建一个准确率达 90% 的分类器。因此，研究者应认识到这一点，要么尽量平衡各组，要么更明确地阐述成功标准。

其次，需要考虑同时设置训练集和测试集。在标准范式中，先在单个训练集上执行特征选择和参数拟合/算法训练，然后在独立的测试集（即先前未以任何方式用于特征选择或算法训练的样本集）上验证预测模型。

如果研究群体较小（例如处理罕见疾病），交叉验证是一种可接受的方法，但每次"折"的交叉验证都应重新在训练子集上执行特征选择和算法训练，然后重新在独立测试子集上检验方法性能。交叉验证方法面临的挑战在于，其结果不是一个单一的分类器，而是一组分类器，其中任何一个都难以与其他分类器进行直接比较。

开发分类器时还需要注意的另一件事是，训练集和测试集都需要有客观真实值；否则，方法无法得到有效训练，其性能也无法得到客观评估。另一种替代方案是反转测试集和训练集范式——使用多个独立的训练集，在每个训练集上学得一个新的分类器，然后在单个测试集上检验多个分类器之间的一致性和稳定性。在许多方面，这提供了一种更好的分类器成功度量指标，因为它说明了无论算法在何处、以何种方式拟合，对同一患者做出相同分类的可能性有多大。

分类器性能度量

在度量分类器准确率时，两个常用指标是敏感度（sensitivity）和特异度（specificity）（见第 7 章）。假设已开发了一个分类器，要在已知样本类别的独立数据集上对其进行检验。我们可以使用分类方法，对测试数据集中的每个样本做出判定，然后检查效果。假设有病例和对照，并希望对病例进行分类，则将病例视为阳性，对照视为阴性。那么，预测可分为真阳性（TP）和真阴性（TN），表示预测结果与真实分类一致；相应的，假阳性（FP）和假阴性（FN）则表示与真实分类不一致。由此可以定义敏感度，即真阳性率（TPR，有时也称命中率或召回率），为检出的 TP 占实际阳性总数的比例：

TPR = TP / (TP + FN)

进而定义特异度，即真阴性率（TNR），为检出的 TN 占实际阴性总数的比例：

TNR = TN / (TN + FP)

第三个有时有用的指标是精确度，或称阳性预测值（PPV），为检出的 TP 占所有判定为阳性总数的比例：

PPV = TP / (TP + FP)

最后，一个整合了这些概念的有用诊断图是受试者工作特征曲线（ROC 曲线，得名于二战时期的雷达探测研究），它将敏感度对特异度（或 FPR）作图。这一图形表示非常有用，因为大多数分类方法都包含可调节敏感度或特异度的参数，理解两者如何相互影响有助于决定如何调整预测模型。随机分类器的 TPR 和 FPR 相等，在 ROC 曲线上表现为对角线（图 10.10）。

特征选择

如前所述，分类器包含两个组成部分：特征集和分类算法。特征选择的目的是利用训练集中样本组之间的比较，识别一组能够区分这些组、并具有足够判别力来分类新样本的基因。

差异表达检验 或许并不令人意外，特征选择最常用的方法之一是使用第 7 步所述的差异表达分析统计学方法。差异表达分析的统计检验会识别在给定数据集中最能区分实验组的基因。虽然直觉上这很合理，但差异表达分析可能识别出大量高度相关的基因，从而对任何下游分类系统产生偏倚。真正需要的是一组在类间具有最大判别能力的特征，这些特征基于可能需要的多重模式来提供完整的类判别。幸运的是，有多种方法可以用于选择此类特征，包括下面介绍的最小冗余最大相关（mRMR）和预后标志物显著性分析（SAPS）。

特征选择面临的挑战在 Venet 等人（2011）的研究中得到了突出强调，他们将已发表的基因集与随机基因集进行比较，检验其将乳腺癌病例分为具有显著生存差异组的能力。随机基因签名能够胜过"显著"基因集这一发现表明，仅凭统计显著性来选择用于分类的基因集是不够的。

最小冗余最大相关（mRMR） 统计学显著的基因集通常包含大量高度相关的基因，因为这些基因共表达，但都代表相似的生物学过程。如果仅使用最显著的基因构建分类器，就有过度采样大型相关基因集的风险，从而错过可能有助于区分不同表型的生物学过程范围。

考虑这样一种情况：比较不同的类别并按基因的显著性排序，选择最显著的基因，然后将与它高度相关的基因从候选中去除，再选择下一个最显著的基因，继续去除相关基因并重复这一过程。这样得到的一组基因应兼具高预测性和相对独立性。正如你可能已认识到的，衡量相似性的方法有很多，包括使用 Pearson 相关系数和欧氏距离。

一种有助于在衡量相似性时平衡不同生物学过程表示的方法称为最小冗余最大相关（mRMR）（Ding and Peng 2005）。mRMR 使用互信息（一种非线性关联方法）同时识别最能区分类别的基因，并减少基因之间表达谱的潜在重叠。R 中 survcomp 包提供了 mRMRe 的实现（De Jay et al. 2013），mRMRe 包的并行版本可在 mRMRe 包中找到。

预后标志物显著性分析（SAPS） Venet 等人（2011）的研究检验了随机基因集预测生存的能力，其发现对许多基于选定基因集发表的预测方法提出了质疑。造成这一发现的部分原因可能是许多已发表分类器使用了相对大量的基因，以及选定特征集中基因之间（以及与基因组中其他基因之间）的相关性等。SAPS（Beck et al. 2013）是一种启发式方法，可以解决前述问题，适用于确定疾病预后或生存的预测因子，尽管一般流程也可以改编用于其他分类问题。

该方法基于三个独立 p 值——P_pure、P_random 和 P_enrichment——为候选基因集计算 SAPS 分数，每个 p 值通过一系列检验来估算：

1. P_pure 的计算方式是：首先使用 k-means 聚类（k = 2）根据选定基因集将患者分为两组，然后计算 log-rank p 值，以估计两组患者样本在生存方面没有差异的概率。
2. 接下来，选择与候选基因集大小相同的随机基因集，并按本列表步骤 1 所述进行检验，以评估它们将人群分为在生存方面存在差异的两组的效果。P_random 是 log-rank p 值至少与 P_pure 一样显著的随机基因集的比例。
3. P_enrichment 检查候选基因集与随机基因集的相对富集程度，以基于一致性指数确定高度预测性基因的相对富集。基因的一致性指数是指，对于数据集中随机选取的一对患者，表达该基因水平较高的患者的肿瘤比表达水平较低者预后更差（或更好）的概率。P_enrichment 的计算方式是：使用预排序基因集富集分析，确定候选基因集中与随机基因集相比，高或低一致性指数基因的富集程度（见本列表步骤 2 中选择的随机基因集）。富集显著性使用置换分析估算。
4. 最后，将这三个分数与候选基因集和预后之间关联的方向（direction，正向关联为"1"，负向关联为"-1"）结合起来，计算 SAPS 分数：

SAPS score = -log₁₀ max(P_pure, P_random, P_enrichment) × direction

SAPS 分数绝对值越大，三个 p 值共同决定的预后关联就越显著。通过置换基因，可以估算 SAPS 分数的统计显著性：生成 SAPS 分数的零分布，计算零分布中与候选基因集的 SAPS 分数绝对值至少同样大的基因集所占的比例。如果评估了多个候选基因集，则每个基因集的原始 SAPS p 值可用于生成相应的 SAPS q 值，即经多重检验校正后的 SAPS p 值。

mRMR 和 SAPS 等方法的价值在于，它们提供了一种在训练和测试算法之前检验基因集质量的方法。使用最优基因集可以大大增加分类方法表现良好的可能性，尽管对基因集加分类算法的组合进行额外验证仍然必不可少。

分类方法

选定候选基因集后，流程中的下一步是选择、训练和验证一个分类算法，用于将新样本分配到研究中的某一表型亚组。为此目的，存在大量从统计学和机器学习领域借鉴的分类方法，包括最近质心法、压缩最近质心法、高斯亚型分类模型、k 近邻法、支持向量机、随机森林、线性判别分析、二次判别分析、偏最小二乘、逻辑回归、神经网络等（Hastie et al. 2001, 2009; Haibe-Kains et al. 2012）（另见第 18 章）。

虽然每种方法的具体细节可能存在很大差异，但每种分类方法都代表一个数学函数，其参数被"拟合"（统计学中称为"估计"，机器学习中称为"学习"），其输入变量是特定样本的基因表达水平，输出是样本的亚组分配。虽然如果能就最佳分类方法提供指导将非常理想，但对文献的调查将表明，就使用何种最佳方法尚无科学共识（尽管许多论文声称某方法优于其他方法）。然而，有一点是公认的：这些方法的性能应按照下述方式得到严格测试和验证。

训练分类器时需要考虑的一个重要问题是 FP 和 FN 的相对成本。大多数方法优化整体性能，但有时存在对某一类成员的过度或不足识别所带来的真实成本。例如，在临床环境中，识别所有患有某种疾病的患者（提高敏感度）可能要好得多，即使这会冒着有一定水平 FP 识别（降低特异度）的风险。根据具体应用，我们可能认为重要的是最小化 FPR、FN 率、PPV 或其他参数。在大多数情况下，以相同的权重优化 FP 和 TP 以及 FN 和 TN。因此，在拟合模型之前，了解这一假设是否适当，或者是否存在某些识别错误比其他错误成本更高，是有用的。这一决定将有助于指导模型拟合和验证的方法选择。

预测模型的验证

预测分类器只有在准确且可重复的情况下才有用。例如，在乳腺癌中，已发表了成千上万个亚型分类器，但临床中使用的不到 10 个。这些预测模型失败的常见原因之一是其性能往往因方法学错误而被高估，这种现象称为过拟合。导致过拟合的原因有很多，下面将考察三种有助于避免过拟合的策略。

使用独立测试集进行群体水平预测验证 大多数分类器的性能基于其将一组样本划分为由其底层生物学或其他指标（如治疗反应或疾病生存）定义的亚组的能力。如果我们从一个有亚组信息的测试数据集开始，可以执行特征选择，然后训练方法并拟合模型参数。如果将同一分类器应用于该测试数据集，应该得到 100% 的分类准确率，但这不是一个公平的检验。恰当的问题是：该分类器应用于真正独立的测试集时，能否提供足够高的准确率而有实际用途。独立验证需要满足以下条件：

1. 从训练数据集出发，利用已知样本类别执行特征选择并拟合模型参数。
2. 使用一个在先前任何研究中都完全未用过的独立测试数据集，且类别标签保持盲法，然后应用所开发的分类模型将样本分配到各亚组。
3. 解盲样本类别标签，然后比较预测类别与真实类别，计算敏感度和特异度等指标，之后再报告结果。

这是一个相对简单且成熟的协议，但在使用时需要注意一些事项。首先，如果独立数据集是真正独立的，它将提供更可靠的答案。如果使用临床样本，尽量让训练样本和测试样本来自不同医院，由不同人员提取 RNA，在不同设施中测量表达量。如果能够生成与训练集以及彼此之间都独立的多个测试集，就更可能准确地估计性能。

转录组谱分析早期常见的一个错误是，将样本混合以进行特征选择，然后再将混合样本分开用于训练和测试。这种方法的问题在于，某些样本中的偏倚可能在特征选择时被捕获，从而潜在地高估最终的性能估计。需要牢记的重要一点是，训练集和测试集在流程的每一步都必须保持真正独立。

使用交叉验证进行群体水平预测验证 独立训练和测试集范式的一种变体，当可用样本数量有限时经常使用，是交叉验证。n 折交叉验证使用单个数据集，将其分成训练和测试集 n 次，然后对每次 n"折"重复训练和测试过程。没有划分初始数据集的绝对正确方法，也没有执行多少折的正确数量，但鉴于这种方法通常用于相对较小的数据集（如罕见疾病），常见的划分比例是训练集 vs. 测试集 90/10，且至少使用 10 折。该过程与独立验证模型一致：

1. 每折将数据集划分为独立的训练集和测试集。
2. 利用训练集中已知的样本类别，执行特征选择并拟合模型参数。
3. 使用拟合的模型对测试集中的样本进行分类。
4. 计算分类器的性能。
5. 将步骤 1-4 重复 n 折，同时汇总分类器及其方法的总体平均性能。

交叉验证是分类器评估中最常被误用的方法之一。一个常见错误是在特征选择中使用整个数据集，然后使用各折来训练模型并应用。问题在于，某些样本中的偏倚可能在特征选择过程中被捕获，然后膨胀整体方法的性能。因此，在每折中必须保持训练集和测试集独立。此外，应多次运行整个过程，以确保特定的分割不会偏倚结果。

交叉验证的问题是，最终没有单一的分类器可以报告。相反，有 n 个特征集和 n 个分类器，每折一个。有些团队取分类基因集的交集（有时是并集），然后将其报告为生物标志物集。有时，这个共识集被用于在整个数据集上训练算法，但没有坚实的理论支持使用这种方法，而且这个新分类器的整体性能在使用真正独立的测试集之前仍然未知。

使用独立训练集进行个体水平分配稳健性验证 使用训练集和一个或多个独立测试集验证分类器，是理解分类方法在群体水平上表现如何的非常好的方式，敏感度和特异度是群体水平性能的良好指标。然而，如果你是诊所里的患者，负责治疗的医生要用某种检测来确定诊断，那么群体水平的性能对你来说可能远不如该检测能否可靠且一致地将你分配到正确的治疗组重要。

在实践中，这种方法相对简单（Haibe-Kains et al. 2012; Beck et al. 2013）。从一个样本类别已知的多个独立训练集和单个独立测试集出发：

1. 对于每个训练集，识别一个能区分类别的基因集；或者，可以基于某些其他标准（如来自已知通路的代表性基因集或从文献中提取的基因集）从候选基因集出发。
2. 对于每个训练集和适当的基因集，拟合分类模型。
3. 每个分类模型独立用于预测测试集中样本的类别。
4. 不同分类器之间的一致性被用来估计方法的稳健性，衡量使用不同训练集构建的分类器对每个样本给出相同分类的频率。

在许多方面，这种稳健性评估符合我们对好的分类模型应该如何表现的理解。无论模型是在伦敦、巴黎、悉尼、东京还是波士顿学得的，都应该无关紧要。如果基因集和预测模型真正可靠，它们应该给出相同的分类。

在实践中，这种方法可能导致与交叉验证相同的困难——即没有单一的分类器。然而，考虑到一个好的分类器无论使用哪个数据集训练都将是高度一致的，人们希望这些分类器或多或少可以互换。无论如何，这种评估方式应被视为使用独立测试集验证的补充，而不是替代。

单细胞测序

尽管单细胞 RNA 测序（scRNA-seq）相对较新，但它特别令人兴奋，因为它允许对单个细胞的转录组进行分析。虽然最初的 scRNA-seq 实验分析的是单个或非常少量的细胞，但现在每次实验可以分析数千个细胞。目前使用的已发表方案包括 SMART-Seq（Ramskold et al. 2012）、CEL-Seq（Hashimshony et al. 2012）和 Drop-Seq（Macosko et al. 2015），以及 10x Chromium 系统等成熟的商业产品。

与传统的 DNA 微阵列和 RNA-seq（分析批量组织样本中的基因表达）不同，scRNA-seq 使用条形码从样本中的每个细胞创建测序文库，产生原始样本中单个细胞的表达谱。因此，scRNA-seq 不仅允许在不同表型之间比较表达，还可以定义细胞群，并研究表达变异和细胞异质性与表型的关联。已发表的 scRNA-seq 应用包括细胞类型及其多样性的识别和探索、基因表达随机性分析，以及跨细胞和细胞类型的基因调控网络估计。

scRNA-seq 中使用的处理步骤与 RNA-seq 类似。通常对序列读数执行质量控制，将质量读数比对到适当的参考序列（例如使用 Salmon（Patro et al. 2017）或 Kallisto（Bray et al. 2016）），并进行比对后质量控制。在这个阶段，通常使用 scRNA-seq 特有的方法进行标准化、基于表达识别亚组、差异表达检验和功能分析。这些方法包括用于标准化的 SCONE（Cole et al. 2018）、用于寻找细胞群的 Seurat（Butler et al. 2018）和 GiniClust/GiniClust2（Jiang et al. 2016; Tsoucas and Yuan 2018），以及提供综合分析的基于网络的处理平台如 Falco（Yang et al. 2017）和 ASAP（Chen et al. 2017; Gardeux et al. 2017）。

scRNA-seq 分析路径与批量 RNA-seq 分析分叉的原因是，存在一些 scRNA-seq 特有的技术干扰因素。例如，作为细胞分离和扩增结果产生的批次效应就是一个例子。但更显著的问题是基因"脱落"或稀疏性。据估计，每个细胞在任何给定时间都有多达 300,000 个 RNA 转录本。然而，在 scRNA-seq 中，每个细胞通常只记录几千（或数万）个读数。脱落发生是由于与 mRNA 计数实验相关的统计特性——我们偶然错过了一些转录本。虽然高表达转录本通常在 scRNA-seq 数据中表现良好，但中低表达转录本可能会被遗漏。"脱落"基因的整体模式产生了一个稀疏的细胞×转录本表达矩阵——观察到的计数值远少于预期。幸运的是，微生物组文献中已探讨了这个问题，scRNA-seq 使用的许多方法可以追溯到微生物组分析。

然而，这些问题意味着，截至本文写作时，在 scRNA-seq 分析的许多方面尚无共识。例如，在标准化方面，全局和缩放方法都被频繁使用，TPM 和 CPM 是常用指标。在比较表达方面，正在开发考虑细胞群产生的基因表达多模态性的方法。然而，许多已发表的研究仍然依赖 t 检验等成熟方法进行表型比较。

scRNA-seq 文献中现在广泛使用的一种新型可视化是 t 分布随机邻域嵌入（t-SNE）图（van der Maaten and Hinton 2008）。类似于 PCA 等方法，t-SNE 降维高维数据，但通过求解一个力图保持相似基因或细胞之间距离的目标函数来实现。虽然 t-SNE 已被证明非常有用，但它需要大量内存和计算资源，运行时间与分析的细胞或基因数量的平方成正比。

第10章 表达分析

10.12 Summary、Internet Resources、Further Reading 与 References

范围：PDF page 329 - PDF page 334；合并 Summary、Internet Resources、Further Reading、References。

---

Summary

本章为基因表达数据分析提供了一张路线图；而且由于这一领域变化极快，这张路线图天然不可能穷尽所有内容。自 20 世纪 90 年代早期通过 expressed sequence tags 测序开展首次全基因组范围表达分析以来，技术平台的变化、分析方法的进步，以及大量辅助数据的爆发式增长——例如众多物种的基因组序列及其基因注释——共同重塑了这一研究领域。

具体选择哪一种软件工具或分析方法，当然始终可以讨论；今天看来最优的方案，明天也可能被更新的方法取代。但良好的实验设计原则与扎实、可靠的分析实践并不会因此改变。

因此，与其把本章当成一本按步骤照做的 cookbook，不如把它理解为一张 roadmap：它为研究者指明一条更有可能获得成功、也更能增强结果可信度的分析路径。本着这样的初衷，我们希望本章概述的方法，能够成为你进入表达分析领域的一份有用导论与实践指南。

---

Internet Resources

---

Further Reading

以下保留原书题录英文原文；仅将题录后的说明性文字译为中文：

Brazma, A., Hingamp, P., Quackenbush, J. et al. (2001). Minimum information about a microarray experiment (MIAME) – toward standards for microarray data. Nat. Genet. 29 (4): 365–371. https://doi.org/10.1038/ng1201-365.

这篇奠基性论文确立了公开共享、且具充分注释的基因组数据集的基本原则，直到今天仍然经得起检验。MIAME 标准中提出的报告要求，在设计和撰写所有大规模实验时都应被认真考虑。

Conesa, A., Madrigal, P., Tarazona, S. et al. (2016). A survey of best practices for RNA-seq data analysis. Genome Biol. 17:13. https://doi.org/10.1186/s13059-016-0881-8.

Conesa, A., Madrigal, P. et al. (2016). Erratum to: A survey of best practices for RNA-seq data analysis. Genome Biol. 17 (1): 181. https://doi.org/10.1186/s13059-016-1047-4.

这篇论文（以及其勘误）很好地综述了 RNA-seq 数据分析方法，并给出了关于 best practices 的指导。尽管该领域仍在快速演进，但所谓“最佳实践”在很大程度上，仍然是建立在 DNA microarray 分析阶段积累下来的那些硬经验之上。

Ching, T., Huang, S., and Garmire, L.X. (2014). Power analysis and sample size estimation for RNA-seq differential expression. RNA 20 (11): 1684–1696. https://doi.org/10.1261/rna.046011.114.

对基因表达分析进行 power calculation 一向非常困难，主要原因在于不同基因之间的表达水平与方差差异都很大。这篇论文为如何估计统计功效，以及 RNA-seq 研究所需的样本量，提供了有价值的方法学指导。

Paulson, J.N., Chen, C.Y., Lopes-Ramos, C.M. et al. (2017). Tissue-aware RNA-seq processing and normalization for heterogeneous and sparse data. BMC Bioinf. 18 (1): 437. https://doi.org/10.1186/s12859-017-1847-x.

这篇论文描述了一条简洁的数据质量控制与标准化流程，可用于检查（部分）样本注释信息，并对来自 heterogeneous samples 的数据进行标准化处理。

Glass, K., Huttenhower, C., Quackenbush, J., and Yuan, G.C. (2013). Passing messages between biological networks to refine predicted interactions. PLoS One 8 (5): e64832. https://doi.org/10.1371/journal.pone.0064832.

基因调控网络推断的方法有很多。作者团队开发这一方法，是基于这样两点认识：其一，transcription factors 会调控基因表达；其二，不同 phenotypes 之间的网络结构本身也会发生变化。

Hung, J.H., Yang, T.H., Hu, Z. et al. (2012). Gene set enrichment analysis: performance evaluation and usage guidelines. Briefings Bioinf. 13 (3): 281–291. https://doi.org/10.1093/bib/bbr049.

尽管 gene set enrichment 方法仍在不断演化，这篇综述仍然较好地比较了多种方法，指出了各自的 strengths 和 weaknesses，并对 best practices 提供了富有启发性的建议。

---

References

以下参考文献题录按原书英文原文保留：

Alexa, A., Rahnenfuhrer, J., and Lengauer, T. (2006). Improved scoring of functional groups from gene expression data by decorrelating GO graph structure. Bioinformatics. 22 (13): 1600–1607. https://doi.org/10.1093/bioinformatics/btl140.

Anders, S. and Huber, W. (2010). Differential expression analysis for sequence count data. Genome Biol. 1 (10): R106. https://doi.org/10.1186/gb-2010-11-10-r106.

Beck, A.H., Knoblauch, N.W., Hefti, M.M. et al. (2013). Significance analysis of prognostic signatures. PLoS Comput. Biol. 9 (1): e1002875. https://doi.org/10.1371/journal.pcbi.1002875.

Benjamini, Y. and Hochberg, Y. (1995). Controlling the false discovery rate: a practical and powerful approach to multiple testing. J. R. Stat. Soc. Series B Methodol. 57 (1): 289–300.

Benjamini, Y. and Yekutieli, D. (2001). The control of the false discovery rate in multiple testing under dependency. Ann. Statist. 29 (4): 1165–1188. https://doi.org/10.1214/aos/1013699998.

Bolstad, B.M., Irizarry, R.A., Astrand, M., and Speed, T.P. (2003). A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics. 19. https://doi.org/10.1093/bioinformatics/19.2.185.

Bolstad, B.M., Collin, F., Simpson, K.M. et al. (2004). Experimental design and low-level analysis of microarray data. Int. Rev. Neurobiol. 60: 25–58.

Bray, N.L., Pimentel, H., Melsted, P., and Pachter, L. (2016). Near-optimal probabilistic RNA-seq quantification. Nat. Biotechnol. 34 (5): 525–527. https://doi.org/10.1038/nbt.3519.

Brettschneider, J., Collin, F., Bolstad, B.M., and Speed, T.P. (2008). Quality assessment for short oligonucleotide microarray data. Technometrics. 50 (3): 241–264.

Butler, A., Hoffman, P., Smibert, P. et al. (2018). Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat. Biotechnol. 36 (5): 411–420. https://doi.org/10.1038/nbt.4096.

Callow, M.J., Dudoit, S., Gong, E.L. et al. (2000). Microarray expression profiling identifies genes with altered expression in HDL-deficient mice. Genome Res. 10 (12): 2022–2029.

Chen, W., Gardeux, V., Meireles-Filho, A., and Deplancke, B. (2017). Profiling of single-cell transcriptomes. Curr. Protoc. Mouse Biol. 7 (3): 145–175. https://doi.org/10.1002/cpmo.30.

Cole, M.B., Risso, D., Wagner, A. et al. (2018). Performance assessment and selection of normalization procedures for single-cell RNA-seq. bioRxiv biorxiv.org/content/early/2018/05/18/235382.abstract.

De Jay, N., Papillon-Cavanagh, S., Olsen, C. et al. (2013). mRMRe: an R package for parallelized mRMR ensemble feature selection. Bioinformatics. 29 (18): 2365–2368. https://doi.org/10.1093/bioinformatics/btt383.

DeRisi, J., Penland, L., Brown, P.O. et al. (1996). Use of a cDNA microarray to analyse gene expression patterns in human cancer. Nat. Genet. 14 (4): 457–460.

Ding, C. and Peng, H. (2005). Minimum redundancy feature selection from microarray gene expression data. J. Bioinform. Comput. Biol. 3 (2): 185–205.

Dobin, A., Davis, C.A., Schlesinger, F. et al. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29 (1): 15–21. https://doi.org/10.1093/bioinformatics/bts635.

Eisen, M.B., Spellman, P.T., Brown, P.O., and Botstein, D. (1998). Cluster analysis and display of genome-wide expression patterns. Proc. Natl. Acad. Sci. USA. 95 (25): 14863–14868.

Gardeux, V., David, F.P.A., Shajkofci, A. et al. (2017). ASAP: a web-based platform for the analysis and interactive visualization of single-cell RNA-seq data. Bioinformatics. 33 (19): 3123–3125. https://doi.org/10.1093/bioinformatics/btx337.

Gautier, L., Cope, L., Bolstad, B.M., and Irizarry, R.A. (2004). affy—analysis of Affymetrix GeneChip data at the probe level. Bioinformatics. 20 (3): 307–315.

Golub, T.R., Slonim, D.K., Tamayo, P. et al. (1999). Molecular classification of cancer: class discovery and class prediction by gene expression monitoring. Science. 286 (5439): 531–537.

Haibe-Kains, B., Desmedt, C., Loi, S. et al. (2012). A three-gene model to robustly identify breast cancer molecular subtypes. J. Natl. Cancer Inst. 104 (4): 311–325. https://doi.org/10.1093/jnci/djr545.

Hashimshony, T., Wagner, F., Sher, N., and Yanai, I. (2012). CEL-Seq: single-cell RNA-seq by multiplexed linear amplification. Cell Rep. 2 (3): 666–673. https://doi.org/10.1016/j.celrep.2012.08.003.

Hastie, T., Tibshirani, R., and Friedman, J.H. (2001). The Elements of Statistical Learning: Data Mining, Inference, and Predictions. New York, NY: Springer.

Hastie, T., Tibshirani, R., and Friedman, J.H. (2009). The Elements of Statistical Learning: Data Mining, Inference, and Prediction, 2e. New York, NY: Springer.

Hegde, P., Qi, R., Abernathy, K. et al. (2000). A concise guide to cDNA microarray analysis. Biotechniques 29 (3): 548–550, 52–44, 56, passim.

da Huang, W., Sherman, B.T., and Lempicki, R.A. (2009a). Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 37 (1): 1–13. https://doi.org/10.1093/nar/gkn923.

da Huang, W., Sherman, B.T., and Lempicki, R.A. (2009b). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4 (1): 44–57. https://doi.org/10.1038/nprot.2008.211.

Ioannidis, J.P., Allison, D.B., Ball, C.A. et al. (2009). Repeatability of published microarray gene expression analyses. Nat. Genet. 41 (2): 149–155.

Irizarry, R.A., Bolstad, B.M., Collin, F. et al. (2003). Summaries of Affymetrix GeneChip probe level data. Nucleic Acids Res. 31 (4): e15.

Irizarry, R.A., Warren, D., Spencer, F. et al. (2005). Multiple-laboratory comparison of microarray platforms. Nat. Methods 2 (5): 345–350. https://doi.org/10.1038/nmeth756.

Ishmael, N., Dunning Hotopp, J.C., Ioannidis, P. et al. (2009). Extensive genomic diversity of closely related Wolbachia strains. Microbiology 155 (Pt 7): 2211–2222.

Jiang, L., Chen, H., Pinello, L., and Yuan, G.C. (2016). GiniClust: detecting rare cell types from single-cell gene expression data with Gini index. Genome Biol. 17 (1): 144. https://doi.org/10.1186/s13059-016-1010-4.

Johnson, W.,.E., Li, C., and Rabinovic, A. (2007). Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics. 8 (1): 118–127. https://doi.org/10.1093/biostatistics/kxj037.

Kahvejian, A., Quackenbush, J., and Thompson, J.F. (2008). What would you do if you could sequence everything? Nat. Biotechnol. 26 (10): 1125–1133. https://doi.org/10.1038/nbt1494.

Konstantinopoulos, P.A., Cannistra, S.A., Fountzilas, H. et al. (2011). Integrated analysis of multiple microarray datasets identifies a reproducible survival predictor in ovarian cancer. PLoS One 6 (3): e18202.

Lander, E.S., Linton, L.M., Birren, B. et al., International Human Genome Sequencing Consortium (2001). Initial sequencing and analysis of the human genome. Nature 409 (6822): 860–921. https://doi.org/10.1038/35057062.

Langmead, B. and Salzberg, S.L. (2012). Fast gapped-read alignment with Bowtie 2. Nat. Methods 9 (4): 357–359. https://doi.org/10.1038/nmeth.1923.

Langmead, B., Trapnell, C., Pop, M., and Salzberg, S.L. (2009). Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 10 (3): R25. https://doi.org/10.1186/gb-2009-10-3-r25.

Larkin, J.E., Frank, B.C., Gavras, H. et al. (2005). Independence and reproducibility across microarray platforms. Nat. Methods. 2 (5): 337–344. https://doi.org/10.1038/nmeth757.

Leek, J.T., Johnson, W.E., Parker, H.S. et al. (2012). The SVA package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. 28 (6): 882–883. https://doi.org/10.1093/bioinformatics/bts034.

Li, H. and Durbin, R. (2009). Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 25 (14): 1754–1760. https://doi.org/10.1093/bioinformatics/btp324.

Li, P., Piao, Y., Shon, H.S., and Ryu, K.H. (2015). Comparing the normalization methods for the differential analysis of Illumina high-throughput RNA-seq data. BMC Bioinf. 16: 347. https://doi.org/10.1186/s12859-015-0778-7.

Love, M.I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15 (12): 550. https://doi.org/10.1186/s13059-014-0550-8.

van der Maaten, L. and Hinton, G.E. (2008). Visualizing high-dimensional data using t-SNE. J. Machine Learn. Res. 9: 2579–2605. prlab.tudelft.nl/sites/default/files/vandermaaten08a.pdf.

Macosko, E.Z., Basu, A., Satija, R. et al. (2015). Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets. Cell. 161 (5): 1202–1204. https://doi.org/10.1016/j.cell.2015.05.002.

Michaels, G.S., Carr, D.B., Askenazi, M. et al. (1998). Cluster analysis and data visualization of large-scale gene expression data. Pac. Symp. Biocomput 1998: 42–53.

Nagalakshmi, U., Wang, Z., Waern, K. et al. (2008). The transcriptional landscape of the yeast genome defined by RNA sequencing. Science 320 (5881): 1344–1349. https://doi.org/10.1126/science.1158441.

Oron, A.P., Jiang, Z., and Gentleman, R. (2008). Gene set enrichment analysis using linear models and diagnostics. Bioinformatics. 24 (22): 2586–2591. https://doi.org/10.1093/bioinformatics/btn465.

Patro, R., Mount, S.M., and Kingsford, C. (2014). Sailfish enables alignment-free isoform quantification from RNA-seq reads using lightweight algorithms. Nat. Biotechnol. 32 (5): 462–464. https://doi.org/10.1038/nbt.2862.

Patro, R., Duggal, G., Love, M.I. et al. (2017). Salmon provides fast and bias-aware quantification of transcript expression. Nat. Methods. 14 (4): 417–419. https://doi.org/10.1038/nmeth.4197.

Paulson, J.N., Chen, C.Y., Lopes-Ramos, C.M. et al. (2017). Tissue-aware RNA-seq processing and normalization for heterogeneous and sparse data. BMC Bioinf. 18 (1): 437. https://doi.org/10.1186/s12859-017-1847-x.

Perou, C.M., Jeffrey, S.S., van de Rijn, M. et al. (1999). Distinctive gene expression patterns in human mammary epithelial cells and breast cancers. Proc. Natl. Acad. Sci. USA. 96 (16): 9212–9217.

Pop, M., Paulson, J.N., Chakraborty, S. et al. (2016). Individual-specific changes in the human gut microbiota after challenge with enterotoxigenic Escherichia coli and subsequent ciprofloxacin treatment. BMC Genomics. 17: 440. https://doi.org/10.1186/s12864-016-2777-0.

Quackenbush, J. (2005). Extracting meaning from functional genomics experiments. Toxicol. Appl. Pharmacol. 207 (2 Suppl): 195–199.

Ramskold, D., Luo, S., Wang, Y.C. et al. (2012). Full-length mRNA-seq from single-cell levels of RNA and individual circulating tumor cells. Nat. Biotechnol. 30 (8): 777–782. https://doi.org/10.1038/nbt.2282.

Robinson, M.D., McCarthy, D.J., and Smyth, G.K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26 (1): 139–140. https://doi.org/10.1093/bioinformatics/btp616.

Schena, M., Shalon, D., Davis, R.W., and Brown, P.O. (1995). Quantitative monitoring of gene expression patterns with a complementary DNA microarray. Science. 270 (5235): 467–470.

Simon, R., Radmacher, M.D., and Dobbin, K. (2002). Design of studies using DNA microarrays. Genet. Epidemiol. 23 (1): 21–36.

Spellman, P.T., Sherlock, G., Zhang, M.Q. et al. (1998). Comprehensive identification of cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization. Mol. Biol. Cell 9 (12): 3273–3297.

Subramanian, A., Tamayo, P., Mootha, V.K. et al. (2005). Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. USA. 102 (43): 15545–15550. https://doi.org/10.1073/pnas.0506580102.

Toker, L., Feng, M., and Pavlidis, P. (2016). Whose sample is it anyway? Widespread misannotation of samples in transcriptomics studies. F1000Res. 5: 2103. https://doi.org/10.12688/f1000research.9471.2.

Tsoucas, D. and Yuan, G.C. (2018). GiniClust2: a cluster-aware, weighted ensemble clustering method for cell-type detection. Genome Biol. 19 (1): 58. https://doi.org/10.1186/s13059-018-1431-3.

Venet, D., Dumont, J.E., and Detours, V. (2011). Most random gene expression signatures are significantly associated with breast cancer outcome. PLoS Comput. Biol. 7 (10): e1002240. https://doi.org/10.1371/journal.pcbi.1002240.

Venter, J.C., Adams, M.D., Myers, E.W. et al. (2001). The sequence of the human genome. Science. 291 (5507): 1304–1351.

Wen, X., Fuhrman, S., Michaels, G.S. et al. (1998). Large-scale temporal gene expression mapping of central nervous system development. Proc. Natl. Acad. Sci. USA. 95 (1): 334–339.

Wilson, C.L. and Miller, C.J. (2005). Simpleaffy: a BioConductor package for Affymetrix Quality Control and data analysis. Bioinformatics. 21 (18): 3683–3685.

Yang, A., Troup, M., Lin, P., and Ho, J.W. (2017). Falco: a quick and flexible single-cell RNA-seq processing framework on the cloud. Bioinformatics. 33 (5): 767–769. https://doi.org/10.1093/bioinformatics/btw732.