差异表达分析
接着上篇文章：https://www.jianshu.com/p/88511070e2dd
继续使用来自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对象中
#怎么获取这些结果以后的推文会讨论
