DESeq2详解
先了解完整的分析流程和工具：http://blog.sina.com.cn/s/blog_9cf2d3640102x9kx.html
这是一个系列文，包括：从标准的workflow开始，到更高级的数据操作和workflow个性化，最后是DESeq2的统计学原理以及常见的question解答
本文介绍在差异表达分析之前的操作步骤，主要是DESeq2导入数据和预处理，DESeq2对导入不同数据类型兼容性很好。
导入数据
Why un-normalized counts? DESeq2要求输入的数据必须是没有标准化的raw count，第i行第j列代表着在样本j分配到featrue(i)
的reads数或fragemnts数（feature可以是基因、转录本、ChIP-seq等NGS测的对应区段bin、定量质谱中的某肽段序列）至于怎么获取count data参考：http://master.bioconductor.org/packages/release/workflows/html/rnaseqGene.html
因为DESeq2模型内部本身会对文库大小做校正，所以不要用已经对文库做过scaled的数据
关于DESeqDataSet对象
这个对象实际上是继承于RangedSummarizedExperiment对象的（来自SummarizedExperiment package）也就是说从assay里获取的每一行(counts)都可以和genomic ranges关系上(比如gene exon/transcript)
这个关联对于下游进一步个性化分析非常有用，比如可以通过DESeq2获取了差异表达基因，对比ChIP-seq数据寻找在基因组坐标上距离该差异基因最近的ChIP-seq peaks
modeling中重要的一步就是design formula,添加适当的变量（DESeq2基于负二项分布线性模型）默认是把control放在前面，感兴趣变量放在后面：design~control+variable
四种构建DESeqDataSet对象的方法
Transcript abundance files and tximport / tximeta
pipeline里推荐的一步是使用快速转录本定量工具获取counts，比如
- Salmon (Patro et al. 2017)
- Sailfish (Patro, Mount, and Kingsford 2014)
- kallisto (Bray et al. 2016)
- RSEM (Li and Dewey 2011)
再用定量数据导入工具tximport导入，直接可以为DESeq2所用
前三者是alignment-free或只是“轻度比对”，具有可以校正在不同样本中基因长度不同（不同isoform)速度快、省内存、省bam文件存储、灵敏度高的优点，参考https://www.bioinfo-scrounger.com/archives/411/
和 https://www.jianshu.com/p/6d4cba26bb60
Salmon软件使用说明：https://combine-lab.github.io/salmon/getting_started/
示例数据演示( gene-level DESeqDataSet object from Salmon quant.sf files)，分析自己的数据时dir更改为/path/to/quant/files，这里tximportData包只是为了提供演示数据而已，做自己的数据分析时不需要
library("tximport")
library("readr")
library("tximportData")
dir <- system.file("extdata", package="tximportData")
samples <- read.table(file.path(dir,"samples.txt"), header=TRUE)
samples$condition <- factor(rep(c("A","B"),each=3))
# pop center assay sample experiment run condition
# 1 TSI UNIGE NA20503.1.M_111124_5 ERS185497 ERX163094 ERR188297 A
# 2 TSI UNIGE NA20504.1.M_111124_7 ERS185242 ERX162972 ERR188088 A
# 3 TSI UNIGE NA20505.1.M_111124_6 ERS185048 ERX163009 ERR188329 A
# 4 TSI UNIGE NA20507.1.M_111124_7 ERS185412 ERX163158 ERR188288 B
# 5 TSI UNIGE NA20508.1.M_111124_2 ERS185362 ERX163159 ERR188021 B
# 6 TSI UNIGE NA20514.1.M_111124_4 ERS185217 ERX163062 ERR188356 B
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
#读取每个run数据
files <- file.path(dir,"salmon", samples$run, "quant.sf.gz")
names(files) <- samples$run
#把转录本和gene关联起来的table
tx2gene <- read_csv(file.path(dir, "tx2gene.gencode.v27.csv"))
#用tximport function导入转录本定量数据
#关于tximport函数的详细说明和怎么构建tx2gene把转录本和基因对应上，参考：http://bioconductor.org/packages/release/bioc/html/tximport.html
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
library("DESeq2")#依赖包非常丰富，其中就包括了GenomicRanges、SummarizedExperiment等
ddsTxi <- DESeqDataSetFromTximport(txi,
colData = samples,
design = ~ condition)
#using counts and average transcript lengths from tximport
#tximeta这个包还可以自动添加annotation metadata(会自动下载gtf文件)，不需要tx2gene，详见https://bioconductor.org/packages/tximeta
coldata <- samples
coldata$files <- files
coldata$names <- coldata$run
library("tximeta")
se <- tximeta(coldata)
ddsTxi <- DESeqDataSet(se, design = ~ condition)
#ddsTxi object here can then be used as dds in the following analysis steps
Count matrix input
关于DESeqDataSetFromMatrix函数, 需要提供的输入对象包括表达矩阵(count-matrix)、样本信息samples (the columns of the count matrix) ，样本信息制作成dataframe，然后就是design formula.
示例数据来自pasilla包http://bioconductor.org/packages/pasilla
library("pasilla")#示例数据来自pasilla包http://bioconductor.org/packages/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和count matrix的rowname是否完全一致
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
rownames(coldata) <- sub("fb", "", rownames(coldata))
all(rownames(coldata) %in% colnames(cts))#修改rownames
## [1] TRUE
all(rownames(coldata) == colnames(cts))#但order不一致
## [1] FALSE
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
#如果有额外的meatdata要添加，可以直接添加到DESeqDataSet，(这里只是为了演示）
featureData <- data.frame(gene=rownames(cts))
mcols(dds) <- DataFrame(mcols(dds), featureData)
mcols(dds)
# DataFrame with 14599 rows and 1 column
# gene
# <factor>
# FBgn0000003 FBgn0000003
# FBgn0000008 FBgn0000008
# FBgn0000014 FBgn0000014
# FBgn0000015 FBgn0000015
# FBgn0000017 FBgn0000017
# ... ...
# FBgn0261571 FBgn0261571
# FBgn0261572 FBgn0261572
# FBgn0261573 FBgn0261573
# FBgn0261574 FBgn0261574
# FBgn0261575 FBgn0261575
htseq-count input
如果是用HTSeq做的转录本定量，可以用DESeqDataSetFromHTSeqCount导入数据，同样这里用的演示数据：
directory <- system.file("extdata", package="pasilla",
mustWork=TRUE)
#用list.files选择要读入的数据，并规定只读取包含"treated"字符串的文件
sampleFiles <- grep("treated",list.files(directory),value=TRUE)
sampleCondition <- sub("(.*treated).*","\\1",sampleFiles)#condition用于制作colData
sampleTable <- data.frame(sampleName = sampleFiles,
fileName = sampleFiles,
condition = sampleCondition)#需要把每一个文件和condition对应起来
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 input
示例数据来自airway，这是一个非常经典的RNA-seq测试数据，导入后实际上是一个RangedSummarizedExperiment 对象，简称为"se"，如果不熟悉请参考：http://bioconductor.org/packages/release/data/experiment/html/airway.html
library("airway")
data("airway")
se <- airway
library("DESeq2")
ddsSE <- DESeqDataSet(se, design = ~ cell + dex)
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
数据预处理
尽管在使用DESeq2函数前过滤低count的gene并不是必须的，但预过滤数据的好处是，去除那些只有很少reads的行以后，可以减少dds的存储，增加程序运行速度。这里演示一下简单的过滤，只保留那些至少有10条reads的基因。至于更严格的过滤方法（ independent filtering ），以后我会在推文里写到。
keep <- rowSums(counts(dds)) >= 10#对行求和函数rowSums()
dds <- dds[keep,]
关于因子水平
和limma里面的contrast 类似的是，DESeq2同样最好说明清楚哪个因子level对应的control，避免后面解释数据遇到麻烦（你不知道到底logFC是以谁做参照的，我们希望是以control作为参照，这样logFC>1就表示实验组表达大于对照组；但是R不知道，R默认按字母顺序排列因子顺序，所以我们要对这个因子顺序set一下：
#factor levels，写在前面的level作为参照
dds$condition <- factor(dds$condition, levels = c("untreated","treated"))
#或者用relevel函数
dds$condition <- relevel(dds$condition, ref = "untreated")
#如果对应的level没有样本，可以这样把对应level删除
dds$condition <- droplevels(dds$condition)
合并技术重复的数据
DESeq2提供了一个collapseReplicates函数来把技术重复的样本数据合并到表达矩阵的一个列中；技术重复就是对同一样本测多次（multiple sequencing runs of the same library）。注意：不能对生物学重复数据用这种方法合并
继续使用来自pasilla包的示例数据
用results函数提取差异分析结果，默认的包括p-value,padj,log2foldchange，而且默认上述值是针对design里最后一个变量的，不过用户也可以自己规定针对什么变量，只要往results里添加参数name或者contrasts.
dds<-DESeq(dds)
res<-results(dds)
res# log2 fold change (MLE): condition treated vs untreated
# Wald test p-value: condition treated vs untreated
# DataFrame with 14599 rows and 6 columns
# baseMean log2FoldChange lfcSE
# <numeric> <numeric> <numeric>
# FBgn0000003 0.171568715207063 1.02601368333522 3.80551160374507
# FBgn0000008 95.1440789963134 0.00215175449141369 0.223883805414937
# FBgn0000014 1.05657219346166 -0.496735199348756 2.16026616643833
# FBgn0000015 0.846723274987709 -1.88276494877056 2.10643527029088
# FBgn0000017 4352.5928987935 -0.240025038615816 0.126024321793908
#如果要明确比较的变量，可以这样做
res<-results(dds,name="condition_treated_vs_untreated")
res<-results(dds,contrast=c("condition","treated","untreated"))
Log fold change shrinkage for visualization and ranking
表达量越低的基因，彼此之间离散程度越大，比如0和10的表达量就会显得方差相对很高；但是10010和10000也是差10但是相对表达量离散程度就很小，为了统计检验准确有效不受count数太大的影响，这里可以做shrink
resultsNames(dds)#查看参数
## [1] "Intercept" "condition_treated_vs_untreated"
#选择要shrink的参数
resLFC<-lfcShrink(dds,coef="condition_treated_vs_untreated",type="apeglm")
# using 'apeglm' for LFC shrinkage. If used in published research, please cite:
# Zhu, A., Ibrahim, J.G., Love, M.I. (2018) Heavy-tailed prior distributions for
# sequence count data: removing the noise and preserving large differences.
# Bioinformatics. https://doi.org/10.1093/bioinformatics/bty895
resLFC#可以用as.data.frame方便查看res
Using parallelization
如果实验的设计更为复杂而且又大量样本，可以考虑并行计算，使用BiocParallel包
Parallelizing DESeq, results, and lfcShrink can be easily accomplished by loading the BiocParallel package, and then setting the following arguments: parallel=TRUE and BPPARAM=MulticoreParam(4)
或者一开始就声明好4核：
library("BiocParallel")register(MulticoreParam(4))
#像上述声明好后，在function里添加parallel=TRUE选项即可
p-values and adjusted p-values
resOrdered<-res[order(res$pvalue),]#按pvalue从小到大排列
summary(res)
# out of 12359 with nonzero total read count
# adjusted p-value < 0.1
# LFC > 0 (up) : 521, 4.2%
# LFC < 0 (down) : 540, 4.4%
# outliers [1] : 1, 0.0081%
# low counts [2] : 4035, 33%
# (mean count < 7)
# [1] see 'cooksCutoff' argument of ?results
# [2] see 'independentFiltering' argument of ?results
sum(res$padj<0.1,na.rm=TRUE)
#results函数参数是很丰富的，建议?results查看，这里把默认pvalue=0.1设置成0.05
res05<-results(dds,alpha=0.05)
summary(res05)
# out of 12359 with nonzero total read count
# adjusted p-value < 0.05
# LFC > 0 (up) : 406, 3.3%
# LFC < 0 (down) : 432, 3.5%
# outliers [1] : 1, 0.0081%
# low counts [2] : 3797, 31%
# (mean count < 5)
# [1] see 'cooksCutoff' argument of ?results
# [2] see 'independentFiltering' argument of ?results
sum(res05$padj<0.05,na.rm=TRUE)
Independent hypothesis weighting
用IHW包对pvalue进行过滤，提高检验功效，通过给假设赋予权重进行多重检验
http://bioconductor.org/packages/release/bioc/html/IHW.html
library("IHW")resIHW<-results(dds,filterFun=ihw)
summary(resIHW)
# out of 12359 with nonzero total read count
# adjusted p-value < 0.1
# LFC > 0 (up) : 498, 4%
# LFC < 0 (down) : 548, 4.4%
# outliers [1] : 1, 0.0081%
# [1] see 'cooksCutoff' argument of ?results
# see metadata(res)$ihwResult on hypothesis weighting
sum(resIHW$padj<0.1,na.rm=TRUE)
metadata(resIHW)$ihwResult
# ihwResult object with 14599 hypothesis tests
# Nominal FDR control level: 0.1
# Split into 8 bins, based on an ordinal covariate
#注意，所有通过DESeq2这个包进行的计算结果都被存储在了DESeqDataSet或者DESeqResults对象中
#怎么获取这些结果以后的推文会讨论
探索、输出结果
MA-plot
每个点是一个gene，横坐标是平均的标准化counts，纵坐标是logFC
plotMA(res, ylim=c(-2,2))
#对log2 fold changes做了shrink以后
plotMA(resLFC, ylim=c(-2,2))
#交互式地判断每个点对应什么基因
idx <- identify(res$baseMean, res$log2FoldChange)
rownames(res)[idx]
如果没有shrink就是这样，低表达量的基因往往离散程度更大
MA-plot.png
做了矫正后：
MA_LFC.png
其他的shrinkage estimators
有时候对于一些数据集之前的shrink会过强，因此也可以考虑其他方法：
通过修改lfcShrink参数type来指定
The options for type are:
normal is the the original DESeq2 shrinkage estimator, an adaptive Normal distribution as prior. This is currently the default, although the default will likely change to apeglm in the October 2018 release given apeglm’s superior performance.
apeglm is the adaptive t prior shrinkage estimator from the apeglm package (Zhu, Ibrahim, and Love 2018).
ashr is the adaptive shrinkage estimator from the ashr package (Stephens 2016). Here DESeq2 uses the ashr option to fit a mixture of Normal distributions to form the prior, with method="shrinkage".
resultsNames(dds)
# because we are interested in treated vs untreated, we set 'coef=2'
resNorm <- lfcShrink(dds, coef=2, type="normal")
resAsh <- lfcShrink(dds, coef=2, type="ashr")
#比较三种shrink方法
par(mfrow=c(1,3), mar=c(4,4,2,1))#修改图的展示方式
xlim <- c(1,1e5); ylim <- c(-3,3)
plotMA(resLFC, xlim=xlim, ylim=ylim, main="apeglm")
plotMA(resNorm, xlim=xlim, ylim=ylim, main="normal")
plotMA(resAsh, xlim=xlim, ylim=ylim, main="ashr")
shrink_methods.png
如果需要校正批次效应，在design里可以声明好batch factor，或者使用其他的包
，比如sva或者RUVseq去捕捉那些潜在会产生异质性数据的无关变量
关于sva是什么，可以参考：https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.0030161
实际上是帮助捕获并模拟那些我们观测或测定不了的会对生物学结果造成干扰的变量
基因表达水平可以受很多因素影响，比如遗传、环境、人群、技术等。除了有些已知因素可以测量以外，很多因素实际上要么是未知的要么是无法测定或者过于复杂以至于不能在单一模型里很好地捕获它们。如果不能把这些因素造成的异质性考虑进去，实际上很有可能对研究结果产生较大的影响。本文介绍的SVA(‘surrogate variable analysis)是这样一种方法，它能够准确地捕获表达信息和任何建模变量间的关系，从而增强生物数据的准确性以及分析的可重复性。
Plot counts
#怎么检查某个基因在不同组里的表达情况？
plotCounts(dds, gene=which.min(res$padj), intgroup="condition")
#还可以用ggplot画
d <- plotCounts(dds, gene=which.min(res$padj), intgroup="condition",
returnData=TRUE)
library("ggplot2")
ggplot(d, aes(x=condition, y=count)) +
geom_point(position=position_jitter(w=0.1,h=0)) +
scale_y_log10(breaks=c(25,100,400))
关于Results
#关于results的详细信息：变量和检验方法的查看
mcols(res)$description
## [1] "mean of normalized counts for all samples"
## [2] "log2 fold change (MLE): condition treated vs untreated"
## [3] "standard error: condition treated vs untreated"
## [4] "Wald statistic: condition treated vs untreated"
## [5] "Wald test p-value: condition treated vs untreated"
## [6] "BH adjusted p-values"
实际上p-value以及一些值可能会输出NA，原因有:
- 如果某基因在所有样本的表达值都是0，将导致pvalue和padj都为NA
- 如果某个基因在某个样本里有异常的count离群值；异常点的检测基于Cook’s distance，不过对异常点的处理是可以自己改的
- 实际上在DESeq2和edgeR中除了library normalization之外，还会进行Independent Filtering。Independent Filtering筛去一些不太可能有生物学差异的基因，而不管它的p值是否显著，以尽量降低假阴性结果。如果某个基因因为这个（可能是low mean count)被自动过滤了，那么虽然有pvalue但padj为NA，详细可参考statQuest课程 http://www.360doc.com/content/18/1007/19/51784026_792756710.shtml
Rich visualization and reporting of results
如果想生成html报告，可以考虑ReportingTools这个包：http://bioconductor.org/packages/release/bioc/html/ReportingTools.html，示例代码在对应vignette的RNA-seq differential expression的文档中可以找到
regionReport、Glimma 、pcaExplorer、iSEE、DEvis等也有交互展示运行结果的功能
Exporting results to CSV files
write.csv(as.data.frame(resOrdered),
file="condition_treated_results.csv")
#还可以通过条件限制需要输出的结果：subset函数
resSig <- subset(resOrdered, padj < 0.1)
resSig
多因子实验设计
这里简单介绍一点，具体的formula设计在以后的推文会讨论。
实际上design的公式可以是非常多样的，还可以考虑因子间交互作用。pasilla包除了condition以外，还有测序类型可以考虑进去：
colData(dds)
# DataFrame with 7 rows and 3 columns
# condition type sizeFactor
# <factor> <factor> <numeric>
# treated1 treated single-read 1.63557509657607
# treated2 treated paired-end 0.761269768042316
# treated3 treated paired-end 0.832652635328833
# untreated1 untreated single-read 1.13826297659084
# untreated2 untreated single-read 1.79300035535039
# untreated3 untreated paired-end 0.649547030603726
# untreated4 untreated paired-end 0.751689223426488
#先拷贝一份dds用来做多因子分析
ddsMF <- dds
levels(ddsMF$type)
## [1] "paired-end" "single-read"
levels(ddsMF$type) <- sub("-.*", "", levels(ddsMF$type))
levels(ddsMF$type)
#注意感兴趣的变量要放到末尾
design(ddsMF) <- formula(~ type + condition)
ddsMF <- DESeq(ddsMF)
resMF <- results(ddsMF)
head(resMF)
实际上我们也可以获得其他变量造成的log2FC、pvalue等参数的结果（这里的其他变量只是测序类型，不是生物学上有意义的变量)
；如果是这样的design:~genotype + condition + genotype:condition，说明我们对由基因型造成的表达差异感兴趣
这里演示只考虑type作为变量的差异分析
resMFType <- results(ddsMF,
contrast=c("type", "single", "paired"))
head(resMFType)
# log2 fold change (MLE): type single vs paired
# Wald test p-value: type single vs paired
# DataFrame with 6 rows and 6 columns
# baseMean log2FoldChange lfcSE
# <numeric> <numeric> <numeric>
# FBgn0000003 0.171568715207063 -1.61158240361812 3.87108289314
# FBgn0000008 95.1440789963134 -0.262254386430198 0.22068580426262
# FBgn0000014 1.05657219346166 3.29058255215038 2.08724091994889
# FBgn0000015 0.846723274987709 -0.581642730889627 2.18165959131568
# FBgn0000017 4352.5928987935 -0.0997652738257474 0.111757030425811
# FBgn0000018 418.614930505122 0.229299212203436 0.1305752643427
# stat pvalue padj
# <numeric> <numeric> <numeric>
# FBgn0000003 -0.416313070038884 0.677180929822828 NA
# FBgn0000008 -1.18836092473855 0.234691244221223 0.543823121250909
# FBgn0000014 1.57652263363587 0.114905406827234 NA
# FBgn0000015 -0.266605630504831 0.789772819994317 NA
# FBgn0000017 -0.892697966701751 0.372018939542609 0.683007220487336
# FBgn0000018 1.7560692935044 0.0790765774697509 0.292463782417282
到这里主要的workflow就结束了，后面会介绍更高阶的操作：数据转换、可视化、个性化分析等
major reference
5人点赞
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
