WGCNA的使用
WGCNA算是一种很经典的分析方法 。可以将临床特征与基因联系起来 ， 但是怎么联系的呢 。这就涉及了相关性分析 ， 也就是聚类 。这里很重要 ， 要将聚类理解成相关性分析 ， 那么生信学起来就会容易很多 。
后面有空我会写一篇关于聚类的文章 ， 就会发现其实生信就是在聚类成的一类一类中去摸索我们想要的 。
【WGCNA的使用】WGCNA将临床特征与基因联系起来的过程中 ， 引入了一个新的东西 ， 叫做模块 。这个模块也就是将基因聚类
那这样我们就有几种研究方法呢 。这就是简单的排列组合而已了 。
可以去研究基因与模块的关系 ， 那就是把基因聚类成模块 。这部分做了吗？做了的 ， 也就是使用软阈值构建模块可以去研究模块与临床特征的关系 。也是做了的 。也就是将模块与临床特征的关系以热图的形式展示出来 ， 而且哈 ， 去研究这两者的关系运用的是cor函数 ， 这个函数特别重要 ， 研究的是皮尔森相关性 。可以去研究基因与临床特征的关系 ， 做了吗 ， 也是做了的 。运用的也是cor函数 。这些研究完还有别的吗 ， 还有！就是将三者放在一起研究 。怎么研究呢？ 可以将基因与模块的关系矩阵与模块与临床特征的关系矩阵放在一起可以将基因与模块的关系矩阵与基因与临床特征的关系矩阵放在一起 ， 这个研究的多一点 ， 因为涉及到基因 。可以将模块与临床特征的关系矩阵与基因与临床特征的关系矩阵放在一起 。
回头看一下 ， 是不是已经有6种研究方法了 ， 其实还不止 ， 稍微将一个变量变一下 ， 又有几种方法出来 。应该反复去阅读思考上面的内容 。
了解完WGCNA能做的事情之后那学起来就会容易很多 。
下面一步一步讲
1.处理数据 。下面是详细的代码
rm(list = ls())library(WGCNA)# 读取基因表达矩阵数据fpkm = read.table("data/fpkm.txt", header = T, row.names = 1, check.names = F)head(fpkm)### 选取基因方法 ##### 第一种 ， 通过标准差选择## 计算每个基因的标准差fpkm_sd = apply(fpkm,1,sd)#1是对每一行 ， 2是对每一列## 使用标准差对基因进行降序排序fpkm_sd_sorted = order(fpkm_sd, decreasing = T)## 选择前5000个标准差较大的基因fpkm_num = fpkm_sd_sorted[1:5000]## 从表达矩阵提取基因fpkm_filter = fpkm[fpkm_num,]## 对表达矩阵进行转置WGCNA_matrix = t(fpkm_filter)#变成行名是样本 ， 列名是基因## 保存过滤后的数据save(WGCNA_matrix, file = "data/Step01-fpkm_sd_filter.Rda")## 第二种 ， 使用绝对中位差选择 ， 推荐使用绝对中位差WGCNA_matrix = t(fpkm[order(apply(fpkm,1,mad), decreasing = T)[1:5000],])#mad代表绝对中位差save(WGCNA_matrix, file = "data/Step01-fpkm_mad_filter.Rda")## 第三种 ， 使用全基因WGCNA_matrix = t(fpkm)save(WGCNA_matrix, file = "data/Step01-fpkm_allgene.Rda")### 去除缺失值较多的基因/样本 ###rm(list = ls())# 加载表达矩阵load(file = "data/Step01-fpkm_mad_filter.Rda")# 加载WGCNA包library(WGCNA)# 判断是否缺失较多的样本和基因datExpr0 = WGCNA_matrixgsg = goodSamplesGenes(datExpr0, verbose = 3)# 是否需要过滤 ， TRUE表示不需要 ， FALSE表示需要gsg$allOK# 当gsg$allOK为TRUE ， 以下代码不会运行 ， 为FALSE时 ， 运行以下代码过滤样本if (!gsg$allOK){# 打印移除的基因名和样本名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 = ", ")));# 提取保留的基因和样本datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]}### 通过样本聚类识别离群样本 ， 去除离群样本 ###sampleTree = hclust(dist(datExpr0), method = "average");#使用hclust函数进行均值聚类# 绘制样本聚类图确定离群样本sizeGrWindow(30,9)pdf(file = "figures/Step01-sampleClustering.pdf", width = 30, 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)# 根据上图判断 ， 需要截取的高度参数habline(h = 120, col = "red")#在120的地方画条线dev.off()# 去除离群得聚类样本 ， cutHeight参数要与上述得h参数值一致clust = cutreeStatic(sampleTree, cutHeight = 120, minSize = 10)table(clust)# clust# 010就是要去除的 ， 1就是要保存的# 15 162# clust 1聚类中包含着我们想要得样本 ， 将其提取出来keepSamples = (clust==1)datExpr = datExpr0[keepSamples, ]# 记录基因和样本数 ， 方便后续可视化nGenes = ncol(datExpr)#基因数nSamples = nrow(datExpr)#样本数save(datExpr, nGenes, nSamples,file = "data/Step01-WGCNA_input.Rda")