摘要
本文，我们将使用Bioconductor包进行端对端基因水平的RNA-seq差异表达分析。我们将从FASTQ文件开始，显示这些数据如何被比对到参考基因组以及如何准备一个计数矩阵（该计数矩阵计算：每个样品中的每个基因内RNA-seq reads/fragments的数量）。我们会进行探索性数据分析（目的是评估测序数据的质量），接着探索不同样本之间的关系，最后进行差异表达基因分析及结果可视化。
前言
Bioconductor中有很多包都支持高通量数据（比如RNA-seq）的分析，我们在此workflow中使用的包主要包括两部分：用于导入和处理原始测序数据并加载基因注释的核心包（由Bioconductor核心团队维护）、对测序数据进行统计分析和可视化的包。Bioconductor项目每半年发布一个计划版本，这便确保每个发布版本中的包都能互相协调。用于本workflow中的包可以通过library功能进行加载，用户可以通过Bioconductor的package安装说明进行安装。
1.1 实验数据
本工作流程中使用的数据存储在airway package中，该软件包总结了一项RNA-seq实验，该实验使用地塞米松对气道平滑肌细胞进行处理，而地塞米松是一种具有抗炎作用的合成糖皮质激素类固醇。
2. 准备计数矩阵
基于count的统计数据（如DESeq2, edgeR，limma（with voom），DSS，EBSeq，baySeq），以及来自RNA-seq或其他高通量测序实验的数据（该数据的值需采用矩阵形式进行整合），都可以作为输入文件。矩阵的第i行和第j列中的值表示在样品j中，有多少reads（或paired-end的RNA-seq的片段）被分配给了基因i。类似地，对于其他类型的分析，矩阵的行可能对应于binding区域（结合了ChIP-seq），细菌种类（使用meta全基因组数据集）或者是多肽序列（结合了定量质谱。）
矩阵中的值应该是测序读数/片段的计数。这对DESeq2的统计模型来说十分重要，因为只有count才能正确评估测量精确度。需要注意的是，输入的count文件绝对不能进行预先标准化（比如基于测序深度/文库大小），因为统计模型在应用于非标准化的count时才更有用，统计模型主要是为了在内部考虑库大小的差异。
2.1 推荐：转录本丰度和tximport管道
在我们演示如何比对并进一步对RNA-seq片段进行计数之前，我们想介绍一个比以往更新更快的可替代性pipeline：那就是使用转录本丰度定量方法（如Salmon，Sailfish，kallisto或RSEM），在不比对reads的情况下估计转录本丰度，然后使用Bioconductor差异基因表达包中的tximport软件包来组装预估的count和偏移矩阵，以供使用。
关于如何使用Salmon软件来量化转录本丰度的教程可以在这里找到。我们推荐使用–gcBias 这个flag，它是RNA-seq数据中最常用的估计系统性偏差的校正因子，除非您能够确定您的数据中不包含这些偏差。在Salmon教程中，你可以使用tximport vignette中的步骤，它将告诉你如何构建一个DESeqDataSet。这是目前我们向您推荐的pipeline，但下文还会包含一些步骤，例如从将reads比对到基因组以及从BAM文件构建count矩阵。
将transcript abundance quantifiers软件与tximport结合起来使用，以产生基因水平计数矩阵和归一化偏移的优点是：这种方法可以纠正样本间基因长度的任何潜在变化（比如使用差异异构体）;与基于比对（alignment）的方法相比，其中一些方法的速度要快得多，需要的内存和磁盘使用量也更少;我们也可以避免丢弃那些可以与同源序列的多个基因进行对齐的片段。请注意，转录本丰度量词会直接跳过生成存储reads比对结果的大文件（而不是生成存储每个转录本的预估丰度、count数及有效长度的小文件。有关更多详细信息，请参阅描述此方法的手稿，以及用于软件详细信息的tximport package vignette。使用转录本量词和tximport后返回的entry point，将是下面分析的重点部分。
2.2 将reads比对到参考基因组
RNA-seq实验的计算分析开始于FASTQ文件（该文件包含每个核苷酸序列的read及其对应位置的质量分数）。这些reads首先要比对到参考基因组或转录组上，或者如上所述，可以在没有比对的情况下估计每个转录本的丰度和预估count数。无论是哪一种情况，都需要预先了解测序数据的情况（是单端测序还是双端测序），如果是双端测序，那么比对软件需要用户提供两个FASTQ文件。该比对步骤通常以SAM/BAM格式输出。
用于将reads比对到参考基因组的软件有很多，我们建议您查看一些讨论每种软件优缺点的文章，这些文章对不同软件的准确性、将reads比对到splice junctions的灵敏度、速度、需要的内存、可用性和许多其他功能进行了对比。在此，我们使用STAR read aligner，将RNA-seq实验产生的reads比对到Ensembl数据库中的人类基因组（第75版）。在下面的代码示例中，假设当前目录中有一个名为files的文件，每行包含每个实验的标识符，并且在fastq子目录中包含所有的FASTQ文件。如果您已从序列读取存档下载FASTQ文件，则标识符将是SRA运行ID，例如，SRR1039520。 对于每个ID，您应该有两个文件用于双端实验，fastq / SRR1039520_1.fastq1和fastq / SRR1039520_2.fastq，它们为双端的fragments提供第一次和第二次的read。如果您运行的是单端实验数据，则每个ID只能有一个文件。 我们还创建了一个叫aligned的子目录，用于STAR软件输出其alignment结果。
for f in `cat files`; do STAR --genomeDir ../STAR/
ENSEMBL.homo_sapiens.release-75 \
--readFilesIn fastq/$f\_1.fastq fastq/$f\_2.fastq \
--runThreadN 12 --outFileNamePrefix aligned/$f.; done
SAMtools用于生成BAM文件。 -@标志可用于分配其他线程。
for f in `cat files`; do samtools view -bS aligned/$f.Aligned.out.sam \
-o aligned/$f.bam; done
然后，可以使用多个测序runs中的BAM文件生成计数矩阵，如以下部分所述。
2.3 对比对文件进行定位
除了我们稍后将使用的计数矩阵之外，airway包还包含八个文件，其中包含实验中reads总数的一小部分。 我们选择了比对到一号染色体上一小段区域的reads。而为了演示，我们选择了该reads的一小部分子集，因为完整的比对文件很大（每个几千兆字节），而且，在对每个sample中的进行计数时，每个sample平均需要花费10-30分钟。 我们将使用这些文件来演示如何利用BAM文件来构造计数矩阵。 然后，我们将加载与所有sample相关的计数矩阵，并将其用于进一步的分析。
我们使用示例数据加载数据包：
library("airway")
R函数system.file可用于查找计算机上已安装package中文件的位置。 在这里，我们要查看extdata目录的完整路径（它是airway包的一部分）。
indir <- system.file("extdata", package="airway", mustWork=TRUE)
在这个目录中，我们找到了八个BAM文件（和一些其他文件）：
list.files(indir)
## [1] "GSE52778_series_matrix.txt" "Homo_sapiens.GRCh37.75_subset.gtf"
## [3] "SRR1039508_subset.bam" "SRR1039509_subset.bam"
## [5] "SRR1039512_subset.bam" "SRR1039513_subset.bam"
## [7] "SRR1039516_subset.bam" "SRR1039517_subset.bam"
## [9] "SRR1039520_subset.bam" "SRR1039521_subset.bam"
## [11] "SraRunInfo_SRP033351.csv" "sample_table.csv"
通常，我们会有一个包含每个sample详细信息的表格，可以将sample关联到其对应的FASTQ和BAM文件。 对于您自己的项目，您可以使用文本编辑器或电子表格软件（如Excel）创建这样的逗号分隔值（CSV）文件。
我们用read.csv加载这样一个CSV文件：
csvfile <- file.path(indir, "sample_table.csv")
sampleTable <- read.csv(csvfile, row.names = 1)
sampleTable
## SampleName cell dex albut Run avgLength Experiment
## SRR1039508 GSM1275862 N61311 untrt untrt SRR1039508 126 SRX384345
## SRR1039509 GSM1275863 N61311 trt untrt SRR1039509 126 SRX384346 ## SRR1039512 GSM1275866 N052611 untrt untrt SRR1039512 126 SRX384349
## SRR1039513 GSM1275867 N052611 trt untrt SRR1039513 87 SRX384350
## SRR1039516 GSM1275870 N080611 untrt untrt SRR1039516 120 SRX384353
## SRR1039517 GSM1275871 N080611 trt untrt SRR1039517 126 SRX384354
## SRR1039520 GSM1275874 N061011 untrt untrt SRR1039520 101 SRX384357
## SRR1039521 GSM1275875 N061011 trt untrt SRR1039521 98 SRX384358
## Sample BioSample
## SRR1039508 SRS508568 SAMN02422669
## SRR1039509 SRS508567 SAMN02422675
## SRR1039512 SRS508571 SAMN02422678
## SRR1039513 SRS508572 SAMN02422670
## SRR1039516 SRS508575 SAMN02422682
## SRR1039517 SRS508576 SAMN02422673
## SRR1039520 SRS508579 SAMN02422683
## SRR1039521 SRS508580 SAMN02422677
一旦reads被比对到基因组上，我们就可以使用许多工具来计算reads/fragments的数量，这些软件需要两个输入文件：一个是SAM/BAM格式的比对文件，另一个是基因组注释文件（GFF3或GTF文件）。
2.4 DESeq2的导入功能
以下工具可用于生成计数矩阵：summarizeOverlaps，featureCounts，tximport，htseq-count。
我们现在使用summarizeOverlaps进行下面的分析。 使用示例表中的Run列，我们构造了要执行计数操作的文件的完整路径：
filenames <- file.path(indir, paste0(sampleTable$Run, "_subset.bam"))
file.exists(filenames)
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
我们在Bioconductor中指出,这些BAM文件是使用Rsamtools包中的BamFileList函数导入的，该函数提供了BAM文件的R接口。 在这里，我们还详细说明了应如何处理BAM文件（例如最大只能同时处理200万个reads）。
library("Rsamtools")
bamfiles <- BamFileList(filenames, yieldSize=2000000)
注意：确保您使用的基因组文件中的染色体名，与用于reads比对的参考注释文件中的染色体名称一致。 否则，由于名称不匹配，脚本可能无法计算任何reads的特征（features）。 请参阅GenomeInfoDb包中的seqlevelsStyle函数以获取帮助。 我们可以检查alignment文件中的染色体名称（这里称为“seqnames”），如下所示：
seqinfo(bamfiles[1])
## Seqinfo object with 84 sequences from an unspecified genome:
## seqnames seqlengths isCircular genome
## 1 249250621 <NA> <NA>
## 10 135534747 <NA> <NA>
## 11 135006516 <NA> <NA>
## 12 133851895 <NA> <NA>
## 13 115169878 <NA> <NA>
## ... ... ... ...
## GL000210.1 27682 <NA> <NA>
## GL000231.1 27386 <NA> <NA>
## GL000229.1 19913 <NA> <NA>
## GL000226.1 15008 <NA> <NA>
## GL000207.1 4262 <NA> <NA>
2.5 定义基因模型（gene model）
接下来，我们要读入gene model，以便对reads/fragments进行计数。我们将使用GenomicFeatures包中的makeTxDbFromGFF从Ensembl的GTF文件中读取gene model。 GTF文件可以从Ensembl的FTP站点或其他数据库下载。 TxDb对象是一个数据库，可用于生成各种基于范围的对象，如外显子，转录本和基因。 我们想要制作一个按基因进行分组的外显子列表，用于计算reads/fragments。
在构建TxDb时，还有一些其他的选项。 对于来自UCSC基因组浏览器的已知基因，我们可以使用其预先构建的转录本数据库：TxDb.Hsapiens.UCSC.hg19.knownGene。 如果AnnotationHub中有相应的注释文件（与Ensembl基因的情况一样），则可以使用makeTxDbFromGRanges导入该GTF文件。
这里我们将演示GTF文件的加载：
library("GenomicFeatures")
我们使用0长度字符向量表明我们的序列（染色体）都不是环形的。
gtffile <- file.path(indir,"Homo_sapiens.GRCh37.75_subset.gtf")
txdb <- makeTxDbFromGFF(gtffile, format = "gtf", circ_seqs = character())
txdb
## TxDb object:
## # Db type: TxDb
## # Supporting package: GenomicFeatures
## # Data source: /home/biocbuild/bbs-3.8-bioc/R/library/airway/extdata/Homo_sapiens.GRCh37.75_subset.gtf
## # Organism: NA
## # Taxonomy ID: NA
## # miRBase build ID: NA
## # Genome: NA
## # transcript_nrow: 65
## # exon_nrow: 279
## # cds_nrow: 158
## # Db created by: GenomicFeatures package from Bioconductor
## # Creation time: 2018-10-31 08:06:55 -0400 (Wed, 31 Oct 2018)
## # GenomicFeatures version at creation time: 1.34.0
## # RSQLite version at creation time: 2.1.1
## # DBSCHEMAVERSION: 1.2
下面这一行产生了由基因进行分组的所有外显子的GRangesList。列表的每个元素都是基因外显子的GRanges对象。
ebg <- exonsBy(txdb, by="gene")
ebg
## GRangesList object of length 20:
## $ENSG00000009724
## GRanges object with 18 ranges and 2 metadata columns:
## seqnames ranges strand | exon_id exon_name
## <Rle> <IRanges> <Rle> | <integer> <character>
## [1] 1 11086580-11087705 - | 98 ENSE00000818830
## [2] 1 11090233-11090307 - | 99 ENSE00000472123
## [3] 1 11090805-11090939 - | 100 ENSE00000743084
## [4] 1 11094885-11094963 - | 101 ENSE00000743085
## [5] 1 11097750-11097868 - | 102 ENSE00003482788
## ... ... ... ... . ... ...
## [14] 1 11106948-11107176 - | 111 ENSE00003467404
## [15] 1 11106948-11107176 - | 112 ENSE00003489217
## [16] 1 11107260-11107280 - | 113 ENSE00001833377
## [17] 1 11107260-11107284 - | 114 ENSE00001472289
## [18] 1 11107260-11107290 - | 115 ENSE00001881401
##
## ...
## <19 more elements>
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengths
2.6 对read进行计数
经过了上述的准备工作，下面的count计数部分实际上是非常容易的，来自GenomicAlignments包中的summarizeOverlaps函数可以完成这个计数工作。这将生成一个SummarizedExperiment对象，其中包含有关实验的各种信息，我们将在下面的步骤中详细描述。
注意：如果您需要使用多个内核执行计数，可以在调用count之前，使用BiocParallel包中的register和MulticoreParam函数或SnowParam函数。 对于具有3000万个比对reads的人类RNA-seq文件来说，summarizeOverlaps函数调用每个RNA-seq文件需要花费至少需要30分钟。 不过，若将文件发送到单独的核心，则可以加快整个计数过程。
library("GenomicAlignments")
library("BiocParallel")
在这里，我们指定使用一个核心，而不是多个核心。 我们也可以跳过这一行，计数步骤将连续运行。
register(SerialParam())
以下调用会创建含有count的SummarizedExperiment对象：
se <- summarizeOverlaps(features=ebg, reads=bamfiles,
mode="Union",
singleEnd=FALSE,
ignore.strand=TRUE,
fragments=TRUE )
除了features和reads以外，我们还指定了许多参数。mode参数描述了哪些重叠的reads将会被计数。图1主要展示了在使用GenomicAlignments包中的summarizeOverlaps函数对reads进行计数的过程中所用到的一些模式。注意，每个fragment只能在同一基因上计数一次，即使它们与基因的多个外显子重叠。将singleEnd参数设置为FALSE表示该实验产生的是单端测序reads，并且将一对reads（一个fragment）在基因上只计数一次。当singleEnd = FALSE时，可以使用fragment参数来指定是否应计算未配对的读数（如果fragment = TRUE，则为yes）。
为了产生正确的计数，重要的是我们要知道该RNA-seq实验是否为链特异性的。如果这个实验不是链特异性的，那我们就把ignore.strand设置为TRUE。但是，有些链特异性的protocol可能会将reads仅仅比对到基因的互补链上。 用户必须检查该实验是否为链特异性的，如果是，那么reads是否应该比对到基因的正向或互补链上。 对于各种对reads进行计数或定量的软件，我们都可以指定以不同方式对正向链或互补链进行计数（目前最容易使用的软件为htseq-count，featureCounts或前面提到的transcript abundance quantifiers）。 在计数完成后，对count矩阵的列数之和进行检查，可以知道该结果是否与我们预期的结果相一致（即比对到基因上的reads数或fragment数是否和预期相匹配）。另外，我们还可以使用基因组可视化工对比对结果进行检查。
2.7 SummarizedExperiment
上图是SummarizedExperiment container的组成部分；这些组分就是SummarizedExperiment的对象。
- assay（粉格）：包含了名为“counts”的矩阵；
- rowRanges（蓝格）：包含了基因组范围内的信息；
- colData（绿格）：包含了样品的信息。
- 每个块中突出显示的行代表第一行（注意colData的第一行与assay的第一列对齐）。
在上面的例子中，我们创建了一个名为“counts”的单一矩阵，其中包含每个基因和sample的片段计数，该矩阵存储于assay对象中。assay对象也可以存储多个矩阵。rowRanges对象存储了用于计数的GRangesList（每一行对应一个外显子的GRanges）。我们使用相同的R函数名称访问SummarizedExperiment的三个组分：assay(assays)，rowRanges和colData。
上面的示例代码实际上只计算了原始实验中的一小部分片段。 尽管如此，我们仍然可以通过查看assay中的计数数据、colData中的样本表型数据（在本例中为空DataFrame）以及rowRanges中的基因组信息来研究得到的SummarizedExperiment。
se
## class: RangedSummarizedExperiment
## dim: 20 8
## metadata(0):
## assays(1): counts
## rownames(20): ENSG00000009724 ENSG00000116649 ... ENSG00000271794
## ENSG00000271895
## rowData names(0):
## colnames(8): SRR1039508_subset.bam SRR1039509_subset.bam ...
## SRR1039520_subset.bam SRR1039521_subset.bam
## colData names(0):
dim(se)
## [1] 20 8
assayNames(se)
## [1] "counts"
head(assay(se), 3)
## SRR1039508_subset.bam SRR1039509_subset.bam
## ENSG00000009724 38 28
## ENSG00000116649 1004 1255
## ENSG00000120942 218 256
## SRR1039512_subset.bam SRR1039513_subset.bam
## ENSG00000009724 66 24
## ENSG00000116649 1122 1313
## ENSG00000120942 233 252
## SRR1039516_subset.bam SRR1039517_subset.bam
## ENSG00000009724 42 41
## ENSG00000116649 1100 1879
## ENSG00000120942 269 465
## SRR1039520_subset.bam SRR1039521_subset.bam
## ENSG00000009724 47 36
## ENSG00000116649 745 1536
## ENSG00000120942 207 400
colSums(assay(se))
## SRR1039508_subset.bam SRR1039509_subset.bam SRR1039512_subset.bam
## 6478 6501 7699
## SRR1039513_subset.bam SRR1039516_subset.bam SRR1039517_subset.bam
## 6801 8009 10849
## SRR1039520_subset.bam SRR1039521_subset.bam
## 5254 9168
rowRanges在打印时只显示第一个GRanges，并告诉我们还有19个元素：
rowRanges(se)
## GRangesList object of length 20:
## $ENSG00000009724
## GRanges object with 18 ranges and 2 metadata columns:
## seqnames ranges strand | exon_id exon_name
## <Rle> <IRanges> <Rle> | <integer> <character>
## [1] 1 11086580-11087705 - | 98 ENSE00000818830
## [2] 1 11090233-11090307 - | 99 ENSE00000472123
## [3] 1 11090805-11090939 - | 100 ENSE00000743084
## [4] 1 11094885-11094963 - | 101 ENSE00000743085
## [5] 1 11097750-11097868 - | 102 ENSE00003482788
## ... ... ... ... . ... ...
## [14] 1 11106948-11107176 - | 111 ENSE00003467404
# [15] 1 11106948-11107176 - | 112 ENSE00003489217
## [16] 1 11107260-11107280 - | 113 ENSE00001833377
## [17] 1 11107260-11107284 - | 114 ENSE00001472289
## [18] 1 11107260-11107290 - | 115 ENSE00001881401
##
## ...
## <19 more elements>
## -------
## seqinfo: 1 sequence from an unspecified genome; no seqlengths
rowRanges还包含了构建gene model时的meta数据。 这里我们使用str函数来显示该meta数据：
str(metadata(rowRanges(se)))
## List of 1
## $ genomeInfo:List of 15
## ..$ Db type : chr "TxDb"
## ..$ Supporting package : chr "GenomicFeatures"
## ..$ Data source : chr "/home/biocbuild/bbs-3.8-bioc/R/library/airway/extdata/Homo_sapiens.GRCh37.75_subset.gtf"
## ..$ Organism : chr NA
## ..$ Taxonomy ID : chr NA
## ..$ miRBase build ID : chr NA
## ..$ Genome : chr NA
## ..$ transcript_nrow : chr "65"
## ..$ exon_nrow : chr "279"
## ..$ cds_nrow : chr "158"
## ..$ Db created by : chr "GenomicFeatures package from Bioconductor"
## ..$ Creation time : chr "2018-10-31 08:06:55 -0400 (Wed, 31 Oct 2018)"
## ..$ GenomicFeatures version at creation time: chr "1.34.0"
## ..$ RSQLite version at creation time : chr "2.1.1"
## ..$ DBSCHEMAVERSION : chr "1.2"
colData：
colData(se)
## DataFrame with 8 rows and 0 columns
虽然colData目前还是空的，但它应该包含所有的meta数据。因为我们使用了一列sampleTable来生成bamfiles向量，所以se的列和sampleTable的行的顺序相同。我们可以将sampleTable变为summarized experiment的colData对象， 方法是使用赋值函数将其转换成数据框（DataFrame）：
colData(se) <- DataFrame(sampleTable)
colData(se)
## DataFrame with 8 rows and 9 columns
## SampleName cell dex albut Run avgLength
## <factor> <factor> <factor> <factor> <factor> <integer>
## SRR1039508 GSM1275862 N61311 untrt untrt SRR1039508 126
## SRR1039509 GSM1275863 N61311 trt untrt SRR1039509 126
## SRR1039512 GSM1275866 N052611 untrt untrt SRR1039512 126
## SRR1039513 GSM1275867 N052611 trt untrt SRR1039513 87
## SRR1039516 GSM1275870 N080611 untrt untrt SRR1039516 120
## SRR1039517 GSM1275871 N080611 trt untrt SRR1039517 126
## SRR1039520 GSM1275874 N061011 untrt untrt SRR1039520 101
## SRR1039521 GSM1275875 N061011 trt untrt SRR1039521 98
## Experiment Sample BioSample
## <factor> <factor> <factor>
## SRR1039508 SRX384345 SRS508568 SAMN02422669
## SRR1039509 SRX384346 SRS508567 SAMN02422675
## SRR1039512 SRX384349 SRS508571 SAMN02422678
## SRR1039513 SRX384350 SRS508572 SAMN02422670
## SRR1039516 SRX384353 SRS508575 SAMN02422682
## SRR1039517 SRX384354 SRS508576 SAMN02422673
## SRR1039520 SRX384357 SRS508579 SAMN02422683
## SRR1039521 SRX384358 SRS508580 SAMN02422677
2.8 分支点
到此为止，我们已经对那些能够比对到指定gene model上的fragment进行了计数。这是一个分支点，接下来，我们便可以使用各种Bioconductor软件包对count结果进行探索和差异表达分析（比如edgeR 、limma、DSS、EBSeq和baySeq）。Schurch等人比较了市面上各种RNA-seq分析软件的利弊，那篇文章可以帮助用户按不同的需求来选择分析软件。在下面的研究中，我们将继续使用DESeq2。 SummarizedExperiment的三个对象是我们开始分析所需的全部内容。 在下一节中，我们将展示如何使用它们来创建DESeq2的数据对象。
3. DESeqDataSet的对象、样本信息和设计公式
Bioconductor软件包一般使用自定义的类（class）来存储数据。此外，Bioconductor有一些通用的类（比如SummarizedExperiment）可以用于在不同的package之间移动数据。核心Bioconductor的类（class）还提供了一些有用的功能：比如，SummarizedExperiment对行或列取子集或重新排序时，会根据相关的rowRanges或colData进行自动取子集或排序，这有助于防止出错。当您使用SummarizedExperiment时，这一切都在后台处理。
在DESeq2软件中，这个自定义类被称为DESeqDataSet。它构建在SummarizedExperiment类之上，我们可以很容易的将SummarizedExperiment对象转换为DESeqDataSet对象。它们之间有两个主要的区别，一是当使用counts accessor函数时，实际上访问的是assay slot；而在DESeqDataSet类中，矩阵中的值必须是非负整数。
第二个区别是DESeqDataSet具有一个关联性的design formula。在分析之前，我们就需要指定实验设计的内容，因为许多DESeq函数将基于实验设计来处理和分析样本（一个例外是size factor estimation，即不同库大小的调整，该因子不取决于design formula）。design formula会告诉样本信息表（colData）哪些列指定了实验设计以及如何在分析中使用这些因子。
用于差异表达的最简单的design formula是~ condition，其中condition是colData（dds）中的一列，用于指定样本属于哪两个组（比如对照组和处理组）或多个组。对于airway实验，我们将指定~ cell + dex，这表示我们想测试地塞米松（dex）在不同细胞中的效果。我们可以在SummarizedExperiment或DESeqDataSet上直接使用$来查看每个列：
se$cell
## [1] N61311 N61311 N052611 N052611 N080611 N080611 N061011 N061011
## Levels: N052611 N061011 N080611 N61311
se$dex
## [1] untrt trt untrt trt untrt trt untrt trt
## Levels: trt untrt
注意：在R中，factor的第一级是对照组（reference level）（比如，对照或未经处理的样本），因此我们可以像这样重新调整dex因子：
library("magrittr")
se$dex %<>% relevel("untrt")
se$dex
## [1] untrt trt untrt trt untrt trt untrt trt
## Levels: untrt trt
%<>%是magrittr包中用于赋值的管道运算符，上面的代码是一种简洁的写法。
se$dex <- relevel(se$dex, "untrt")
在运行DESeq2时，您可以使用R中的formula notation来实现各种实验设计。请注意，DESeq2将使用和一些基础R包（比如lm函数）相同的formula notation。如果研究目的是为了确定在不同的处理组之间的差异基因，那我们可以使用如下的交互命令：~ group + treatment + group:treatment。有关更多示例，请参见手册页以获取结果。在下面的分析中，我们将采取交互命令来测试不同condition下的变化。
接下来，我们将从两个方面开始构造DESeqDataSet：
- 来自一个SummarizedExperiment对象
- 来自一个计数矩阵和一个样本信息表
有关使用HTSeq Python包进行读取计数的完整示例，请参阅pasilla vignette。而使用htseq-count生成的计数文件来进一步生成DESeqDataSet的示例，请参阅DESeq2 vignette。
3.1 从SummarizedExperiment开始
我们现在使用一段R代码来加载一个已经准备好的SummarizedExperiment（是由Himes等人公开的测序数据文件生成的）。我们用于生成此对象的步骤与您在之前所做的步骤相同，只是这次我们用了测序数据的全部reads和所有基因（上面的步骤只使用一小部分数据来示范）。有关用于创建此对象的确切步骤的更多详细信息，请在R窗口中键入vignette(“airway”)。
data("airway")
se <- airway
同样，我们也要指定untrt是dex变量的参考级别（reference level）：
se$dex %<>% relevel("untrt")
se$dex
## [1] untrt trt untrt trt untrt trt untrt trt
## Levels: untrt trt
我们可以对数百万个比对岛基因上的fragments进行快速检查（round的第二个参数表示需要保留的小数点位数）
round( colSums(assay(se)) / 1e6, 1 )
## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520
## 20.6 18.8 25.3 15.2 24.4 30.8 19.1
## SRR1039521
## 21.2
假设我们用上一节描述的方法构建了一个SummarizedExperiment，我们现在需要确保该对象包含样本的所有必要信息，比如，存储在colData slot计数矩阵中的meta数据表：
colData(se)
## DataFrame with 8 rows and 9 columns
## SampleName cell dex albut Run avgLength
## <factor> <factor> <factor> <factor> <factor> <integer>
## SRR1039508 GSM1275862 N61311 untrt untrt SRR1039508 126
## SRR1039509 GSM1275863 N61311 trt untrt SRR1039509 126
## SRR1039512 GSM1275866 N052611 untrt untrt SRR1039512 126
## SRR1039513 GSM1275867 N052611 trt untrt SRR1039513 87
## SRR1039516 GSM1275870 N080611 untrt untrt SRR1039516 120
## SRR1039517 GSM1275871 N080611 trt untrt SRR1039517 126
## SRR1039520 GSM1275874 N061011 untrt untrt SRR1039520 101
## SRR1039521 GSM1275875 N061011 trt untrt SRR1039521 98
## Experiment Sample BioSample
## <factor> <factor> <factor>
## SRR1039508 SRX384345 SRS508568 SAMN02422669
## SRR1039509 SRX384346 SRS508567 SAMN02422675
## SRR1039512 SRX384349 SRS508571 SAMN02422678
## SRR1039513 SRX384350 SRS508572 SAMN02422670
## SRR1039516 SRX384353 SRS508575 SAMN02422682
## SRR1039517 SRX384354 SRS508576 SAMN02422673
## SRR1039520 SRX384357 SRS508579 SAMN02422683
## SRR1039521 SRX384358 SRS508580 SAMN02422677
在这里，我们看到这个object已经包含了一个信息全面的colData slot（因为我们已经提前为您准备好了），但是，当您使用自己的数据时，您必须在此阶段添加与实验相关的所有信息（如上图所示）。我们强烈建议您将实验信息保存为csv（以逗号分隔）或TSV格式（以制表符分割），该文件可以从Excel中导出，并可以进一步导入到SummarizedExperiment中（确保csv中的行对应于SummarizedExperiment中的列）。我们之前通过使用示例表中的列信息来指定BAM文件，从而确保了这一对应关系。
library("DESeq2")
dds <- DESeqDataSet(se, design = ~ cell + dex)
3.2 从计数矩阵开始
在本节中，我们将展示如何构建一个DESeqDataSet（假设我们只有一个计数矩阵和一个sample信息表）。
注意：如果您已经准备好了一个SummarizedExperiment，则应跳过本节。虽然在上一节，我们从SummarizedExperiment中构造了DESeqDataSet，但在本节中，为了重新构建回新的对象（仅用于演示目的），我们将先从SummarizedExperiment中提取两个独立的对象（计数矩阵和样本信息）。在应用中，计数矩阵要么从文件中读入，要么由R函数生成（比如Rsubread包中的featureCounts函数）。
我们可以通过accessor 函数来访问SummarizedExperiment对象中的信息。举例来说，为了查看实际数据，比如fragment计数，我们可以使用assay函数来完成（head函数将输出限制为前几行）：
countdata <- assay(se)
head(countdata, 3)
## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516
## ENSG00000000003 679 448 873 408 1138
## ENSG00000000005 0 0 0 0 0
## ENSG00000000419 467 515 621 365 587
## SRR1039517 SRR1039520 SRR1039521
## ENSG00000000003 1047 770 572
## ENSG00000000005 0 0 0
## ENSG00000000419 799 417 508
在上面这个计数矩阵中，每行代表一个Ensembl基因，每一列代表一个测序的RNA文库，并且在每个文库中，fragment的数量仅对应到有且一个基因上（体现为value值）。我们还有每个样本的信息（计数矩阵的列）。如果您使用其他的软件对reads进行计数，则必须要检查计数矩阵的列是否与样本信息表中的行相对应。
coldata <- colData(se)
我们现在具备了可供分析的全部要素：
- countdata：fragment计数的表
- coldata：样本信息表
现在，要想从计数矩阵和样本信息表中构建DESeqDataSet对象，我们使用一下命令：
ddsMat <- DESeqDataSetFromMatrix(countData = countdata,
colData = coldata,
design = ~ cell + dex)
演示完毕，接下来，我们将继续使用从SummarizedExperiment章节中生成的对象。
4. 探索性分析和可视化
此工作流程中有两个单独的路径（path）：第一部分，我们会对计数进行转换，以便可视化地探索不同sample之间的关系。第二部分，我们将回到原始计数，对数据进行统计测试。这是十分重要的，因为通过统计测试，我们可以在测量精度方面对原始的计数数据进行评估。
4.1 预过滤数据集
在DESeqDataSet计数矩阵中，有许多行全是零，还有许多行总共只有几个片段。为了减小对象的大小，并提高函数的速度，我们可以删除没有或几乎没有关于基因表达量的信息的行。在这里，我们应用最小的过滤规则：删除没有计数的DESeqDataSet行，或者只删除所有样本中的单个计数。
nrow(dds)
## [1] 64102
dds <- dds[ rowSums(counts(dds)) > 1, ]
nrow(dds)
## [1] 29391
4.2 方差稳定转换和rlog
很多统计方法都可以用于多维数据探索性分析（比如聚类和主成分分析（PCA））。对于RNA-seq计数，预期方差会随平均值的增加而增加。举例来说，如果直接在计数矩阵或标准化计数上进行PCA（例如校正测序深度的差异），则得到的可视化图通常取决于具有最高计数的基因，因为它们显示样品之间的最大绝对差异。我们可以采用对count值进行归一化处理的方法来避免这种情况;但是，在归一化以后，此时具有最低计数的基因又会对可视化图产生大量噪声，因为采用小计数的对数实际上会夸大它们的方差。我们可以使用一些模拟数据快速显示计数的这个属性（这里，泊松计数的范围为0.1到100的lambda）。
我们绘制每行（基因）与平均值的标准差：
lambda <- 10^seq(from = -1, to = 2, length = 1000)
cts <- matrix(rpois(1000*100, lambda), ncol = 100)
library("vsn")
meanSdPlot(cts, ranks = FALSE)
对于对数转换计数：
log.cts.one <- log2(cts + 1)
meanSdPlot(log.cts.one, ranks = FALSE)
当值接近0时，具有小对数值的基因会放大整张图的差异，而这种低信噪比的低计数基因将严重的影响PCA作图的精确性。在此，DESeq2为计数数据提供了两种转换，以稳定均值的方差：第一个是variance stabilizing transformation（VST），其具有色散均值趋势，可以在 vst 函数中实现；第二个则是 regularized-logarithm transformation 或 rlog。
对于具有高计数的基因，VST和rlog都能给出与普通log2转换类似的结果。然而，对于计数较低的基因，这些值会缩小到中间值。然后，VST或rlog转换的数据变得近似同一性（在meanSdPlot中可以显示出一个更平坦的趋势），并且其值直接用于计算样本之间的距离、制作PCA图或用于下游分析。
该选择哪种转换？ VST的计算速度要快得多，并且对高计数异常值的敏感度低于rlog。rlog倾向于在小数据集（n <30）上很好地工作，当样本中存在大范围的测序深度（一个数量级的差异）时，其可能优于VST。因此，我们建议将VST用于中到大型数据集（n> 30）。您可以执行这两种转换并比较生成的meanSdPlot或PCA图，如下所述。
请注意，DESeq2提供的两个转换是为除差异测试以外的应用程序提供的。若仅仅想进行差异测试，我们建议您将DESeq函数应用于原始计数（而不是转换后的数据，如本分析的后续流程）。
vst和rlog都返回一个基于SummarizedExperiment类的 DESeqTransform对象。转化的值不再进行计数，而是直接存储在assay中。附加到dds的colData仍然可以访问：
首先使用vsd方法：
vsd <- vst(dds, blind = FALSE)
head(assay(vsd), 3)
## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516
## ENSG00000000003 9.742074 9.430420 9.867627 9.645845 10.183143
## ENSG00000000419 9.333669 9.581707 9.486145 9.523093 9.427283
## ENSG00000000457 8.765274 8.698449 8.651473 8.732426 8.592787
## SRR1039517 SRR1039520 SRR1039521
## ENSG00000000003 9.880416 10.010366 9.639782
## ENSG00000000419 9.574860 9.325999 9.509246
## ENSG00000000457 8.702674 8.761945 8.724101
colData(vsd)
## DataFrame with 8 rows and 10 columns
## SampleName cell dex albut Run avgLength
## <factor> <factor> <factor> <factor> <factor> <integer>
## SRR1039508 GSM1275862 N61311 untrt untrt SRR1039508 126
## SRR1039509 GSM1275863 N61311 trt untrt SRR1039509 126
## SRR1039512 GSM1275866 N052611 untrt untrt SRR1039512 126
## SRR1039513 GSM1275867 N052611 trt untrt SRR1039513 87
## SRR1039516 GSM1275870 N080611 untrt untrt SRR1039516 120
## SRR1039517 GSM1275871 N080611 trt untrt SRR1039517 126
## SRR1039520 GSM1275874 N061011 untrt untrt SRR1039520 101
## SRR1039521 GSM1275875 N061011 trt untrt SRR1039521 98
## Experiment Sample BioSample sizeFactor
## <factor> <factor> <factor> <numeric>
## SRR1039508 SRX384345 SRS508568 SAMN02422669 1.02364764926706
## SRR1039509 SRX384346 SRS508567 SAMN02422675 0.896166721793923
## SRR1039512 SRX384349 SRS508571 SAMN02422678 1.17948612081678
## SRR1039513 SRX384350 SRS508572 SAMN02422670 0.670053829048202
## SRR1039516 SRX384353 SRS508575 SAMN02422682 1.17767135405022
## SRR1039517 SRX384354 SRS508576 SAMN02422673 1.39903646915342
## SRR1039520 SRX384357 SRS508579 SAMN02422683 0.920778683328161
## SRR1039521 SRX384358 SRS508580 SAMN02422677 0.944514089340919
接着使用rlog方法：
rld <- rlog(dds, blind = FALSE)
head(assay(rld), 3)
## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516
## ENSG00000000003 9.385681 9.052599 9.516877 9.285335 9.839093
## ENSG00000000419 8.869611 9.138274 9.036117 9.075296 8.972125
## ENSG00000000457 7.961373 7.881385 7.824075 7.921493 7.751690
## SRR1039517 SRR1039520 SRR1039521
## ENSG00000000003 9.530313 9.663260 9.277695
## ENSG00000000419 9.131828 8.861529 9.060906
## ENSG00000000457 7.886441 7.957126 7.912125
在上面的函数调用中，我们指定了blind = FALSE，这意味着细胞系和处理之间的差异（设计中的变量）将不会影响实验的预期方差-均值趋势。实验设计不直接用于转换，仅用于估计计数中的全局变异量。对于完全无监督（unsupervised）的转换，可以设置blind = TRUE（这是默认值）。
为了显示转换的效果，在下图中我们绘制了两个样本，第一个样本简单地使用log2函数（添加1后，以避免记录为零），第二个样本使用VST和rlog对数值进行转换。对于log2方法，我们需要首先对size factors进行评估，以考察测序深度，然后指定normalized = TRUE。而vst和rlog则会自动进行排序深度校正。
library("dplyr")
library("ggplot2")
dds <- estimateSizeFactors(dds)
df <- bind_rows(
as_data_frame(log2(counts(dds, normalized=TRUE)[, 1:2]+1)) %>%
mutate(transformation = "log2(x + 1)"),
as_data_frame(assay(vsd)[, 1:2]) %>% mutate(transformation = "vst"),
as_data_frame(assay(rld)[, 1:2]) %>% mutate(transformation = "rlog"))
colnames(df)[1:2] <- c("x", "y")
ggplot(df, aes(x = x, y = y)) + geom_hex(bins = 80) +
coord_fixed() + facet_grid( . ~ transformation)
两个样本的转换计数的散点图。显示的是使用标准化计数的log2变换（左），使用VST（中）和使用rlog（右）的散点图。虽然rlog的大小与log2计数大致相同，但VST对较小的值有向上的变化。这是样本之间的差异（在这些散点图中偏离y = x），这将有助于距离计算和PCA图。我们可以看到具有低计数的基因（左下角）在普通对数标度上似乎过度变化，而VST和rlog压缩低计数基因的差异，其中数据提供关于差异表达的很少信息。
4.3 样本距离
在RNA-seq分析中，第一步一般是评估各个sample之间的整体相似度：哪些sample之间是类似的，哪些是不同的？这些sample的数据符合实验设计的期望值吗？
我们使用R函数dist来计算样本之间的欧几里德距离。为了确保所有基因的贡献值相等，我们使用了VST数据。我们需要t函数来对矩阵进行转置（行-列转换），因为dist函数要求行为不同的sample，而列为不同的维度（比如样本中的基因）。
sampleDists <- dist(t(assay(vsd)))
sampleDists
## SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517
## SRR1039509 44.08478
## SRR1039512 36.82952 51.50926
## SRR1039513 59.23174 43.11327 47.21641
## SRR1039516 41.08424 54.79377 40.18856 59.72262
## SRR1039517 59.08409 47.46042 54.13979 46.36635 44.74738
## SRR1039520 37.27578 54.00938 34.56797 55.44564 42.81024 58.22766
## SRR1039521 59.26716 42.84573 54.13120 35.24746 60.64298 47.80456
## SRR1039520
## SRR1039509
## SRR1039512
## SRR1039513
## SRR1039516
## SRR1039517
## SRR1039520
## SRR1039521 48.24754
我们使用pheatmap函数（来自pheatmap包）对Sample distances进行可视化：
library("pheatmap")
library("RColorBrewer")
为了使用距离矩阵中按distances进行排序的行/列绘制样本距离矩阵，我们手动将sampleDists提供给pheatmap函数的clustering_distance参数。否则，pheatmap函数会假设距离矩阵包含数据值本身，并且会对距离矩阵中的行/列之间的距离进行计算，这不是我们所期望的。我们还使用了RColorBrewer包中的colorRampPalette函数将颜色设置为蓝色。
sampleDistMatrix <- as.matrix( sampleDists )
rownames(sampleDistMatrix) <- paste( vsd$dex, vsd$cell, sep = " - " )
colnames(sampleDistMatrix) <- NULL
colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
clustering_distance_rows = sampleDists,
clustering_distance_cols = sampleDists,
col = colors)
使用rlog转换值的样本到样本距离的热图。
请注意，我们更改了距离矩阵的行名称以包含治疗类型和患者编号（而不是样本ID），以便在查看热图时我们可以查看所有这些信息。
计算样本距离的另一个选择是使用在PoiClaClu包中的泊松距离。在计算样本之间的距离时，这种衡量计数之间差异的方法也考虑了计数的固有方差结构。PoissonDistance函数采用原始计数矩阵（未规范化），样本为行而不是列，因此我们需要在dds中对计数进行转置。
library("PoiClaClu")
poisd <- PoissonDistance(t(counts(dds)))
我们绘制了以下的热图：
samplePoisDistMatrix <- as.matrix( poisd$dd )
rownames(samplePoisDistMatrix) <- paste( dds$dex, dds$cell, sep=" - " )
colnames(samplePoisDistMatrix) <- NULL
pheatmap(samplePoisDistMatrix,
clustering_distance_rows = poisd$dd,
clustering_distance_cols = poisd$dd,
col = colors)
4.4 PCA图
另一种可视化样本到样本距离的方法是主成分分析（PCA）。在这种排序方法中，数据点（这里是样本）被映射到2D平面上，x轴是分离数据点最多的方向，此方向上的样本值写为PC1；y轴是将数据分成第二大的方向（必须与第一方向正交）。此方向上的样本值写为PC2。方向中包含的总方差的百分比将显示在轴标签中。请注意，这些百分比不会增加到100％，因为有更多维度包含剩余的差异（尽管这些剩余维度中的每一个都将解释少于我们看到的两个维度）。
plotPCA(vsd, intgroup = c("dex", "cell"))
使用VST数据的PCA图。每种独特的治疗组合和细胞系都有自己的颜色。
在这里，我们使用了 DESeq2 附带的函数plotPCA。intgroup指定的两个术语用于标记样本的两个不同组;它们将告诉函数使用该参数来选择颜色。我们还可以使用ggplot2包从头开始构建PCA图。（有关使用ggplot的更多详细信息，请参阅ggplot2文档）
pcaData <- plotPCA(vsd, intgroup = c( "dex", "cell"), returnData = TRUE)
pcaData
## PC1 PC2 group dex cell name
## SRR1039508 -16.228435 -2.3110688 untrt:N61311 untrt N61311 SRR1039508
## SRR1039509 8.618864 0.1055443 trt:N61311 trt N61311 SRR1039509
## SRR1039512 -10.025968 -4.9857836 untrt:N052611 untrt N052611 SRR1039512
## SRR1039513 16.424047 -4.0524552 trt:N052611 trt N052611 SRR1039513
## SRR1039516 -13.635225 13.0513834 untrt:N080611 untrt N080611 SRR1039516
## SRR1039517 10.153386 17.2437947 trt:N080611 trt N080611 SRR1039517
## SRR1039520 -11.948262 -10.4259788 untrt:N061011 untrt N061011 SRR1039520
## SRR1039521 16.641592 -8.6254361 trt:N061011 trt N061011 SRR1039521
percentVar <- round(100 * attr(pcaData, "percentVar"))
然后我们可以使用这些数据在下图中创建第二个图，我们将点的颜色映射到地塞米松处理，而将形状映射到不同的细胞系。
ggplot(pcaData, aes(x = PC1, y = PC2, color = dex, shape = cell)) +
geom_point(size =3) +
xlab(paste0("PC1: ", percentVar[1], "% variance")) +
ylab(paste0("PC2: ", percentVar[2], "% variance")) +
coord_fixed()
使用VST值和自定义ggplot2代码的PCA图。在这里我们指定细胞系（形状符号）和地塞米松处理（颜色）。
从PCA图中，我们看到细胞之间的差异（不同的形状）是相当大的，但地塞米松处理的差异更大（红色与蓝色）。这表明为什么通过使用配对设计（“配对”，因为每个dex处理的样品与来自相同细胞系的一个未处理的样品配对）在差异测试中考虑这一点是重要的。我们已经为此设计设置了早期分配公式~cell + dex。
4.5 MDS图
我们可以使用R中的MDS函数来绘制一个与PCA图非常相似的图。当我们只有距离矩阵，而没有数据矩阵时，MDS值就显得很有用。在这里，我们根据VST数据来计算距离的MDS值，并在下图中进行绘制：
mds <- as.data.frame(colData(vsd)) %>%
cbind(cmdscale(sampleDistMatrix))
ggplot(mds, aes(x = `1`, y = `2`, color = dex, shape = cell)) +
geom_point(size = 3) + coord_fixed()
使用VST数据的MDS图。
在下图中，我们显示了PoissonDistance的MDS图：
mdsPois <- as.data.frame(colData(dds)) %>%
cbind(cmdscale(samplePoisDistMatrix))
ggplot(mdsPois, aes(x = `1`, y = `2`, color = dex, shape = cell)) +
geom_point(size = 3) + coord_fixed()
使用泊松距离的MDS图。
5. 差异表达分析
5.1 运行差异表达的pipeline
由于我们在创建DESeqDataSet时已经指定了实验设计，因此我们可以调用DESeq函数来对原始计数进行差异表达分析：
dds <- DESeq(dds)
此功能将print执行的各个步骤。这些在DESeq的手册页中有更详细的描述，可以通过键入？DESeq来访问它们。
该函数会返回DESeqDataSet中所有符合拟合参数的结果，下面我们将描述如何从该object中提取出我们感兴趣的结果。
5.2 构建结果表
如果我们不带任何参数对results函数进行调用，那DESeq将提取design formula中最后一个变量的log2倍数变化和p值。如果此变量含有两个level（比如trt和untrt），则results函数将提取结果表，并对两个水平进行比较。这种comparison将print在输出结果的顶部：dex trt vs untrt。
res <- results(dds)
res
## log2 fold change (MLE): dex trt vs untrt
## Wald test p-value: dex trt vs untrt
## DataFrame with 29391 rows and 6 columns
## baseMean log2FoldChange lfcSE
## <numeric> <numeric> <numeric>
## ENSG00000000003 708.602169691234 -0.381253887429316 0.100654430187038
## ENSG00000000419 520.297900552084 0.206812715390391 0.112218674572537
## ENSG00000000457 237.163036796015 0.0379205923945955 0.143444716340181
## ENSG00000000460 57.9326331250967 -0.088167696263787 0.287141995230743
## ENSG00000000938 0.318098378392895 -1.378233965407 3.49987526565862
## ... ... ... ...
## ENSG00000273485 1.28644765243289 -0.127098087720873 1.60042118107398
## ENSG00000273486 15.4525365439045 -0.150988513292619 0.486499370223246
## ENSG00000273487 8.1632349843654 1.04640806907172 0.698967253900779
## ENSG00000273488 8.58447903624707 0.107831134309983 0.638099483385526
## ENSG00000273489 0.275899382507797 1.48372584344306 3.51394515550545
## stat pvalue padj
## <numeric> <numeric> <numeric>
## ENSG00000000003 -3.78775069036569 0.000152017272634607 0.00128121723830778
## ENSG00000000419 1.84294384315429 0.0653372100766794 0.196207903558039
## ENSG00000000457 0.264356843264037 0.791504963002102 0.911194055629064
## ENSG00000000460 -0.307052600205472 0.758803335537914 0.894634203883451
## ENSG00000000938 -0.39379516719652 0.693732272741941 NA
## ... ... ... ...
## ENSG00000273485 -0.0794153996609706 0.936702220169112 NA
## ENSG00000273486 -0.310357058064295 0.756289445616055 0.893661169189342
## ENSG00000273487 1.49707738557415 0.134373123316847 0.329056883636818
## ENSG00000273488 0.168987966794566 0.865806106383284 0.945182646119088
## ENSG00000273489 0.422239328669782 0.672850337762335 NA
我们可以使用以下更具体的命令等效地生成此结果表。因为dex是D实验设计中的最后一个变量，所以我们可以选择不使用contrast参数来提取dex的两个级别的comparison。
res <- results(dds, contrast=c("dex","trt","untrt"))
由于res是一个DataFrame对象，它携带的meda数据包含了列信息：
mcols(res, use.names = TRUE)
## DataFrame with 6 rows and 2 columns
## type description
## <character> <character>
## baseMean intermediate mean of normalized counts for all samples
## log2FoldChange results log2 fold change (MLE): dex trt vs untrt
## lfcSE results standard error: dex trt vs untrt
## stat results Wald statistic: dex trt vs untrt
## pvalue results Wald test p-value: dex trt vs untrt
## padj results BH adjusted p-values
第一列baseMean是标准化计数值的平均值除以size factors（取自DESeqDataSet中的所有样本)。剩余的四列指的是特定的contrast，即因子变量dex的trt水平与untrt水平的比较。我们将在下面的分析中阐述如何获得其他contrast。
log2FoldChange列是effect size的估计值。它告诉我们，与未经处理的样品相比，由于地塞米松处理，基因的表达会有多少变化。
当然，log2FoldChange列的估计值具有不确定性，其可在列lfcSE中获得，该列是log2倍变化估计的标准误差估计。作为统计检验的结果，我们还可以表示特定效应大小估计的不确定性。差异表达测试的目的是检测数据是否能提供足够的证据来支撑其差异性。DESeq2将对每个基因进行假设检验，以确定证据是否足以判断零假设（即不同处理方式对基因的影响为零，治疗和对照之间观察到的差异仅仅是由实验的变量所引起的）。与统计中常见的情况一样，此测试的结果将以p值呈现，可在pvalue列中找到。请记住，p值表示在零假设描述的情况下，倍数变化与观察到的变化一样强或甚至更强的概率。
我们还可以使用以下代码行汇总结果，这些代码会report一些其他信息，这些信息将在后面的章节中介绍。
summary(res)
##
## out of 29391 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up) : 2607, 8.9%
## LFC < 0 (down) : 2218, 7.5%
## outliers [1] : 0, 0%
## low counts [2] : 11397, 39%
## (mean count < 5)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results
注意，由于地塞米松处理的FDR水平为10％，有许多基因具有差异表达。这是可信的，因为已知气道的平滑肌细胞会对糖皮质激素类固醇起反应。就目前而言，有两种更严格的确定差异基因的方法：
-
降低错误发现率（FDR）阈值（结果表中padj上的阈值）
-
使用results函数中的lfcThreshold参数升高log2倍数变化阈值为0的
如果我们降低错误发现率（FDR）阈值，我们还应该inform results()函数，以便函数可以使用此阈值来执行最佳筛选：
res.05 <- results(dds, alpha = 0.05)
table(res.05$padj < 0.05)
##
## FALSE TRUE
## 12829 4026
如果我们想提高log2倍数变化的阈值，以便我们筛选由于治疗而显示出差异变化的基因，我们只需提供log2量表的值。例如，通过指定lfcThreshold = 1，我们将筛选出treatment对基因计数的显著影响超过一倍或小于一半的基因，因为\（2 ^ 1 = 2 \）。
resLFC1 <- results(dds, lfcThreshold=1)
table(resLFC1$padj < 0.1)
##
## FALSE TRUE
## 20034 240
有时，res中p值的一个子集是NA（即“not available”）。这是DESeq报告该基因的所有计数均为零的一种方式。此外，如果基因被排除在分析之外，则可以将p值指定为NA，因为它包含极端计数异常值。有关更多信息，请参阅DESeq2 vignette的异常值检测部分。
如果您在已发布的研究中使用R分析包的结果，您可以通过键入citation（“pkgName”）找到该软件的正确引用，您可以在其中用pkgName替换包的名称。引用方法论文有助于支持和奖励那些将时间投入开源软件进行基因组数据分析的人。
5.3 其他比较
通常，可以使用results函数的contrast参数来提取变量的任何两个水平的比较结果。用户应指定三个值：变量名，分子的level名以及分母的level名。在此，我们提取一个细胞系相对于另一个细胞系的倍数变化的log2的结果：
results(dds, contrast = c("cell", "N061011", "N61311"))
## log2 fold change (MLE): cell N061011 vs N61311
## Wald test p-value: cell N061011 vs N61311
## DataFrame with 29391 rows and 6 columns
## baseMean log2FoldChange lfcSE
## <numeric> <numeric> <numeric>
## ENSG00000000003 708.602169691234 0.306326362739978 0.143530866244254
## ENSG00000000419 520.297900552084 -0.0540467035087632 0.159716090785549
## ENSG00000000457 237.163036796015 0.0163084038114845 0.203025956979501
## ENSG00000000460 57.9326331250967 0.279127044592999 0.400667397006204
## ENSG00000000938 0.318098378392895 0 4.9975798104491
## ... ... ... ...
## ENSG00000273485 1.28644765243289 -2.32141260766577 2.35291638140544
## ENSG00000273486 15.4525365439045 -0.0727157623832606 0.711890912073144
## ENSG00000273487 8.1632349843654 -0.0269346906402055 1.00507452886325
## ENSG00000273488 8.58447903624707 1.09584759066511 0.898948246150507
## ENSG00000273489 0.275899382507797 0 4.9975798104491
## stat pvalue padj
## <numeric> <numeric> <numeric>
## ENSG00000000003 2.1342194244037 0.0328248233823354 0.214730145872613
## ENSG00000000419 -0.338392351346314 0.735067537267964 0.935493501438637
## ENSG00000000457 0.0803266934637877 0.935977428491719 0.989390053166383
## ENSG00000000460 0.696655247416293 0.486018572065076 0.84616382331195
## ENSG00000000938 0 1 NA
## ... ... ... ...
## ENSG00000273485 -0.986610755066083 0.323833495704888 NA
## ENSG00000273486 -0.102144529660451 0.918641956272752 0.985420987200879
## ENSG00000273487 -0.0267986998642467 0.978620290204837 NA
## ENSG00000273488 1.21903301481234 0.222831662001977 NA
## ENSG00000273489 0 1 NA
在DESeq 运行结束后，我们还提供了一些其他方法来比较结果表。如果需要交互项的结果，则应使用results函数中的name参数。有关构建结果表的其他方法的详细信息，请参阅results函数的帮助页面，results帮助页面的示例部分给出了一些相关的示例。
5.4 多重测验
在高通量生物学中，我们要注意：不能直接使用p值作为反对null的证据，而是要用多重测验进行矫正。如果我们简单地将p值的阈值设置为较低值（比如0.05），这会发生什么？在29391个基因中，有5677个基因的p值低于0.05：
sum(res$pvalue < 0.05, na.rm=TRUE)
## [1] 5677
sum(!is.na(res$pvalue))
## [1] 29391
现在，假设所有基因的零假设都是正确的，即没有基因受到地塞米松处理的影响。然后，通过p值的定义，我们预计高达5％的基因的p值低于0.05。这相当于，在29391个基因中，有1470个基因的p值低于0.05。如果我们只将p值低于0.05的基因选为差异表达基因，那么该列表将包含高达26％的假阳性（1470/5677）。
DESeq2中的p.adjust函数使用的是Benjamini-Hochberg（BH）调整;简而言之，这种方法为每个基因计算一个调整后的p值，该值可以回答以下问题：如果该基因的adj.P value小于或等于P value，那么结果的错误发现率（FDR）是多少？上述这些值称为BH调整过的p值，存在于res对象的padj这一列中。
对于许多高通量实验而言，FDR是一个有用的统计数据。如果我们认为10％假阳性是可以接受的，那我们可以判定：adj.P value低于0.1的所有基因都是显著的。那么会有多少这样的基因呢？
sum(res$padj < 0.1, na.rm=TRUE)
## [1] 4825
我们将结果表分配到这些基因中，然后通过log2倍变化对其进行排序，以获得具有最大下调表达的显著基因：
resSig <- subset(res, padj < 0.1)
head(resSig[ order(resSig$log2FoldChange), ])
## log2 fold change (MLE): dex trt vs untrt
## Wald test p-value: dex trt vs untrt
## DataFrame with 6 rows and 6 columns
## baseMean log2FoldChange lfcSE
## <numeric> <numeric> <numeric>
## ENSG00000128285 6.62474139908411 -5.32590514755406 1.25781646244798
## ENSG00000267339 26.2335725682209 -4.61155039131948 0.673094468700224
## ENSG00000019186 14.0876049307722 -4.32590664693342 0.857767767163718
## ENSG00000183454 5.80417114112773 -4.26407747502168 1.16687734166995
## ENSG00000146006 46.8075965458274 -4.21185023752769 0.528851947804938
## ENSG00000141469 53.4365283575179 -4.12478447879794 1.1297976633824
## stat pvalue padj
## <numeric> <numeric> <numeric>
## ENSG00000128285 -4.23424665406963 2.29319143028288e-05 0.000237557205506679
## ENSG00000267339 -6.85126769831372 7.31983899180854e-12 2.0484165290607e-10
## ENSG00000019186 -5.04321427376246 4.5777613182452e-07 6.60462531203733e-06
## ENSG00000183454 -3.65426366829546 0.000257921111014126 0.00204722065415641
## ENSG00000146006 -7.96413864978559 1.66378621897198e-15 7.14247557421908e-14
## ENSG00000141469 -3.65090547846338 0.000261317404823323 0.00206960624224951
……以及最大上调表达：
head(resSig[ order(resSig$log2FoldChange, decreasing = TRUE), ])
## log2 fold change (MLE): dex trt vs untrt
## Wald test p-value: dex trt vs untrt
## DataFrame with 6 rows and 6 columns
## baseMean log2FoldChange lfcSE
## <numeric> <numeric> <numeric>
## ENSG00000179593 67.2430476358842 9.50597241247416 1.05450220396067
## ENSG00000109906 385.071028862516 7.35262793318718 0.536390220039323
## ENSG00000250978 56.3181939282226 6.3273837397078 0.677797382672819
## ENSG00000132518 5.65465445433956 5.88511165460783 1.32404324478883
## ENSG00000127954 286.384119280247 5.20716032804828 0.493082760580173
## ENSG00000249364 8.83906075961539 5.09810727071616 1.15961389777208
## stat pvalue padj
## <numeric> <numeric> <numeric>
## ENSG00000179593 9.01465390662066 1.9749305904616e-19 1.25129933256218e-17
## ENSG00000109906 13.7076099796304 9.14198752722158e-43 2.25343730910719e-40
## ENSG00000250978 9.33521418267575 1.00787266739018e-20 7.19669078453133e-19
## ENSG00000132518 4.44480320244104 8.79723586811634e-06 9.98722159059214e-05
## ENSG00000127954 10.5604185429671 4.54630170772714e-26 5.04976252647174e-24
## ENSG00000249364 4.39638338287508 1.10069452799683e-05 0.000122258625535648
6. 绘图结果
6.1 Counts plot
可视化特定基因计数的一种快速方法是使用plotCounts函数，该函数将DESeqDataSet、基因名称和分组作为参数（如下图所示）。
topGene <- rownames(res)[which.min(res$padj)]
plotCounts(dds, gene = topGene, intgroup=c("dex"))
Normalized counts for a single gene over treatment group.
我们还可以使用 ggplot2 包中的 ggplot 函数进行绘图（如下图所示）。
library("ggbeeswarm")
geneCounts <- plotCounts(dds, gene = topGene, intgroup = c("dex","cell"),
returnData = TRUE)
ggplot(geneCounts, aes(x = dex, y = count, color = cell)) +
scale_y_log10() + geom_beeswarm(cex = 3)
ggplot(geneCounts, aes(x = dex, y = count, color = cell, group = cell)) +
scale_y_log10() + geom_point(size = 3) + geom_line()
Normalized counts with lines connecting cell lines.Note that the DESeq test actually takes into account the cell line effect, so this figure more closely depicts the difference being tested.
6.2 MA-plot
MA图主要展示了模型中估计系数的分布情况。在y轴上，“M”代表“减” - 对数值的减法等于比率的对数，而在x轴上，“A”代表“平均”。你可能会听到有人也称为MA图为均值-差异图或Bland-Altman图。
在绘制MA图之前，我们使用lfcShrink函数来缩小log2 fold changes ，以便更好地比较dex处理样本与未处理样本。在DESeq2中，有三种类型的shrinkage estimators，这里我们指定收缩系数的apeglm方法，这有利于收缩嘈杂的LFC估计，同时给出真正大差异的低偏差LFC估计。因为系数出现在resultsNames（dds）中，如果我们要使用apeglm，那我们需要通过命名或数字指定模型中要收缩的系数。
library("apeglm")
resultsNames(dds)
## [1] "Intercept" "cell_N061011_vs_N052611"
## [3] "cell_N080611_vs_N052611" "cell_N61311_vs_N052611"
## [5] "dex_trt_vs_untrt"
res <- lfcShrink(dds, coef="dex_trt_vs_untrt", type="apeglm")
plotMA(res, ylim = c(-5, 5))
如果需要指定对比度（其未出现在resultsNames（dds）中），则可以使用其他两种收缩方法中的任何一种，或者在某些情况下，重新分解相关变量并运行nbinomWaldTest后跟lfcShrink就足够了。有关详细信息，请参阅DESeq2 vignette。
由treatment造成的MA变化图。特定比较的log2倍数变化绘制在y轴上，并且通过尺寸因子归一化的计数的平均值显示在x轴上。每个基因用点表示。adj.p低于阈值的基因（此处为0.1，默认值）以红色显示。
DESeq2程序包使用贝叶斯程序来缓和（或“缩小”）具有非常低计数和高度可变计数的基因的log2倍变化，这可以从MA-左侧点的垂直扩展变窄看出。如上所示，lfcShrink函数执行此操作。有关缓和倍数变化的基本原理的详细说明，请参阅DESeq2论文。
如果我们没有使用统计调节来收缩嘈杂的log2倍变化，我们会看到以下情况：
我们也可以在MA图上标记各个点。这里我们使用with 函数为结果对象的选定行绘制圆还有文本信息。在with函数中，仅使用所选res行的baseMean和log2FoldChange值。
plotMA(res, ylim = c(-5,5))
topGene <- rownames(res)[which.min(res$padj)]
with(res[topGene, ], {
points(baseMean, log2FoldChange, col="dodgerblue", cex=2, lwd=2)
text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")
})
另一个有用的diagnostic图是p值的直方图（如下图所示）。该图最好通过排除计数非常小的基因而形成，否则会在直方图中产生峰值。
hist(res$pvalue[res$baseMean > 1], breaks = 0:20/20,
col = "grey50", border = "white")
Histogram of p values for genes with mean normalized count larger than 1.
6.3 基因聚类
在先前绘制的样本距离热图中，侧面的树状图显示了样本的层次聚类。我们也可以只对基因进行这种聚类。由于聚类仅与实际携带信号的基因相关，因此通常只聚类最高度可变基因的子集。在这里，为了演示，我们选择样本中方差最大的20个基因。我们将使用VST数据。
library("genefilter")
topVarGenes <- head(order(rowVars(assay(vsd)), decreasing = TRUE), 20)
如果我们不考虑基因的绝对表达强度，而是根据基因偏离平均值的量来绘制热图，那我会得到一个很有趣的结果。因此，我们将每个基因的值集中在样本上，并绘制热图（如下图所示）。我们提供了一个data.frame，它指示pheatmap函数如何标记列。
mat <- assay(vsd)[ topVarGenes, ]
mat <- mat - rowMeans(mat)
anno <- as.data.frame(colData(vsd)[, c("cell","dex")])
pheatmap(mat, annotation_col = anno)
样本间的相对的VST变换值的热图。治疗状态和细胞系信息显示在热图顶部的彩色条。每一个方格为不同患者的不同基因。请注意，热图顶部的一组基因将N061011细胞系与其他细胞系分开。在热图的中心，我们看到一组经地塞米松处理的样本，其基因具有更高的基因表达。
6.4 独立过滤
MA图突出了RNA-seq数据的重要特性。对于弱表达的基因，我们没有机会看到差异表达，因为低reads计数会遭受非常高的泊松噪声，任何生物效应都会有一定的概率从取样的不确定性中淹没。我们可以通过检查归一化计数后的基因的小p值（例如，小于0.05）的比率来证明这一点。我们将使用经过阈值处理的结果表来显示：当小p值很少的情况下，这看起来是什么样的。
qs <- c(0, quantile(resLFC1$baseMean[resLFC1$baseMean > 0], 0:6/6))
bins <- cut(resLFC1$baseMean, qs)
levels(bins) <- paste0("~", round(signif((qs[-1] + qs[-length(qs)])/2, 2)))
fractionSig <- tapply(resLFC1$pvalue, bins, function(p)
mean(p < .05, na.rm = TRUE))
barplot(fractionSig, xlab = "mean normalized count",
ylab = "fraction of small p values")
The ratio of small p values for genes binned by mean normalized count.p值来自log2倍数变化大于1或小于-1的测验。该图表明，平均计数非常低的基因很少或没有power，最好从将它们从test中排除。
乍一看，过滤这些基因似乎没什么好处。但是呢，这些基因会对多重测验的调整产生影响，如果去除这些基因，多重测验的结果会好一些。通过从FDR程序的输入中去除低计数基因，我们可以发现更多基因在我们保留的基因中具有重要意义，从而提高了我们测验的power。这种方法称为独立过滤。
DESeq2软件会自动执行独立过滤，使得adj.p的于临界值的基因数量最大化（默认情况下，alpha设置为0.1）。这种自动独立过滤由results函数执行。
独立一词突出了一个重要的警告。只有当我们过滤的统计量（这里是所有样本的归一化计数的平均值）与零假设下的实际检验统计量（p值）无关时，才允许这种过滤。否则，过滤将使测验无效，从而使BH过程的假设无效。DESeq2内部使用的独立过滤软件来自genefilter软件包。
7. 注释并导出结果文件
到目前为止，我们的结果表中仅包含Ensembl基因ID，但使用一些其他的基因名称则可以提供更好的注释信息。Bioconductor的注释包可以将各种ID相互mapping。首先，我们加载AnnotationDbi包和注释包org.Hs.eg.db：
library("AnnotationDbi")
library("org.Hs.eg.db")
这是Homo sapiens（“Hs”）的organ注释包（“org”），组织为AnnotationDbi数据库包（“db”），使用Entrez Gene ID（“eg”）作为主键。要获取所有可用键类型的列表，请使用：
columns(org.Hs.eg.db)
## [1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT"
## [5] "ENSEMBLTRANS" "ENTREZID" "ENZYME" "EVIDENCE"
## [9] "EVIDENCEALL" "GENENAME" "GO" "GOALL"
## [13] "IPI" "MAP" "OMIM" "ONTOLOGY"
## [17] "ONTOLOGYALL" "PATH" "PFAM" "PMID"
## [21] "PROSITE" "REFSEQ" "SYMBOL" "UCSCKG"
## [25] "UNIGENE" "UNIPROT"
我们可以使用mapIds函数将各列添加到结果表中。我们将结果表的行名称作为键，并指定keytype = ENSEMBL。column参数告诉mapIds函数我们想要哪些信息，multiVals参数告诉函数如果单个输入值有多个可能值，该怎么办。在这里，我们要求只返回数据库中出现的第一个。要添加基因符号和Entrez ID，我们调用mapIds两次。
res$symbol <- mapIds(org.Hs.eg.db,
keys=row.names(res),
column="SYMBOL",
keytype="ENSEMBL",
multiVals="first")
res$entrez <- mapIds(org.Hs.eg.db,
keys=row.names(res),
column="ENTREZID",
keytype="ENSEMBL",
multiVals="first")
现在结果具有所需的外部基因ID：
resOrdered <- res[order(res$pvalue),]
head(resOrdered)
## log2 fold change (MAP): dex trt vs untrt
## Wald test p-value: dex trt vs untrt
## DataFrame with 6 rows and 7 columns
## baseMean log2FoldChange lfcSE
## <numeric> <numeric> <numeric>
## ENSG00000152583 997.439773207048 4.56334706476054 0.1856198978357
## ENSG00000165995 495.092906698546 3.28035231082466 0.133330866932756
## ENSG00000120129 3409.02937530956 2.93799351420539 0.122027770833094
## ENSG00000101347 12703.3870621821 3.75351791746965 0.157054428118604
## ENSG00000189221 2341.76725275591 3.34218643467603 0.14309394729784
## ENSG00000211445 12285.6151498691 3.71695196968584 0.167889677028122
## pvalue padj symbol
## <numeric> <numeric> <character>
## ENSG00000152583 2.22231974839155e-136 3.99884215525576e-132 SPARCL1
## ENSG00000165995 7.83975720860937e-135 7.05342956058585e-131 CACNB2
## ENSG00000120129 3.66732005915504e-130 2.19965857148119e-126 DUSP1
## ENSG00000101347 9.58233815136126e-130 4.31061481738986e-126 SAMHD1
## ENSG00000189221 1.09937074995478e-123 3.95641545493725e-120 MAOA
## ENSG00000211445 4.6181378193729e-112 1.38497953202993e-108 GPX3
## entrez
## <character>
## ENSG00000152583 8404
## ENSG00000165995 783
## ENSG00000120129 1843
## ENSG00000101347 25939
## ENSG00000189221 4128
## ENSG00000211445 2878
7.1 导出结果
我们可以轻松地将结果表以CSV文件的形式进行保存，然后我们可以使用Excel等电子表格程序共享或加载该文件。我们需要调用as.data.frame才能将DataFrame对象（ IRanges 包）转换为可由write.csv处理的data.frame对象。在这里，我们只采用前100个基因进行演示。
resOrderedDF <- as.data.frame(resOrdered)[1:100, ]
write.csv(resOrderedDF, file = "results.csv")
一种更复杂的导出方式是使用Bioconductor包ReportingTools。 ReportingTools将自动生成动态HTML文档，包括使用基因标识符（gene identifiers）链接到外部数据库，以及总结各组标准化计数的boxplot图。有关完整详细信息，请参阅ReportingTools *vignettes。我以下代码是创建ReportingTools*报告的最简单版本：
library("ReportingTools")
htmlRep <- HTMLReport(shortName="report", title="My report",
reportDirectory="./report")
publish(resOrderedDF, htmlRep)
url <- finish(htmlRep)
browseURL(url)
7.2 绘制基因组空间中的倍数变化
如果我们使用summarizeOverlaps函数来计算读数（count），那么我们的DESeqDataSet对象将建立在指定基因的基因组坐标的即用型Bioconductor对象之上。因此，我们可以轻松地在基因组空间中绘制差异表达结果。虽然results函数默认返回DataFrame，但使用format参数，我们也可以输出GRanges或GRangesList。lfcShrink还没有输出GRanges的选项，所以我们手动添加收缩系数列。
resGR <- results(dds, name="dex_trt_vs_untrt", format="GRanges")
resGR$log2FoldChange <- res$log2FoldChange
resGR
## GRanges object with 29391 ranges and 6 metadata columns:
## seqnames ranges strand | baseMean
## <Rle> <IRanges> <Rle> | <numeric>
## ENSG00000000003 X 99883667-99894988 - | 708.602169691234
## ENSG00000000419 20 49551404-49575092 - | 520.297900552084
## ENSG00000000457 1 169818772-169863408 - | 237.163036796015
## ENSG00000000460 1 169631245-169823221 + | 57.9326331250967
## ENSG00000000938 1 27938575-27961788 - | 0.318098378392895
## ... ... ... ... . ...
## ENSG00000273485 10 105209953-105210609 + | 1.28644765243289
## ENSG00000273486 3 136556180-136557863 - | 15.4525365439045
## ENSG00000273487 1 92654794-92656264 + | 8.1632349843654
## ENSG00000273488 3 100080031-100080481 + | 8.58447903624707
## ENSG00000273489 7 131178723-131182453 - | 0.275899382507797
## log2FoldChange lfcSE stat
## <numeric> <numeric> <numeric>
## ENSG00000000003 -0.3584541167481 0.100654430187038 -3.78775069036569
## ENSG00000000419 0.186072582398103 0.112218674572537 1.84294384315429
## ENSG00000000457 0.031346823539879 0.143444716340181 0.264356843264037
## ENSG00000000460 -0.0472622279238664 0.287141995230743 -0.307052600205472
## ENSG00000000938 -0.0123430867945908 3.49987526565862 -0.39379516719652
## ... ... ... ...
## ENSG00000273485 -0.00445615139561941 1.60042118107398 -0.0794153996609706
## ENSG00000273486 -0.0430928313534234 0.486499370223246 -0.310357058064295
## ENSG00000273487 0.191163979777062 0.698967253900779 1.49707738557415
## ENSG00000273488 0.0219401136602299 0.638099483385526 0.168987966794566
## ENSG00000273489 0.0121980340179439 3.51394515550545 0.422239328669782
## pvalue padj
## <numeric> <numeric>
## ENSG00000000003 0.000152017272634607 0.00128121723830778
## ENSG00000000419 0.0653372100766794 0.196207903558039
## ENSG00000000457 0.791504963002102 0.911194055629064
## ENSG00000000460 0.758803335537914 0.894634203883451
## ENSG00000000938 0.693732272741941 <NA>
## ... ... ...
## ENSG00000273485 0.936702220169112 <NA>
## ENSG00000273486 0.756289445616055 0.893661169189342
## ENSG00000273487 0.134373123316847 0.329056883636818
## ENSG00000273488 0.865806106383284 0.945182646119088
## ENSG00000273489 0.672850337762335 <NA>
## -------
## seqinfo: 722 sequences (1 circular) from an unspecified genome
我们需要再次添加符号来标记图上的基因：
resGR$symbol <- mapIds(org.Hs.eg.db, names(resGR), "SYMBOL", "ENSEMBL")
我们将使用Gviz软件包绘制GRanges和相关元数据：由于地塞米松处理所导致的log倍数变化。
library("Gviz")
以下代码块指定具有最小p值的基因上游和下游的100万个碱基对的窗口。我们为窗口内的基因创建了完整结果的子集。如果gene symbol存在且我们的子集中没有重复，我们将gene symbol添加为一个name。
window <- resGR[topGene] + 1e6
strand(window) <- "*"
resGRsub <- resGR[resGR %over% window]
naOrDup <- is.na(resGRsub$symbol) | duplicated(resGRsub$symbol)
resGRsub$group <- ifelse(naOrDup, names(resGRsub), resGRsub$symbol)
我们创建一个向量，指定该子集中的基因是否具有较低的padj值。
status <- factor(ifelse(resGRsub$padj < 0.1 & !is.na(resGRsub$padj),
"sig", "notsig"))
然后我们可以使用Gviz函数绘制结果（如下图所示）。我们创建一个轴轨道，指定我们在基因组中的位置，一条轨道将显示基因及其名称，按重要性着色，数据轨道将绘制垂直条，显示DESeq2产生的log倍数变化（我们知道，当计数中的信息很好地支持这种效果时，它们才会很大）。
options(ucscChromosomeNames = FALSE)
g <- GenomeAxisTrack()
a <- AnnotationTrack(resGRsub, name = "gene ranges", feature = status)
d <- DataTrack(resGRsub, data = "log2FoldChange", baseline = 0,
type = "h", name = "log2 fold change", strand = "+")
plotTracks(list(g, d, a), groupAnnotation = "group",
notsig = "grey", sig = "hotpink")
log2 fold changes in genomic region surrounding the gene with smallest adjusted p value. Genes highlighted in pink have adjusted p value less than 0.1.
8. 删除隐藏的batch effect
假设我们不知道实验中涉及不同的细胞系，只知道它们都被用地塞米松处理。那么这些隐藏的变异，则会对细胞系的计数结果产生影响，进而影响数据集（dataset）中的部分或全部基因。我们可以使用来自sva包用于RNA-seq的统计方法或Bioconductor中的RUVSeq包来检测样本的这种分组；为了进一步解释这些变异，我们可以将它们添加到DESeqDataSet design中。
对于我们的分析，SVA软件包使用surrogate variables这个术语来表示我们想要考虑的估计变量，而RUV软件包使用factors of unwanted variation术语和首字母缩写“删除不需要的变化”来解释软件包标题。我们首先使用SVA查找隐藏的batch effect，然后使用RUV。
8.1 将SVA与DESeq2一起使用
library("sva")
在下面的步骤中，将我们获得一个标准化计数矩阵，其中样本的平均计数大于1。如上所述，我们试图恢复任何隐藏的batch effect（假设我们不知道细胞系信息）。因此，我们使用带有dex变量的完整模型矩阵，以及仅具有截距项的简化或零模型矩阵。最后，我们指定我们想要估计2个代理变量。有关更多信息，请通过键入 ？svaseq 来阅读svaseq函数的手册页。
dat <- counts(dds, normalized = TRUE)
idx <- rowMeans(dat) > 1
dat <- dat[idx, ]
mod <- model.matrix(~ dex, colData(dds))
mod0 <- model.matrix(~ 1, colData(dds))
svseq <- svaseq(dat, mod, mod0, n.sv = 2)
## Number of significant surrogate variables is: 2
## Iteration (out of 5 ):1 2 3 4 5
svseq$sv
## [,1] [,2]
## [1,] 0.2481108 -0.52600157
## [2,] 0.2629867 -0.58115433
## [3,] 0.1502704 0.27428267
## [4,] 0.2023800 0.38419545
## [5,] -0.6086586 -0.07854931
## [6,] -0.6101210 -0.02923693
## [7,] 0.1788509 0.25708985
## [8,] 0.1761807 0.29937417
因为我们预先知道这些细胞系的信息，我们可以看到SVA方法在恢复这些变量方面的表现如何（见下图）。
par(mfrow = c(2, 1), mar = c(3,5,3,1))
for (i in 1:2) {
stripchart(svseq$sv[, i] ~ dds$cell, vertical = TRUE, main = paste0("SV", i))
abline(h = 0)
}
Surrogate variables 1 and 2 plotted over cell line. Here, we know the hidden source of variation (cell line), and therefore can see how the SVA procedure is able to identify a source of variation which is correlated with cell line.
最后，为了使用SVA从我们的代理变量中删除对计数的任何影响，我们只需将这两个代理变量作为列添加到DESeqDataSet中，然后将它们添加到design中：
ddssva <- dds
ddssva$SV1 <- svseq$sv[,1]
ddssva$SV2 <- svseq$sv[,2]
design(ddssva) <- ~ SV1 + SV2 + dex
然后，我们可以通过使用新的design来运行DESeq，以生成控制代理变量的结果。
8.2 将RUV与DESeq2一起使用
我们还可以使用RUVSeq包中的RUV方法来检测隐藏的batch effect。
library("RUVSeq")
我们可以使用RUVg函数来估计不需要的变化因子，类似于SVA的替代变量。与上述 SVA 程序相比,RUVg函数的差异在于，我们首先运行DESeq并得到结果以获得分析的p值，而不知道batches，例如只是~dex。假设我们有这个结果表res，我们然后通过查看不具有小p值的基因来提取一组empirical control genes。
set <- newSeqExpressionSet(counts(dds))
idx <- rowSums(counts(set) > 5) >= 2
set <- set[idx, ]
set <- betweenLaneNormalization(set, which="upper")
not.sig <- rownames(res)[which(res$pvalue > .1)]
empirical <- rownames(set)[ rownames(set) %in% not.sig ]
set <- RUVg(set, empirical, k=2)
pData(set)
## W_1 W_2
## SRR1039508 -0.17841946 0.49379636
## SRR1039509 -0.20749027 0.61226276
## SRR1039512 -0.05113443 -0.07266363
## SRR1039513 -0.14974859 -0.07839422
## SRR1039516 0.58209686 -0.05446432
## SRR1039517 0.61285776 -0.04648514
## SRR1039520 -0.28519844 -0.42231712
## SRR1039521 -0.32296344 -0.43173467
我们可以绘制RUV估计的因子：
par(mfrow = c(2, 1), mar = c(3,5,3,1))
for (i in 1:2) {
stripchart(pData(set)[, i] ~ dds$cell, vertical = TRUE, main = paste0("W", i))
abline(h = 0)
}
Factors of unwanted variation plotted over cell line.
和上面一样，如果我们想控制这些因素，我们只需将它们添加到DESeqDataSet和design中：
ddsruv <- dds
ddsruv$W1 <- set$W_1
ddsruv$W2 <- set$W_2
design(ddsruv) <- ~ W1 + W2 + dex
然后我们将使用新design运行DESeq以重新估计参数和结果。
9. Time course实验
DESeq2可用于分析Time course实验，例如，与一组基线样本相比，找到那些随时间变化，以条件特异性方式进行响应的基因。在这里，我们展示了fission数据包中的Time course分析，其中包含裂殖酵母的RNA-seq数据（随时间变化的基因count）。酵母暴露于氧化应激，并且一半样品含有atf21基因的缺失。我们使用设计公式来模拟以下三种情况：时间为0时的strain差异，随时间变化而产生的差异以及随时间变化产生的特异性的strain差异（相互作用项：strain:minute）。
library("fission")
data("fission")
ddsTC <- DESeqDataSet(fission, ~ strain + minute + strain:minute)
以下代码块将执行似然比检验，其中，我们会随时间的变化，而移除特异性的strain差异。来自该测试的具有小p值的基因是在时间0之后的一个或多个时间点显示具有菌株特异性效应的基因。因此，请注意，对于在两种菌株中以相同方式随时间上移或下移的基因，这不会给出小的p值。
ddsTC <- DESeq(ddsTC, test="LRT", reduced = ~ strain + minute)
resTC <- results(ddsTC)
resTC$symbol <- mcols(ddsTC)$symbol
head(resTC[order(resTC$padj),], 4)
## log2 fold change (MLE): strainmut.minute180
## LRT p-value: '~ strain + minute + strain:minute' vs '~ strain + minute'
## DataFrame with 4 rows and 7 columns
## baseMean log2FoldChange lfcSE
## <numeric> <numeric> <numeric>
## SPBC2F12.09c 174.671161802578 -2.656719531675 0.752261272511001
## SPAC1002.18 444.504950375698 -0.0509321395109172 0.204299484857659
## SPAC1002.19 336.373206596558 -0.392748982380898 0.573494009172242
## SPAC1002.17c 261.773132733438 -1.13876476912795 0.606128756826878
## stat pvalue padj
## <numeric> <numeric> <numeric>
## SPBC2F12.09c 97.2833864058632 1.97415107530233e-19 1.33452612690438e-15
## SPAC1002.18 56.9535986216169 5.16955159517375e-11 1.74730843916873e-07
## SPAC1002.19 43.5339085771997 2.87980357331442e-08 6.48915738520182e-05
## SPAC1002.17c 39.315835535309 2.0513707836516e-07 0.000346681662437121
## symbol
## <character>
## SPBC2F12.09c atf21
## SPAC1002.18 urg3
## SPAC1002.19 urg1
## SPAC1002.17c urg2
这只是可应用于时间序列数据的测试之一。另一种选择是将计数建模为时间的平滑函数，并将条件的相互作用项包括在平滑函数中。可以使用R中的spline basis函数来构建这样的模型，而另一种更现代的方法是使用Gaussian processes。
我们可以使用ggplot2绘制group随时间变化的计数、绘制具有最小adj.p值的基因的计数，也可以绘制与条件相关的时间曲线，甚至是绘制count图来评估时间为0时的差异（下图）。请记住，交互条件（interaction terms）是在考虑时间0时的差异后，在给定时间两组之间的差异。
fiss <- plotCounts(ddsTC, which.min(resTC$padj),
intgroup = c("minute","strain"), returnData = TRUE)
ggplot(fiss,
aes(x = as.numeric(minute), y = count, color = strain, group = strain)) +
geom_point() + geom_smooth(se = FALSE, method = "loess") + scale_y_log10()
Normalized counts for a gene with condition-specific changes over time.
可以使用results函数的test参数来研究各个时间点的log2倍数变化的Wald测验：
resultsNames(ddsTC)
## [1] "Intercept" "strain_mut_vs_wt" "minute_15_vs_0"
## [4] "minute_30_vs_0" "minute_60_vs_0" "minute_120_vs_0"
## [7] "minute_180_vs_0" "strainmut.minute15" "strainmut.minute30"
## [10] "strainmut.minute60" "strainmut.minute120" "strainmut.minute180"
res30 <- results(ddsTC, name="strainmut.minute30", test="Wald")
res30[which.min(resTC$padj),]
## log2 fold change (MLE): strainmut.minute30
## Wald test p-value: strainmut.minute30
## DataFrame with 1 row and 6 columns
## baseMean log2FoldChange lfcSE
## <numeric> <numeric> <numeric>
## SPBC2F12.09c 174.671161802578 -2.6004690287543 0.634342916344087
## stat pvalue padj
## <numeric> <numeric> <numeric>
## SPBC2F12.09c -4.09946885470339 4.14099389609047e-05 0.279931187375715
我们还可以通过它们的表达谱来聚类一些重要的基因。我们使用coef函数提取缩小的log倍数变化的矩阵：
betas <- coef(ddsTC)
colnames(betas)
## [1] "Intercept" "strain_mut_vs_wt" "minute_15_vs_0"
## [4] "minute_30_vs_0" "minute_60_vs_0" "minute_120_vs_0"
## [7] "minute_180_vs_0" "strainmut.minute15" "strainmut.minute30"
## [10] "strainmut.minute60" "strainmut.minute120" "strainmut.minute180"
我们现在可以在热图中绘制log2倍数变化（如下图所示）。
topGenes <- head(order(resTC$padj),20)
mat <- betas[topGenes, -c(1,2)]
thr <- 3
mat[mat < -thr] <- -thr
mat[mat > thr] <- thr
pheatmap(mat, breaks=seq(from=-thr, to=thr, length=101),
cluster_col=FALSE)
Heatmap of log2 fold changes for genes with smallest adjusted p value. The bottom set of genes show strong induction of expression for the baseline samples in minutes 15-60 (red boxes in the bottom left corner), but then have slight differences for the mutant strain (shown in the boxes in the bottom right corner).
