WGCNA：加权基因共表达网络分析，简而言之，就是将基因划分为若干个模块，探究与表型数据与基因模块之间的相关关系，并找到模块中的核心基因。
适用于复杂的数据模式，推荐5组(或者15个样品)以上的数据。一般可应用的研究方向有：不同器官或组织类型发育调控、同一组织不同发育调控、非生物胁迫不同时间点应答、病原菌侵染后不同时间点应答。
更多原理还有专业术语可以看生信技能树的帖子：https://cloud.tencent.com/developer/article/1516749
1.R包下载
BiocManager::install("WGCNA")
2.数据导入和清洗
这是任何网络分析的第一步。在这里展示了如何加载典型的表达数据，将其预处理成适合网络分析的格式，并通过去除明显的离群样本以及缺失条目数量过多的基因和样本来清洗数据。
2.1 导入表达数据
#设置工作目录
> setwd("D:/RNA-seq/WGCNA/fpkm")
#载入WGCNA包
> library('WGCNA')
#这个设置很重要，不要忽略，允许R语言以最大线程运行
> options(stringsAsFactors = FALSE)
> allowWGCNAThreads()
Allowing multi-threading with up to 4 threads.
#导入表达数据
> fpkm <- read.csv(file="fpkm.csv")
#简单查看数据
> dim(fpkm)
[1] 138749 61
请注意，每一行对应一个基因，而每一列对应一个样本或辅助信息。需要移除辅助信息，因为WGCNA针对的是基因进行聚类，而一般我们的聚类是针对样本，所以这个时候需对表达数据进行转置以进行进一步分析。
因为我的基因数目比较多，所以在这里选取差异最大的前30%进行模块构建。
> rownames(fpkm) <- fpkm[,1]
> fpkm<- fpkm[,-1]
> WGCNA_matrix = t(fpkm[order(apply(fpkm,1,mad), decreasing = T)[1:41625],]) #按mad进行排序，取前30%的基因进行后续分析
> datExpr0 = as.data.frame(WGCNA_matrix)
#全部基因
> datExpr_all = as.data.frame(t(fpkm))
2.2 检查缺失值和离群样本
如果是全部基因，那首先检查基因和样本是否有过多的缺失值。如果是Top xx%选择后的基因，那通常 gsg$allOK会返回TRUE，因为前期已经进行了选择。
#全部基因
> gsg_all = goodSamplesGenes(datExpr_all, verbose = 3)
> gsg_all$allOK
[1] FALSE
#筛选后基因
> gsg = goodSamplesGenes(datExpr0, verbose = 3)
> gsg$allOK
[1] TRUE
如果返回值是TURE，证明所有的基因都符合要求，如果返回值是FALSE，需要去除离群值。
> if (!gsg_all$allOK)
+ {
+ # Optionally, print the gene and sample names that were removed:
+ if (sum(!gsg_all$goodGenes)>0)
+ printFlush(paste("Removing genes:", paste(names(datExpr_all)[!gsg_all$goodGenes], collapse = ", ")));
+ if (sum(!gsg_all$goodSamples)>0)
+ printFlush(paste("Removing samples:", paste(rownames(datExpr_all)[!gsg_all$goodSamples], collapse = ", ")));
+ # Remove the offending genes and samples from the data:
+ datExpr_all = datExpr_all[gsg_all$goodSamples, gsg_all$goodGenes]
+ }
#查看新的datExpr_all
> dim(datExpr_all)
[1] 60 122752
#可以看到基因数目变少
再次检测
> gsg_all = goodSamplesGenes(datExpr_all, verbose = 3)
> gsg_all$allOK
[1] TRUE
#这次返回值是TRUE，可以继续后续计算
由于基因数目太多，后续计算都是使用datExpr0数据集。
对样本进行聚类
> sampleTree = hclust(dist(datExpr0), method = "average")
> sizeGrWindow(12,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)
这里的样本有明显离群值，可以通过设定阈值去掉离群值。
Fig.1 Clustering dendrogram of samples based on their Euclidean distance
# 画一条阈值线，根据实际情况人工设定，这里是15
> abline(h = 15, col = "red")
# 确定阈值线下的集群
> clust = cutreeStatic(sampleTree, cutHeight = 15, minSize = 10)
> table(clust)
# clust 1包含想要留下的样本
> keepSamples = (clust==1)
> datExpr = datExpr0[keepSamples, ]
现在变量datExpr包含用于网络分析的表达量数据
2.3 导入表型数据
> datTraits <- read.csv("phe_blue.csv",header=TRUE)
> dim(datTraits)
[1] 60 5
#去掉第一列样本名
> rownames(datTraits) <- datTraits[,1]
> datTraits <- datTraits[,-1]
> dim(datTraits)
[1] 60 4
现在表型数据和基因表达量数据都已经成功导入，在继续网络构建之前，可以可视化样本和表型之间的关系。
# 将表型值转换为颜色表示:白色表示低，红色表示高，灰色表示缺失
> traitColors = numbers2colors(datTraits, signed = FALSE)
#在样本聚类图的基础上，增加表型值热图
> plotDendroAndColors(sampleTree, traitColors,
+ groupLabels = names(datTraits),
+ main = "Sample dendrogram and trait heatmap")
Figure 2: Clustering dendrogram of samples based on their Euclidean distance.
2.4 保存数据
最后一步是保存处理好的表达量数据和表型数据，以便在本教程的后续步骤中使用。
save(datExpr, datTraits, file = "WGCNA0.3-dataInput.RData")
网友评论