代谢网络图四部曲（二）--网络图数据准备实战演练

中秋假期已经结束，国庆节还会远么？在这双节交替的节奏里，麦特绘谱愿意助您快速学习新技能，画最美的图，过最舒心的假。本期推送内容：代谢网络图四部曲（二）--网络图数据准备实战演练。

1. R中环境配置；

2. 创建模拟的代谢物数据集以及代谢物Class分类数据集；

3. 计算数据集中各个代谢物之间的Spearman相关性，获得相关性系数r矩阵、p值矩阵；

4. 将p、r矩阵整理为网络图所需“边文件（edge file）”，并生成附加的用于表示边的颜色、粗细的列；

5. 准备网络图中所需“节点文件（node file）”,利用各个代谢物的分类生成控制颜色的变量，利用节点在网络中的连接数生成权重变量。

##### 载入相应的包 ###############################

if( "pacman" %in% installed.packages() == F){ ### 检查包管理的R包package manger是否已经安装.

install.packages("pacman") ### 如果需要的话，就安装

 cat (" 注意 pacman 包需要R升级到 3.5及以上版本 ,注意更新.")

} 

library ( "pacman")

p_load("Hmisc") ## 提供rcorr函数用于相关性计算

p_load("dplyr") ## 提供管道操作形如“%>%”,用于逻辑更清晰的代码书写

p_load(reshape2) ## 提供数据形式变换的melt函数

###############################

set.seed(20191001) ### 设定随机数种子，保障结果可重复。这里将随机数种子设定为十一，祝大家看完本推送都像十一放假，拥有一个好心情。

2.1  生成原始代谢物浓度示例数据：30个样本、100个代谢物

my_data <-  matrix(rnorm(3000),nrow=30) ### 创建一个30行、100列的矩阵

colnames(my_data) <- paste0("Meta",1:ncol(my_data))   ## 生成列名

rownames(my_data) <- paste0("Sample",1:nrow(my_data)) ## 生成行名

my_data[1:6,1:6] ##查看数据的局部情况

2.2  生成模拟的代谢物Class数据(用于生成节点文件)

node_df <- my_data %>% colnames() %>% data.frame(node = .,stringsAsFactors = F)

nrow(node_df) ### 查看节点文件的行数

## 根据节点文件的行数生成一个模拟的代谢物class矩阵

node_df$class <-  rep(c("Organic nitrogen compounds","Lipids and lipid-like molecules","Nucleosides, nucleotides, and analogues"),c(30,30,40)) ## 采用rep函数生成模拟的代谢物class

## 查看代谢物class矩阵

head(node_df)

### 这里采用的是Hmisc包中的rcorr函数, 确保输入两个矩阵都是 数值型（numeric)的matrix

### 这里使用type="spearman" 控制进行spearman秩相关分析。也可以使用type="pearson" 进行pearson线性相关分析。

result_df <-  Hmisc::rcorr(my_data,type="spearman") ## 写成Hmisc::rcorr是为了防止rcorr函数重名。

r_df <-  result_df$r # 获取r矩阵

p_df <-  result_df$P # 获取p矩阵

###务必注意result_df$P这里的P是大写的，如果用$p进行提取，将得不到正确结果。

## 将结果保存到文件中

write.csv(r_df,"01-1 raw r.csv") #保存r结果

write.csv(p_df,"01-2 raw p.csv") #保存p结果

4.1  读取数据

### 由于许多分析公司、软件都是以这种m*n矩阵形式返回相关性结果，这里我们特意展示了如何将数据读取并进行准备。如果有现成的相关性数据可以忽略前面的步骤，从此处开始

### 获取默认的工作路径，

getwd() 

### 读取文件时，若仅提供文件名，不提供文件夹路径，R会根据文件名在工作路径中查找文件。因此如果我们想要读取自己的数据，将要读取的相关性r矩阵、p矩阵结果文件放入getwd()返回的文件夹中，将文件名修改为自己的文件名，就可以用R进行读取了。上图所示例子返回“E:/”说明工作路径为E盘，我们确认将r矩阵、p矩阵放到E盘即可进行读取。

r_filename <-  "01-1 raw r.csv" ## 此处输入r矩阵的文件名

p_filename <-  "01-2 raw p.csv" ## 此处输入p矩阵的文件名

