• wgcna 原文复现 小胶质细胞亚群在脑发育时髓鞘形成的作用 microglial

    手把手10分文章WGCNA复现：小胶质细胞亚群在脑发育时髓鞘形成的作用

    Original 生信技能树 生信技能树 2020-01-09 14:2

    Hi大家好，我是老米！生信入门一个月，感谢生信技能树平台。这是我的第3篇Markdown。不知道多少人还记得我的前两个投稿：

    现在继续完成生信技能树的作业题：重复一篇WGCNA分析的文章。WGCNA学习教程：一文学会WGCNA分析。

    复现的文献：A novel microglial subset plays a key role in myelinogenesis in developing brain。EMBO J，好杂志。

    该文章对五个处理组，共17个老鼠：

    • orange represents neonatal CD11c+ microglia (n = 4),

    • green neonatal CD11c microglia (n = 4),

    • blue EAE CD11c+ microglia (n = 3),

    • purple EAE CD11c microglia (n = 3),

    • black adult microglia (n = 3).

    其实就是 neonatal(新生的) 和 EAE的Microglia，还有CD11C阳性和阴性，然后和成年小鼠的Microglia进行比较。这些样本进行了RNASeq测序，数据在GEO可供下载：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE78809。当然我偷了个懒，该文章Supplemental material提供了整理之后的csv矩阵，大概1万3千个基因。

    作者进行了WGCNA分析，以此证明小胶质细胞的作用及寻找关键信号通路及基因（嗯，估计这个实验室后续还有大文章出来）。

    组学及生信分析的篇幅占了整个文章的一半以上。在此，仅重复文章的Figure3。

           其实我跌跌撞撞学了WGCNA半个月，大概了解了这个原理，现在还一知半解，对很多细节摸索了很久，头很大。下面进入正题，请大家一起跟着我头大。

    1.数据处理

    1. rm(list = ls())
    2. dev.off()
    3. library(WGCNA)# 偷懒，读取作者给出的附件的表达矩阵
    4. norm.counts <- read.csv("Dataset_EV2.csv",header = T,sep = ",") ##original counts
    5. rownames(norm.counts)<-norm.counts[,1]
    6. norm.counts <- norm.counts[,2:18]
    7. save(norm.counts,file = "norm.counts.Rdata")
    8.  
    1. ## 处理表型信息，就是老鼠的分组信息
    2. a <- colnames(norm.counts)
    3. library(stringr)
    4. str_tmp=as.data.frame(str_split(a,"[_]",simplify = T))
    5. rownames(str_tmp) <- colnames(norm.counts)
    6. colnames(str_tmp) <- c("group","sampleRep")
    1. datTraits <- str_tmp
    2. datTraits$groupNo <- c(rep(1,3),rep(2,4),rep(3,4),rep(4,3),rep(5,3)) 
    1. ###数据初步处理完成
    2.  
    3. ## 筛选中位绝对偏差前75%的基因，取MAD大于0.3
    4. ## 筛选后会降低运算量，也会失去部分信息
    5. ## 也可不做筛选，使MAD大于0即可
    6. norm.counts <- log2(norm.counts)  ##注意我在这里进行了log2，因为我发现如果不log的话，出现了个很有意思的现象，请您往下看
    7. m.mad <- apply(norm.counts,1,mad)
    8. dataExprMad <- norm.counts[which(m.mad > max(quantile(m.mad, probs=seq(0, 1, 0.25))[2],0.3)),] #25%  = 0.22
    9. datExpr0 <- as.data.frame(t(dataExprMad))
    1. ##最终取到8223个基因。
    2.  
    3. ###########start handling the data
    4. ##看看有没有为空的值，需要剔除。这里下载的都是作者处理好的，基本没问题。
    5. gsg <- goodSamplesGenes(datExpr0,verbose = 3)
    6. gsg$allOK
    7. if (!gsg$allOK){
    8.   # Optionally, print the gene and sample names that were removed:
    9.   if (sum(!gsg$goodGenes)>0)
    10.     printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes],
    11.                                               collapse = ", ")));
    12.   if (sum(!gsg$goodSamples)>0)
    13.     printFlush(paste("Removing samples:",
    14.                      paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", ")));
    15.   # Remove the offending genes and samples from the data:
    16.   datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
    17. }
    18. gsg <- goodSamplesGenes(datExpr0,verbose = 3)
    19. gsg$allOK
    20.  
    21. ###画个树看看样本有没有离群值！！！
    22. if(T){
    23.   # Plot a line to show the cut
    24.   sampleTree <- hclust(dist(datExpr0),method = "average")
    25.   par(cex=.6)
    26.   par(mar=c(0,4,2,0))
    27.   plot(sampleTree, main = "Sample clustering to detect outliers", sub = "",
    28.        xlab = "",cex.lab = 1.5, cex.axis =1.5, cex.main=2
    29.   )
    30. }
    31.  
    32. ##然后保存数据
    33. dim(datExpr0) ##17，8222
    34. datExpr <- datExpr0
    35. nGenes <- ncol(datExpr)
    36. nSamples <- nrow(datExpr)
    37. save(nGenes,nSamples,datExpr,datTraits,file="input.Rdata")

    结果如图：

           这结果挺好的，符合实验设计，都聚类好了。但是！但是如果不对原始数值进行log2处理，会出现一个outlier老鼠，如图：

           看到没，那个新生的NEO1号鼠有点调皮！他超脱于所有老鼠之上，统领其他所有老鼠！这要么是问题少年，要么是智商超群，比如科大少年班的那种！而我大概深究了一下，控制这个离群值的大概是500个基因，哈哈哈。。。

    2. MDS分析

           专业名词叫多维标度法（Multidimensional Scaling）MDS，是一种维数缩减方法，把高维的数据点映射到一个低维的流形上；同时也是一种可视化方法，实践中通常利用2D或3D的MDS 结果观察（投影后）点的分布和聚集来研究数据的性质。实际上你也可以把它看做是PCA分析，反正都是看样本间的相关性计算距离远近。

    说人话叫做：通过基因表达log2组间差异，看看这17个老鼠是否能够分到一个类别，比如成年鼠归到一类，EAE归到一类，新生鼠归到一类。和上面的有点类似，又不完全一样。

           参考学习了stat的视频，也参考了他提供的代码。如下：

    1. rm(list = ls())
    2. load(file = 'input.Rdata')
    3. log2.data.matrix <- as.data.frame(t(datExpr))
    4. ## now create an empty distance matrix
    5. log2.distance.matrix <- matrix(0,
    6.                                nrow=ncol(log2.data.matrix),
    7.                                ncol=ncol(log2.data.matrix),
    8.                                dimnames=list(colnames(log2.data.matrix),
    9.                                              colnames(log2.data.matrix)))
    10. ## now compute the distance matrix using avg(absolute value(log2(FC)))
    11. for(i in 1:ncol(log2.distance.matrix)) {
    12.   for(j in 1:i) {
    13.     log2.distance.matrix[i, j] <-
    14.       mean(abs(log2.data.matrix[,i] - log2.data.matrix[,j]))
    15.   }
    16. }
    17. ## do the MDS math (this is basically eigen value decomposition)
    18. ## cmdscale() is the function for "Classical Multi-Dimensional Scalign"
    19. mds.stuff <- cmdscale(as.dist(log2.distance.matrix),
    20.                       eig=TRUE,
    21.                       x.ret=TRUE)
    22. save(mds.stuff,file = "step7_MDS.Rdata")
    23. ## calculate the percentage of variation that each MDS axis accounts for...
    24. mds.var.per <- round(mds.stuff$eig/sum(mds.stuff$eig)*100, 1)
    25. ## now make a fancy looking plot that shows the MDS axes and the variation:
    26. mds.values <- mds.stuff$points
    27. mds.data <- data.frame(Sample=rownames(mds.values),
    28.                        X=mds.values[,1],
    29.                        Y=mds.values[,2])
    30. ## 然后画图
    31. library(ggpubr)
    32. library(ggplot2)
    33. groups <- datTraits$group
    34. png("figures/step7_MDS_logfc.png",width = 800,height=600)
    35. ggplot(mds.data, aes(x=X, y=Y, label=Sample,col=groups)) + 
    36.   geom_point(size = 10, alpha = 0.8) +
    37.   theme(panel.grid.minor = element_blank()) +
    38.   coord_fixed() + theme_bw()+
    39.   ggtitle("MDS plot using avg(logFC) as the distance")+
    40.   xlab(paste("Leading logFC dim1 ", mds.var.per[1], "%", sep="")) +
    41.   ylab(paste("Leading logFC dim2 ", mds.var.per[2], "%", sep="")) +
    42.   ggtitle("MDS plot using avg(logFC) as the distance")
    43. dev.off()

    结果如下：

           嗯，，，结果表示完美！别问我为什么要先做这个图，因为它好看，我忍不住！

    点评：老米在这里开始秀技术了，实际上大家随便跑一下PCA分析即可，暂时无需深入写这么长代码！

           下面正式进入WGCNA分析。

    3. 计算beta值

      这个没啥好讲的，代码如下：

    1. rm(list = ls())
    2. library(WGCNA)
    3. load(file = "input.Rdata")
    4. dim(datExpr)
    5. #################################################
    6. if(T){
    7.         powers = c(c(1:10), seq(from = 12, to=30, by=2))
    8.         # Call the network topology analysis function
    9.         sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
    10.         png("figures/step2-beta-value.png",width = 800,height = 600)
    11.         # Plot the results:
    12.         ##sizeGrWindow(9, 5)
    13.         par(mfrow = c(1,2));
    14.         cex1 = 0.9;
    15.         # Scale-free topology fit index as a function of the soft-thresholding power
    16.         plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
    17.              xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
    18.              main = paste("Scale independence"));
    19.         text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
    20.              labels=powers,cex=cex1,col="red");
    21.         # this line corresponds to using an R^2 cut-off of h
    22.         abline(h=0.9,col="red")
    23.         # Mean connectivity as a function of the soft-thresholding power
    24.         plot(sft$fitIndices[,1], sft$fitIndices[,5],
    25.              xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
    26.              main = paste("Mean connectivity"))
    27.         text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
    28.         dev.off()
    29. }
    30. sft$powerEstimate  ## beta=22 SCI文章里面用了20
    31. save(sft,file = "step2_beta_value.Rdata")

    4. 采用一步法构建加权共表达网络

    (weight co-expression network)

    1. rm(list = ls())
    2. library(WGCNA)
    3. load(file = "input.Rdata")
    4. load(file = "step2_beta_value.Rdata")
    5. enableWGCNAThreads()
    6. if(T){
    7.   net = blockwiseModules(
    8.     datExpr,
    9.     power = sft$powerEstimate,
    10.     maxBlockSize = nGenes,
    11.     TOMType = "unsigned", minModuleSize = 30,
    12.     reassignThreshold = 0, mergeCutHeight = 0.28, ## 注意这个0.28
    13.     numericLabels = TRUE, pamRespectsDendro = FALSE,
    14.     saveTOMs = F,
    15.     verbose = 3
    16.   )
    17.   table(net$colors) 
    18. }
    19.  
    20. ##模块可视化，画那个传说中的树
    21. if(T){
    22.   # Convert labels to colors for plotting
    23.   moduleColors=labels2colors(net$colors)
    24.   table(moduleColors)
    25.   # Plot the dendrogram and the module colors underneath
    26.   png("figures/step3-genes-modules.png",width = 800,height = 600)
    27.   plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]],
    28.                       "Module colors",
    29.                       dendroLabels = FALSE, hang = 0.03,
    30.                       addGuide = TRUE, guideHang = 0.05)
    31.   dev.off()
    32. }
    33. save(net,moduleColors,file = "step3_genes_modules.Rdata")

           

    结果如下：

            这里需要停顿一下：SCI文章只取到8个模块，图注说用了12691个基因（被我发现全部加起来也才9999，没错！9999个基因）。事实上，能够分这么少模块，我不知道作者参数细节。有两个参数特别影响模块大小：

    • 一个是minModuleSize，表示每个模块里面最少放多少个基因，这很好理解，设定越大，模块越少；

    • 另一个是mergeCutHeight，这个参数表示你在哪里砍树。值设定越小，树枝越多，通常是0.25。

    google一下，发现这么一个说法：

    I am not aware of a principle from which

          这里，我的mergeCutHeight = 0.28，最终取到13个模块！

    5. 模块与基因与表型

    5.1 模块和老鼠表型对应

    这是关键，其实就是开始解读这些关联的基因群对老鼠表型的影响，比如敲除CDC11之后的老鼠，哪些信号通路激活或抑制，EAE的老鼠哪些基因激活？等等等。

    1. rm(list = ls())
    2. library(WGCNA)
    3. load(file = "input.Rdata")
    4. load(file = "step2_beta_value.Rdata")
    5. load(file = "step3_genes_modules.Rdata")
    6.  
    7. if(T){
    8.   nGenes = ncol(datExpr)
    9.   nSamples = nrow(datExpr)
    10.   design <- model.matrix(~0+datTraits$group)
    11.   colnames(design)= levels(datTraits$group) ## get the group
    12.   MES0 <- moduleEigengenes(datExpr,moduleColors)$eigengenes
    13.   MEs = orderMEs(MES0)
    14.   moduleTraitCor <- cor(MEs,design,use = "p")
    15.   moduleTraitPvalue <- corPvalueStudent(moduleTraitCor,nSamples)
    16.   textMatrix = paste(signif(moduleTraitCor,2),"\n(",
    17.                      signif(moduleTraitPvalue,1),")",sep = "")
    18.   dim(textMatrix)=dim(moduleTraitCor)
    19.  
    20.   png("figures/step4-Module-trait-relationship.png",width = 800,height = 1200,res = 120)
    21.   par(mar=c(6, 8.5, 3, 3))
    22.   labeledHeatmap(Matrix = moduleTraitCor,
    23.                  xLabels = colnames(design),
    24.                  yLabels = names(MEs),
    25.                  ySymbols = names(MEs),
    26.                  colorLabels = FALSE,
    27.                  colors = blueWhiteRed(50),
    28.                  textMatrix = textMatrix,
    29.                  setStdMargins = FALSE,
    30.                  cex.text = 0.5,
    31.                  zlim = c(-1,1),
    32.                  main = paste("Module-trait relationships"))
    33.   dev.off()
    34.   save(design,file = "step4_design.Rdata")
    35. }

     

    通过模块与各种表型的相关系数，可发现自己感兴趣的模块。如图：

    可以看到有些模块与成年鼠关联，有些与EAE相关性很高。和原SCI的图有点像。

    5.2 各模块与表型相关性bar图

    就是上面的图转化成bar图。

    1. if(T){
    2.   mes_group <- merge(MEs,datTraits,by="row.names") # 小技巧，可以通过rowname进行merge
    3.   ##写了个画图的function
    4.   library(gplots)
    5.   library(ggplot2)
    6.   library(ggpubr)
    7.   library(grid)
    8.   library(gridExtra) 
    9.   draw_ggboxplot <- function(data,gene="P53",group="group"){
    10.     #print(gene)
    11.     ggboxplot(data,x=group, y=gene,
    12.               ylab = sprintf("Expression of %s",gene),
    13.               xlab = group,
    14.               fill = group,
    15.               palette = "jco",
    16.               #add="jitter",
    17.               legend = "") +stat_compare_means()
    18.   }
    19.   ###开始批量画boxplot
    20.   colorNames = names(MEs)
    21.   png("figures/step4-expression-group.png",width = 800,height=2000,res = 120)
    22.   #par(mfrow=c(ceiling(length(colorNames)/2),2))
    23.   p <- lapply(colorNames,function(x) {
    24.     draw_ggboxplot(mes_group,gene= x,group="group")
    25.   })
    26.   do.call(grid.arrange,c(p,ncol=2))
    27.   dev.off()
    28. }

    mes_group <- merge(MEs,datTraits,by="row.names") # 小技巧，可以通过rowname进行merge 

    部分结果如下：

    5.3 模块与基因的相关性

    1. if(T){
    2.   modNames = substring(names(MEs), 3)
    3.   geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"))
    4.   ## 算出每个模块跟基因的皮尔森相关系数矩阵
    5.   ## MEs是每个模块在每个样本里面的值
    6.   ## datExpr是每个基因在每个样本的表达量
    7.   MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples))
    8.   names(geneModuleMembership) = paste("MM", modNames, sep="")
    9.   names(MMPvalue) = paste("p.MM", modNames, sep="")
    10.  
    11.   geneTraitSignificance <- as.data.frame(cor(datExpr,datTraits$groupNo,use = "p"))
    12.   GSPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance),nSamples))
    13.   names(geneTraitSignificance)<- paste("GS.",names(datTraits$group),sep = "")
    14.   names(GSPvalue)<-paste("GS.",names(datTraits$group),sep = "")
    15.  
    16.   #selectModule<-c("blue","green","purple","grey")  ##可以选择自己喜欢的模块
    17.   selectModule <- modNames  ## 也可以批量作图
    18.   png("figures/step4-Module-trait-significance.png",width = 800,height=2000,res = 120)
    19.   par(mfrow=c(ceiling(length(selectModule)/2),2)) #批量作图开始
    20.   for(module in selectModule){
    21.     column <- match(module,selectModule)
    22.     print(module)
    23.     moduleGenes <- moduleColors==module
    24.     verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
    25.                        abs(geneTraitSignificance[moduleGenes, 1]),
    26.                        xlab = paste("Module Membership in", module, "module"),
    27.                        ylab = paste("Gene significance for", module, "module"),
    28.                        main = paste("Module membership vs. gene significance\n"),
    29.                        cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
    30.   }
    31.   dev.off()
    32. }

     

     

     

    截取部分如下：

    5.4 关注感兴趣的模块，导出基因进行GO分析

     

     

     

     

    1. rm(list = ls())
    2. library(WGCNA)
    3. load(file = 'input.Rdata')
    4. load(file = "step2_beta_value.Rdata")
    5. load(file = "step3_genes_modules.Rdata")
    6. load(file = "step4_design.Rdata")
    7. load(file="step4_datTraits.Rdata")
    8.  
    9. table(moduleColors)
    10. group_g <- data.frame(gene=colnames(datExpr),
    11.                       group=moduleColors)
    12. write.csv(group_g,file = "figures/group_g.csv") ## 导出了对应模块所有基因
    13.  
    14. # 选择mouse的基因组进行注释及ID转化啥的，如果是人的，另有R包
    15. if(!require("clusterProfiler")) BiocManager::install("clusterProfiler",ask = F,update = F)
    16. library(clusterProfiler)
    17. if(!require("org.Mm.eg.db")) BiocManager::install("org.Mm.eg.db",ask = F,update = F)
    18. library(org.Mm.eg.db)
    19. tmp <- bitr(group_g$gene,fromType = "ENSEMBL",
    20.             toType = "ENTREZID",
    21.             OrgDb = "org.Mm.eg.db")
    22. de_gene_cluster <- merge(tmp,group_g,by.x="ENSEMBL",by.y="gene")
    23. table(de_gene_cluster$group)
    24.  
    25. ###run go analysis
    26. formula_res <- compareCluster(
    27.   ENTREZID~group,
    28.   data = de_gene_cluster,
    29.   fun = "enrichGO",
    30.   OrgDb = "org.Mm.eg.db",
    31.   style="margin: 0px; padding: 0px; font-size: inherit; line-height: inherit; color: rgb(80, 161, 79); overflow-wrap: inherit !important; word-break: inherit !important;">"BP",
    32.   pAdjustMethod = "BH",
    33.   pvalueCutoff = 0.01,
    34.   qvalueCutoff = 0.05
    35. )
    36.  
    37. lineage1_ego <-simplify(
    38.   formula_res,
    39.   cutoff=0.5,
    40.   by="p.adjust",
    41.   select_fun=min
    42. )
    43. save(group_g,formula_res,lineage1_ego,file="step5_GOananlysis.Rdata")
    44. #出图
    45. pdf("figures/step5_Microglia_GO_term_DE.pdf",width = 15,height = 15)
    46. dotplot(lineage1_ego,showCategory=10)
    47. dev.off()
    48. write.csv(lineage1_ego@compareClusterResult,
    49.           file="figures/Microglia_GO_term_DE.csv")

     

     

    GO结果分析如下：

           

    通过与上面的几个图表配合，发现EAE鼠与炎症通路/免疫相关，而新生鼠的胶质细胞则是cell cycle/胶质细胞分化/神经细胞发育/轴树突形成通路相关。与原SCI分析相符，与预期相符。

    5.5 画WGCNA的标配热图

    1. rm(list = ls())
    2. library(WGCNA)
    3. load(file = 'input.Rdata')
    4. load(file = "step2_beta_value.Rdata")
    5. load(file = "step3_genes_modules.Rdata")
    6. load(file = "step4_design.Rdata")
    7. load(file="step4_datTraits.Rdata")
    8. load(file="step5_GOananlysis.Rdata")
    9. # 主要是可视化 TOM矩阵，WGCNA的标准配图
    10. # 然后可视化不同 模块 的相关性 热图
    11. # 不同模块的层次聚类图
    12. # 还有模块诊断，主要是 intramodular connectivity
    13. if(T){
    14.   #geneTree = net$dendrograms[[1]]
    15.   TOM=TOMsimilarityFromExpr(datExpr,power=20)
    16.   dissTOM=1-TOM
    17.   #plotTOM = dissTOM^7
    18.   #diag(plotTOM)=NA
    19.   #TOMplot(plotTOM,geneTree,moduleColors,main="Network heapmap plot of all genes")
    20.   ### 我这里只取了1000个基因哈，我试了一下全部基因，结果跑了半个小时没跑完，被我强行退出！
    21.   nSelect =1000
    22.   set.seed(20)
    23.   select=sample(nGenes,size = nSelect)
    24.   selectTOM = dissTOM[select,select]
    25.   selectTree = hclust(as.dist(selectTOM),method = "average")
    26.   selectColors = moduleColors[select]
    27.   plotDiss=selectTOM^7
    28.   diag(plotDiss)=NA
    29.   png("figures/step6_select_Network-heatmap.png",width = 800,height=600)
    30.   TOMplot(plotDiss,selectTree,selectColors,main="Network heapmap of selected gene")
    31.   dev.off()
    32. }

    如下：

    5.6 提取指定模块的基因并做热图

    1. if(T){
    2.   module="turquoise"
    3.   which.module=module
    4.   dat=datExpr[,moduleColors==which.module]
    5.   library(pheatmap)
    6.   pheatmap(dat,show_colnames = F,show_rownames = F)
    7.   n=scale(t(dat+1)) # 'scale'可以对log-ratio数值进行归一化
    8.   n[n>2]=2 
    9.   n[n< -2]= -2
    10.   n[1:4,1:4]
    11.   pheatmap(n,show_colnames =F,show_rownames = F)
    12.   group_list=datTraits$group
    13.   ac=data.frame(g=group_list)
    14.   rownames(ac)=colnames(n)
    15.   png("figures/step6-moduleGene-heatmap.png",width = 800,height = 600)
    16.   pheatmap(n,show_colnames =F,show_rownames = F,annotation_col =ac )
    17.   dev.off()
    18.  
    19. }

    都挺好的。

    点评：这个热图有问题，具体代码是  n=scale(t(dat+1)) 里面少了一个转置！

    5.7 性状与模块的关系

    1. if(T){
    2.   ## 连续性的变量，如体重等才好和模块进行比较。
    3.   MEs=moduleEigengenes(datExpr,moduleColors)$eigengenes
    4.   MET = orderMEs(cbind(MEs,datTraits$groupNo))
    5.   par(cex = 0.9)
    6.   png("figures/step6-Eigengene-dendrogram.png",width = 800,height = 600)
    7.   plotEigengeneNetworks(MET, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle
    8.                         = 90)
    9.   dev.off()
    10. }

    如图：

    6. 模块导出

    感兴趣的模块导出到cytoscape，VisANT等软件进一步进行可视化分析。

    1. #######gene export for VisANT or cytoscape
    2. if(T){
    3.  
    4.   module="turquoise"
    5.   probes = colnames(datExpr)
    6.   inModule = (moduleColors==module)
    7.   modProbes=probes[inModule]
    8.   head(modProbes)
    9.   modTOM = TOM[inModule,inModule]
    10.   dimnames(modTOM)=list(modProbes,modProbes)
    11.   ### 这里只选了top100的基因
    12.   nTop=100
    13.   IMConn = softConnectivity(datExpr[,modProbes])
    14.   top=(rank(-IMConn)<=nTop)
    15.   filterTOM=modTOM[top,top]
    16.   # for visANT
    17.   vis = exportNetworkToVisANT(filterTOM,file = paste("figures/visANTinput-",module,".txt",sep = ""),
    18.                               weighted = T,threshold = 0)
    19.  
    20.   # for cytoscape
    21.   cyt = exportNetworkToCytoscape(filterTOM,
    22.                                  edgeFile = paste("figures/CytoscapeInput-edges-", paste(module, collapse="-"), ".txt", sep=""),
    23.                                  nodeFile = paste("figures/CytoscapeInput-nodes-", paste(module, collapse="-"), ".txt", sep=""),
    24.                                  weighted = TRUE,
    25.                                  threshold = 0.02,
    26.                                  nodeNames = modProbes[top], 
    27.                                  nodeAttr = moduleColors[inModule][top])
    28. }

    总结

           至此，复现结束。如果你仔细的看到了这里，说明你是想学WGCNA的了。回到课题本身，是否对小胶质细胞的认识更近了一步？提个小小的科学问题：在PD中，是否有一个胶质细胞亚群可对MPP+诱导的黑质神经细胞死亡有抵抗和保护作用？哈哈，别说我是搞肿瘤的了。。。

           这些内容，换个思路应该够一个硕士毕业了。。。

    1. rm(list = ls())
    2. dev.off()
    3. library(WGCNA)
    4. # 偷懒，读取作者给出的附件的表达矩阵
    5. getwd()
    6. file="G:/papers_for_ARDS/r_scripts_for_ARDS_From_zhongda/lncRNA_mRNA_cicleRNA/microglial cells/EMBJ-36-3292-s002/embj201696056-sup-0003-DatasetEV1"
    7. dir.create(file)
    8. setwd(file)
    9. getwd()
    10. ??read.csv
    11. norm.counts <- openxlsx::read.xlsx("Dataset_EV1.xlsx",
    12.                         sheet = 2) ##original counts
    13. head(norm.counts)
    14. rownames(norm.counts)<-norm.counts[,1]
    15. norm.counts <- norm.counts[,2:18]
    16. head(norm.counts)
    17. save(norm.counts,file = "norm.counts.Rdata")
    18. load("norm.counts.Rdata")
    19.  
    20. ## 处理表型信息，就是老鼠的分组信息
    21. a <- colnames(norm.counts)
    22. library(stringr)
    23. str_split(a,"[,]",simplify = T)
    24. strsplit("UK, USA, Germany", ",(?=[^,]+$)", perl=TRUE)
    25. str_split(a, ",(?=[^,]+$)",simplify = T)
    26.  
    27. #以最后一个逗号为分割符号 切割字符串
    28. str_tmp=as.data.frame(str_split(a, ",(?=[^,]+$)",simplify = T))
    29. rownames(str_tmp) <- colnames(norm.counts)
    30. colnames(str_tmp) <- c("group","sampleRep")
    31. str_tmp
    32. datTraits <- str_tmp
    33. datTraits$groupNo <- c(rep(1,3),rep(2,4),rep(3,4),rep(4,3),rep(5,3)) 
    34. datTraits
    35. ###数据初步处理完成
    36.  
    37. ## 筛选中位绝对偏差前75%的基因，取MAD大于0.3
    38. ## 筛选后会降低运算量，也会失去部分信息
    39. ## 也可不做筛选，使MAD大于0即可
    40. head(norm.counts)
    41. norm.counts <- log2(norm.counts)  ##注意我在这里进行了log2，因为我发现如果不log的话，出现了个很有意思的现象，请您往下看
    42. m.mad <- apply(norm.counts,1,mad)
    43. dataExprMad <- norm.counts[which(m.mad > max(quantile(m.mad, probs=seq(0, 1, 0.25))[2],0.3)),] #25%  = 0.22
    44. datExpr0 <- as.data.frame(t(dataExprMad))
    45. dim(datExpr0)
    46. head(datExpr0)[1:10,1:10]
    47. ##最终取到8222个基因。
    48.  
    49.  
    50.  
    51. ###########start handling the data
    52. ##看看有没有为空的值，需要剔除。这里下载的都是作者处理好的，基本没问题。
    53. gsg <- goodSamplesGenes(datExpr0,verbose = 3)
    54. gsg$allOK
    55. if (!gsg$allOK){
    56.   # Optionally, print the gene and sample names that were removed:
    57.   if (sum(!gsg$goodGenes)>0)
    58.     printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes],
    59.                                               collapse = ", ")));
    60.   if (sum(!gsg$goodSamples)>0)
    61.     printFlush(paste("Removing samples:",
    62.                      paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", ")));
    63.   # Remove the offending genes and samples from the data:
    64.   datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
    65. }
    66. gsg <- goodSamplesGenes(datExpr0,verbose = 3)
    67. gsg$allOK
    68.  
    69.  
    70. ###画个树看看样本有没有离群值！！！
    71. if(T){
    72.   # Plot a line to show the cut
    73.   sampleTree <- hclust(dist(datExpr0),method = "average")
    74.   par(cex=.6)
    75.   par(mar=c(0,4,2,0))
    76.   plot(sampleTree, main = "Sample clustering to detect outliers", sub = "",
    77.        xlab = "",cex.lab = 1.5, cex.axis =1.5, cex.main=2
    78.   )
    79. }
    80.  
    81. ##然后保存数据
    82. dim(datExpr0) ##17，8222
    83. datExpr <- datExpr0
    84. nGenes <- ncol(datExpr)
    85. nSamples <- nrow(datExpr)
    86. save(nGenes,nSamples,datExpr,datTraits,file="input.Rdata")
    87.  
    88.  
    89.  
    90.  
    91. rm(list = ls())
    92. library(WGCNA)
    93. load(file = "input.Rdata")
    94. dim(datExpr)
    95. head(datTraits)
    96. #################################################
    97. if(T){
    98.   powers = c(c(1:10), seq(from = 12, to=30, by=2))
    99.   # Call the network topology analysis function
    100.   sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
    101.   
    102.   dev.off()
    103.   getwd()
    104.   dir.create("./figures")
    105.   png("figures/step2-beta-value.png",width = 800,height = 600)
    106.   # Plot the results:
    107.   ##sizeGrWindow(9, 5)
    108.   par(mfrow = c(1,2));
    109.   cex1 = 0.9;
    110.   # Scale-free topology fit index as a function of the soft-thresholding power
    111.   plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
    112.        xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
    113.        main = paste("Scale independence"));
    114.   text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
    115.        labels=powers,cex=cex1,col="red");
    116.   # this line corresponds to using an R^2 cut-off of h
    117.   abline(h=0.9,col="red")
    118.   # Mean connectivity as a function of the soft-thresholding power
    119.   plot(sft$fitIndices[,1], sft$fitIndices[,5],
    120.        xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
    121.        main = paste("Mean connectivity"))
    122.   text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
    123.   dev.off()
    124. }
    125. sft$powerEstimate  ## beta=22 SCI文章里面用了20
    126. save(sft,file = "step2_beta_value.Rdata")
    127.  
    128.  
    129.  
    130.  
    131.  
    132. rm(list = ls())
    133. library(WGCNA)
    134. load(file = "input.Rdata")
    135. load(file = "step2_beta_value.Rdata")
    136. enableWGCNAThreads()
    137. if(T){
    138.   net = blockwiseModules(
    139.     datExpr,
    140.     power = sft$powerEstimate,
    141.     maxBlockSize = nGenes,
    142.     TOMType = "unsigned", minModuleSize = 30,
    143.     reassignThreshold = 0, mergeCutHeight = 0.28, ## 注意这个0.28
    144.     numericLabels = TRUE, pamRespectsDendro = FALSE,
    145.     saveTOMs = F,
    146.     verbose = 3
    147.   )
    148.   table(net$colors) 
    149. }
    150.  
    151. ##模块可视化，画那个传说中的树
    152. if(T){
    153.   # Convert labels to colors for plotting
    154.   moduleColors=labels2colors(net$colors)
    155.   table(moduleColors)
    156.   # Plot the dendrogram and the module colors underneath
    157.   png("figures/step3-genes-modules.png",width = 800,height = 600)
    158.   plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]],
    159.                       "Module colors",
    160.                       dendroLabels = FALSE, hang = 0.03,
    161.                       addGuide = TRUE, guideHang = 0.05)
    162.   dev.off()
    163. }
    164. save(net,moduleColors,file = "step3_genes_modules.Rdata")
    165.  
    166.  
    167.  
    168. rm(list = ls())
    169. library(WGCNA)
    170. load(file = "input.Rdata")
    171. load(file = "step2_beta_value.Rdata")
    172. load(file = "step3_genes_modules.Rdata")
    173.  
    174. if(T){
    175.   nGenes = ncol(datExpr)
    176.   nSamples = nrow(datExpr)
    177.   datTraits
    178.   design <- model.matrix(~0+datTraits$group)
    179.   colnames(design)= levels(as.factor(datTraits$group)) ## get the group
    180.   design
    181.   MES0 <- moduleEigengenes(datExpr,moduleColors)$eigengenes
    182.   head(MES0)
    183.   MEs = orderMEs(MES0)
    184.   head(MEs)
    185.   head(design)
    186.   moduleTraitCor <- cor(MEs,design,use = "p")
    187.   head(moduleTraitCor)
    188.   moduleTraitPvalue <- corPvalueStudent(moduleTraitCor,nSamples)
    189.   textMatrix = paste(signif(moduleTraitCor,2),"\n(",
    190.                      signif(moduleTraitPvalue,1),")",sep = "")
    191.   dim(textMatrix)=dim(moduleTraitCor)
    192.   
    193.   textMatrix
    194.   
    195.   png("figures/step4-Module-trait-relationship.png",width = 800,height = 1200,res = 120)
    196.   par(mar=c(6, 8.5, 3, 3))
    197.   labeledHeatmap(Matrix = moduleTraitCor,
    198.                  xLabels = colnames(design),
    199.                  yLabels = names(MEs),
    200.                  ySymbols = names(MEs),
    201.                  colorLabels = FALSE,
    202.                  colors = blueWhiteRed(50),
    203.                  textMatrix = textMatrix,
    204.                  setStdMargins = FALSE,
    205.                  cex.text = 0.5,
    206.                  zlim = c(-1,1),
    207.                  main = paste("Module-trait relationships"))
    208.   dev.off()
    209.   save(design,file = "step4_design.Rdata")
    210. }
    211.  
    212.  
    213.  
    214. if(T){
    215.   mes_group <- merge(MEs,datTraits,by="row.names") # 小技巧，可以通过rowname进行merge
    216.   ##写了个画图的function
    217.   library(gplots)
    218.   library(ggplot2)
    219.   library(ggpubr)
    220.   library(grid)
    221.   library(gridExtra) 
    222.   draw_ggboxplot <- function(data,gene="P53",group="group"){
    223.     #print(gene)
    224.     ggboxplot(data,x=group, y=gene,
    225.               ylab = sprintf("Expression of %s",gene),
    226.               xlab = group,
    227.               fill = group,
    228.               palette = "jco",
    229.               #add="jitter",
    230.               legend = "") +stat_compare_means()
    231.   }
    232.   ###开始批量画boxplot
    233.   colorNames = names(MEs)
    234.   png("figures/step4-expression-group.png",width = 800,height=2000,res = 120)
    235.   #par(mfrow=c(ceiling(length(colorNames)/2),2))
    236.   p <- lapply(colorNames,function(x) {
    237.     draw_ggboxplot(mes_group,gene= x,group="group")
    238.   })
    239.   do.call(grid.arrange,c(p,ncol=2))
    240.   dev.off()
    241. }
    242.  
    243.  
    244.  
    245. if(T){
    246.   modNames = substring(names(MEs), 3)
    247.   modNames
    248.   head(MEs)
    249.   head(datExpr)[1:6:1:9]
    250.   geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"))
    251.   ## 算出每个模块跟基因的皮尔森相关系数矩阵
    252.   ## MEs是每个模块在每个样本里面的值
    253.   ## datExpr是每个基因在每个样本的表达量
    254.   MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples))
    255.   names(geneModuleMembership) = paste("MM", modNames, sep="")
    256.   names(MMPvalue) = paste("p.MM", modNames, sep="")
    257.   head(MMPvalue)
    258.   
    259.   geneTraitSignificance <- as.data.frame(cor(datExpr,datTraits$groupNo,use = "p"))
    260.   head(geneTraitSignificance)
    261.   GSPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance),nSamples))
    262.   head(GSPvalue)
    263.   names(geneTraitSignificance)<- paste("GS.",names(datTraits$group),sep = "")
    264.   names(GSPvalue)<-paste("GS.",names(datTraits$group),sep = "")
    265.   head(GSPvalue)
    266.   
    267.   #selectModule<-c("blue","green","purple","grey")  ##可以选择自己喜欢的模块
    268.   selectModule <- modNames  ## 也可以批量作图
    269.   png("figures/step4-Module-trait-significance.png",width = 800,height=2000,res = 120)
    270.   par(mfrow=c(ceiling(length(selectModule)/2),2)) #批量作图开始
    271.   for(module in selectModule){
    272.     column <- match(module,selectModule)
    273.     print(module)
    274.     moduleGenes <- moduleColors==module
    275.     verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
    276.                        abs(geneTraitSignificance[moduleGenes, 1]),
    277.                        xlab = paste("Module Membership in", module, "module"),
    278.                        ylab = paste("Gene significance for", module, "module"),
    279.                        main = paste("Module membership vs. gene significance\n"),
    280.                        cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
    281.   }
    282.   dev.off()
    283. }
    284.  
    285.  
    286.  
    287.  
    288.  
    289.  
    290.  
    291.  
    292.  
    293.  
    294.  
    295.  
    296.  
    297.  
    298.  
    299.  
    300.

  • 相关阅读:
    （原创）[C#] GDI+ 之鼠标交互：原理、示例、一步步深入、性能优化
    Greenplum外表gpfdist加载数据
    6.6 - Windows网络相关命令
    最长回文串-力扣409-C++
    基于YOLOv8/YOLOv7/YOLOv6/YOLOv5的零售柜商品检测软件（Python+PySide6界面+训练代码）
    【Java】继承练习
    Windows上运行Redis
    vue如何使用深度选择器？
    想做扫码看图效果，你需要学会这一招
    【Flutter】使用Android Studio 创建第一个flutter应用。
  • 原文地址：https://blog.csdn.net/qq_52813185/article/details/127681071