接着前两篇文章：
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就是这样，低表达量的基因往往离散程度更大
做了矫正后：
其他的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")
如果需要校正批次效应，在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就结束了，后面会介绍更高阶的操作：数据转换、可视化、个性化分析等