### 如果r和p值是我们自己进行计算的，可以跳过前面写入csv和下面读取csv的操作，直接利用r_df和p_df两个矩阵进行后续分析。

r_df <-  read.csv(r_filename,stringsAsFactors = F,encoding = "UTF-8") ### stringsAsFactors 决定读取数据时是否将非数字的字符串变成因子。默认为是，这里设定为否。encoding = "UTF-8" 是控制编码方式，可以防止出现乱码。

head(r_df)

p_df <-  read.csv(p_filename,stringsAsFactors = F,encoding = "UTF-8") 

r_df[1:8,1:3] ## 查看r_df 文件

p_df[1:8,1:3] ## 查看p_df 文件

4.2  数据形式转换（“宽数据”转换为“长数据”）

### 数据的“长”、“宽”转换是数据前处理中的重要环节。如下图所示，左侧为“宽”数据，表现形式通常为每行一个观测（observe），每列一个变量（variable），首行为列名（变量名），首列（去除最左上角格子）为观测名，除掉首行首列外，每个格子的值（value）代表行对应的观测在列对应的变量中的表现，这种“宽数据”也是通常我们见到和准备的数据形式。下图右侧为“长”数据，这种数据的特点为：前n-2列（n为长数据的总列数）为使各个观测能够唯一确定的标识变量，如样本ID，受试者编号、动物编号+测试日期（对于多次重复测量，可能需要两个变量组合才能将各个观测彼此区分）。第n-1列为标识变量外的其余各个列的列名，指示该行代表的是何种变量，第n列为该行前n-2行确定的观测在第n-1行确定的变量的值。

## 由于宽数据往往不设循环难以对所有变量进行操作，而长数据往往比宽数据格式更规整，更便于计算机处理。储存某些数据时看起来更规整易读（比如网络图的边文件和节点文件数据），所以一些工具会要求数据的数据是长数据格式的。比如ggplot绘制分组图/分面图时就会要求我们将宽数据整理为长数据。

###我们要想把计算好的r和p值数据转换为edge_file所需要的形式，就要完成相应的长、宽转换。R中有reshape2包中的melt函数可以专门进行这种处理变换。melt完成的工作如下图红色箭头所示。可以将我们前面得到的r、p矩阵由宽数据转换为长数据。

### 提取边文件信息，对边文件信息进行处理

### 设定r 绝对值的cutoff值和 p 的cutoff值

r_abs_cutoff <- 0.30 

p_cutoff <- 0.05

### 使用melt函数将宽数据转化为长数据。data参数是要进行转化的数据集。id.vars 后面输入要进行转换时的依据。

length_r_df <-  melt(data = r_df, id.vars = c("X")) ## 转化r矩阵

length_r_df %>% head() ## 查看转换后的数据

colnames(length_r_df) <-  c("node_1", "node_2", "r") ## 修改转化后数据的列名

length_p_df <-  melt(data = p_df, id.vars =  c("X")) ## 转化p矩阵

colnames(length_p_df) <-  c("node_1", "node_2", "p") ## 修改转化后数据的列名

edge_df <-  cbind(length_r_df,p=length_p_df[,3]) ## 将r、p数据进行合并，得到边文件

4.3  剔除两端连接同一个节点的边

head(edge_df) ## 查看合并后的数据，可以看到数据被转换成了一个四列的数据框

### 可以观察到，第四列的p值中存在NA值，这是由于左右两个节点是一样的，函数在计算完全相同的两个向量的相关系数p值时会返回NA导致的。事实上，这种边也是我们要删除的。因为在相关性网络图中，两端连接同一个节点的边（节点自己连接自己）是没有意义的，无法传递有用的信息。

## 这种节点的自连接会在网络图中显示为某个节点中发出一条弧形的边，绕了一圈回来又连接自己（如下图A中绿色箭头所示）。这种边在相关性网络图中是没有意义的，因此我们要删除这种边数据，再绘制网络图（如下图B中所示）

IS_node_differ <-  edge_df$node_1 != edge_df$node_2 ### 选出左右不相等的行

edge_df <-  edge_df[IS_node_differ,] ### 提取出左右不相等的行

head(edge_df) ### 查看数据的前n行，可以发现左右不相等的行已经删除了。