DESeq2详解系列(3)

接着前两篇文章：
DESeq2详解系列(1)：https://www.jianshu.com/p/88511070e2dd
DESeq2详解系列(2)：https://www.jianshu.com/p/ffc16dacc711

探索、输出结果

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

http://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#exploring-and-exporting-results

最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
平台声明：文章内容（如有图片或视频亦包括在内）由作者上传并发布，文章内容仅代表作者本人观点，简书系信息发布平台，仅提供信息存储服务。

推荐阅读更多精彩内容

  • 2019-12-26 复习表观调控相关内容（图1-3）
    复习资料： 表观调控13张图之一证明基因干扰有效性 表观调控13张图之二相关性热图看不同样本相关性 表观调控13张...
    程凉皮儿阅读 4,077评论 0赞 4
  • R 记录
    library(pheatmap) data<- read.table("new 1.txt",header = ...
    Weiyx阅读 3,848评论 0赞 0
  • DESeq2详解系列(2)
    差异表达分析 接着上篇文章：https://www.jianshu.com/p/88511070e2dd继续使用来...
    Shaoqian_Ma阅读 5,426评论 2赞 5
  • DESeq2的建模原理及简单用法
    写在前面的废话 不研究不知道，一研究吓一跳，原来DESeq2这么复杂，这10000多的引用量真不是吹的…… 废话超...
    鹿无为阅读 26,977评论 11赞 62
  • 我喜欢你
    我喜欢 喜欢春天的阳光 喜欢夏天的花 喜欢秋天的落叶 喜欢冬天的雪 和 喜欢爱可爱的你 暖暖的样子 贴心的偎依 完美四季
    野派阅读 1,508评论 0赞 0