WGCNA新手入门笔记2(含代码和数据)

WGCNA新手入门笔记2（含代码和数据）

上次我们介绍了WGCNA的入门（WGCNA新手入门笔记（含代码和数据）），大家在安装WGCNA包的时候，可能会遇到GO.db这个包安装不了的问题。主要问题应该是出在电脑的防火墙，安装时请关闭防火墙。

如果还有问题，请先单独安装AnnotationDbi这个包，biocLite("AnnotationDbi")

再安装GO.db，并尝试从本地文件安装该包。

如果还有问题，请使用管理员身份运行R语言，尝试上述步骤。

另外如果大家问题解决了请在留言处留个言，告知大家是在哪一步解决了问题，谢谢！因为本人没有进行单因素实验，不知道到底是哪个因素改变了实验结果。。。

今天给大家过一遍代码。网盘中有代码和数据。

链接：/s/1bpvu9Dt

密码：w7g4

##导入数据##

library(WGCNA)options(stringsAsFactors =

FALSE)enableWGCNAThreads()

enableWGCNAThreads()是指允许R语言程序最大线程运行，像我这电脑是4核CPU的，那么就能用上3核：

当然如果当前电脑没别的事，也可以满负荷运作

samples=read.csv( 'Sam_info.txt',sep = 't',s =

1)expro=read.csv( 'ExpData.txt',sep = 't',s =

1)dim(expro)

这部分代码是为了让R语言读取外部数据。当然了在读取数据之前首先改变一下工作目录，这一点在周二的文章中提过了。R语言读取外部数据的方式常用的有read.table和read.csv，这里用的是read.csv，想要查看某一函数的具体参数，可以用？函数名查看，比如：

大家可以注意到read.table和read.csv中header参数的默认值是不同的，header=true表示第一行是标题，第二行才是数据，header=false则表示第一行就是数据，没有标题。##筛选方差前25%的基因##

m.vars=apply(expro,

1,var)expro.upper=expro[which(m.vars>quantile(m.vars, probs = seq( 0, 1, 0.25))[ 4]),]dim(expro.upper)datExpr= as.data.frame(t(expro.upper));nGenes =

ncol(datExpr)nSamples = nrow(datExpr)

这一步是为了减少运算量，因为一个测序数据可能会有好几

万个探针，而可能其中很多基因在各个样本中的表达情况并没有什么太大变化，为了减少运算量，这里我们筛选方差前25%的基因。

当然了，以下几种情况，你可以忽略此步，转而执行下列代码

datExpr= as.data.frame(t(expro));nGenes =

ncol(datExpr)nSamples = nrow(datExpr)

情况1：所用的数据是比较老的芯片数据，探针数量较少；情况2：你的电脑足够强大，不必减少运算；

情况3：你第一步导入的数据是差异基因的数据，已经作过初筛，比如这一文章的套路（PMID：27133569）（个人不建议这么做）。

##样本聚类检查离群值##

gsg = goodSamplesGenes(datExpr, verbose =

3);gsg$allOKsampleTree = hclust(dist(datExpr), method = "average")plot(sampleTree, main = "Sample clustering to detect outliers", sub= "", xlab= "")save(datExpr, file = "FPKM-01-dataInput.RData")

执行这一段代码我们会得到下面这张图：

从结果上来，我们分析的样本没啥离群值，所以代码里说不

作处理。在网上的一个案例中，离群的样本就比较明显了。

（https://www.shengxin.ren/article/88）

如果需要去除离群样本，则执行下列代码，其中cutHeight=多少就看你自己了。

clust = cutreeStatic(sampleTree, cutHeight = 20000, minSize = 10)table(clust)keepSamples = (clust==

1)datExpr = datExpr[keepSamples, ]nGenes =

ncol(datExpr)nSamples = nrow(datExpr)save(datExpr, file = "FPKM-01-dataInput.RData")

执行上述代码的话，就会去掉8个样本

##软阈值筛选##

powers = c(c( 1: 10), seq( from= 12, to= 20, by= 2))sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)par(mfrow = c( 1, 2));cex1 =

0.9;plot(sft$fitIndices[, 1], -sign(sft$fitIndices[,

3])*sft$fitIndices[, 2], xlab= "Soft Threshold (power)",ylab= "Scale Free Topology Model Fit,signed R^2",type= "n", main = paste( "Scale independence"));text(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3])*sft$fitIndices[, 2],

labels=powers,cex=cex1,col= "red");abline(h= 0.90,col=