下一代测序 (NGS)/RNA

本指南旨在为没有分析下一代数据的经验的人提供一个易于理解的 RNA-seq 数据分析指南。但是，假设对 R、下一代测序程序和使用 UNIX shell 有基本的了解。这里描述的大多数步骤都在评论文章中概述。^[1]

• 它主要由 Matthew Young ([email protected]) 编写，目前仍在开发中。
• 病原体示例由 B. Usadel 提供，并使用了一组不同的工具。

RNA

转录组学简介。

对于运行 RNA-seq 的每个样本，您通常会收到一个包含数百万个短 (25-300bp) DNA 序列（称为读）和指示每个碱基调用置信度的质量分数的文件。但是，有一些重要的常见变化，这些变化取决于使用的平台和协议。这些包括但不限于

• 碱基空间或颜色空间
• 双端/配对或单端/非配对
• 链特异性或非链特异性

以下部分将详细介绍每种类型。

碱基空间与颜色空间

用于第二代 RNA 测序的两个主要平台由 Illumina 和 ABI Solid 生产。虽然两者都产生数百万个短读，但它们是以略微不同的方式测序和报告的。Solid 平台使用一种测序技术，该技术通过一次连接两个碱基对来生成读信息。读中的每个碱基对都通过两个（可能不同的）二核苷酸测序两次。为了利用这种测序化学，Solid 报告其读不是作为核苷酸序列，而是作为四种颜色的序列，其中每种颜色代表碱基之间的转换，称为“颜色空间编码”。^[2] 另一方面，Illumina 平台一次读取片段上的一个碱基，直到达到所需的读长度。然后，观察到的碱基被报告为输出。

结果是，分析 RNA-seq 数据的工具取决于用于生成短读数据的平台。虽然可以将颜色空间读转换为碱基空间（反之亦然），但这样做会导致数据中出现严重的偏差，应不惜一切代价避免。因此，读应保留其本机格式，应使用适当的工具分析它们。

双端/配对读

标准 RNA-seq 协议涉及逆转录 mRNA (cDNA) 的随机剪切，然后从片段的一端测序一个短“读”。这意味着只有片段的前 25-300bp 是已知的（取决于读的长度），而其余片段仍然未测序。事实上，由于片段化是随机的，每个片段的长度也是未知的，尽管通常会应用大小选择步骤。

虽然化学物质的精确度不足以对整个片段进行测序，但可以应用一个巧妙的技巧，即从片段的两端获取一个短读，从而产生一对短读，一个来自片段的一端，一个来自另一端。执行此操作的读称为双端读或配对读。双端读允许推断有关间隔序列的更多信息，并且对于从头转录组构建和检测结构变异（如插入缺失）特别有用。

链特异性读

RNA-seq 数据可以分为链特异性和非链特异性两种。如果数据是非链特异性的，则无法直接从序列中识别片段转录的链。此外，由于 RNA-seq 协议通常涉及形成双链 cDNA 以便于测序，因此返回的序列与来源 DNA 序列的反向互补序列的可能性与原始 DNA 序列的可能性相同。在实际应用中，这意味着虽然有一半的 reads 映射到正向链，一半映射到反向链，但这种映射不包含任何关于 RNA 从哪条链转录的信息。

另一方面，链特异性 RNA-seq 数据保留了链信息，从而可以识别 RNA 从哪条链转录。

典型工作流程

TopHat

Tophat 是一款工具，可以作为 RNA-Seq reads 的快速拼接连接映射器。Tophat 使用 Bowtie 将 RNA-Seq reads 映射到哺乳动物基因组，Bowtie 是一款超高通量短读比对器，利用 Burrows-Wheeler 索引。在 Bowtie 比对 reads 后，Tophat 会分析映射结果，以识别外显子之间的拼接连接。Tophat 是为数不多的几种能够有效处理内含子区域的工具之一，内含子区域在长度大于 50 bp 的 reads 中更为普遍。我们可以使用 HISAT 代替 tophat

参考数据集

人类

基因组参考联盟提供了一个人类参考基因组，该基因组由来自不同人群的几个个体构建而成。并且，通过修正现有的组装错误（包括但不限于在所使用个体之间错误混合单倍型结构）不断提高参考的质量。

CCDS

一致性 CDS (CCDS) 项目是一项合作努力，旨在识别一组核心的人类和小鼠蛋白质编码区域，这些区域是一致注释的且质量很高。注释是 NCBI 的 RefSeq、EBI 的 Ensembl 和桑格研究所的 Havana 之间的共识。

GENCODE

GENCODE 是 ENCODE 扩展项目的一个子项目，旨在注释整个人类基因组中所有基于证据的基因特征。

HAVANA

桑格研究所的 HAVANA 小组手动注释人类基因组，为基因位点的全部复杂性和可能不被自动化注释系统很好地满足的特征提供全面的注释。

RefSeq

RefSeq 是 NCBI 提供的一种参考序列注释服务。

UCSC 基因

UCSC 已知基因是通过完全自动化的过程构建的，该过程基于来自 Swiss-Prot/TrEMBL (UniProt) 的蛋白质数据以及来自 Genbank 的相关 mRNA 数据。

SEQwiki 内容转储

RNASeq 数据可以做什么

为了分析你的 RNASeq 数据，不同的目标需要不同的预处理步骤。

寻找差异表达基因

如果你想寻找差异表达基因，许多 BioConductor 软件包要求你不要将你的读取计数转换为例如 rpkm，而要使用原始读取计数。原因是统计模型使用原始读取计数数据。（读取次数越少，短噪声越多，而较低的 rpkm 值仍可能与许多读取相关联，例如在非常长的基因情况下）。请注意，即使你使用 RNASeq 数据，你仍然需要生物学重复。

