用R自己写一个多组DEGs计算的代码

函数主体

DEGs_find=function(data, design){
    ####安装需要的包
    if (!requireNamespace("BiocManager", quietly = TRUE))
        install.packages("BiocManager")

    if(!require("limma")) BiocManager::install("limma")
    if(!require("data.table")) install.packages("data.table")
    if(!require("dplyr")) install.packages("dplyr")
    if(!require("splitstackshape")) install.packages("splitstackshape")

    library(limma);library(data.table);library(dplyr);library(splitstackshape)

    #####开始计算
    de=design
    design <- model.matrix(~0+factor(de))
    #colnames(design)=paste("ii",levels(factor(de)),sep="")
    colnames(design)=levels(factor(de))
    design<<-design
    
    if(length(table(de))==2){
        contrast.matrix<-makeContrasts(
            paste(colnames(design)[1], 
            colnames(design)[2], sep="-"), 
            levels = design)
    }
    
    if(length(table(de)) > 2){
        
        makeContrasts_generate=function(de){
            lis=matrix(1:length(levels(factor(de)))^2, ncol=length(levels(factor(de))))
            colnames(lis)=rownames(lis)=colnames(design)

            lis[upper.tri(lis), drop=FALSE]

                locate_group = function(x){
                    lis_all=(1:dim(lis)[1])*dim(lis)[1]
                    col=length(which(lis_all<x))+1
                    row=x-length(which(lis_all<x))*dim(lis)[1]
                    return(paste(rownames(lis)[row], colnames(lis)[col], sep="-"))
                    
                }
            return(sapply(lis[upper.tri(lis), drop=FALSE], locate_group))
        }
            
        lis=matrix(1:length(levels(factor(de)))^2, ncol=length(levels(factor(de))))
        colnames(lis)=rownames(lis)=colnames(design)
        shan2=matrix(ncol=length(lis[upper.tri(lis), drop=FALSE]), nrow=length(colnames(design)))
        colnames(shan2)=makeContrasts_generate(de)
        rownames(shan2)=colnames(design)

        shan=colnames(shan2)
        shan1=rownames(shan2)
        library(data.table);library(splitstackshape)
        shan_1=as.data.table(shan);colnames(shan_1)="V1"
        shan_1=cSplit(shan_1, "V1", "-")

        ####定位函数
        shan2_locate=function(x, y){
        grep(x, as.character(as.matrix(y)))
        }

        ####横坐标定位_正定位
        out_x1=alist()
        for (i in 1:length(shan1)){
            out_x1[[i]]=shan2_locate(x=shan1[i], y=shan_1[,1])
        }

        for (i in 1:length(out_x1)){
            if (length(out_x1[[i]])==0){"skip"}
            if (length(out_x1[[i]])>0){
                for (j in 1:length(out_x1[[i]])){
                    shan2[i,out_x1[[i]][j]]=1
                }
            }
        }
        ####横坐标定位_负定位
        out_x2=alist()
        for (i in 1:length(shan1)){
            out_x2[[i]]=shan2_locate(x=shan1[i], y=shan_1[,2])
        }

        for (i in 1:length(out_x2)){
            if (length(out_x2[[i]])==0){"skip"}
            if (length(out_x2[[i]])>0){
                for (j in 1:length(out_x2[[i]])){
                    shan2[i,out_x2[[i]][j]]=-1
                }
            }
        }

        shan2[which(is.na(shan2))]=0
        names(dimnames(shan2))=c("Levels", "Contrasts") 
        contrast.matrix=shan2
}


    #####开始计算
                
    fit <- lmFit(data,design)
    fit2 <- contrasts.fit(fit, contrast.matrix) 
    fit2 <- eBayes(fit2) 
    
    if(length(table(de))==2){
    DEG<-topTable(fit2, n=nrow(fit2))
    return(DEG)
    }
    
    if(length(table(de)) > 2){
    out=alist()
    out[[1]]=fit2
    out[[2]]=decideTests(fit2)
    out[[3]]=topTableF(fit2, number=dim(data)[1])
    ####自动返回韦恩图
    out[[4]]=vennDiagram(out[[2]])
    
    l=length(colnames(contrast.matrix))
    colnames(out[[3]])[1:l]=colnames(contrast.matrix)

    return(out)
    }
}

生成测试样本

####规定生成多少基因、多少样本
n_gene = 1000; n_sample = 60
exp = runif(n_gene * n_sample, 1,2)


####生成矩阵
gene_matrix = matrix(exp, n_gene, n_sample)


rownames(gene_matrix) = paste0("gene_", 1:n_gene)
colnames(gene_matrix) = paste0("sam_", 1:n_sample)


####生成DEGs
gene_matrix[1:round(n_gene/3), 1:round(n_sample/3)] = gene_matrix[1:round(n_gene/3), 1:round(n_sample/3)] + 1
gene_matrix[(round(n_gene/3)*2):dim(gene_matrix)[1], (round(n_sample/3)*2):dim(gene_matrix)[2]] = gene_matrix[(round(n_gene/3)*2):dim(gene_matrix)[1], (round(n_sample/3)*2):dim(gene_matrix)[2]] + 1


####生成分组信息
group=vector(length=n_sample)
group[1:round(n_sample/3)]="g1"
group[(round(n_sample/3)*2):dim(gene_matrix)[2]]="g3"
group[which(group=="FALSE")]="g2"

计算

####计算3组的DEGs
deg=DEGs_find(data=gene_matrix, design=group)
####查看DEGs结果
head(deg[[3]])

####计算2组的DEGs
deg=DEGs_find(data=gene_matrix[,1:39], design=group[1:39])
####查看DEGs结果
head(deg)