#安装需要的R包BiocManager::install("DESeq2")#载入需要的R包library(DESeq2)

#读取基因count矩阵文件expr <- read.table("input_counts.txt",sep = "t",header = T,check.names = F,stringsAsFactors = F,row.names = 1)# 读取样本分组信息文件subt <- read.table("input_subtype.txt", sep = "t", check.names = F, stringsAsFactors = F, header = T, row.names = 1)# 亚型名称n.sub.label <- unique(subt$TCGA_Subtype) # 亚型个数n.sub <- length(table(subt$TCGA_Subtype)) #创建配对比较的列表信息 group <- subt$TCGA_Subtypenames(group) <- rownames(subt)# 创建需要配对比较的列表函数,创建了三个分组。createList <- function(group=NULL) {  tumorsam <- names(group)  sampleList = list()  treatsamList =list()  treatnameList <- c()  ctrlnameList <- c()  #A-1: 类1 vs 其他  sampleList[[1]] = tumorsam  treatsamList[[1]] = intersect(tumorsam, names(group[group=="immune"])) # 亚型名称需要根据情况修改  treatnameList[1] <- "immune" # 该亚型的命名  ctrlnameList[1] <- "Others" # 其他亚型的命名  #A-2: 类2 vs 其他  sampleList[[2]] = tumorsam  treatsamList[[2]] = intersect(tumorsam, names(group[group=="keratin"]))  treatnameList[2] <- "keratin"  ctrlnameList[2] <- "Others"  #A-3: 类3 vs 其他  sampleList[[3]] = tumorsam  treatsamList[[3]] = intersect(tumorsam, names(group[group=="MITF-low"]))  treatnameList[3] <- "MITF-low"  ctrlnameList[3] <- "Others"  #如果有更多类，按以上规律继续写  return(list(sampleList, treatsamList, treatnameList, ctrlnameList))}complist <- createList(group=group)# 配对DESeq2函数twoclassDESeq2 <- function(res.path=NULL, countsTable=NULL, prefix=NULL, complist=NULL, overwt=FALSE) {  sampleList <- complist[[1]]  treatsamList <- complist[[2]]  treatnameList <- complist[[3]]  ctrlnameList <- complist[[4]]  allsamples <- colnames(countsTable)  options(warn=1)  for (k in 1:length(sampleList)) { # 循环读取每一次比较的内容    samples <- sampleList[[k]]    treatsam <- treatsamList[[k]]    treatname <- treatnameList[k]    ctrlname <- ctrlnameList[k]    compname <- paste(treatname, "_vs_", ctrlname, sep="") # 生成最终文件名    tmp = rep("others", times=length(allsamples))    names(tmp) <- allsamples    tmp[samples]="control"    tmp[treatsam]="treatment"    outfile <- file.path( res.path, paste(prefix, "_deseq2_test_result.", compname, ".txt", sep="") )    if (file.exists(outfile) & (overwt==FALSE)) { # 因为差异表达分析较慢，因此如果文件存在，在不覆盖的情况下（overwt=F）不再次计算差异表达      cat(k, ":", compname, "exists and skipped;n")      next    }    saminfo <- data.frame("Type"=tmp[samples],"SampleID"=samples,stringsAsFactors = F)    cts <- countsTable[,samples]    coldata <- saminfo[samples,]    # 差异表达过程，具体参数细节及输出结果解释，请参阅相关document    dds <- DESeqDataSetFromMatrix(countData = cts,                                  colData = coldata,                                  design = as.formula("~ Type")) # 设计矩阵仅包含亚型信息    dds$Type <- relevel(dds$Type,ref = "control")    dds <- DESeq(dds)    res <- results(dds, contrast=c("Type","treatment","control"))    #将分析结果转化为数据框resData <- as.data.frame(res[order(res$padj),])    #将行名作为id列resData$id <- rownames(resData)    #提取想要的列数据resData <- resData[,c("id","baseMean","log2FoldChange","lfcSE","stat","pvalue","padj")]    #修改列名colnames(resData) <- c("id","baseMean","log2FC","lfcSE","stat","PValue","FDR")    #输出到文件    write.table(resData, file=outfile, row.names=F, col.names=T, sep="t", quote=F)    cat(k, ",")  }  options(warn=0)}# 差异表达分析过程比较慢请耐心等待twoclassDESeq2(res.path = ".", #所有配对差异表达结果都会输出在res.path路径下               countsTable = expr[,intersect(colnames(expr),rownames(subt))],               prefix = "SKCM", #文件名以SKCM开头               complist = complist,               overwt = F)