美文网首页
2021-03-09

2021-03-09

作者: 世界很大_我想去体验 | 来源:发表于2021-03-09 23:13 被阅读0次

    WGCNA与免疫浸润相关分析

    rm(list=ls())
    library("WGCNA")            #引用WGCNA包
    library("limma")            #引用limma包
    expFile="uniq.symbol.txt"     #输入文件名字
    # normalCount=49              #正常样品数目
    # tumorCount=502              #肿瘤样品数目
    
    
    #读取文件,并对输入文件整理
    rt=read.table(expFile,sep="\t",header=T,check.names=F)
    rt=as.matrix(rt)
    rownames(rt)=rt[,1]
    exp=rt[,2:ncol(rt)]
    dimnames=list(rownames(exp),colnames(exp))
    data=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
    data=avereps(data)
    data=log2(data+1)
    data=data[apply(data,1,sd)>0.5,]
    datExpr0=t(data)
    group_1=sapply(strsplit(rownames(datExpr0),"\\-"),"[",4)
    group_2=sapply(strsplit(group_1,""),"[",1)
    group_2=gsub("2","1",group_2)
    group_2=as.data.frame(group_2)
    datExpr0=cbind(group_2,datExpr0)
    datExpr0 <- datExpr0[which(datExpr0$group_2 %in% 0),]
    dim(datExpr0)[2]
    datExpr0 <- datExpr0[,-1]
    
    ###检查缺失值
    gsg = goodSamplesGenes(datExpr0, verbose = 3)
    if (!gsg$allOK)
    {
      # Optionally, print the gene and sample names that were removed:
      if (sum(!gsg$goodGenes)>0)
        printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ", ")))
      if (sum(!gsg$goodSamples)>0)
        printFlush(paste("Removing samples:", paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", ")))
      # Remove the offending genes and samples from the data:
      datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
    }
    
    
    # 
    # #读取文件,并对输入文件整理
    # exp=read.table(expFile,sep="\t",header=T,row.names = 1,check.names = F)
    # exp=avereps(exp)
    # data=exp
    # data[1:4,1:4]
    # data=log2(data+1)
    # data=data[apply(data,1,sd)>0.5,]
    # datExpr0=t(data)
    # group_1=sapply(strsplit(rownames(datExpr0),"\\-"),"[",4)
    # group_2=sapply(strsplit(group_1,""),"[",1)
    # group_2=gsub("2","1",group_2)
    # group_2=as.data.frame(group_2)
    # datExpr0=cbind(group_2,datExpr0)
    # datExpr0=datExpr0[which(datExpr0$group_2==0),-1]
    # rownames(datExpr0)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*?)\\-.*",
    #                         "\\1\\-\\2\\-\\3",rownames(datExpr0))
    # 
    # ###检查缺失值
    # gsg = goodSamplesGenes(datExpr0, verbose = 3)
    # if (!gsg$allOK)
    # {
    #   # Optionally, print the gene and sample names that were removed:
    #   if (sum(!gsg$goodGenes)>0)
    #     printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ", ")))
    #   if (sum(!gsg$goodSamples)>0)
    #     printFlush(paste("Removing samples:", paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", ")))
    #   # Remove the offending genes and samples from the data:
    #   datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
    # }
    
    ###样品聚类
    sampleTree = hclust(dist(datExpr0), method = "average")
    pdf(file = "1_sample_cluster.pdf", width = 12, height = 9)
    par(cex = 0.6)
    par(mar = c(0,4,2,0))
    plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5, cex.axis = 1.5, cex.main = 2)
    ###剪切线
    abline(h = 19000, col = "red")
    dev.off()
    
    ###删除剪切线以下的样品
    clust = cutreeStatic(sampleTree, cutHeight = 19000, minSize = 10)
    table(clust)
    keepSamples = (clust==1)
    datExpr0 = datExpr0[keepSamples, ]
    
    
    ###准备临床数据
    clinical=read.table(file = "CIBERSORT-Results.txt",sep="\t",header = T,row.names = 1,check.names = F)
    clinical=na.omit(clinical)
    group_1=sapply(strsplit(rownames(clinical),"\\-"),"[",4)
    group_2=sapply(strsplit(group_1,""),"[",1)
    group_2=gsub("2","1",group_2)
    group_2=as.data.frame(group_2)
    clinical=cbind(group_2,clinical)
    clinical <- clinical[which(clinical$group_2 %in% 0),]
    dim(clinical)[2]
    clinical <- clinical[,-c(1,24,25,26)]
    flag <- apply(clinical,2,function(x) 
      sum(x == 0) < dim(clinical)[1]*4/5)
    clinical<- clinical[,which(flag)] 
    
    sample=as.character(rownames(datExpr0))
    traitsRows=match(sample,rownames(clinical))
    clinical=clinical[traitsRows,]
    
    
    ###样品聚类
    sampleTree2 = hclust(dist(datExpr0), method = "average")
    traitColors = numbers2colors(clinical, signed = FALSE)
    pdf(file="2_sample_heatmap.pdf",width=15,height=12)
    plotDendroAndColors(sampleTree2, traitColors,
                        groupLabels = names(clinical),
                        main = "Sample dendrogram and trait heatmap")
    dev.off()
    
    ###power值散点图
    enableWGCNAThreads()   #多线程工作
    powers = c(1:20)       #幂指数范围1:20
    sft = pickSoftThreshold(datExpr0, powerVector = powers, verbose = 5)
    pdf(file="3_scale_independence.pdf",width=9,height=5)
    par(mfrow = c(1,2))
    cex1 = 0.9
    ###拟合指数与power值散点图
    plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
         xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
         main = paste("Scale independence"));
    text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
         labels=powers,cex=cex1,col="red");
    abline(h=0.90,col="red") #可以修改
    ###平均连通性与power值散点图
    plot(sft$fitIndices[,1], sft$fitIndices[,5],
         xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
         main = paste("Mean connectivity"))
    text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
    dev.off()
    
    
    ###邻接矩阵转换
    sft #查看最佳power值
    softPower =sft$powerEstimate #最佳power值
    softPower
    adjacency = adjacency(datExpr0, power = softPower)
    
    ###TOM矩阵
    TOM = TOMsimilarity(adjacency);
    dissTOM = 1-TOM
    
    ###基因聚类
    geneTree = hclust(as.dist(dissTOM), method = "average");
    pdf(file="4_gene_clustering.pdf",width=12,height=9)
    plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
         labels = FALSE, hang = 0.04)
    dev.off()
    
    
    ###动态剪切模块识别
    minModuleSize = 50      #模块基因数目
    dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
                                deepSplit = 2, pamRespectsDendro = FALSE,
                                minClusterSize = minModuleSize);
    table(dynamicMods)
    dynamicColors = labels2colors(dynamicMods)
    table(dynamicColors)
    pdf(file="5_Dynamic_Tree.pdf",width=8,height=6)
    plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut",
                        dendroLabels = FALSE, hang = 0.03,
                        addGuide = TRUE, guideHang = 0.05,
                        main = "Gene dendrogram and module colors")
    dev.off()
    
    
    ###相似模块聚类
    MEList = moduleEigengenes(datExpr0, colors = dynamicColors)
    MEs = MEList$eigengenes
    MEDiss = 1-cor(MEs);
    METree = hclust(as.dist(MEDiss), method = "average")
    pdf(file="6_Clustering_module.pdf",width=7,height=6)
    plot(METree, main = "Clustering of module eigengenes",
         xlab = "", sub = "")
    MEDissThres = 0.25 #剪切高度可修改
    abline(h=MEDissThres, col = "red")
    dev.off()
    
    
    ###相似模块合并
    merge = mergeCloseModules(datExpr0, dynamicColors, cutHeight = MEDissThres, verbose = 3)
    mergedColors = merge$colors
    mergedMEs = merge$newMEs
    pdf(file="7_merged_dynamic.pdf", width = 9, height = 6)
    plotDendroAndColors(geneTree, mergedColors,"Dynamic Tree Cut",
                        dendroLabels = FALSE, hang = 0.03,
                        addGuide = TRUE, guideHang = 0.05,
                        main = "Gene dendrogram and module colors(TCGA)")
    dev.off()
    moduleColors = mergedColors
    table(moduleColors)
    colorOrder = c("grey", standardColors(50))
    moduleLabels = match(moduleColors, colorOrder)-1
    MEs = mergedMEs
    
    a=as.matrix(table(moduleColors))
    pdf("barplot.pdf",width = 21,height = 8)
    barplot(table(moduleColors),col = rownames(a))
    dev.off()
    ###模块与性状数据热图
    datTraits=clinical
    nGenes = ncol(datExpr0)
    nSamples = nrow(datExpr0)
    moduleTraitCor = cor(MEs, datTraits, use = "p")
    moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
    pdf(file="8_Module_trait.pdf",width=11,height=20)
    textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                       signif(moduleTraitPvalue, 1), ")", sep = "")
    dim(textMatrix) = dim(moduleTraitCor)
    par(mar = c(10, 10, 3, 3))
    labeledHeatmap(Matrix = moduleTraitCor,
                   xLabels = names(datTraits),
                   yLabels = names(MEs),
                   ySymbols = names(MEs),
                   colorLabels = FALSE,
                   colors = blueWhiteRed(50),
                   textMatrix = textMatrix,
                   setStdMargins = FALSE,
                   cex.text = 0.5,
                   zlim = c(-1,1),
                   main = paste("Module-trait relationships(TCGA)"))
    dev.off()
    
    
    ###计算MM和GS值
    modNames = substring(names(MEs), 3)
    geneModuleMembership = as.data.frame(cor(datExpr0, MEs, use = "p"))
    MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples))
    names(geneModuleMembership) = paste("MM", modNames, sep="")
    names(MMPvalue) = paste("p.MM", modNames, sep="")
    traitNames=names(datTraits)
    geneTraitSignificance = as.data.frame(cor(datExpr0, datTraits, use = "p"))
    GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))
    names(geneTraitSignificance) = paste("GS.", traitNames, sep="")
    names(GSPvalue) = paste("p.GS.", traitNames, sep="")
    
    
    ###批量输出性状和模块散点图
    for (trait in traitNames){
      traitColumn=match(trait,traitNames)  
      for (module in modNames){
        column = match(module, modNames)
        moduleGenes = moduleColors==module
        if (nrow(geneModuleMembership[moduleGenes,]) > 1){
          outPdf=paste("9_", trait, "_", module,".pdf",sep="")
          pdf(file=outPdf,width=7,height=7)
          par(mfrow = c(1,1))
          verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
                             abs(geneTraitSignificance[moduleGenes, traitColumn]),
                             xlab = paste("Module Membership in", module, "module"),
                             ylab = paste("Gene significance for ",trait),
                             main = paste("Module membership vs. gene significance\n"),
                             cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
          abline(v=0.8,h=0.5,col="red")
          dev.off()
        }
      }
    }
    
    
    ###输出GS_MM数据
    probes = colnames(datExpr0)
    geneInfo0 = data.frame(probes= probes,
                           moduleColor = moduleColors)
    for (Tra in 1:ncol(geneTraitSignificance))
    {
      oldNames = names(geneInfo0)
      geneInfo0 = data.frame(geneInfo0, geneTraitSignificance[,Tra],
                             GSPvalue[, Tra])
      names(geneInfo0) = c(oldNames,names(geneTraitSignificance)[Tra],
                           names(GSPvalue)[Tra])
    }
    
    for (mod in 1:ncol(geneModuleMembership))
    {
      oldNames = names(geneInfo0)
      geneInfo0 = data.frame(geneInfo0, geneModuleMembership[,mod],
                             MMPvalue[, mod])
      names(geneInfo0) = c(oldNames,names(geneModuleMembership)[mod],
                           names(MMPvalue)[mod])
    }
    geneOrder =order(geneInfo0$moduleColor)
    geneInfo = geneInfo0[geneOrder, ]
    write.table(geneInfo, file = "GS_MM.xls",sep="\t",row.names=F)
    
    
    ###输出每个模块的基因
    for (mod in 1:nrow(table(moduleColors)))
    {  
      modules = names(table(moduleColors))[mod]
      probes = colnames(datExpr0)
      inModule = (moduleColors == modules)
      modGenes = probes[inModule]
      write.table(modGenes, file =paste0("TCGA_",modules,".txt"),sep="\t",row.names=F,col.names=F,quote=F)
    }
    
    

    相关文章

      网友评论

          本文标题:2021-03-09

          本文链接:https://www.haomeiwen.com/subject/endiqltx.html