RNA-seq 差异分析的点点滴滴（2）
## 引言
[本系列](https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html "Source")将开展全新的转录组分析专栏，主要针对使用 `DESeq2`时可能出现的问题和方法进行展开。
## Tximeta：自动导入并附加元数据
Bioconductor 家族中的 tximeta 包，在 tximport 的基础上进行了扩展，不仅保留了原有功能，还增加了一项新特性：能够自动为常用的转录组数据（包括人类和小鼠的 GENCODE、Ensembl、RefSeq）添加注释元数据。该包生成的 SummarizedExperiment 对象可以便捷地通过 DESeqDataSet 函数导入到 DESeq2 中，如下所示：
```R
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` 输出中提取这些数据对象。
```R
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)
```
检查计数矩阵和列数据，看看它们在样本顺序方面是否一致。
```R
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 treatedpaired-end
## treated3fb treatedpaired-end
## untreated1fb untreated single-read
## untreated2fb untreated single-read
## untreated3fb untreatedpaired-end
## untreated4fb untreatedpaired-end
```
请留意，这些数据的顺序与样本的顺序并不一致！
非常重要的一点是，计数矩阵的列顺序和样本信息（列数据的行）必须匹配。DESeq2 不会自动推断计数矩阵的哪一列对应于列数据的哪一行，这些信息在提供给 DESeq2 时必须是一致排序的。
由于它们没有按照正确的顺序排列，需要对其中一个进行重新排序，以确保它们在样本顺序上是一致的（如果不这样做，后续的操作将会出现错误）。此外，还需要将 coldata 的行名中的 "fb" 删除，以保持命名的一致性。
```R
rownames(coldata) <- sub("fb", "", rownames(coldata))
all(rownames(coldata) %in% colnames(cts))
## TRUE
all(rownames(coldata) == colnames(cts))
## FALSE
cts <- cts[, rownames(coldata)]
all(rownames(coldata) == colnames(cts))
## TRUE
```
如果您之前使用了 Rsubread 包中的 featureCounts 函数（Liao, Smyth, 和 Shi 2013），可以直接从该函数输出的列表中的 "counts" 项获取读数计数矩阵。通常情况下，计数矩阵和样本信息可以通过 R 基础函数如 read.csv 或 read.delim 从文本文件中导入。对于 htseq-count 文件，请参阅下面的专门输入函数。
拥有计数矩阵 `cts` 和样本信息 `coldata` 后，就可以构建一个 DESeqDataSet 对象：
```R
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 的行名存在了。）
```R
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 输出文件的目录。
```R
directory <- "/path/to/your/files/"
```
但是，仅出于演示目的，以下代码行指向 pasilla 包的演示 htseq-count 输出文件包的目录。
```R
directory <- system.file("extdata", package="pasilla",
mustWork=TRUE)
```
通过 `list.files` 函数来指定需要导入的文件，并利用 `grep` 函数筛选出包含 "treated" 字符串的文件。接着，使用 `sub` 函数对样本文件名进行拆分，以获取样本的条件状态；或者，您也可以选择使用 `read.table` 函数直接导入一个包含表型信息的表格。
```R
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：
```R
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
```
页:
[1]