刘小泽写于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 -
利用
releveldds$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
欢迎关注我们的公众号～_～
我们是两个农转生信的小硕，打造生信星球，想让它成为一个不拽术语、通俗易懂的生信知识平台。需要帮助或提出意见请后台留言或发送邮件到jieandze1314@gmail.com
Welcome to our bioinfoplanet!
网友评论