DeSeq2：基于负二项分布的差异表达分析
描述：估计来自高通量测序分析的计数（count）数据中的方差-均值相关性，并基于使用负二项分布的模型来检验差异表达
说明：DeSeq2包用于高维度count数据的标准化，可视化以及差异分析。通过经验贝叶斯算法估计logFC以及离散度先验，计算这些量的后验估计。（先验与后验从原因到结果的论证称为“先验的”，而从结果到原因的论证称为“后验的”。）
主要函数：
- DESeqDataSet-构建数据集，请参阅tximeta和tximport包以准备输入
- DESeq-执行差异分析
- Results-构建结果表
- lfcShrink-使用apeglm和ASHR包估计缩小的LFC(后验估计)
- vst-应用方差稳定化变换，例如用于PCA或样本聚类
-
Plots，例如：plotPCA、plotMA、plotCounts
image.png
其他函数
- coef(object, SE = FALSE, ...)
- object，由DeSeq返回的一个DESeqDataSet。
- SE（standard error），除了系数意外是否给出标准误。
提取模型系数矩阵，通常使用results函数一次性查看结果表，该函数仅用于想要观察立即所有模型系数。
suppressMessages(library(DESeq2)) #加载包
dds <- makeExampleDESeqDataSet(m=4) #构建数据集
dds <- DESeq(dds) #差异分析
coef(dds)[1,] #查看系数
coef(dds, SE=TRUE)[1,] ##查看系数标准误
- collapseReplicates函数
- 通过在分组因子GROUP BY的级别内求和来折叠对象中的列。此函数的目的是汇总技术复制的读取计数，以便为每个样本创建一个具有单列读取计数的对象。注：技术复制指的是同一文库的多次测序运行，而生物复制是指从不同的生物单位制备多个文库。可以选择使用分组系数的级别重命名返回对象的列。
dds <- makeExampleDESeqDataSet(m=12) #1000个基因（行），12个样本（列）
# make data with two technical replicates for three samples
dds$sample <- factor(sample(paste0("sample",rep(1:9, c(2,1,1,2,1,1,2,1,1))))) #设定分组因子（数目与数据集列一致）
dds$run <- paste0("run",1:12)
ddsColl <- collapseReplicates(dds, dds$sample, dds$run) #按照因子等级求和将数据折叠
# examine the colData and column names of the collapsed data
colData(ddsColl)
colnames(ddsColl)
###对该函数结果验证
matchFirstLevel <- dds$sample == levels(dds$sample)[1] #返回判断sample1所得到的逻辑值
#sample1对应dds第2，6列，对两列数据每行基因求和后，将其与折叠后的数据ddscoll第一列sample1对比。
stopifnot(all(rowSums(counts(dds[,matchFirstLevel])) == counts(ddsColl[,1])))
- counts函数
- 输入： DESeqDataSet
- 输出：行为基因，列为样本的矩阵
dds <- makeExampleDESeqDataSet(m=4) #构建示例数据集
class(counts(dds)) #查看数据类型，counts(dds)输出为表达矩阵
head(counts(dds)) #查看样本基因count矩阵前6行
dds<- estimateSizeFactors(dds)
head(counts(dds, normalized=TRUE))
normalized 标准化这里大小因子或归一化因子是什么？回答：使用标准化因子(scale factor)来解释库深度的差异，简单说就是校正测序深度。
mark一个原理说明：DESeq2对原始reads进行建模，使用标准化因子(scale factor)来解释库深度的差异。然后，DESeq2估计基因的离散度，并缩小这些估计值以生成更准确的离散度估计，从而对reads count进行建模。最后，DESeq2拟合负二项分布的模型，并使用Wald检验或似然比检验进行假设检验。
DESeq2的建模原理及简单用法 - 简书 (jianshu.com)
- DESeq函数
- 基于负二项分布的差异表达分析
- 函数默认分析步骤：1.大小因子的估计：estimateSizeFactors；2. 离散程度估计：estimateDispersions；3. 负二项GLM拟合与Wald统计量：nbinomWaldTest
- DESeq函数返回结果表（log2FC以及pvalues）可以用results函数查看；可以使用lfcShrink函数生成收缩的LFC。
DESeq(
object,
test = c("Wald", "LRT"),
fitType = c("parametric", "local", "mean", "glmGamPoi"),
sfType = c("ratio", "poscounts", "iterate"),
betaPrior,
full = design(object),
reduced,
quiet = FALSE,
minReplicatesForReplace = 7,
modelMatrixType,
useT = FALSE,
minmu = if (fitType == "glmGamPoi") 1e-06 else 0.5,
parallel = FALSE,
BPPARAM = bpparam()
) -
差异分析使用广义线性模型：
image.png
看例子
#
# 来自转录组测序的count数据
cnts <- matrix(rnbinom(n=1000, mu=100, size=1/0.5), ncol=10)
cond <- factor(rep(1:2, each=5))
# 构建DESeqDataSet数据集
# 输入表达矩阵cnts，样本分组因子矩阵DataFrame(cond)，以及模型矩阵公式~ cond
# 输出DESeqDataSet数据集
dds <- DESeqDataSetFromMatrix(cnts, DataFrame(cond), ~ cond)
# 标准差异分析过程
dds <- DESeq(dds) #差异分析
res <- results(dds) #查看差异分析结果并赋值给res
# baseMean log2FoldChange lfcSE stat pvalue padj
# moderated log2 fold changes
# baseMean表示所有样本经过归一化系数矫正的read counts（counts/sizeFactor）的均值。log2Foldchange表示该基因的表达发生了多大的变化，是估计的效应大小effect size。对差异表达的倍数取以2为底的对数，变化倍数=2^log2Foldchange。（log2FC反映的是不同分组间表达量的差异，这个差异由两部分构成，一种是样本间本身的差异，比如生物学重复样本间基因的表达量就有一定程度的差异，另外一部分就是我们真正感兴趣的，由于分组不同或者实验条件不同造成的差异。用归一化之后的数值直接计算出的log2FC包含了以上两种差异，而我们真正关注的只有分组不同造成的差异，DESeq2在差异分析的过程中已经考虑到了样本本身的差异，其最终提供的log2FC只包含了分组间的差异，所以会与手动计算的不同）。
#
# lfcSE(logfoldchange Standard Error)是对于log2Foldchange估计的标准误差估计，效应大小估计有不确定性。
#
# stat是Wald统计量，它是由log2Foldchange除以标准差所得。
# pvalue和padj分别代表原始的p值以及经过校正后的p值。
resultsNames(dds)
resLFC <- lfcShrink(dds, coef=2, type="apeglm") #logFC缩放
BiocManager::install("apeglm")
# an alternate analysis: likelihood ratio test
ddsLRT <- DESeq(dds, test="LRT", reduced= ~ 1)
resLRT <- results(ddsLRT)
网友评论