此外，你将希望使用为 RNAseq（即计数数据）设计的软件包，例如DESeq，bayseq，NBPSeq 或 edgeR，它们都使用负二项式模型。

一个典型的流程可能如下所示

1. 对你的数据进行质量控制
2. 读取比对（参见说明）
3. 计算每个特征的读取次数
4. 统计分析

读取计数/汇总

在将读取比对到你的基因组后，你需要汇总这些读取。这可能是一个关键步骤，因为有许多不同的汇总方法。一种方法可以是十分严格，只统计完全明确的读取，但是在下游分析中需要考虑到这一点。可以使用例如 HTSeq count 或 Bioconductor Iranges 模块来汇总数据

统计分析

在测试差异表达时，你应该使用一个可以对计数数据进行建模的软件包。在 Bioconductor 中，可以选择 DESeq, edgeR 和 bayseq。

差异外显子使用

除了基因表达的总体变化外，还可能存在等位基因或特定外显子的差异丰度。可以使用 DEXSeq 来测试差异丰度外显子，另请参见 Genome Research 论文。

类别富集（例如 GO 富集）

与所有富集分析一样，你发现的类别可能会受到检测能力的影响。使用 RNASeq，检测差异表达的能力（在相同类型的错误下，例如相同的调整后的 p 值）对于计数更多的基因来说更高。例如，如果所有光合作用基因都高度表达，你更有可能发现其中任何一个基因发生差异表达，而不是其他基因。现在，如果你发现光合作用基因在你的差异表达基因列表中富集，这可能仅仅是因为它们更容易被发现。

Bioconductor 中的Goseq 软件包的作者指出，更长的基因往往会产生更多的计数，并且他们提供了进行 GO 类别富集分析的软件，该软件可以根据基因长度或直接根据表达强度进行调整。

估计表达高度

一个典型的流程可能如下所示

1. 读取比对
2. 计算每个特征的读取次数
3. 根据基因长度进行归一化

质量控制 (QC)

有很多方法可以对测序数据进行质量检查。一些流行的选择是：FastQC， DNAA 或者你可以在 R 中进行操作或使用 HTSeq

尽管可以使用更复杂的质量指标，但一个基本的检查是，序列组成在读取的长度上没有太大变化，并且质量得分没有过低，这是一个良好的起点。这些检查（以及其他一些检查）可以通过将 fastq 文件加载到 fastQC 程序中来执行。

比对读取

除非你的目标是进行从头转录组组装，否则任何分析的第一步都是将你的数百万条短读取比对到某种参考。这通常是提取 RNA 的物种的参考基因组。有许多不同的比对器可以执行短读取比对到参考（列表可在 http://en.wikipedia.org/wiki/List_of_sequence_alignment_software 和 http://seqanswers.com/wiki/Special:BrowseData/Bioinformatics%20application?Bioinformatics_method=Mapping 中找到）。每一个都有特定的优点和缺点，通常涉及速度和敏感性之间的权衡。在本指南中，我们将使用比对器 BOWTIE，因为它正在积极开发、速度快、使用广泛，并支持基本空间和颜色空间读取。

作为起点，我假设你拥有以下文件

对于单端读取

sample.fa 或 sample.cfa（cfa 用于颜色空间，fa 用于基本空间）

对于双端读取

sample_pair1.fa 和 sample_pair2.fa 或 sample_pair1.cfa 和 sample_pair2.cfa（cfa 用于颜色空间，fa 用于基本空间）

获取参考

为了将短序列比对到参考基因组，该基因组必须被转换为BOWTIE可以使用的索引。一些常见的基因组的预构建索引可以从BOWTIE网站下载。如果您的基因组不可用，您将需要从包含参考基因组的fasta文件中自己构建索引，这些文件可以从UCSC 通过FTP获取。可以使用以下命令完成此操作：

   bowtie-build reference.fa reference_name

参数“reference_name”是用于从现在开始引用此参考基因组的唯一标识符。如果您的序列为颜色空间，则通过添加-C命令构建颜色空间索引：

   bowtie-build -C reference.fa reference_name

此命令将输出6个文件，分别命名为reference_name.1.ebwt、reference_name.2.ebwt、reference_name.3.ebwt、reference_name.4.ebwt、reference_name.rev.1.ebwt和reference_name.rev.2.ebwt。

无论您如何获取索引，为了使BOWTIE能够使用它，上述六个文件都需要放在BOWTIE索引目录中。如果您不确定您的索引目录是什么，它由环境变量BOWTIE_INDEXES指向，因此：

   echo $BOWTIE_INDEXES

将显示您应该放置索引文件的路径。

如果您希望使用除了BOWTIE之外的其他比对工具，您也需要从参考基因组构建索引。有关更多信息，请参阅您首选比对工具的文档。

将序列比对到参考基因组

将参考基因组构建成BOWTIE可以使用的索引后，我们现在想将我们的数据比对到这个索引。BOWTIE提供了丰富的命令行选项，可用于调整比对算法及其处理输入/输出的方式。这些命令行选项在BOWTIE手册中进行了详细描述，但有一些标志在RNA-seq分析中经常使用，值得在这里提及。

--sam或-S告诉BOWTIE以SAM格式而不是BOWTIE格式输出比对结果。SAM正迅速成为报告短序列比对的标准，并得到许多下游分析工具的支持。除非有充分的理由，否则应始终指定此标志。

--best标志告诉BOWTIE保证它报告的比对在所有找到的匹配中与参考基因组的错配数量最少。它还执行一些其他理想的功能，例如消除链偏向性。 这些好处的代价是--best稍微慢一些，但在几乎所有情况下，这种差异都可以忽略不计。请注意，--best不适用于双端序列。除非速度是一个主要考虑因素，否则应启用此标志。

