差异分析流程（四）DESeq2

参考链接：https://www.plob.org/article/9971.html

#加载包和数据
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
if (!require("DESeq2", quietly = TRUE))
  BiocManager::install("DESeq2")
options(stringsAsFactors = F)
load("expr.RData")

expr.RData储存了经过处理后的表达矩阵和分组信息，详见差异分析流程（一）数据预处理

1.构建dds对象

colData <- data.frame(row.names=colnames(expr), group_list=group_list)

#需要保证输入的表达矩阵数值为整数
dds <- DESeqDataSetFromMatrix(countData = round(expr,0),colData = colData,design = ~ group_list)

2.标准化

量化因子 (size factor, SF)： 首先计算每个基因在所有样品中表达的几何平均值。每个细胞的量化因子(size factor)是所有基因与其在所有样品中的表达值的几何平均值的比值的中位数。由于几何平均值的使用，只有在所有样品中表达都不为0的基因才能用来计算。这一方法又被称为 RLE (relative log expression)

参考：DESeq2的标准化方法

• 估计SF和离散度及负二项式GLM拟合和Wald统计
  这一步其实已经出来差异结果了，如果不需要可视化数据分布的话可以略过后面的步骤

dds1 <- DESeq(dds)

• log2转换

rld <- rlogTransformation(dds1)
expr_new=assay(rld)
head(expr_new[,1:5])

• 可视化标准化结果

n.sample=ncol(expr)
par(cex = 0.9)
cols <- rainbow(n.sample*1.2)
par(mfrow=c(1,2))
boxplot(expr, col = cols,main="expression value",las=2)
boxplot(expr_new, col = cols,main="expression value",las=2)

3.提取差异分析结果

res <-  results(dds1, contrast=c("group_list","untrt","trt")) 
res <- as.data.frame(res)

res <- na.omit(res)
padj = 0.05
foldChange= 1
diff_signif = res[(res$padj< padj & abs(res$log2FoldChange)>foldChange),]                    
diff_signif = diff_signif[order(diff_signif$log2FoldChange),]
save(diff_signif, file = 'limma_diff.Rdata')