limma，edgeR，DESeq2三大包基本是做转录组差异分析的金标准，大多数转录组的文章都是用这三个R包进行差异分析。
edgeR差异分析速度快，得到的基因数目比较多，假阳性高（实际不差异结果差异）。
DESeq2差异分析速度慢，得到的基因数目比较少，假阴性高（实际差异结果不差异）。
需要注意的是制作分组信息的因子向量时，因子水平的前后顺序，在R的很多模型中，默认将因子向量的第一个水平看作对照组。
logFC值:
limma使用对数化的表达矩阵作为输入，所以将零假设指定在平均log值（对数几何平均数）。
edgeR和DESeq2使用原始的count矩阵作为输入，所以将原假设指定在count的平均值上（算术平均数）。
logFC一般是在1.2到2之间，筛选到差异基因的数目在500-1000左右为宜•
根据logFC的统计指标确定阈值的方法是计算logFC绝对值的平均数与2倍标准差之和。（正态分布）
DEG表示差异表达基因，differential gene express
0. 准备
判断式安装R包：如果该R包存在，可以顺带加载该R包，不需要再次library。
if(!require(stringr))install.packages("stringr") # 常规R包
if(!require(ggplotify))install.packages("ggplotify") # 热图转换为ggplot的类型
if(!require(patchwork))install.packages("patchwork") # 拼图R包
if(!require(cowplot))install.packages("cowplot") # 拼图R包
if(!require(DESeq2))BiocManager::install("DESeq2") # 差异分析R包
if(!require(edgeR))BiocManager::install("edgeR") # 差异分析R包
if(!require(limma))BiocManager::install("limma") # 差异分析R包
rm(list = ls())
load("TCGA-CHOLgdc.Rdata")
table(group_list)
差异分析输入数据：表达矩阵exp，分组信息group_list。
1. DESeq2
差异分析
library(DESeq2)
library(stringr)
colData <- data.frame(row.names = colnames(exp),
condition = group_list) # 列出每个样品是肿瘤样品还是正常样品
dds <- DESeqDataSetFromMatrix(countData = exp,
colData = colData,
design = ~ condition) %>%
DESeq() # 将数据框转为DESeq2的数据集类型，然后用DESeq函数做差异分析
提取结果，两两比较
res <- results(dds,
contrast = c("condition",rev(levels(group_list)))) # 按设置的比较水平提取dds的数据（从DESeq分析中提取结果表，给出样本的基本均值，logFC及其标准误，检验统计量，p值和矫正后的p值）。rev表示调换顺序，levels(group_list)是因子水平
DEG <- res[order(res$pvalue),] %>%
as.data.frame() # 将res按照P值从小到大排序，并转换回数据框。order函数获取向量每个元素从小到大排序后的第几位
head(DEG)
添加change列标记基因上调下调
logFC_cutoff <- with(DEG,mean(abs(log2FoldChange)) + 2*sd(abs(log2FoldChange))) # 设置logFC的阈值，可以计算得出，也可以设置固定值。abs函数计算log2FoldChange的绝对值。均数+2倍标准差可以包含log2FoldChange的95%置信区间 k1 = (DEG$pvalue < 0.05)&(DEG$log2FoldChange < -logFC_cutoff) # 判断下调基因 k2 = (DEG$pvalue < 0.05)&(DEG$log2FoldChange > logFC_cutoff) # 判断上调基因 DEG$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT")) # 标记上、下调基因 table(DEG$change) # 查看上下调基因的数量 head(DEG) DESeq2_DEG <- DEG # 保存DESeq2差异分析的结果
2. edgeR
差异分析
library(edgeR) library(stringr) # 准备count表达矩阵和分组信息，构建edgeR的DGEList对象 dge <- DGEList(counts = exp, group = group_list) # 过滤low count数据，保留在至少在两个样本中cpm值（Counts per million）大约1的基因。 # keep <- rowSums(cpm(dge) > 1 ) >= 2 # table(keep) # dge <- dge[keep, , keep.lib.sizes = FALSE] # 由于在上一步整理数据的过程中已经过滤了一遍数据了，所以这里不再过滤。 dge$samples$lib.size <- colSums(dge$counts) # 标准化基因表达分布，以TMM标准化为例。 dge <- calcNormFactors(dge) # calcNormFactors函数并不会标准化数据，只是计算标准化因子 # 假设数据符合正态分布，构建线性模型。 design <- model.matrix( ~ 0 + group_list) # 0代表x线性模型的截距为0，需要修改 rownames(design) <- colnames(dge) colnames(design) <- levels(group_list) # 计算线性模型的参数 dge <- estimateGLMCommonDisp(dge,design) %>% estimateGLMTrendedDisp(design) %>% estimateGLMTagwiseDisp(design) fit <- glmFit(dge, design) %>% # 拟合线性模型 glmLRT(contrast = c(-1,1)) # 进行差异分析。（1,-1）意味着前比后
提取结果
DEG = topTags(fit, n=nrow(exp)) %>% as.data.frame() # 提取差异分析结果 # 添加change列标记基因上调下调 logFC_cutoff <- with(DEG,mean(abs(logFC)) + 2*sd(abs(logFC))) # 设置logFC的阈值 # logFC_cutoff <- 2 k1 = (DEG$PValue < 0.05)&(DEG$logFC < -logFC_cutoff) # 判断下调基因 k2 = (DEG$PValue < 0.05)&(DEG$logFC > logFC_cutoff) # 判断上调基因 DEG$change = ifelse(k1, "DOWN", ifelse(k2, "UP", "NOT")) # 标记上、下调基因 head(DEG) table(DEG$change) edgeR_DEG <- DEG
3. limma
差异分析
library(limma) library(stringr) # 创建设计矩阵和对比。假设数据符合正态分布，构建线性模型 design <- model.matrix( ~ 0 + group_list) colnames(design) = levels(group_list) rownames(design) = colnames(exp) dge <- DGEList(counts=exp) %>% # 将count表达矩阵转换为edgeR的DGEList对象 calcNormFactors() # 进行标准化和表达量计算。calcNormFactors函数并不会标准化数据，只是计算标准化因子。 # 标准化、线性模型拟合 fit <- voom(dge, design, normalize = "quantile") %>% # voom是limma中的一种标准化方法 lmFit(design) # 线性模型拟合，出图 # 设置需要进行对比的分组 comp = paste(rev(levels(group_list)), collapse = "-") cont.matrix <- makeContrasts(contrasts = comp, levels = design) # 比较每个基因 fit2 <- contrasts.fit(fit, cont.matrix) %>% # 根据对比模型进行差值计算 eBayes() # 贝叶斯检验 # 使用plotSA绘制了log2残差标准差与log-CPM均值的关系。平均log2残差标准差由水平蓝线标出 plotSA(fit2, main = "Final model: Mean-variance trend")
提取结果
DEG = topTable(fit2, coef = comp, n=Inf) %>% # 提取差异分析结果 na.omit() # 去除差异分析结果中包含NA值的行 length(which(DEG$adj.P.Val < 0.05)) # p值<0.05的基因有多少个 # 添加change列标记基因上调下调 logFC_cutoff <- with(DEG,mean(abs(logFC)) + 2*sd(abs(logFC))) # 设置logFC的阈值 k1 = (DEG$P.Value < 0.05)&(DEG$logFC < -logFC_cutoff) # 判断下调基因 k2 = (DEG$P.Value < 0.05)&(DEG$logFC > logFC_cutoff) # 判断上调基因 DEG$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT")) # 标记上、下调基因 table(DEG$change) head(DEG) limma_voom_DEG <- DEG
4. 数据检查
整合三大R包差异分析的数据：
tj = data.frame(deseq2 = as.integer(table(DESeq2_DEG$change)),
edgeR = as.integer(table(edgeR_DEG$change)),
limma_voom = as.integer(table(limma_voom_DEG$change)),
row.names = c("down","not","up")
);tj
验证数据：
验证logFC是否由tumor/normal计算而来：
boxplot(as.integer(exp[rownames(DESeq2_DEG)[1],]) ~ group_list) # 取exp的第一行基因表达量，以group_list为分组依据，做箱线图。修改“DESeq2_DEG”即可分别分析三种方法。 DESeq2_DEG$log2FoldChange[1] # 得到DESeq2_DEG的第一个logFC值。如果箱线图中normal大tumor小，则logFC应为负值，反之为正值。 # 如果检查发现logFC的正负值反了，则修改为负值即可：DESeq2_DEG$log2FoldChange = -DESeq2_DEG$log2FoldChange
保存数据
save(DESeq2_DEG, edgeR_DEG, limma_voom_DEG, group_list, tj, file = paste0(cancer_type,"DEG.Rdata"))
5.三大R包差异分析可视化（热图、火山图）
rm(list = ls())
setwd("D:/WorkingDirectory/01.mRNA_Expdata")
load("TCGA-CHOL_01.mRNA_Expdata_count.Rdata") # 加载表达数据整理结果
setwd("D:/WorkingDirectory/03.mRNA_Diff")
load("TCGA-CHOL_DEG.Rdata") # 加载差异表达基因结果
if(!require(tinyarray))devtools::install_local("tinyarray-master.zip", upgrade = F)
library(ggplot2)
library(tinyarray)
exp[1:4,1:4] # 原始count表达矩阵数据
dat = log2(exp+1)
pca.plot = draw_pca(dat, group_list); pca.plot # 作PCA图
save(pca.plot, file = paste0(cancer_type, "_PCA.Plot.Rdata")) # 保存PCA图数据以便后续拼图
# 提取差异基因
cg1 = rownames(DESeq2_DEG)[DESeq2_DEG$change != "NOT"]
cg2 = rownames(edgeR_DEG)[edgeR_DEG$change != "NOT"]
cg3 = rownames(limma_voom_DEG)[limma_voom_DEG$change != "NOT"]
# 画热图
h1 = draw_heatmap(dat[cg1,], group_list, scale_before = T)
h2 = draw_heatmap(dat[cg2,], group_list, scale_before = T)
h3 = draw_heatmap(dat[cg3,], group_list, scale_before = T)
# 构建“均数+2倍标准差”函数
m2d = function(x){
mean(abs(x))+2*sd(abs(x))
}
# 画火山图
v1 = draw_volcano(DESeq2_DEG, pkg = 1, logFC_cutoff = m2d(DESeq2_DEG$log2FoldChange))
v2 = draw_volcano(edgeR_DEG, pkg = 2, logFC_cutoff = m2d(edgeR_DEG$logFC))
v3 = draw_volcano(limma_voom_DEG, pkg = 3, logFC_cutoff = m2d(limma_voom_DEG$logFC))
# 拼图
library(patchwork)
(h1 + h2 + h3) / (v1 + v2 + v3) + plot_layout(guides = "collect") # & theme(legend.position = "none")
# +表示横向，/表示纵向，guides表示将图例统一放在右侧。
# &theme(legend.position = "none")表示不显示图例。
ggsave(paste0(cancer_type, "_Heat.Vo.png"),width = 15, height = 10)
6.三大R包差异基因对比（热图、韦恩图、PCA图）
rm(list = ls())
setwd("D:/WorkingDirectory/01.mRNA_Expdata")
load("TCGA-CHOL_01.mRNA_Expdata_count.Rdata")
setwd("D:/WorkingDirectory/03.mRNA_Diff")
load("TCGA-CHOL_DEG.Rdata")
load("TCGA-CHOL_PCA.Plot.Rdata")
if(!require(tinyarray))devtools::install_local("tinyarray-master.zip", upgrade = F)
library(ggplot2)
library(tinyarray)
# 构建提取上调基因的函数：
UP = function(df){
rownames(df)[df$change == "UP"]
}
# 构建提取下调基因的函数：
DOWN = function(df){
rownames(df)[df$change == "DOWN"]
}
# 取三大R包上调基因交集：
up = intersect(intersect(UP(DESeq2_DEG),
UP(edgeR_DEG)),
UP(limma_voom_DEG)) # intersect函数表示取两个向量的交集。
# 取三大R包下调基因交集：
down = intersect(intersect(DOWN(DESeq2_DEG),
DOWN(edgeR_DEG)),
DOWN(limma_voom_DEG))
# 取三大R包差异基因交集画热图：
dat = log2(exp+1)
hp = draw_heatmap(dat[c(up, down), ], group_list, scale_before = T)
#上调、下调基因分别画维恩图：
up_genes = list(Deseq2 = UP(DESeq2_DEG),
edgeR = UP(edgeR_DEG),
limma = UP(limma_voom_DEG))
down_genes = list(Deseq2 = DOWN(DESeq2_DEG),
edgeR = DOWN(edgeR_DEG),
limma = DOWN(limma_voom_DEG))
up.plot <- draw_venn(up_genes, "UPgene")
down.plot <- draw_venn(down_genes, "DOWNgene")
#维恩图拼图
library(patchwork)
pca.plot + hp + up.plot + down.plot + plot_layout(guides = "collect")
ggsave(paste0(cancer_type, "_Heat.Ve.PCA.png"), width = 15, height = 10)