每个碱基都有一个与之相关的质量得分，该质量得分以PHRED量表 (http://en.wikipedia.org/wiki/Phred_quality_score) 报告，其中较低的得分意味着对碱基调用准确性的信心较低。对于每个候选比对，BOWTIE会将与参考基因组不匹配的碱基的质量得分加起来。任何错配质量得分总和大于-e的匹配位置都被认为是无效的，不会被报告。-e的默认值为70，是在序列长度较短（~30bp）时优化的，但不适用于当今常用的较长的序列长度。此外，任何来自参考基因组的真正生物学变异（例如SNP）理论上都应该具有很高的质量得分。因此，具有SNP的序列在错配质量得分总和指标上的得分将非常低。出于所有这些原因，建议用户将-e提高到默认值以上，除非已知参考基因组是RNA生物学来源的极佳代表，并且序列中的错误数量很少。

-p标志设置要使用的同时线程数。由于短序列比对实际上是无限可并行的，因此将其设置为可用的CPU核心数。

考虑到这些因素，我们使用以下命令将短序列数据映射到参考基因组：

   bowtie -p 8 --sam --best reference_name sample.fa aligned_reads.sam

如果使用颜色空间，则需要添加-C标志。

   bowtie -p 8 --sam -C --best reference_name sample.cfa aligned_reads.sam

如果您的数据是双端的，则需要使用-1和-2标志指定两个配对文件。

   bowtie -p 8 --sam reference_name -1 sample_pair1.fa -2 sample_pair2.fa aligned_reads.sam

更复杂的比对

正在测序的cDNA片段不是来自基因组，而是来自转录组。转录组是通过组合来自基因组的外显子形成的，这就是为什么将序列映射到基因组是一个很好的近似。但是，这样做会丢失任何跨越外显子-外显子边界的序列的映射能力。序列越长，这个问题就越严重，因为序列更有可能跨越边界而无法映射。根据所需的后续分析，这种在外显子边界处的覆盖丢失可能是或可能不是一个问题。

外显子连接库

可以从已知的注释外显子构建外显子连接库，然后尝试将那些无法映射到该序列的序列映射到库。这种方法只能捕获跨越已知注释的序列，限制了其用途，并使结果偏向于注释良好的基因/基因组。为此，需要构建一个新的参考基因组，其中包含参考基因组以及新的外显子-外显子连接序列。即使您只将无法比对到基因组的序列映射到库，也要确保原始基因组参考仍然包含在参考基因组中，以便在确定比对时，序列可以在所有可能的映射位置之间竞争。您还需要决定从外显子-外显子边界两侧取多少序列。从每个外显子取的序列量必须小于序列长度，否则任何靠近外显子边界但没有超过边界的序列都将映射到两个位置，即外显子连接库和基因组本身。这意味着，如果您打算修剪您的序列，那么从外显子连接两侧取的序列量必须小于修剪后的序列的长度。最后，您需要决定要考虑哪些外显子-外显子连接组合。您只考虑基因内、染色体内的剪接，而不考虑外显子重新定位？所有这些情况在不同的样本中以不同的频率发生，您考虑的可能性越多，计算复杂度就越大。

要构建连接库，您首先需要一个基因组注释。这些可以从UCSC表格浏览器下载。然后，您可以编写自己的工具来创建一个包含所有外显子连接的fasta文件，或者使用现有的工具，例如USeq软件包中包含的“Make Splice Junction Fasta”应用程序。 一旦您拥有包含外显子-外显子连接库的fasta文件，您就需要将其与参考基因组的fasta文件合并。

   cat reference.fa junctions.fa >reference_and_junctions.fa

然后，如上所述构建新的bowtie索引：

   bowtie-build reference_and_junctions.fa junction_lib

最后，您以与参考基因组相同的方式将序列映射到库。例如：

   bowtie -p 8 --sam --best junction_lib sample.fa aligned_reads.sam

其他选项

如果这种方法仍然无法映射合理的序列数量，那么还有其他替代方法可以探索，即使计算成本更高。例如，您可以尝试使用“从头”剪接连接查找器从数据本身估计剪接连接。有很多工具尝试这样做，一些例子包括TopHat, SplitSeek, PerM, 和SOAPsplice（以前称为SOAPals）。 转录组的完全从头组装也是一种可能性，尽管为此工作良好需要非常高的覆盖率。双端数据对此任务非常有益，因为片段的任一端映射到不同的外显子，是外显子连接的非常强烈的证据。在综述文章中，您可以找到此类工具的列表（见表格1）。

最近发表在基因组生物学上的论文进一步提出了一种自动框架，可以根据数据集的关键特征（例如参考基因组、SNP 率、读长等）轻松地对比和选择比对器。Teaser 可作为在线工具使用，或者您可以下载并自定义自己的版本。

差异表达

表达分析的一个常见用途是在两个实验条件（例如野生型与敲除）之间寻找基因或其他感兴趣对象的表达水平差异。为此，我们需要将数据从基因组坐标的读序列列表转换为计数表。我们在这里采用的策略是使用 Rsamtools 包将短读序列加载到 R 中，然后计算与某些注释对象重叠的读序列数量，该对象通常类似于从 UCSC 下载的基因集合。转换后，可以执行测试以寻找表达水平的统计学显著差异。

读序列的汇总

使用 SAMtools 压缩比对的读序列

为了能够将数百万个比对的读序列加载到 R 中的内存中，我们需要创建人类可读的 SAM 输出的二进制压缩版本。该文件将包含所有相同的信息，但内存占用更小，并且可以快速搜索。要创建此类文件，我们首先需要安装 SAMtools。 接下来，我们需要使用 fasta 文件构建参考索引。执行以下操作：

   samtools faidx reference.fa

这将创建一个文件 reference.fa.fai。接下来，我们将 SAM 转换为 BAM。BAM 文件包含与 SAM 文件相同的所有信息，但经过压缩，使其更加节省空间且可搜索。

   samtools import reference.fa.fai aligned_reads.sam aligned_reads.bam

最后，我们需要创建读序列索引，以便可以快速搜索它们。为此，我们首先需要对 BAM 文件进行排序。

   samtools sort aligned_reads.bam aligned_reads_sorted

这将创建 aligned_reads_sorted.bam，我们现在对其进行索引。

   samtools index aligned_reads_sorted.bam

这将创建索引文件 aligned_reads_sorted.bam.bai。

如果您没有参考的原始 fasta 文件，因为您从 BOWTIE 网站下载了预构建的索引（或者因为您在制作自己的索引后丢失了 fasta 文件），您可以通过运行以下命令重建源 fasta 文件。

   bowtie-inspect reference_name>reference.fa

在 R 中使用 BAM 文件

我们的目标是使用 R 根据每个样本的基因对读序列进行汇总。为此，我们将使用我们新创建的短读序列的压缩表示（已排序、已索引的 BAM 文件）。

获取基因信息

我们要做的第一件事是定义基因在染色体坐标中的位置。为此，我们使用 GenomicFeatures 包。 该包允许我们使用以下命令从 UCSC 基因组浏览器下载基因信息

   library(GenomicFeatures)
   txdb=makeTranscriptDbFromUCSC(genome='hg19',tablename='ensGene')

可以使用各种基因组和基因 ID，但作为示例，我们将使用最新的 human 基因组和 ENSEMBL 基因 ID。变量“txdb”现在包含了我们所需的所有信息，但为了使用它，我们需要进行一些处理。

   tx_by_gene=transcriptsBy(txdb,'gene')

这将生成一个 GRangesList 对象，它是一个 GRanges 对象列表，其中每个 GRanges 对象都是一个基因，条目是其转录本的基因组坐标。我们将使用 countOverlaps 函数将此对象与包含短读序列的对象重叠，从而计算出哪些读序列与哪些基因重叠。

虽然我们选择通过包含落入基因内的所有读序列来进行汇总，但该过程对于任何 GRanges 或 GRangesList 对象都将相同。例如，如果您希望只包含与外显子重叠的读序列，您可以创建一个不同的 GRangesList 对象。

   ex_by_gene=exonsBy(txdb,'gene')

在 R 中汇总读序列

首先，我们将 aligned_reads_sorted.bam 文件加载到 R 中。

   library(Rsamtools)
   reads=readBamGappedAlignments("aligned_reads_sorted.bam")

检查注释和读序列的兼容性

在我们尝试将读序列与注释进行比较之前，我们需要进行一些检查，以确保一切正常。如果您的 RNA-seq 数据没有附带链信息，那么我们无法知道读序列是从哪个链转录的。但是，比对过程会将其比对到一个链或另一个链（理论上，两者可能性相同），因此读序列将被分配一个人工链。当我们计算与基因（或其他特征）重叠的读序列数量时，只有比对到基因所在的链的读序列才会被计数，大约一半的读序列会丢失。为了避免这种情况，我们需要将读序列的链值设置为“*”（未知）。

此外，比对软件使用的染色体名称（最终由参考基因组 fasta 文件中的染色体名称决定）通常与注释中给出的染色体名称不同。为了使比较函数正常工作，这些名称需要转换为相同的命名约定。通过列出读序列和注释的全部染色体，可以轻松检查名称是否匹配。请记住，我们的注释数据存储在“tx_by_gene”中，我们的短读序列存储在“reads”中。

   #The annotations have chromosomes called
   names(seqlengths(tx_by_gene))
   #The reads have chromosomes called
   as.character(unique(rname(reads)))

如果染色体名称相同（或对于您关心的染色体名称相同），则不需要进行名称转换。另一方面，如果它们不同，我们需要更改读序列或注释的命名约定。事实证明，更改读序列的名称通常更容易，但无论哪种情况，过程都是相同的。以下示例可以很好地说明这一点。假设您有以下情况

   #The annotations have chromosomes called
   > names(seqlengths(tx_by_gene))
   [1] "chr1" "chr10" "chr11" "chr12" "chr13"
   [6] "chr13_random" "chr14" "chr15" "chr16" "chr17"
   [11] "chr17_random" "chr18" "chr19" "chr1_random" "chr2"
   [16] "chr3" "chr3_random" "chr4" "chr4_random" "chr5"
   [21] "chr5_random" "chr6" "chr7" "chr7_random" "chr8"
   [26] "chr8_random" "chr9" "chr9_random" "chrM" "chrUn_random"
   [31] "chrX" "chrX_random" "chrY" "chrY_random"

   #The reads have chromosomes called
   >as.character(unique(rname(reads)))
   [1] "10.1-129993255" "11.1-121843856" "1.1-197195432" "12.1-121257530"
   [5] "13.1-120284312" "14.1-125194864" "15.1-103494974" "16.1-98319150"
   [9] "17.1-95272651" "18.1-90772031" "19.1-61342430" "2.1-181748087"
   [13] "3.1-159599783" "4.1-155630120" "5.1-152537259" "6.1-149517037"
   [17] "7.1-152524553" "8.1-131738871" "9.1-124076172" "MT.1-16299"
   [21] "X.1-166650296" "Y.1-15902555"

因此，我们需要将读序列染色体名称“NO.1-Length”转换为注释名称“chrNO”。

   new_read_chr_names=gsub("(.*)[T]*\\..*","chr\\1",rname(reads))

如果您不熟悉正则表达式，请参阅 R 中 gsub 的帮助文件。new_read_chr_names 现在将包含读序列染色体名称，已转换为与注释对象 tx_by_genes 相同的格式。

现在，我们可以通过从每个读序列对象构建 GenomicRanges 对象来同时修复染色体名称和链问题。如果您有无链 RNA-seq 数据，并且需要转换染色体，我们会运行以下命令

   reads=GRanges(seqnames=new_read_chr_names,ranges=IRanges(start=start(reads),end=end(reads)), 
     strand=rep("*",length(reads)))

如果我们只是想转换染色体名称。

       reads=GRanges(seqnames=new_read_chr_names,ranges=IRanges(start=start(reads),end=end(reads)), 
     strand=strand(reads))

如果我们只是想使每个读序列在链方面都模棱两可。

   reads=GRanges(seqnames=rname(reads),ranges=IRanges(start=start(reads),end=end(reads)),
   strand=rep("*",length(reads)))

请注意，如果您将读序列比对到外显子连接库，则每个外显子连接将具有自己的“染色体”。如果您希望这些连接读序列包含在汇总中，则必须将它们中的每一个都转换为基因组坐标。由于每个读序列来自两个不同的基因组位置（外显子连接的任一侧），因此您必须决定如何为每个读序列分配基因组坐标。

计算读序列数量

最后，我们获得了与每个基因（或您感兴趣的其他任何内容）重叠的读序列数量。

   counts=countOverlaps(tx_by_gene,reads)

“counts”现在将包含一个数值向量，其中第 i 个条目是与“tx_by_gene”中第 i 个基因重叠的读序列数量。

差异表达测试

由于我们的目标是比较条件，因此我们将拥有多个读序列泳道，每个条件可能有多个泳道。使用上一节中概述的程序，我们可以计算与每个实验条件下每个重复的感兴趣特征（例如基因）重叠的读序列数量。接下来，我们将它们组合成一个计数表。

   toc=data.frame(condition1_rep1=counts1.1,condition1_rep2=counts1.2,
     condition2_rep1=counts2.1,condition2_rep2=counts2.2,stringsAsFactors=FALSE)
   rownames(toc)=names(tx_by_gene)

等等，只要有条件和重复，就重复多少次。这里的约定是 countsn.m 是一个向量，包含条件 n 中重复 m 中与 tx_by_gene 给出的基因重叠的读取次数。

归一化

已经证明，少数高表达基因会消耗大量的总序列。由于这种情况会在泳道和实验条件之间发生变化，以及文库大小，因此在测试差异表达时，需要进行某种样本间归一化。归一化的选择与用于确定任何基因在条件之间是否存在显著差异表达 (DE) 的检验无关。例如，分位数归一化会产生非整数计数，使得基于计数数据的检验（例如广泛使用的泊松模型或负二项式模型）变得不适用。我们选择使用比例因子归一化方法，因为它保留了数据的计数性质，并且已被证明是提高 DE 检测的有效方法。

为了进行归一化和差异表达检验，我们将使用 R 包 edgeR， 尽管还有其他选择。我们现在可以使用 TMM 方法计算归一化因子。

   library(edgeR)
   norm_factors=calcNormFactors(as.matrix(toc))

统计检验

接下来，我们需要创建一个 edgeR 使用的 DGE 对象。使用 TMM 方法计算的比例因子通过按归一化因子加权文库大小（然后在统计模型中用作偏移量）来合并到统计检验中。

   DGE=DGEList(toc,lib.size=norm_factors*colSums(toc),group=rep(c("Condition1","Condition2"),c(2,2)))

组变量标识内容表中的哪些列来自哪个实验条件或“组”。为了进行显著性统计检验，我们首先估计公共离散参数

   disp=estimateCommonDisp(DGE)

最后，我们计算基因是 DE 的 p 值

   tested=exactTest(disp)

基因集检验 (GO)

为了使用 RNA-seq 数据准确地检验 DE 基因中基因集的过度表达，我们需要使用一种方法来考虑这种技术的特定偏差。GOseq 包 是一种用于在执行 GO（和其他基于基因集的检验）分析时解释某些 RNA-seq 特定偏差的方法。

首先，我们必须格式化 edgeR 的输出以供 goseq 读取。我们将任何 Benjamini-Hochberg FDR 小于 .05 的基因称为 DE。

   library(goseq)
   genes = as.integer(p.adjust(tested$table$p.value, method = "BH") < 0.05)
   names(genes) = row.names(tested$table)

接下来，我们计算概率加权函数，校正长度偏差，这是所有形式的 RNA-seq 数据中存在的技术偏差。

   pwf=nullp(genes,'hg19','ensGene')

如果我们想校正总读取计数偏差，我们会使用每个基因的计数数量来计算 pwf，如下所示

   pwf=nullp(genes,bias.data=rowsum(counts[match(names(genes),rownames(counts))]))

最后，我们计算每个 GO 类别在 DE 基因中过度表达的 p 值。

   GO.pvals=goseq(pwf,'hg19','ensGene')

示例 1：差异表达和 GO 术语分析：李前列腺癌数据集

本节提供了一个易于理解的示例，用于说明上述分析流程。

描述

该数据集比较了经过和未经过类固醇激素雄激素处理的前列腺癌 LNcap 细胞系。 测序使用 Illumina GA I 完成，产生 36 bp 的单端非链读取。机器的输出是 7 个文件（每个来自不同的测序泳道）

untreated1.fa

untreated2.fa

untreated3.fa

untreated4.fa

treated1.fa

treated2.fa

treated3.fa

请注意，此数据集有点不寻常，因为读取中的质量得分缺失。因此，我们在进行分析时需要牢记这一点。

质量控制

这是已发布的数据，所以这是一个相当不错的质量控制，人们希望……

序列比对

构建参考

流程中的第一步是将所有读取比对到参考序列。由于该数据取自人类 LNcap 细胞，因此人类基因组的最新版本是一个显而易见的选择。我们已使用所有标准选项安装了 BOWTIE（版本 0.10.0）。从 BOWTIE 网站可以获得人类索引的预构建副本，但是，为了说明从头开始构建基因组，我们改为从 UCSC 下载基因组的 .fa 文件。我们创建一个包含 7 个 RNA-seq 数据文件和从 加州大学圣克鲁兹分校 下载的 chromFA.tar.gz 文件的工作目录。文件 chromFA.tar.gz 包含所有人类染色体的序列，包括未分配的重叠群。为了制作 BOWTIE 索引，我们需要将它们连接到一个文件中。我们将从我们的 fasta 文件中排除所有重叠群。

   tar -zxvf chromFA.tar.gz

我们只想要 chr1-22.fa、chrX.fa、chrY.fa 和 chrM.fa，所以删除其他所有内容

   rm chr*_*.fa

现在，我们将所需的文件连接起来（最好将参考序列分别保存在每个染色体一个文件和每个基因组一个文件格式中，尽管 BOWTIE 索引可以从任何一种格式中生成）。

   cat chr*.fa>hg19.fa

我们构建 BOWTIE 索引并将其命名为 hg19。

   bowtie-build hg19.fa hg19

并将 BOWTIE 索引文件移动到 BOWTIE 可以找到它们的适当位置。

   mv *.ebwt $BOWTIE_INDEXES

比对读取

构建完人类基因组的 BOWTIE 索引后，我们现在继续将每个泳道的读取映射到索引。我们想使用 --best 和 --sam 标记。此时，我们回想起我们的数据缺少读取的质量信息。因此，我们使用 -v 3 选项，该选项在将读取比对到基因组时忽略质量得分。

   for src_fastafile in *treated*.fa
   do
     bowtie -v 3 -p 8 --best --sam hg19 ${src_fastafile} ${src_fastafile%%.fa}.sam
   done

for 循环仅用于方便，等同于编写

   bowtie -v 3 -p 8 --best --sam hg19 untreated1.fa untreated1.sam
   bowtie -v 3 -p 8 --best --sam hg19 untreated2.fa untreated2.sam
   bowtie -v 3 -p 8 --best --sam hg19 untreated3.fa untreated3.sam
   bowtie -v 3 -p 8 --best --sam hg19 untreated4.fa untreated4.sam
   bowtie -v 3 -p 8 --best --sam hg19 treated1.fa treated1.sam
   bowtie -v 3 -p 8 --best --sam hg19 treated2.fa treated2.sam
   bowtie -v 3 -p 8 --best --sam hg19 treated3.fa treated3.sam

BOWTIE 将输出 7 个 SAM 文件，其中包含比对的读取。为了统计比对读取的比例，我们在 shell 中运行以下命令（此信息也由 BOWTIE 直接报告，但能够自己计算它很有用）。

   awk '$3!="*"' untreated1.fa|wc -l
   wc -l untreated1.fa

第一个命令打印 BOWTIE 输出文件中映射的读取数量，第二个命令输出输入文件中读取数量的 4 倍（因为 fasta 格式是每个读取 4 行）。

我们将忽略那些跨越外显子-外显子边界的读取，并继续进行分析。

读序列的汇总

转换为 BAM

为了将比对后的reads汇总成基因，我们首先需要将BOWTIE输出的SAM文件转换为压缩的BAM索引格式。首先，我们需要为人类基因组创建samtools格式的索引。为此，我们需要参考序列的fasta文件（应该与用于创建比对reads索引的fasta文件相同），我们已经有了，但我们会从BOWTIE索引重建它。

   bowtie-inspect hg19>hg19.fa

现在我们构建samtools索引。以下命令会生成一个.fai文件，我们可以用它将SAM文件转换为BAM文件。

   samtools faidx hg19.fa

接下来，我们将所有7个SAM文件转换为BAM文件。

   for SAM in *.sam
   do
     #Convert to BAM format
     samtools import hg19.fa.fai ${SAM} ${SAM%%.sam}.bam
     #Sort everything
     samtools sort ${SAM%%.sam}.bam ${SAM%%.sam}_sorted
     #Create an index for fast searching
     samtools index ${SAM%%.sam}_sorted.bam
     #Delete temporary files
     rm ${SAM%%.sam}.bam
   done

同样，bash for循环只是为了方便，上面等同于对所有7个文件运行以下命令。

   samtools import hg19.fa.fai untreated1.sam untreated1.bam
   samtools sort untreated1.bam untreated1_sorted
   samtools index untreated1_sorted.bam
   rm untreated1.bam

在R中进行处理

接下来我们需要将排序后的、索引的BAM文件加载到R中。由于我们使用的是非链特异性RNA-seq，我们需要将每个read的链指示器设置为模糊（通过将其设置为"*"）。在启动R后，我们运行：

   library(Rsamtools)
   #Create a list of bam file object containing the short reads
   bamlist=list()
   src_files=list.files(pattern="*_sorted.bam$")
   for(filename in src_files){
     #Since we do not know which strand the reads were originally transcribed,
     #so set the strand to be ambiguous
     tmp=readBamGappedAlignments(filename)
     bamlist[[length(bamlist)+1]]=GRanges(seqnames=rname(tmp),
       ranges=IRanges(start=start(tmp),end=end(tmp)),
       strand=rep("*",length(tmp)))
   }
   names(bamlist)=src_files

将文件加载到R中后，我们需要创建一个注释对象。由于我们使用的是hg19，我们可以使用GenomicFeatures包从UCSC轻松下载一个。我们选择使用ENSEMBL基因注释。

   library(GenomicFeatures)
   txdb=makeTranscriptDbFromUCSC(genome="hg19",tablename="ensGene")

我们想比较基因的差异表达，因此我们将按基因进行汇总，并选择将所有落在基因体内（包括内含子）的reads计入基因的计数。

   tx_by_gene=transcriptsBy(txdb,"gene")

最后，我们计算每个lane中落在每个基因的reads数量，并将结果记录在一个计数表中。

   #Initialize table of counts
   toc=data.frame(rep(NA,length(tx_by_gene)))
   for(i in 1:length(bamlist)){
     toc[,i]=countOverlaps(tx_by_gene,bamlist[[i]])
   }
   #Fix up labels
   rownames(toc)=names(tx_by_gene)
   colnames(toc)=names(bamlist)

差异表达测试

最终获得一个计数表后，我们现在想比较处理组和未处理组，并寻找每个基因的计数数量是否存在任何统计学上的显著差异。我们将使用edgeR使用的负二项式模型来做到这一点。

归一化

我们使用TMM方法（以第一个lane为参考）计算合适的缩放因子进行归一化。

   library(edgeR)
   norm_factors=calcNormFactors(as.matrix(toc))

计数本身并没有改变，而是将这些缩放因子用作负二项式模型中的偏移量。这被包含在edgeR需要的DGE列表对象中。

   DGE=DGEList(toc,lib.size=norm_factors*colSums(toc),group=gsub("[0-9].*","",colnames(toc)))

统计测试

接下来我们计算一个共同的离散参数，它表示数据中额外的泊松可变性。

   disp=estimateCommonDisp(DGE)

这使我们能够计算基因差异表达的p值。

   tested=exactTest(disp)

基因本体测试

为了测试DE基因中过度表达的GO类别，我们首先需要选择一个阈值，用于在应用多重假设校正后将基因称为差异表达。我们选择流行的显著性阈值为0.05。

   library(goseq)
   #Apply benjamini hochberg correction for multiple testing
   #choose genes with a p-value less than .05
   genes=as.integer(p.adjust(tested$table$p.value,method="BH") <.05)
   names(genes)=row.names(tested$table)

现在我们计算概率加权函数，它量化了长度偏差效应。

   pwf=nullp(genes,"hg19","ensGene")

最后，我们计算每个GO类别在DE基因中过度表达的p值。

   GO.pvals=goseq(pwf,"hg19","ensGene")

示例2 差异表达：拟南芥病原体数据

 Donwnload the six fastq.gz files from here: http://www.ebi.ac.uk/ena/data/view/SRP004047&display=html

这些文件来自一项拟南芥研究，使用受感染植物和模拟感染植物的三个重复样本。这是NBPSeq R包的基础数据集。

质量控制

 Download FastQC and open the files in FastQC one by one.
 You can open the files by using File->Open. 

截至版本0.94，如果您使用的是Windows，请使用Linux版本并双击run_fastqc.bat。

 You will see a very wiggly line for the first library. If you just look at the peaks and note the sequence you will see the pattern
 AAGAGCTCGTATGC starting at the green plateau towards the right. This is an illumina adapter sequence, which you will also see in the
 overrepresented counts tab.

文件：FastQC.png

读取比对

从这里下载拟南芥基因组版本。

 ftp://ftp.arabidopsis.org/home/tair/Sequences/blast_datasets/other_datasets/CURRENT/arabidopsis.seq.gz
 It doesn't include mitochondrial or chloroplast sequences but is good enough for the tutorial purpose

 unzip all read files 

使用bowtie需要构建索引。

 Unix systems
 bowtie-build arabidopsis.seq arabidopsis

 Windows
 C:\Prog\bowtie-0.12.7\bowtie-build.exe arabidopsis.seq arabidopsis

或者，如果您使用的是支持的生物体之一，您可以从bowtie网站下载提供的索引文件。

现在您可以使用bowtie比对reads。 Bowtie 具有许多选项，您最好检查它们。这里我们告诉它使用两个处理器（-p 2）来报告基于SAM的比对结果 -S，所有（-a）比对最多只有一个错配（-v 1）。我们进一步限制了这一点，因为只有每个reads提供多个有效比对结果才应该被丢弃（-m 1）。

 bowtie -p 2 -S -a -m 1 -v 1 arabidopsis SRR074262.fastq aligned_074262
 bowtie -p 2 -S -a -m 1 -v 1 arabidopsis SRR074263.fastq aligned_074263
 bowtie -p 2 -S -a -m 1 -v 1 arabidopsis SRR074264.fastq aligned_074264
 bowtie -p 2 -S -a -m 1 -v 1 arabidopsis SRR074284.fastq aligned_074284
 bowtie -p 2 -S -a -m 1 -v 1 arabidopsis SRR074285.fastq aligned_074285
 bowtie -p 2 -S -a -m 1 -v 1 arabidopsis SRR074286.fastq aligned_074286

 C:\Prog\bowtie-0.12.7\bowtie.exe -p 2 -S -a -m 1 -v 1 arabidopsis SRR074262.fastq aligned_074262
 

这需要对所有6个库重复。

输出SRR074262.fastq

1. 已处理的reads：9619406
2. 至少有一个报告的比对结果的reads：5034111（52.33%）
3. 无法比对的reads：4132070（42.96%）
4. 由于-m被抑制的比对结果的reads：453225（4.71%）

报告了5034111个比对结果到1个输出流中。

输出SRR074263.fastq

1. 已处理的reads：4199495
2. 至少有一个报告的比对结果的reads：2298490（54.73%）
3. 无法比对的reads：1701822（40.52%）
4. 由于-m被抑制的比对结果的reads：199183（4.74%）

报告了2298490个比对结果到1个输出流中。

输出SRR074264.fastq

1. 已处理的reads：4810892
2. 至少有一个报告的比对结果的reads：2730242（56.75%）
3. 无法比对的reads：1852040（38.50%）
4. 由于-m被抑制的比对结果的reads：228610（4.75%）

报告了2730242个比对结果到1个输出流中。

输出SRR074284.fastq

1. 已处理的reads：4763394
2. 至少有一个报告的比对结果的reads：2622118（55.05%）
3. 无法比对的reads：1908441（40.06%）
4. 由于-m被抑制的比对结果的reads：232835（4.89%）

报告了2622118个比对结果到1个输出流中。

输出SRR074285.fastq

1. 已处理的reads：8217239
2. 至少有一个报告的比对结果的reads：4515790（54.96%）
3. 无法比对的reads：3206578（39.02%）
4. 由于-m被抑制的比对结果的reads：494871（6.02%）

报告了4515790个比对结果到1个输出流中。

输出SRR074286.fastq

1. 已处理的reads：3776979
2. 至少有一个报告的比对结果的reads：2014498（53.34%）
3. 无法比对的reads：1506142（39.88%）
4. 由于-m被抑制的比对结果的reads：256339（6.79%）

报告了2014498个比对结果到1个输出流中。

汇总Reads

如果您想使用HTSeq-count 并且使用的是Fedora或CentOS，您将需要付出一些额外的努力，因为HTSeq使用的是python 2.6，而Fedora和CentOS只安装了python 2.4。（在执行以下任何操作之前，请咨询您的系统管理员是否允许这样做）

这里有一些建议 http://stackoverflow.com/questions/1465036/install-python-2-6-in-centos 我最喜欢的是使用EPEL（见下文），但YMMV。

 rpm -Uvh http://download.fedora.redhat.com/pub/epel/5/i386/epel-release-5-4.noarch.rpm
 yum install python26
 yum install python26-numpy
 yum install python26-numpy-devel

然后按照网站上的详细说明继续，但键入

 python26 setup.py build
 python26 setup.py install

而不是

 python setup.py build
 python setup.py install

---

现在我们准备好下载gff文件了： ftp://ftp.arabidopsis.org/home/tair/Genes/TAIR10_genome_release/TAIR10_gff3/TAIR10_GFF3_genes.gff

该文件不包含RNA基因和转座因子。

现在，父ID可能适用于HTSeq的一些用例。如果您需要一个真正的GTF文件，您可以使用提到的perl脚本 这里 底部。如果我们将其命名为celeste_script.pl，这里将是执行的内容。（我们不能使用 -i Parent 的原因是，否则外显子可能属于两个剪接变体，例如 AT1G01040.1 和 AT1G01040.2，而不会被计入，因为我们想使用交集严格模式）

 cat TAIR10_GFF3_genes.gff |perl celeste_script.pl > TAIR10_GTF.gtf 

现在终于

 htseq-count -m intersection-strict -s no aligned_074262 TAIR10_GTF.gff > counts_074262
 htseq-count -m intersection-strict -s no aligned_074263 TAIR10_GTF.gtf > counts_074263
 htseq-count -m intersection-strict -s no aligned_074264 TAIR10_GTF.gtf > counts_074264
 htseq-count -m intersection-strict -s no aligned_074284 TAIR10_GTF.gtf > counts_074284
 htseq-count -m intersection-strict -s no aligned_074285 TAIR10_GTF.gtf > counts_074285
 htseq-count -m intersection-strict -s no aligned_074286 TAIR10_GTF.gtf > counts_074286

现在我们开始得到一些我们可以使用的東西。

 head counts_074262

 AT1G01010       75
 AT1G01020       73
 AT1G01030       33
 AT1G01040       109
 AT1G01046       0
 AT1G01050       41
 AT1G01060       10
 AT1G01070       2
 AT1G01073       0
 AT1G01080       281

 head counts_074263

 AT1G01010       35
 AT1G01020       34
 AT1G01030       18
 AT1G01040       52
 AT1G01046       0
 AT1G01050       51
 AT1G01060       5
 AT1G01070       19
 AT1G01073       0
 AT1G01080       212

 head counts_074264

 AT1G01010       80
 AT1G01020       46
 AT1G01030       45
 AT1G01040       50
 AT1G01046       0
 AT1G01050       69
 AT1G01060       13
 AT1G01070       35
 AT1G01073       0
 AT1G01080       226

 head counts_074284

 AT1G01010       45
 AT1G01020       50
 AT1G01030       33
 AT1G01040       65
 AT1G01046       0
 AT1G01050       50
 AT1G01060       0
 AT1G01070       7
 AT1G01073       0
 AT1G01080       178

 head counts_074285

 AT1G01010       35
 AT1G01020       43
 AT1G01030       26
 AT1G01040       72
 AT1G01046       0
 AT1G01050       51
 AT1G01060       3
 AT1G01070       6
 AT1G01073       0
 AT1G01080       531

 head counts_074286

 AT1G01010       61
 AT1G01020       27
 AT1G01030       36
 AT1G01040       25
 AT1G01046       0
 AT1G01050       25
 AT1G01060       16
 AT1G01070       20
 AT1G01073       0
 AT1G01080       128

差异表达基因的分析

启动R

使用DESeq（现在参考 vignette）

 R
 

在R中

 library(DESeq)
 setwd("where the data is")
 c1<-read.table("counts_074284",row.names=1) 
 c2<-read.table("counts_074286",row.names=1)
 c3<-read.table("counts_074262",row.names=1)
 c4<-read.table("counts_074263",row.names=1)
 c5<-read.table("counts_074264",row.names=1)
 c6<-read.table("counts_074285",row.names=1)
 counts<-cbind(c1,c2,c3,c4,c5,c6)
 counts<-counts[-c(32679:32683),]  #remove the more general lines
 colnames(counts)<-c("P1","P2","P3","M1","M2","M3")
 design <- rep (c("P","Mo"),each=3)
 de  <-  newCountDataSet(counts, design)
 de  <-  estimateSizeFactors( de)
 # de  <-  estimateVarianceFunctions(  de  ) # This function has been removed. Use 'estimateDispersions' instead.  
 de <- estimateDispersions( de )
 res  <-  nbinomTest(  de,  "P",  "Mo")

有多少基因是显著的？

 sum(na.omit(res$padj<0.05))