• GEO生信数据挖掘（九）肺结核数据-差异分析-WGCNA分析（900行代码整理注释更新版本）

    第六节，我们使用结核病基因数据，做了一个数据预处理的实操案例。例子中结核类型，包括结核，潜隐进展，对照和潜隐，四个类别。第七节延续上个数据，进行了差异分析。 第八节对差异基因进行富集分析。本节进行WGCNA分析。

    WGCNA分析 分段代码（附运行效果图）请查看上节

    运行后效果

    1. rm(list = ls()) ######清除环境数据
    2.  
    3.  
    4. #============================================================================
    5. #======================================================================
    6. #+========step0数据预处理和检查,已经做过step0==========================
    7. #+========================================
    8. #+=============================
    9. """
    10. ##############设置工作路径###################
    11. workingDir = "C:/Users/Desktop/GSE152532"############工作路径，可以修改，可以设置为数据存放路径
    12. setwd(workingDir)
    13. getwd()
    15. ################载入R包########################
    16. library(WGCNA)
    17. library(data.table)
    19. #############################导入数据##########################
    20. # The following setting is important, do not omit.
    21. options(stringsAsFactors = FALSE)
    22. #Read in the female liver data set
    23. fpkm = fread("Gene_expression.csv",header=T)##############数据文件名，根据实际修改，如果工作路径不是实际数据路径，需要添加正确的数据路径
    24. # Take a quick look at what is in the data set
    25. dim(fpkm)
    26. names(fpkm)
    28. ####################导入平台数据##########################
    29. library(idmap2)
    30. ids=get_soft_IDs('GPL10558')
    31. head(ids)
    33. #####################将探针ID改为基因ID##########################
    34. fpkm<-merge(fpkm,ids,by='ID')#merge()函数将dat1的探针id与芯片平台探针id相匹配，合并到dat1
    35. library(limma)
    36. fpkm<-avereps(fpkm[,-c(1,99)],ID=fpkm$symbol)#多个探针检测一个基因，合并一起，取其平均值
    37. fpkm<-as.data.frame(fpkm)#将矩阵转换为表格
    38. write.table(fpkm, file="FPKM_genesymbol.csv",row.names=T, col.names=T,quote=FALSE,sep=",")
    39. ###结束后查看文件，进行修改！！！
    44. # 加载自己的数据
    46. # load( "group_data_TB_LTBI.Rdata")
    48. load("exprSet_clean_mean_filter_log1.RData")  #exprSet_clean
    50. load( "dataset_TB_LTBI.Rdata")
    51. exprSet_clean = dataset_TB_LTBI
    52. gene_var <- apply(exprSet_clean, 1, var)##### 计算基因的方差
    53. keep_genes <- gene_var >= quantile(gene_var, 0.75)##### 筛选方差较大的基因,选择方差前25%的基因
    54. exprSet_clean <- exprSet_clean[keep_genes,]##### 保留筛选后的基因
    55. dim(exprSet_clean)
    56. save (exprSet_clean,file="方差前25per_TB_LTBI.Rdata")
    60. #######################基于方差筛选基因#################################
    61. fpkm_var <- read.csv("FPKM_genesymbol.csv", header = TRUE, row.names = 1)##### 读入表达矩阵,矩阵的行是基因，列是样本
    62. gene_var <- apply(fpkm_var, 1, var)##### 计算基因的方差
    63. keep_genes <- gene_var >= quantile(gene_var, 0.75)##### 筛选方差较大的基因,选择方差前25%的基因
    64. fpkm_var <- fpkm_var[keep_genes,]##### 保留筛选后的基因
    66. write.table(fpkm_var, file="FPKM_var.csv",row.names=T, col.names=T,quote=FALSE,sep=",")
    67. ###结束后查看文件，进行修改！！！
    69. ##################重新载入数据########################
    70. # The following setting is important, do not omit.
    71. options(stringsAsFactors = FALSE)
    72. #Read in the female liver data set
    73. fpkm = fread("FPKM_var_filter.csv",header=T)##############数据文件名，根据实际修改，如果工作路径不是实际数据路径，需要添加正确的数据路径
    74. # Take a quick look at what is in the data set
    75. dim(fpkm)
    76. names(fpkm)
    78. datExpr0 = as.data.frame(t(fpkm[,-1]))
    79. names(datExpr0) = fpkm$ID;##########如果第一行是ID命名，就写成fpkm$ID
    80. rownames(datExpr0) = names(fpkm[,-1])
    85. ##################check missing value and filter ####################
    86. load("方差前25per_TB_LTBI.Rdata")
    87. datExpr0  = exprSet_clean
    89. ##check missing value
    90. library(WGCNA)
    91. gsg = goodSamplesGenes(datExpr0, verbose = 3)
    92. gsg$allOK
    94. if (!gsg$allOK)
    95. {
    96.   # Optionally, print the gene and sample names that were removed:
    97.   if (sum(!gsg$goodGenes)>0)
    98.     printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ", ")))
    99.   if (sum(!gsg$goodSamples)>0)
    100.     printFlush(paste("Removing samples:", paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", ")))
    101.   # Remove the offending genes and samples from the data:
    102.   datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
    103. }
    105. ##filter
    106. #meanFPKM=0.5  ####过滤标准，可以修改
    107. #n=nrow(datExpr0)
    108. #datExpr0[n+1,]=apply(datExpr0[c(1:nrow(datExpr0)),],2,mean)
    109. #datExpr0=datExpr0[1:n,datExpr0[n+1,] > meanFPKM]  # for meanFpkm in row n+1 and it must be above what you set--select meanFpkm>opt$meanFpkm(by rp)
    112. filtered_fpkm=t(datExpr0)  #行 样本 列 基因
    113. filtered_fpkm=data.frame(rownames(filtered_fpkm),filtered_fpkm)
    114. names(filtered_fpkm)[1]="sample"
    115. write.table(filtered_fpkm, file="FPKM_filter.csv",row.names=F, col.names=T,quote=FALSE,sep="\t")
    117. """
    118.  
    119.  
    120.  
    121.  
    122. #&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    123. #+&&&&&&&&&&&&&&&&&&加载数据&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    124. #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    125. load('DEG_TB_LTBI_step13.Rdata')  # DEG,res,all_diff,limma_clean_res,dataset_TB_LTBI_DEG,
    126. #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    127. #+&&&&&&&&&&&&&&&&&&加载数据&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    128. #&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    129.  
    130.  
    131.  
    132.  
    133. library(WGCNA)
    134. #读取目录名称，方便复制粘贴
    135. dir()
    136.  
    137. #============================================================================
    138. #======================================================================
    139. #+========step1样品聚类step1=================================
    140. #+========================================
    141. #+=============================
    142.  
    143.  
    144. ################################样品聚类#################### 
    145. #这里行是样品名，列为基因名，做转置处理
    146. datExpr = t(dataset_TB_LTBI_DEG)
    147. #初次聚类
    148. sampleTree = hclust(dist(datExpr), method = "average")
    149. # Plot the sample tree: Open a graphic output window of size 20 by 15 inches
    150. # The user should change the dimensions if the window is too large or too small.
    151. #设置绘图窗口
    152. sizeGrWindow(12,9)
    153. pdf(file='1_sampleCluestering_初次聚类检查偏离样本.pdf',width = 12,height = 9)
    154. par(cex=0.6)
    155. par(mar=c(0,4,2,0))
    156. plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
    157.      cex.axis = 1.5, cex.main = 2)
    158.  
    159.  
    160. dev.off()
    161.  
    162.  
    163. #============================================================================
    164. #======================================================================
    165. #+========step2切割离群样本=================================
    166. #+========================================
    167. #+=============================
    168.  
    169.  
    170. pdf(file='2_sampleCluestering_初次聚类进行切割删除样本.pdf',width = 12,height = 9)
    171. par(cex=0.6)
    172. par(mar=c(0,4,2,0))
    173. plot(sampleTree, main = "Sample clustering to detect outliers ", sub="", xlab="", cex.lab = 1.5,
    174.      cex.axis = 1.5, cex.main = 2)
    175.  
    176. ### 测试画线，可以多次尝试
    177. ##############剪切高度问题,这个根据实际设置后可用
    178. abline(h = 87, col = "red")##剪切高度不确定，故无红线
    179. dev.off()
    180.  
    181. ### Determine cluster under the line
    182. clust = cutreeStatic(sampleTree, cutHeight = 87, minSize = 10)
    183. table(clust)
    184. #clust
    185. #0  1  2 
    186. #5 57 40
    187. #由于本人案例，一刀切出三段，需要保留两段，用了’或‘逻辑运算符号
    188. ### 需要保留哪个，就传如要保留clust编号
    189. keepSamples = (clust==1|clust==2)
    190. #剔除离群样本
    191. datExpr0 = datExpr[keepSamples, ]
    192. #观察新表达矩阵
    193. dim(datExpr0) #[1]   97 2813
    194.  
    195.  
    196. #&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    197. #+&&&&&&&&&&&&&&&&&&数据保存&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    198. #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    199. save(datExpr0,file='3.聚类后剔除离群样本datExpr0.Rdata')#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    200. #+&&&&&&&&&&&&&&&&&&数据保存&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    201. #&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    202.  
    203.  
    204. load('datExpr0_cluster_filter.Rdata')
    205.  
    206.  
    207. #============================================================================
    208. #======================================================================
    209. #+========step3临床性状数据整理，与新表达矩阵保持一致=================================
    210. #+========================================
    211. #+=============================
    212.  
    213.  
    214. #加载自己的临床性状数据
    215. load('design_TB_LTBI.Rdata')
    216. traitData=design
    217.  
    218. dim(traitData)
    219.  
    220. # Form a data frame analogous to expression data that will hold the clinical traits.
    221. fpkmSamples = rownames(datExpr0)
    222. traitSamples =rownames(traitData)
    223. #匹配样本名称，性状数据与表达数据保证一致
    224. traitRows = match(fpkmSamples, traitSamples)
    225. datTraits = traitData[traitRows,]
    226. rownames(datTraits) 
    227. collectGarbage()
    228.  
    229. #&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    230. #+&&&&&&&&&&&&&&&&&&数据保存&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    231. #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    232. save(datTraits,file='4.剔除离群样本的临床性状数据datTraits.Rdata')#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    233. #+&&&&&&&&&&&&&&&&&&数据保存&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    234. #&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    235.  
    236.  
    237.  
    238.  
    239.  
    240. #============================================================================
    241. #======================================================================
    242. #+========step4 增加临床性状数据后再次聚类=======================
    243. #+========================================
    244. #+=============================
    245. # Re-cluster samples
    246. sampleTree2 = hclust(dist(datExpr0), method = "average")
    247. # Convert traits to a color representation: white means low, red means high, grey means missing entry
    248. traitColors = numbers2colors(datTraits, signed = FALSE)
    249. # Plot the sample dendrogram and the colors underneath.
    250.  
    251. #sizeGrWindow(20,20)
    252.  
    253. pdf(file="5_Sample cluster dendrogram and trait heatmap.pdf",width=12,height=12)
    254. plotDendroAndColors(sampleTree2, traitColors,
    255.                     groupLabels = names(datTraits),
    256.                     main = "Sample dendrogram and trait heatmap")
    257.  
    258. #Error in .plotOrderedColorSubplot(order = order, colors = colors, rowLabels = rowLabels,  : 
    259. #                                    Length of colors vector not compatible with number of objects in 'order'.
    260.                                   
    261. dev.off()
    262.  
    263.  
    264.  
    265.  
    266. #============================================================================
    267. #======================================================================
    268. #+========step5 构建WGCNA网络=======================
    269. #+========================================
    270. #+=============================
    271.  
    272. # Allow multi-threading within WGCNA. At present this call is necessary.
    273. # Any error here may be ignored but you may want to update WGCNA if you see one.
    274. # Caution: skip this line if you run RStudio or other third-party R environments.
    275. # See note above.
    276. #检查环境，能开几个线程
    277. enableWGCNAThreads()
    278.  
    279. # Choose a set of soft-thresholding powers
    280. #设置阈值范围，WGCNA是无标度网络（scale free）,节点连结数服从幂次定律分布。（连接数越多核心节点越少）
    281. powers = c(1:15)
    282.  
    283. # Call the network topology analysis function
    284. #网络拓扑分析
    285. sft = pickSoftThreshold(datExpr0, powerVector = powers, verbose = 5)
    286.  
    287. # Plot the results:
    288. sizeGrWindow(15, 9)
    289. pdf(file="6_Scale independence选阈值测试过程.pdf",width=9,height=5)
    290. #pdf(file="Rplot03.pdf",width=9,height=5)
    291. par(mfrow = c(1,2))
    292. cex1 = 0.9
    293. # Scale-free topology fit index as a function of the soft-thresholding power
    294. #无标度拓扑拟合指标作为软阈值能力的函数，根据下图结果，挑选合适阈值
    295. plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
    296.      xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
    297.      main = paste("Scale independence"));
    298. text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
    299.      labels=powers,cex=cex1,col="red");
    300. # this line corresponds to using an R^2 cut-off of h
    301. abline(h=0.90,col="red")
    302. # Mean connectivity as a function of the soft-thresholding power
    303. plot(sft$fitIndices[,1], sft$fitIndices[,5],
    304.      xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
    305.      main = paste("Mean connectivity"))
    306. text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")
    307. dev.off()
    308.  
    309.  
    310.  
    311. ######chose the softPower
    312. #选择阈值
    313. softPower =sft$powerEstimate
    314. adjacency = adjacency(datExpr0, power = softPower)
    315.  
    316. ##### Turn adjacency into topological overlap
    317. #将邻接转换为拓扑重叠
    318. TOM = TOMsimilarity(adjacency);
    319. dissTOM = 1-TOM
    320.  
    321. # Call the hierarchical clustering function
    322. #无标度网络阈值参数确定后，调用分层聚类函数
    323. #基于TOM的不相似性基因聚类
    324. geneTree = hclust(as.dist(dissTOM), method = "average");
    325. # Plot the resulting clustering tree (dendrogram)
    326.  
    327. #sizeGrWindow(12,9)
    328. pdf(file="7_Gene clustering on TOM-based dissimilarity基因聚类图.pdf",width=12,height=9)
    329. plot(geneTree, xlab="", sub="", main = "Gene clustering on TOM-based dissimilarity",
    330.      labels = FALSE, hang = 0.04)
    331. dev.off()
    332.  
    333.  
    334.  
    335.  
    336. #聚类模块，最小的基因数量
    337. # We like large modules, so we set the minimum module size relatively high:
    338. minModuleSize = 30
    339.  
    340. # Module identification using dynamic tree cut:
    341. #使用dynamic tree cut进行模块识别
    342. dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM,
    343.                             deepSplit = 2, pamRespectsDendro = FALSE,
    344.                             minClusterSize = minModuleSize);
    345. table(dynamicMods)
    346.  
    347. # Convert numeric lables into colors
    348. #给不同模块分配颜色
    349. dynamicColors = labels2colors(dynamicMods)
    350. table(dynamicColors)
    351. # Plot the dendrogram and colors underneath
    352. #sizeGrWindow(8,6)
    353. pdf(file="8_带颜色标识的聚类树Dynamic Tree Cut.pdf",width=8,height=6)
    354. plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut",
    355.                     dendroLabels = FALSE, hang = 0.03,
    356.                     addGuide = TRUE, guideHang = 0.05,
    357.                     main = "Gene dendrogram and module colors")
    358. dev.off()
    359.  
    360.  
    361. # Calculate eigengenes
    362. MEList = moduleEigengenes(datExpr0, colors = dynamicColors)
    363. MEs = MEList$eigengenes
    364. # Calculate dissimilarity of module eigengenes
    365. MEDiss = 1-cor(MEs);
    366. # Cluster module eigengenes
    367. METree = hclust(as.dist(MEDiss), method = "average")
    368. # Plot the result
    369. #sizeGrWindow(7, 6)
    370. pdf(file="9_Clustering of module eigengenes.pdf",width=7,height=6)
    371. plot(METree, main = "Clustering of module eigengenes",
    372.      xlab = "", sub = "")
    373. MEDissThres = 0.25######剪切高度可修改
    374. # Plot the cut line into the dendrogram
    375. abline(h=MEDissThres, col = "red")
    376. dev.off()
    377.  
    378.  
    379.  
    380.  
    381. # Call an automatic merging function
    382. #根据前面设置的剪切高度，对模块进行合并
    383. merge = mergeCloseModules(datExpr0, dynamicColors, cutHeight = MEDissThres, verbose = 3)
    384. # The merged module colors
    385. mergedColors = merge$colors
    386. # Eigengenes of the new merged modules:
    387. mergedMEs = merge$newMEs
    388.  
    389. #sizeGrWindow(12, 9)
    390. pdf(file="10_合并模块后的聚类树merged dynamic.pdf", width = 9, height = 6)
    391. plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors),
    392.                     c("Dynamic Tree Cut", "Merged dynamic"),
    393.                     dendroLabels = FALSE, hang = 0.03,
    394.                     addGuide = TRUE, guideHang = 0.05)
    395. dev.off()
    396.  
    397. # Rename to moduleColors
    398. moduleColors = mergedColors
    399. # Construct numerical labels corresponding to the colors
    400. #构建相应颜色的数字标签
    401. colorOrder = c("grey", standardColors(50))
    402. moduleLabels = match(moduleColors, colorOrder)-1
    403. MEs = mergedMEs
    404.  
    405.  
    406.  
    407.  
    408. # Save module colors and labels for use in subsequent parts
    409. #&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    410. #+&&&&&&&&&&&&&&&&&&数据保存&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    411. #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    412. save(datExpr0,datTraits,MEs, TOM, dissTOM,  moduleLabels, moduleColors, geneTree, sft, file = "11_networkConstruction-stepByStep.RData")
    413. #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    414. #+&&&&&&&&&&&&&&&&&&数据保存&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    415. #&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    416.  
    417.  
    418.  
    419.  
    420.  
    421.  
    422. load("11_networkConstruction-stepByStep.RData")
    423.  
    424.  
    425.  
    426.  
    427.  
    428.  
    429.  
    430. #=====================================================================================
    431. #===============================================================================
    432. #+========step6 计算模块和临床性状相关系数（核心挑选色块）==============
    433. #+========================================
    434. #+=============================
    435. ##############################relate modules to external clinical triats
    436.  
    437. # Define numbers of genes and samples
    438. nGenes = ncol(datExpr0)
    439. nSamples = nrow(datExpr0)
    440.  
    441. #可以修改参数 p值pvalue 更换 
    442. moduleTraitCor = cor(MEs, datTraits, use = "p") 
    443. moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
    444.  
    445. #sizeGrWindow(10,6)
    446. pdf(file="12_模块和临床形状关系图Module-trait relationships.pdf",width=10,height=6)
    447. # Will display correlations and their p-values
    448. textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
    449.                    signif(moduleTraitPvalue, 1), ")", sep = "")
    450.  
    451. dim(textMatrix) = dim(moduleTraitCor)
    452. par(mar = c(6, 8.5, 3, 3))
    453.  
    454. # Display the correlation values within a heatmap plot #修改性状类型 data.frame
    455. labeledHeatmap(Matrix = moduleTraitCor,
    456.                xLabels = names(data.frame(datTraits)),
    457.                yLabels = names(MEs),
    458.                ySymbols = names(MEs),
    459.                colorLabels = FALSE,
    460.                colors = greenWhiteRed(50),
    461.                textMatrix = textMatrix,
    462.                setStdMargins = FALSE,
    463.                cex.text = 0.5,
    464.                zlim = c(-1,1),
    465.                main = paste("Module-trait relationships"))
    466. dev.off()
    467.  
    468. #色块 red相关度 0.75
    469.  
    470.  
    471.  
    472.  
    473.  
    474.  
    475.  
    476.  
    477. #=====================================================================================
    478. #===============================================================================
    479. #+========step7 定义包含所有datTraits列的可变权重（MM and GS）==============
    480. #+========================================
    481. #+=============================
    482.  
    483.  
    484. #定义包含所有datTraits列的可变权重
    485.  
    486. ######## Define variable weight containing all column of datTraits
    487.  
    488. ###MM（gene Module Membership） and GS（gene Trait Significance）
    489.  
    490.  
    491. # names (colors) of the modules
    492. modNames = substring(names(MEs), 3)
    493.  
    494. geneModuleMembership = as.data.frame(cor(datExpr0, MEs, use = "p"))
    495. MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples))
    496.  
    497. names(geneModuleMembership) = paste("MM", modNames, sep="")
    498. names(MMPvalue) = paste("p.MM", modNames, sep="")
    499.  
    500. #names of those trait
    501. traitNames=names(data.frame(datTraits))
    502. class(datTraits)
    503.  
    504. geneTraitSignificance = as.data.frame(cor(datExpr0, datTraits, use = "p"))
    505. GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples))
    506.  
    507. names(geneTraitSignificance) = paste("GS.", traitNames, sep="")
    508. names(GSPvalue) = paste("p.GS.", traitNames, sep="")
    509.  
    510.  
    511. ####plot MM vs GS for each trait vs each module
    512.  
    513.  
    514. ##########example:royalblue and CK
    515. module="red"
    516. column = match(module, modNames)
    517. moduleGenes = moduleColors==module
    518.  
    519. trait="TB"
    520. traitColumn=match(trait,traitNames)
    521.  
    522.  
    523. sizeGrWindow(7, 7)
    524.  
    525. #par(mfrow = c(1,1))
    526. verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
    527. abs(geneTraitSignificance[moduleGenes, traitColumn]),
    528. xlab = paste("Module Membership in", module, "module"),
    529. ylab = paste("Gene significance for ",trait),
    530. main = paste("Module membership vs. gene significance\n"),
    531. cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
    532. ######
    533.  
    534. for (trait in traitNames){
    535.   traitColumn=match(trait,traitNames)
    536.   
    537.   for (module in modNames){
    538.     column = match(module, modNames)
    539.     moduleGenes = moduleColors==module
    540.     
    541.     if (nrow(geneModuleMembership[moduleGenes,]) > 1){####进行这部分计算必须每个模块内基因数量大于2，由于前面设置了最小数量是30，这里可以不做这个判断，但是grey有可能会出现1个gene,它会导致代码运行的时候中断，故设置这一步
    542.       
    543.       #sizeGrWindow(7, 7)
    544.       pdf(file=paste("13_", trait, "_", module,"_Module membership vs gene significance.pdf",sep=""),width=7,height=7)
    545.       par(mfrow = c(1,1))
    546.       verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
    547.                          abs(geneTraitSignificance[moduleGenes, traitColumn]),
    548.                          xlab = paste("Module Membership in", module, "module"),
    549.                          ylab = paste("Gene significance for ",trait),
    550.                          main = paste("Module membership vs. gene significance\n"),
    551.                          cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
    552.       dev.off()
    553.     }
    554.   }
    555. }
    556.  
    557. #####
    558. names(datExpr0)
    559. probes = names(data.frame(datExpr0))
    560.  
    561. #=====================================================================================
    562. #===============================================================================
    563. #+========step8 导出计算完毕的（MM and GS）==============
    564. #+========================================
    565. #+=============================
    566. #################export GS and MM############### 
    567.  
    568. geneInfo0 = data.frame(probes= probes,
    569.                        moduleColor = moduleColors)
    570.  
    571. for (Tra in 1:ncol(geneTraitSignificance))
    572. {
    573.   oldNames = names(geneInfo0)
    574.   geneInfo0 = data.frame(geneInfo0, geneTraitSignificance[,Tra],
    575.                          GSPvalue[, Tra])
    576.   names(geneInfo0) = c(oldNames,names(geneTraitSignificance)[Tra],
    577.                        names(GSPvalue)[Tra])
    578. }
    579.  
    580. for (mod in 1:ncol(geneModuleMembership))
    581. {
    582.   oldNames = names(geneInfo0)
    583.   geneInfo0 = data.frame(geneInfo0, geneModuleMembership[,mod],
    584.                          MMPvalue[, mod])
    585.   names(geneInfo0) = c(oldNames,names(geneModuleMembership)[mod],
    586.                        names(MMPvalue)[mod])
    587. }
    588. geneOrder =order(geneInfo0$moduleColor)
    589. geneInfo = geneInfo0[geneOrder, ]
    590.  
    591. write.table(geneInfo, file = "14_GS_and_MM.xls",sep="\t",row.names=F)
    592.  
    593.  
    594. #=====================================================================================
    595. #===============================================================================
    596. #+========step9 基因网络热图进行可视化（非常耗费内存资源）==============
    597. #+========================================
    598. #+=============================
    599.  
    600. nGenes = ncol(datExpr0)
    601. nSamples = nrow(datExpr0)
    602.  
    603.  
    604. # Transform dissTOM with a power to make moderately strong connections more visible in the heatmap
    605. plotTOM = dissTOM^7
    606. # Set diagonal to NA for a nicer plot
    607. diag(plotTOM) = NA
    608.  
    609.  
    610.  
    611. # Call the plot function
    612.  
    613. sizeGrWindow(9,9)  #这个耗电脑内存
    614. pdf(file="15_所有基因数量太多Network heatmap plot_all gene.pdf",width=9, height=9)
    615. TOMplot(plotTOM, geneTree, moduleColors, main = "Network heatmap plot, all genes")
    616. dev.off()
    617.  
    618.  
    619. nSelect = 400
    620. # For reproducibility, we set the random seed
    621. set.seed(10)
    622. select = sample(nGenes, size = nSelect)
    623. selectTOM = dissTOM[select, select]
    624. # There's no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.
    625. selectTree = hclust(as.dist(selectTOM), method = "average")
    626. selectColors = moduleColors[select]
    628. # Open a graphical window
    629. #sizeGrWindow(9,9)
    630. # Taking the dissimilarity to a power, say 10, makes the plot more informative by effectively changing
    631. # the color palette; setting the diagonal to NA also improves the clarity of the plot
    632. plotDiss = selectTOM^7
    633. diag(plotDiss) = NA
    635. pdf(file="16_400个基因试试Network heatmap plot_selected genes.pdf",width=9, height=9)
    636. TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")
    637. dev.off()
    641. #=====================================================================================
    642. #===============================================================================
    643. #+========step10 新模块和临床性状热图 合并和拆分两个版本==============
    644. #+========================================
    645. #+=============================
    647. #sizeGrWindow(5,7.5)
    648. pdf(file="17新模块和临床性状热图_Eigengene dendrogram and Eigengene adjacency heatmap.pdf", width=5, height=7.5)
    649. par(cex = 0.9)
    650. plotEigengeneNetworks(MEs, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle= 90)
    651. dev.off()
    653. #or devide into two parts
    654. # Plot the dendrogram
    655. #sizeGrWindow(6,6);
    656. pdf(file="18_Eigengene dendrogram_2.pdf",width=6, height=6)
    657. par(cex = 1.0)
    658. plotEigengeneNetworks(MEs, "Eigengene dendrogram", marDendro = c(0,4,2,0), plotHeatmaps = FALSE)
    659. dev.off()
    661. pdf(file="19_Eigengene adjacency heatmap_2.pdf",width=6, height=6)
    662. # Plot the heatmap matrix (note: this plot will overwrite the dendrogram plot)
    663. par(cex = 1.0)
    664. plotEigengeneNetworks(MEs, "Eigengene adjacency heatmap", marHeatmap = c(3,4,2,2), plotDendrograms = FALSE, xLabelsAngle = 90)
    665. dev.off()
    670. ###########################Exporting to Cytoscape all one by one ##########################
    673. #=====================================================================================
    674. #===============================================================================
    675. #+========step11 导出每个模块的边和节点关系（Cytoscape 绘图所需）==============
    676. #+========================================
    677. #+=============================
    679. # Select each module
    680. '''
    681. Error in exportNetworkToCytoscape(modTOM, edgeFile = paste("CytoscapeInput-edges-",  : 
    682.   Cannot determine node names: nodeNames is NULL and adjMat has no dimnames.
    683.   datExpr0 格式需要dataframe
    684. '''
    685. modules =module  #module="red"
    686. for (mod in 1:nrow(table(moduleColors)))
    687. {
    688.   
    689.   modules = names(table(moduleColors))[mod]
    690.   # Select module probes
    691.   probes = names(data.frame(datExpr0))  # 
    692.   inModule = (moduleColors == modules)
    693.   modProbes = probes[inModule]
    694.   modGenes = modProbes
    695.   # Select the corresponding Topological Overlap
    696.   modTOM = TOM[inModule, inModule]
    697.   
    698.   dimnames(modTOM) = list(modProbes, modProbes)
    699.   # Export the network into edge and node list files Cytoscape can read
    700.   cyt = exportNetworkToCytoscape(modTOM,
    701.                                  edgeFile = paste("20_CytoscapeInput-edges-", modules , ".txt", sep=""),
    702.                                  nodeFile = paste("20_CytoscapeInput-nodes-", modules, ".txt", sep=""),
    703.                                  weighted = TRUE,
    704.                                  threshold = 0.02,
    705.                                  nodeNames = modProbes,
    706.                                  altNodeNames = modGenes,
    707.                                  nodeAttr = moduleColors[inModule])
    708. }

    WGCNA关系网络的构建完毕，绘图找核心基因，Cytoscape 到底怎么玩？

  • 相关阅读:
    AT800(3000) +昇腾300V 之 第一个例子图片分类
    平移矩阵中的数学思考
    hyperf 十六 session
    为什么PDF打开没有密码，但是不能编辑？
    力扣511. 游戏玩法分析 I
    密码保护工具的编写
    b2b2c o2o 多商家入驻商城 直播带货商城 电子商务
    mysql编译安装教程
    Linux驱动开发：掌握SPI通信机制
    Linux系统无法连接网络的情况
  • 原文地址：https://blog.csdn.net/zzh1464501547/article/details/133915414