108-如何利用DESeq2分析转录组数据？

Apr 20, 2019 12 min read cnposts

刘小泽写于19.4.20+4.25 来自于：https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#how-do-i-use-vst-or-rlog-data-for-differential-testing 这个教程还比较新，当前版本2019-4-12，基于DESeq2 的1.23.10版本，有兴趣可以自行探索英文说明书。我目前使用的版本是1.22.2
另外还有一个 An RNA-seq workflow 教程，但是其中包含的内容有点复杂

前言

做转录组RNA-seq的一个重要目的就是找到差异基因，一般我们会从公司那里或者自己上游分析得到的原始表达矩阵再进行下游分析，其中count值就是每个样本中比对到每个基因的reads数。能做差异基因筛选这件事的包除了DESeq2还有 edgeR, limma, DSS, EBSeq和 baySeq.

DESeq2是利用负二项分布广义线性模型( negative binomial generalized linear models），同时，还利用了离散型估计、logFoldChange

具体的模型含义可以搜索，但首先要明白的是负二项分布是一个离散分布，符合测序reads分布

快速了解

现有一个大体印象，假设现在有一个表达矩阵cts，样本信息coldata，其中的design 表示怎样设计样本的模型(这里表示既考虑样本条件condition，又考虑了批次效应batch，并且这两项都要是coldata 的因子型变量)

# 这里演示的就是普通的矩阵导入
dds <- DESeqDataSetFromMatrix(countData = cts,
                              colData = coldata,
                              design= ~ batch + condition)
dds <- DESeq(dds)
resultsNames(dds) # lists the coefficients
res <- results(dds, name="condition_trt_vs_untrt")
# or to shrink log fold changes association with condition:
res <- lfcShrink(dds, coef="condition_trt_vs_untrt", type="apeglm")

DESeq2可以和一些上游的定量软件兼容，比如：

如果用Salmon、Sailfish、kallisto 得到表达矩阵，那么就可以用DESeqDataSetFromTximport 导入
如果有htseq得到的，就利用DESeqDataSetFromHTSeq
如果有RangedSummarizedExperiment 得到的，利用DESeqDataSet， RangedSummarizedExperiment

输入数据要注意：

首先=》为什么用原始数据(没有标准化的数据)？

DESeq2规定：不论是RNA-seq还是其他高通量数据，都要用原始的整数型矩阵(matrix of integer values)，比如i行xj列的矩阵中取出来第i行，第j列，表示的就是：在sample j中有多少reads比对到了gene i 上；RNA-seq这里的行表示的是基因，同样类比到另一种高通量数据，比如ChIP-seq，矩阵的行就可能表示结合区域(binding regions)

官方解释不使用标准化数据原因如下：
The DESeq2 model internally corrects for library size, so transformed or normalized values such as counts scaled by library size should not be used as input.
就是为了利用自身算法校正文库差异

然后=》关于DESeqDataSet对象

DESeq2包会使用DESeqDataSet对象存储原始表达量以及中间的估算值，在代码中常常用dds表示。这个对象可以看做是 SummarizedExperiment 包的RangedSummarizedExperiment对象扩展，那么其中的Ranged表示含义就是表达矩阵的每一行与基因组区间(基因的外显子)相联系，这样可以方便下游探索(比如：找到差异表达基因的ChIP-seq最临近的峰)

另外，这个对象必须有一个design表达式，用于估算数据分布以及后期log2FC值，形式一般是~ variable1 + variable2，并且control一般在前，处理在后。

the formula expresses how the counts for each gene depend on the variables in colData

比如：用group list作为design表达式，那么group list的样子一般就是：

> group_list
[1] untrt trt   untrt trt   untrt trt   untrt trt  
Levels: trt untrt

有4种构建DESeqDataSet的方法，取决于上游得到counts矩阵的方法

利用tximport

支持Salmon、Sailfish、kallisto、RSEM的表达量数据，这些方法产生的数据特点如下：

校正了样本间基因长度的差异 (Trapnell et al. 2013)
不基于比对的方法节省内存与空间
避免将比对到多个基因同源区段的fragments丢掉，提高敏感度 (Robert and Watson 2015)

在进行上游Salmon分析时，这里推荐使用—gcBias参数，可以进行RNAseq数据的系统误差校正(Love, Hogenesch, and Irizarry 2016; Patro et al. 2017)，除非很明确数据中没有这种bias。这里为了测试选择了tximportData包中的Salmon定量数据quant.sf

注意这里寻找文件的方法，这样不用自己去找路径，利用system.file和file.path 很方便。另外一般样本数据中会构建好了condition信息，但是这里为了演示，手动构建一下。

首先读入样本信息

library("tximport")
library("readr")
library("tximportData")
dir <- system.file("extdata", package="tximportData")
samples <- read.table(file.path(dir,"samples.txt"), header=TRUE)
# 构建condition
samples$condition <- factor(rep(c("A","B"),each=3))
rownames(samples) <- samples$run
samples[,c("pop","center","run","condition")]

##           pop center       run condition
## ERR188297 TSI  UNIGE ERR188297         A
## ERR188088 TSI  UNIGE ERR188088         A
## ERR188329 TSI  UNIGE ERR188329         A
## ERR188288 TSI  UNIGE ERR188288         B
## ERR188021 TSI  UNIGE ERR188021         B
## ERR188356 TSI  UNIGE ERR188356         B

然后准备表达量信息，需要表达矩阵文件以及转录本、基因关系文件

# 找到各个样本对应的表达文件，并更改列名为样本名
files <- file.path(dir,"salmon", samples$run, "quant.sf.gz")
names(files) <- samples$run
# 之后读入一个转录本对应基因的文件
tx2gene <- read_csv(file.path(dir, "tx2gene.gencode.v27.csv"))
head(tx2gene)

# A tibble: 6 x 2
  TXNAME            GENEID           
  <chr>             <chr>            
1 ENST00000456328.2 ENSG00000223972.5
2 ENST00000450305.2 ENSG00000223972.5
3 ENST00000473358.1 ENSG00000243485.5
4 ENST00000469289.1 ENSG00000243485.5
5 ENST00000607096.1 ENSG00000284332.1
6 ENST00000606857.1 ENSG00000268020.3

接着就是导入了

> txi <- tximport(files, type="salmon", tx2gene=tx2gene)
reading in files with read_tsv
1 2 3 4 5 6 
summarizing abundance
summarizing counts
summarizing length

最后，利用txi对象和样本信息构建DESeqDataSet

library("DESeq2")
ddsTxi <- DESeqDataSetFromTximport(txi,
                                   colData = samples,
                                   design = ~ condition)
# ddsTxi对象就可以继续进行下游分析了

利用count matrix

如果之前已经有了表达矩阵，那么就可以主要使用 DESeqDataSetFromMatrix函数。主要需要三个东西：原始表达矩阵(例如featureCounts产生的)、样本信息(就是表达矩阵的列名，要求是一个数据框)、design表达式

这里利用 pasilla 包的示例数据进行演示，其中表达矩阵是cts，样本信息是coldata

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")]
# 检查一下表达量矩阵是否和样本信息的顺序相同
> 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是不会智能区分的，否则后面运行就会报错；除了让它们的顺序保持一致以外，还要让coldata的行名中fb字符去掉

# 先去掉coldata的行名中fb字符
> rownames(coldata) <- sub("fb", "", rownames(coldata))
> all(rownames(coldata) %in% colnames(cts))
[1] TRUE
# 然后判断现在表达矩阵与样本信息是否一致
> all(rownames(coldata) == colnames(cts))
[1] FALSE
# 重新按照coldata的行名排序即可
> cts <- cts[, rownames(coldata)]
> all(rownames(coldata) == colnames(cts))
[1] TRUE

接下来直接导入即可

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

另外如果有额外的基因feature信息，也可以使用S4的DataFrame函数添加到dds的metadata中

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

对于htseq-count

主要利用DESeqDataSetFromHTSeqCount 导入Htseq数据(虽然目前来讲featureCounts速度要比Htseq快得多，但是htseq依然许多人在用)

首先还是先指定数据位置

# 实际分析时可以这样
directory <- "/path/to/your/files/"
# 这里测试还是用pasilla包
directory <- system.file("extdata", package="pasilla",
                         mustWork=TRUE)
# 使用list.files()可以列出当前文件夹下有哪些文件
# 并找到文件名中包含有treated字符的
sampleFiles <- grep("treated",list.files(directory),value=TRUE)
# 直接在文件名基础上提取condition的信息，这样保证了表达量矩阵和样本的一一对应
sampleCondition <- sub("(.*treated).*","\\1",sampleFiles) #其中\\1代表括号中匹配到的内容
> sampleTable
        sampleName         fileName condition
1   treated1fb.txt   treated1fb.txt   treated
2   treated2fb.txt   treated2fb.txt   treated
3   treated3fb.txt   treated3fb.txt   treated
4 untreated1fb.txt untreated1fb.txt untreated
5 untreated2fb.txt untreated2fb.txt untreated
6 untreated3fb.txt untreated3fb.txt untreated
7 untreated4fb.txt untreated4fb.txt untreated

最后

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

利用SummarizedExperiment

library("airway")
data("airway")
se <- airway
library("DESeq2")
ddsSE <- DESeqDataSet(se, design = ~ cell + dex)
# design表达式：~ group + condition
> ddsSE
class: DESeqDataSet 
dim: 64102 8 
metadata(2): '' version
assays(1): counts
rownames(64102): ENSG00000000003 ENSG00000000005 ... LRG_98 LRG_99
rowData names(0):
colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
colData names(9): SampleName cell ... Sample BioSample

刘小泽写于19.4.26 前面导入了DESeq2需要的数据 https://www.jianshu.com/p/4fc6d40b0a09 那么接下来应该怎么再进行分析呢？

导入数据后首先初步过滤(Pre-filtering)

这不是必须，只是为了减少后面计算量，主要就是去除表达量非常少的行，比如设置阈值为每行表达量为10

keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

告诉DESeq2你想怎么比较

默认情况下，R会根据字母表顺序排列因子型变量。如果不告诉DESeq2哪组作为对照组的话，它会自动根据字母表顺序设置，这里有两种解决方案：

明确指出contrast 是怎么设置(下面进行详细探讨)

在design中设置factor level，必须在分析之前就设置好，不能分析进行中或完成后再设置，它的设置方法也有两种

利用factor：

dds$condition <- factor(dds$condition, levels = c("untreated","treated"))
    
# 看下前后对比
> dds$condition # 之前不设置时，默认control组(untreated)排在后面
[1] treated   treated   treated   untreated untreated untreated untreated
Levels: treated untreated
> factor(dds$condition, levels = c("untreated","treated")) # 设置后，control组(untreated)排在了前面
[1] treated   treated   treated   untreated untreated untreated untreated
Levels: untreated treated

利用relevel

dds$condition <- relevel(dds$condition, ref = "untreated")

另外，这个因子型的比较公式不是一成不变的，如果我们从原来的表达矩阵中去除了某些样本，那么比较公式中也应该去除，利用droplevels()可以很方便去掉没有样本的factor level

# 比如这里去掉了第一个样本，那么最后的比较公式就是
> droplevels(dds[,-1]$condition)
[1] treated   treated   untreated untreated untreated untreated
Levels: untreated treated

合并技术重复到同一个表达矩阵中[可选]

技术重复technical replicate就是指一个文库的多个测序的run

提供了collapseReplicates函数可以做，但是不要用这个去合并生物学重复

差异表达分析

标准的差异分析流程被写进了单独的函数DESeq ，然后结果利用results函数生成，提取出log2FC、P值、校正后的p值等6列。

如果没有额外的参数，log2FC和P值是默认对design公式中的最后一个变量或者最后一个因子与参考因子进行比较，举个例子

# 之前构建的时候，选择的design是condition这个因子型变量
dds <- DESeqDataSetFromMatrix(countData = cts,
                              colData = coldata,
                              design = ~ condition)
# condition又被手动设置成了参考在前，比较在后的样子
> dds$condition
[1] treated   treated   treated   untreated untreated untreated untreated
Levels: untreated treated
# 这里results利用的就是最后一个因子(treated)与参考因子(untreated)进行比较，算出的log2FC也就是log2(treated/untreated)
## !!!这里不能反，因为即使反过来也能得到结果，但与事实就是完全相反的
dds <- DESeq(dds)
res <- results(dds)
# 截取头信息看一看，的确是treated vs untreated 
> res
log2 fold change (MLE): condition treated vs untreated 
Wald test p-value: condition treated vs untreated 
DataFrame with 9921 rows and 6 columns
# 如果还不放心或者进行更正，可以手动设置(方法二选一)
res <- results(dds, name="condition_treated_vs_untreated")
res <- results(dds, contrast=c("condition","treated","untreated"))

log2FC(LFC) vs lfcShrink

The shrunken fold changes are useful for ranking genes by effect size and for visualization. –Michael Love (the developer of deseq2)

关于lfcShrink： New function lfcShrink() in DESeq2 这个函数2年前就已经开发出来了

**为何采用lfcShrink？**log2FC estimates do not account for the large dispersion we observe with low read counts.因此，两种数据特别需要：低表达量占比高的；数据特别分散的

As with the shrinkage of dispersion estimates, LFC shrinkage uses information from all genes to generate more accurate estimates.

举一个来自 DESeq2 paper的例子

左图中相同两组（6J和2J）中绿色和紫色基因的平均值是一样的，但是绿色基因的方差更小，离散程度更小，因此右图中unshrunken LFC estimate (绿色实线) 与 shrunken LFC estimate (绿色虚线)是基本一致的；但是紫色基因由于高离散度导致这两者差别较大。另外可以看出，即使两个基因的标准化count值相似，但也可能存在不同的LFC shrinkage

注意：Shrinking the log2 fold changes不会改变显著差异的基因总数
例如，如果要根据LFC值提取差异基因，需要shrunken values
另外，进行功能分析例如GSEA时，需要提供shrunken values

作者还是非常推荐使用lfcShrink的https://support.bioconductor.org/p/95695/

DESeq2的老版本1.16.0之前是默认进行shrink的，但是并不是所有数据都是需要shrink的，因此作者把这个功能设成了默认关闭。把原来DESeq中自带的 log2 fold change shrinkage移到了lfcShrink中，可以手动添加[最新版的DESeq2是1.22]
如果使用的是旧版本的，并且不需要shrink，可以这样关闭：
DESeq(dds, betaPrior=FALSE)
# 这个结果和下面直接计算的结果是一样的
log2 (normalized_counts_group1 / normalized_counts_group2)

# 对于新版本而言，可以这样手动操作
resultsNames(dds)
## [1] "Intercept"                      "condition_treated_vs_untreated"
resLFC <- lfcShrink(dds, coef="condition_treated_vs_untreated", type="apeglm")
# 选择的apeglm参数(Zhu, Ibrahim, and Love 2018)进行effect size shrinkage,which improves on the previous estimator

# resLFC相比res数据更加紧凑

# 发现shrink前后少了一列stat
> names(resLFC)
[1] "baseMean"       "log2FoldChange" "lfcSE"          "pvalue"        
[5] "padj"          
> names(res)
[1] "baseMean"       "log2FoldChange" "lfcSE"          "stat"          
[5] "pvalue"         "padj"

开启多线程

上面的分析步骤大约只需要30s，但是对于复杂的design公式、大样本数据，就可以采用多线程加速。DESeq、results、lfcShrink 都可以通过加载BiocParallel实现

library("BiocParallel")
# 先预定4个核，等需要的时候直接使用parallel=TRUE来调用
register(MulticoreParam(4))

探索p值与adjusted p值

作者给出的建议是：Need to filter on adjusted p-values, not p-values, to obtain FDR control. 10% FDR is common because RNA-seq experiments are often exploratory and having 90% true positives in the gene set is ok

先从小到大排个序：

resOrdered <- res[order(res$pvalue),]
# order函数是给出从小到大排序后的位置，默认decreasing = FALSE

利用summary可以探索一些基本的统计量

> summary(res)

out of 9921 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 518, 5.2%
LFC < 0 (down)     : 536, 5.4%
outliers [1]       : 1, 0.01%
low counts [2]     : 1539, 16%
(mean count < 6)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

看看有多少p值小于0.1的？

> sum(res$padj < 0.1, na.rm=TRUE) # 有多少p值小于0.1的？
[1] 1054

独立假设加权Independent hypothesis weighting

主要在过滤p值时使用，与Benjamini和Hochberg的方法相比，它可以为每个假设分配数据驱动的权重来进行多重测试，从而增强统计效果。IHW的输入是p值和协变量的两列表，其中协变量可以是任何连续值或者分类变量，它可以为每个假设检验的统计特性提供信息，而与零假设下的p值无关

利用IHW包，结果存储在metadata中

Ignatiadis, N., Klaus, B., Zaugg, J.B., Huber, W. (2016) Data-driven hypothesis weighting increases detection power in genome-scale multiple testing. Nature Methods, 13:7. 10.1038/nmeth.3885

library("IHW")
resIHW <- results(dds, filterFun=ihw)
> summary(resIHW)

out of 9921 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)       : 502, 5.1%
LFC < 0 (down)     : 538, 5.4%
outliers [1]       : 1, 0.01%
[1] see 'cooksCutoff' argument of ?results
see metadata(res)$ihwResult on hypothesis weighting
> sum(resIHW$padj < 0.1, na.rm=TRUE)
[1] 1040

探索结果

MA-plot

利用plotMA函数画出纵坐标为log2FC，横坐标为mean of normalized counts的散点图，如果adjusted p值小于0.1就标记颜色。可以设置临界值，超过的部分就标记为三角形

plotMA(res, ylim=c(-2,2))

对于shrunken的LFC值，MA-plot就更有用了。它直接将与LFC相关的噪音从低表达的基因中去除，而不需要人工设置阈值

plotMA(resLFC, ylim=c(-2,2))

RNA-Seq DESeq2