RNA-seq 差异分析的点点滴滴（2）

引言

本系列^[1]将开展全新的转录组分析专栏，主要针对使用DESeq2时可能出现的问题和方法进行展开。

Tximeta：自动导入并附加元数据

Bioconductor 家族中的 tximeta 包，在 tximport 的基础上进行了扩展，不仅保留了原有功能，还增加了一项新特性：能够自动为常用的转录组数据（包括人类和小鼠的 GENCODE、Ensembl、RefSeq）添加注释元数据。该包生成的 SummarizedExperiment 对象可以便捷地通过 DESeqDataSet 函数导入到 DESeq2 中，如下所示：

coldata <- samples
coldata$files <- files
coldata$names <- coldata$run

library("tximeta")
se <- tximeta(coldata)
ddsTxi <- DESeqDataSet(se, design = ~ condition)

这个 ddsTxi 对象接下来可以在分析流程中作为 dds 对象使用。如果 tximeta 能够识别出参考转录组是预设的几种之一，并且具有预先计算好的哈希校验和，那么 dds 对象中的 rowRanges 将会被自动填充。

计数矩阵输入

另外，如果你已经有了从其他来源准备好的读数计数矩阵，可以使用 DESeqDataSetFromMatrix 函数。快速从比对文件生成计数矩阵的另一种方法是使用 Rsubread 包中的 featureCounts 函数。使用 DESeqDataSetFromMatrix 时，用户需要提供计数矩阵、样本信息（计数矩阵的列）以 DataFrame 或 data.frame 的形式，以及设计公式。

为了展示如何使用 DESeqDataSetFromMatrix，将从 pasilla 包中导入计数数据。导入一个计数矩阵，并将其命名为 cts，同时导入样本信息表，并将其命名为 coldata。在后续部分，会描述如何从例如 featureCounts 输出中提取这些数据对象。

library("pasilla")
pasCts <- system.file("extdata",
                      "pasilla_gene_counts.tsv",
                      package="pasilla", mustWork=TRUE)
pasAnno <- system.file("extdata",
                       "pasilla_sample_annotation.csv",
                       package="pasilla", mustWork=TRUE)
cts <- as.matrix(read.csv(pasCts,sep="\t",row.names="gene_id"))
coldata <- read.csv(pasAnno, row.names=1)
coldata <- coldata[,c("condition","type")]
coldata$condition <- factor(coldata$condition)
coldata$type <- factor(coldata$type)

检查计数矩阵和列数据，看看它们在样本顺序方面是否一致。

head(cts,2)

##             untreated1 untreated2 untreated3 untreated4 treated1 treated2
## FBgn0000003          0          0          0          0        0        0
## FBgn0000008         92        161         76         70      140       88
##             treated3
## FBgn0000003        1
## FBgn0000008       70

coldata

##              condition        type
## treated1fb     treated single-read
## treated2fb     treated  paired-end
## treated3fb     treated  paired-end
## untreated1fb untreated single-read
## untreated2fb untreated single-read
## untreated3fb untreated  paired-end
## untreated4fb untreated  paired-end

请留意，这些数据的顺序与样本的顺序并不一致！

非常重要的一点是，计数矩阵的列顺序和样本信息（列数据的行）必须匹配。DESeq2 不会自动推断计数矩阵的哪一列对应于列数据的哪一行，这些信息在提供给 DESeq2 时必须是一致排序的。

由于它们没有按照正确的顺序排列，需要对其中一个进行重新排序，以确保它们在样本顺序上是一致的（如果不这样做，后续的操作将会出现错误）。此外，还需要将 coldata 的行名中的 "fb" 删除，以保持命名的一致性。

rownames(coldata) <- sub("fb", "", rownames(coldata))
all(rownames(coldata) %in% colnames(cts))

## [1] TRUE

all(rownames(coldata) == colnames(cts))

## [1] FALSE

cts <- cts[, rownames(coldata)]
all(rownames(coldata) == colnames(cts))

## [1] TRUE

如果您之前使用了 Rsubread 包中的 featureCounts 函数（Liao, Smyth, 和 Shi 2013），可以直接从该函数输出的列表中的 "counts" 项获取读数计数矩阵。通常情况下，计数矩阵和样本信息可以通过 R 基础函数如 read.csv 或 read.delim 从文本文件中导入。对于 htseq-count 文件，请参阅下面的专门输入函数。

拥有计数矩阵 cts 和样本信息 coldata 后，就可以构建一个 DESeqDataSet 对象：

library("DESeq2")
dds <- DESeqDataSetFromMatrix(countData = cts,
                              colData = coldata,
                              design = ~ condition)
dds

## class: DESeqDataSet 
## dim: 14599 7 
## metadata(1): version
## assays(1): counts
## rownames(14599): FBgn0000003 FBgn0000008 ... FBgn0261574 FBgn0261575
## rowData names(0):
## colnames(7): treated1 treated2 ... untreated3 untreated4
## colData names(2): condition type

如果您拥有额外的特征数据，可以通过将这些数据添加到新创建对象的元数据列中，进而将它们整合到 DESeqDataSet 中。（此处为了演示目的添加了一些重复的数据，实际上基因名称已经作为 dds 的行名存在了。）

featureData <- data.frame(gene=rownames(cts))
mcols(dds) <- DataFrame(mcols(dds), featureData)
mcols(dds)

## DataFrame with 14599 rows and 1 column
##                    gene
##             <character>
## FBgn0000003 FBgn0000003
## FBgn0000008 FBgn0000008
## FBgn0000014 FBgn0000014
## FBgn0000015 FBgn0000015
## FBgn0000017 FBgn0000017
## ...                 ...
## FBgn0261571 FBgn0261571
## FBgn0261572 FBgn0261572
## FBgn0261573 FBgn0261573
## FBgn0261574 FBgn0261574
## FBgn0261575 FBgn0261575

htseq-count 数据输入

如果您之前使用了 HTSeq python 包中的 htseq-count 工具，那么可以通过 DESeqDataSetFromHTSeqCount 函数来处理数据。首先，您需要设置一个变量，指向存放 htseq-count 输出文件的目录。

directory <- "/path/to/your/files/"

但是，仅出于演示目的，以下代码行指向 pasilla 包的演示 htseq-count 输出文件包的目录。

directory <- system.file("extdata", package="pasilla",
                         mustWork=TRUE)

通过 list.files 函数来指定需要导入的文件，并利用 grep 函数筛选出包含 "treated" 字符串的文件。接着，使用 sub 函数对样本文件名进行拆分，以获取样本的条件状态；或者，您也可以选择使用 read.table 函数直接导入一个包含表型信息的表格。

sampleFiles <- grep("treated",list.files(directory),value=TRUE)
sampleCondition <- sub("(.*treated).*","\\1",sampleFiles)
sampleTable <- data.frame(sampleName = sampleFiles,
                          fileName = sampleFiles,
                          condition = sampleCondition)
sampleTable$condition <- factor(sampleTable$condition)

然后使用以下函数构建 DESeqDataSet：

library("DESeq2")
ddsHTSeq <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable,
                                       directory = directory,
                                       design= ~ condition)
ddsHTSeq

## class: DESeqDataSet 
## dim: 70463 7 
## metadata(1): version
## assays(1): counts
## rownames(70463): FBgn0000003:001 FBgn0000008:001 ... FBgn0261575:001
##   FBgn0261575:002
## rowData names(0):
## colnames(7): treated1fb.txt treated2fb.txt ... untreated3fb.txt
##   untreated4fb.txt
## colData names(1): condition