加权基因共表达网络构建
- 安装WGCNA包
- 准备基因表达矩阵
- 离群样本检测
安装WGCNA包
首先最好下载R4.1.3版本,然后按照这位博主的的方法安装好WGCNA包:WGCNA包安装 当你输入library(WGCNA)没有提示错误时,说明安装成功
准备基因表达矩阵
在做xx(后面以水稻为例)加权基因共表达网络之前,我们需要先获得若干水稻样本中各基因的表达情况,即基因表达矩阵,类似如下图(后期我会开放下载渠道,记得提醒我)
第一列表示基因的ID号,第一行表示样本的名称(一共有十个样本),矩阵中的数字表示该基因在该样本中的表达量但是最开始时,我得到的表达矩阵是下图这样的
可以看出是把每行所有的数据放在一列了,按如下操作(选择所有行->数据->分列->完成)
先把第二行删掉,因为有count字符会影响后续的操作但并非所有基因都在稻瘟病致病过程中有表达,不仅表达量差异较大,每个样本中的基因表达量也有部分差异,本文通过计算每个基因在10个样本中表达量的中值绝对偏差,作为后续过滤与排序的标准。
我决定用excel去做,在L2中输入下图所示公式(即计算中值绝对偏差)
公式:=MEDIAN(ABS(B2-MEDIAN(B2:K2)),ABS(C2-MEDIAN(B2:K2)),ABS(D2-MEDIAN(B2:K2)),ABS(E2-MEDIAN(B2:K2)),ABS(F2-MEDIAN(B2:K2)),ABS(G2-MEDIAN(B2:K2)),ABS(H2-MEDIAN(B2:K2)),ABS(I2-MEDIAN(B2:K2)),ABS(J2-MEDIAN(B2:K2)),ABS(K2-MEDIAN(B2:K2)))
然后再拖到所有列中
最后选中L列,按降序进行排序
得到下图所示矩阵
再复制其中前5000个数据,新建excel粘贴即可得到用于处理的表达矩阵,同时删除最后一列(因为最后一列存放的是中值绝对偏差),按下图操作,然后点击另存为,选择其他格式中的csv。千万不要直接改后缀,否则R语言读取时会出错的。
离群样本检测
library(WGCNA)
options(stringsAsFactors =FALSE)
enableWGCNAThreads() # 开多线程,为了速度更快
datExpr=read.csv("D:\\desktop\\DiffAnalysis\\top5000.csv") # 读取文件并用对象datExpr表示
因为要用hclust()进行离群样本检测,所以我们进行聚类的对象是样本,而将每个样本对应的5000个基因表达量作为特征向量。所以我们需要对矩阵进行转置
matrix=t(datExpr) # 用matrix保存转置后的datExpr
但是你按下面代码查看matrix矩阵的一部分,会发现
matrix[1:6,1:6] # 显示1到6行以及1到6列
发现第一行(gene这行)不是表达量形式,你如果不删除这行,直接聚类的话会出错
所以运行下面代码删除第一行
matrix=matrix[-1,] # 删除第一行
删除成功
对10个样本进行层次聚类
# hclust是聚类函数,dist()计算距离,method确定距离的方式
geneTree = hclust(dist(matrix), method = "average")
# 若需要分类,添加聚类分类矩形,如分为3类
rect.hclust(geneTree,k=3)
plot(geneTree) # 画出聚类图
聚类树不存在高度明显偏移现象,因此本文所使用的10个水稻样本未出现离群样本。
# 确定模型软阈值