差异分析
• limma
• edgeR
• DESeq2
title: "三大R包差异分析"
output: html_document
editor_options:
chunk_output_type: console
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
knitr::opts_chunk$set(fig.width = 8,fig.height = 6,collapse = TRUE)
knitr::opts_chunk$set(message = FALSE)
三大R包差异分析
if(!require(stringr))install.packages('stringr')
if(!require(ggplotify))install.packages("ggplotify")
if(!require(patchwork))install.packages("patchwork")
if(!require(cowplot))install.packages("cowplot")
if(!require(DESeq2))BiocManager::install('DESeq2')
if(!require(edgeR))BiocManager::install('edgeR')
if(!require(limma))BiocManager::install('limma')
rm(list = ls())
load("TCGA-CHOL_gdc.Rdata")
table(Group)
deseq2----
library(DESeq2)
colData <- data.frame(row.names =colnames(exp),
condition=Group)
if(!file.exists(paste0(cancer_type,"_dd.Rdata"))){ #第一次运行需要将数据转换为deseq2所需要的类型，并保存_dd.Rdata格式，再次运行将跳过此步骤
dds <- DESeqDataSetFromMatrix(
countData = exp,
colData = colData,
design = ~ condition)
dds <- DESeq(dds)
save(dds,file = paste0(cancer_type,"_dd.Rdata"))
}
load(file = paste0(cancer_type,"_dd.Rdata"))
# 两两比较
res <- results(dds, contrast = c("condition",rev(levels(Group)))) #三包都需要提供分组信息，默认因子normol在前，tumor在后，加上rev函数逆转顺序
resOrdered <- res[order(res$pvalue),] # 按照P值排序
DEG <- as.data.frame(resOrdered) #转换为数据框
DEG = na.omit(DEG) #当取log及其他操作时，数据太小无法识别则返回NA，去除NA操作
head(DEG)
#添加change列标记基因上调下调
logFC_cutoff <- with(DEG,mean(abs(log2FoldChange)) + 2*sd(abs(log2FoldChange)) )#计算logFC_cutoff 95%置信区间； #with的作用是替代$将列变为局部变量进行操作
#logFC_cutoff <- 2/以2作为判断标准和前面置信区间计算出来的值都可以
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
edgeR----
library(edgeR)
dge <- DGEList(counts=exp,group=Group) #提供表达数据及分组信息
dge$samples$lib.size <- colSums(dge$counts)
dge <- calcNormFactors(dge)
design <- model.matrix(~0+Group)
rownames(design)<-colnames(dge)
colnames(design)<-levels(Group)
dge <- estimateGLMCommonDisp(dge, design)
dge <- estimateGLMTrendedDisp(dge, design)
dge <- estimateGLMTagwiseDisp(dge, design)
fit <- glmFit(dge, design)
fit2 <- glmLRT(fit, contrast=c(-1,1)) #之前内容都可看为一个函数，fit2是我们想要的东西
DEG=topTags(fit2, n=nrow(exp)) #赋值给DEG
DEG=as.data.frame(DEG)
logFC_cutoff <- with(DEG,mean(abs(logFC)) + 2*sd(abs(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
limma----
library(limma)
design <- model.matrix(~0+Group) #0是为了给后面的比较矩阵cont.matrix留下位置
colnames(design)=levels(Group)
rownames(design)=colnames(exp)
dge <- DGEList(counts=exp)
dge <- calcNormFactors(dge)
v <- voom(dge,design, normalize="quantile")
fit <- lmFit(v, design)
constrasts = paste(rev(levels(Group)),collapse = "-")
cont.matrix <- makeContrasts(contrasts=constrasts,levels = design)
fit2=contrasts.fit(fit,cont.matrix)
fit2=eBayes(fit2)
DEG = topTable(fit2, coef=constrasts, n=Inf)
DEG = na.omit(DEG)
logFC_cutoff <- with(DEG,mean(abs(logFC)) + 2*sd(abs(logFC)) )
#logFC_cutoff <- 2
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)
检验我们的上下调基因是否比反
exp["ENSG00000042781.11",] #看到第一行基因不同样本表达量情况
as.numeric(exp["ENSG00000042781.11",])
boxplot(as.numeric(exp["ENSG00000042781.11",])~Group) #显示下调 #根据基因表达情况画箱线图，对比基因上调下调情况
DEG$logFC[1] #结果返回值也是负数，表示下调，操作正确
不同方式算出来基因表达数量
limma_voom_DEG <- DEG
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 注意稳定的基因占大多数才是正常情况
save(DESeq2_DEG,edgeR_DEG,limma_voom_DEG,Group,tj,file = paste0(cancer_type,"_DEG.Rdata"))
image.png
具体代码见下图
image.png
