最近好几个小伙伴都在用我们之前的一个bulk RNA多组差异分析函数(重启之普通R转录组分析(3)：写一个通用的Deseq2多组差异分析函数)，有些小问题，或者说一些功能不完善。所以我们这次进行了升级和优化。首先将bulk差异分析3大包DEseq2 edgeR limma都纳入进来了。第二，让您的分析更加简单，不用担心复杂的代码和分组设置(很多时候还搞错谁vs谁)。第三, 不论你的两组样本重复数是否相同，都可以进行(无重复样本分析没有，主要是我认为你算法再优化，也清除不了个体差异随机性)。第四，多组样本的分析！轻松化解！ 相信随着目前bulk的普及和便宜化，我们这个函数还是很有用的！
参考：
https://www.bilibili.com/video/BV1n2421A7oJ/?spm_id_from=333.999.0.0&vd_source=05b5479545ba945a8f5d7b2e7160ea34
首先我们看一下函数的参数，也是很简洁：
然后我们测试一下：首先是两组等重复样本，很简单！
setwd('D:\KS项目\公众号文章\Bulk多组差异基因分析函数')
#一、两组分析（样本数相等）-------------------------------------------------------------------
df1 <- read.csv("Two_group.csv", header = T, row.names = 1)
colnames(df1)
# [1] "Cancer1" "Cancer2" "Cancer3" "Health1" "Health2" "Health3"
meta1 <- data.frame(Cancer=c("Cancer1" ,"Cancer2" ,"Cancer3"),
Health=c("Health1", "Health2", "Health3"))
deg1_Deseq2 <- KS_bulkRNA_MultiGroup_DEGs(exprSet = df1,
meta = meta1,
methods = "DESeq2",
test = "Cancer",
control = "Health",
repNum1 = 3,
repNum2 = 3)
deg1_edgeR <- KS_bulkRNA_MultiGroup_DEGs(exprSet = df1,
meta = meta1,
methods = "edgeR",
test = "Cancer",
control = "Health",
repNum1 = 3,
repNum2 = 3)
deg1_limma <- KS_bulkRNA_MultiGroup_DEGs(exprSet = df1,
meta = meta1,
methods = "limma",
test = "Cancer",
control = "Health",
repNum1 = 3,
repNum2 = 3)
重复数不等的两组样本分析也是如此，只要设置好test和control即可：
#二、两组分析（样本数不同）-------------------------------------------------------------------
df4 <- read.csv("two_diff_group.csv", header = T, row.names = 1)
colnames(df4)
# [1] "HC_1" "HC_2" "HC_3" "HC_4" "PPC_1" "PPC_2" "PPC_3" "PPC_4" "PPC_5" "PPC_6"
meta4 <- data.frame(HC=c("HC_1","HC_2","HC_3","HC_4","",""),
PPC=c("PPC_1","PPC_2","PPC_3","PPC_4","PPC_5","PPC_6"))
deg4_Deseq2 <- KS_bulkRNA_MultiGroup_DEGs(exprSet = df4,
meta = meta4,
methods = "DESeq2",
test = "PPC",
control = "HC",
repNum1 = 6,
repNum2 = 4)
deg4_edgeR <- KS_bulkRNA_MultiGroup_DEGs(exprSet = df4,
meta = meta4,
methods = "edgeR",
test = "PPC",
control = "HC",
repNum1 = 6,
repNum2 = 4)
deg4_limma <- KS_bulkRNA_MultiGroup_DEGs(exprSet = df4,
meta = meta4,
methods = "limma",
test = "PPC",
control = "HC",
repNum1 = 6,
repNum2 = 4)
我们可以对比下三种方法差异基因结果，发现肯定不是完全一样，那是因为算法不一样，不用纠结。在我们这个数据里面可以看出，DEseq2和edgeR的结果重合还是挺多的。
#韦恩图
deg4_Deseq2_sig <- deg4_Deseq2[which(abs(deg4_Deseq2$log2FoldChange)>0 & deg4_Deseq2$pvalue <=0.05),]
deg4_edgeR_sig <- deg4_edgeR[which(abs(deg4_edgeR$logFC)>0 & deg4_edgeR$PValue <=0.05),]
deg4_limma_sig <- deg4_limma[which(abs(deg4_limma$logFC)>0 & deg4_limma$P.Value <=0.05),]
library(ggvenn)
library(tidyverse)
Venn_list <- list(deg4_Deseq2_sig=rownames(deg4_Deseq2_sig),
deg4_edgeR_sig=rownames(deg4_edgeR_sig),
deg4_limma_sig=rownames(deg4_limma_sig))
#使用list_to_data_frame将list转化为data.frame
data_veen = list_to_data_frame(Venn_list)
ggvenn(Venn_list,
show_percentage = T,
show_elements = F,
text_size=3,
digits = 1,
set_name_size=4,
stroke_color = "grey30",
fill_color = c("#FF8C00","#4DAF4A","#B64E89"),
set_name_color = c("#FF8C00","#4DAF4A","#B64E89"))
多组样本的分析就更简单了：大多设置都设置好了。这里需要强调一点，那就是对于多组样本，在进行meta设置的时候，test组靠前，control靠后！
#四、多组test------------------------------------------------------------------
df3 <- read.csv("count_gene.csv", header = T, row.names = 1, check.names = F)
colnames(df3)
# [1] "Fbrain 1" "Fbrain 2" "Fbrain 3" "Fbrain 4" "Fhom 1" "Fhom 2" "Fhom 3"
# [8] "Mbrain 1" "Mbrain 2" "Mbrain 3" "Mbrain 4" "Mbrain 5" "Mhom 1" "Mhom 2"
# [15] "Mhom 3" "Mhom 4"
meta3 <- data.frame(Mhom = c("Mhom 1","Mhom 2", "Mhom 3","Mhom 4",""),
Fhom = c("Fhom 1","Fhom 2","Fhom 3","",""),
Mbrain=c("Mbrain 1","Mbrain 2","Mbrain 3","Mbrain 4","Mbrain 5"),
Fbrain=c("Fbrain 1","Fbrain 2","Fbrain 3","Fbrain 4",""))
deg3_Deseq2 <- KS_bulkRNA_MultiGroup_DEGs(exprSet = df3,
meta = meta3,
methods = "DESeq2",
separator=" ")
deg3_edgeR <- KS_bulkRNA_MultiGroup_DEGs(exprSet = df3,
meta = meta3,
methods = "edgeR",
separator=" ")
deg3_limma <- KS_bulkRNA_MultiGroup_DEGs(exprSet = df3,
meta = meta3,
methods = "limma",
separator=" ")
好了，这就是所有的内容了，希望对你有帮助！