GEO数据挖掘系列，第20篇学习笔记：WGCNA代码实操之准备阶段
0.数据预整理
本例用GEO数据库下载的数据，GSE199335
rm(list = ls())
library(tinyarray)
gse = "GSE199335"
geo = geo_download(gse,destdir=".",colon_remove = T)
geo$gpl
做WGCNA必须有15及以上的样本，否则不符合要求，没有做的意义
图0
#View(geo$pd)
library(stringr)
Group = paste(geo$pd$genotype,geo$pd$age,sep="_") %>%
str_remove(" months of age| weeks of age") %>%
str_remove(" type") %>%
str_replace("/",".")
table(Group)
Group = factor(Group,levels = c("wild_6", "wild_9", "R6.1_6", "R6.2_9"))
ids $gpl,destdir = tempdir(),type = "soft")
ids = na.omit(ids)#将NA行去除
exp = trans_array(geo$exp,ids,from = "ID")#把探针表达矩阵转为基因表达矩阵
pd = geo$pd
Group
save(exp,Group,pd,file = "Dat.Rdata")
以上整理数据的代码，尤其是分组信息和探针注释，因数据而异。请观察本例的结果，当遇到其他数据时，需要变动以上代码使结果和本例的尽可能相同，才能便于后续代码的运行
后续代码主要来自WGCNA的官方文档，部分进行了一些改动
1.表达矩阵数据整理
①需要有至少15个样本
②需要行是样本，列是基因，所以需要转置
③不推荐使用全部基因，因为计算量太大
④不推荐使用差异分析所得的差异基因，包作者说的
⑤比较推荐的方法是按照方差或者mad去取前多少个基因，例如3500，5000，8000个，或者保留基因总数的前1/4。无标准答案（有选择困难的可以选5000，出现问题后再试试其他数值）
rm(list = ls())
library(WGCNA)
library(tinyarray)
load("Dat.Rdata")
#exp = log2(geo$exp+1)
png("1.exp.png", width = 2000, height = 1200,res = 300)
boxplot(exp)#检查是否有离群的样本
dev.off()
从图1可以看出数据是取过log的，且没有异常样本，可以用
图1 exp
#方法一：前5000个
datExpr0 = t(exp[order(apply(exp, 1, mad), decreasing = T)[1:5000],])
#方法二：基因总数的前1/4
# datExpr0 = t(exp[order(apply(exp, 1, var), decreasing = T)[1:round(0.25*nrow(exp))],])
datExpr0[1:4,1:4]
1.1基因过滤
主要看缺失值，GEO的芯片数据大多数没缺失值
gsg = goodSamplesGenes(datExpr0, verbose = 3)
gsg$allOK # 返回TRUE则继续，若不是TRUE，以下代码将进行处理
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]
}
1.2样本过滤
sampleTree = hclust(dist(datExpr0), method = "average")
png(file = "2.sampleClustering.png", width = 2000, height = 2000,res = 300)
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)
dev.off()
看有没有自己单独一个分支的样本，如果有异常值就需要去掉
图2 sampleClustering
# 如果有离群值，以下代码将去掉该离群值
# 根据聚类图自己设置 cutHeight 参数的值
if(F){
clust = cutreeStatic(sampleTree, cutHeight = 170, minSize = 10)
table(clust) # 0代表切除的，1代表保留的
keepSamples = (clust!=0)
datExpr = datExpr0[keepSamples, ]
}
#没有异常样本就不需要去除
datExpr = datExpr0
2.表型信息整理
即临床信息的整理，这个信息来自芯片数据的pData，也就是上面的pd
