| 撰文：莫北

在之前《如何用Gephi繪製漂亮的網路圖》一文中，已經為大家介紹過如何使用Gephi繪製下圖這樣的共發生網路圖（Co-occurrence network），該文主要是介紹透過準備“

邊檔案

”和“

點檔案

”來建立網路圖。

但有時候繪圖比較容易，難的是準備資料。

前文也提過，最常見的共發生網路分析過程是先計算相關係數矩陣，接著使用igraph生成gml格式的網路圖檔案，最後匯入Gephi對網路圖進行個性化繪製。

當然，如果你不願使用R語言計算，也可以使用OmicShare Tools中的組內相關性分析工具計算得到相關係數矩陣，然後整理成“邊檔案”，使用Gephi或Cytoscape做個性化繪圖。

工具網址：

https：//www。

omicshare。com/tools/

但如果資料量較大，手工整理顯然非常費時費力，有OmicShare使用者在OmicShare Forum中留言：“資料的準備太難了！”

因此，下面我將從最初的OTU絕對丰度表開始，介紹如何使用R語言進行網路圖資料的準備、網路圖檔案生成和網路圖繪製。

本文的重點內容是

資料的過濾、相關係數的計算和“graphml”格式網路圖檔案的生成。

資料匯入與過濾

為了便於大家練習使用，本文所用到的範例資料已上傳到OmicShare論壇上，大家可前去下載。

連結：

https：//www。

omicshare。com/forum/thr

ead-6116-1-1。html

#載入所需R包；

library（igraph）

library（Hmisc）

#工作目錄設定；

setwd（“C：\\Users\\MHY\\Desktop\\Co-occurrence_Network\\example”）

#讀入OTU絕對丰度表；

otu=read。table（“

otu_data。xls

” ，header=T，row。names = 1，sep = “\t”）

#轉成矩陣，之後的相關性計算需要矩陣物件；

otu<-as。matrix（otu）

dim（otu）

head（otu）

#將丰度值大於1的值替換為1，便於計算不同otu的檢測率；

dt<-otu

dt［dt>1］<-1

#將樣本發現率低於20%的otu過濾掉；

no<-which（rowSums（dt）/ncol（dt）>0。2）

length（no）

otu<-otu［no，］

相關性係數計算

主要利用Hmisc包的

rcorr（）函式

計算spearman相關係數，當然，也可計算

pearson相關係數

。rcorr（）函式預設是對列進行兩兩計算，故須轉置，計算otu間的相關性。

#計算相關性係數；

sp。cor<-< span=“”>rcorr（t（otu），type=“spearman”）

#提取r、p值矩陣；

r。cor<-< span=“”>sp。cor$r

p。cor<-< span=“”>sp。cor$P

#使用Benjamini-Hochberg（“FDR-BH”）法進行多重檢驗校正；

p。adj <-< span=“”> p。adjust（p。cor， method=“BH”）

#指定閾值；

r。cutoff=0。6

p。cutoff=0。001

#在只考慮“正相關”情況下；

r。matrix<-< span=“”>r。cor

p<-< span=“”>p。adj

#將矩陣中不符合條件的r值替換為0；

r。matrix［which（r。cor <=< span=“”> r。cutoff）］=0

r。matrix［which（p。adj>p。cutoff）］=0

#刪掉相關係數矩陣資料全都為0（

對角線

處的1不計）的行和列；

r。matrix<-< span=“”>r。matrix［which（rowSums（r。matrix）！=1），］

r。matrix<-< span=“”>r。matrix［，which（colSums（r。matrix）！=0）］

#檢視過濾後的矩陣；

dim（r。matrix）

r。matrix［1：7，1：7］

生成網路圖

從結構來看，對稱矩陣是展示網路圖非常合適的

資料結構

，比如相關係數矩陣是兩個OTU的丰度相關係數，而網路圖恰是Source和Target兩個OTU結點的連線，相關係數可作為權重。如果是n維矩陣（即n個OTU），除去對角線，共有n（n-1）/2個關係對，當然用組合公式也可算得。

#使用

鄰接矩陣

（即相關係數矩陣）建立網路；

g1<-< span=“”>graph。adjacency（r。matrix，weight=T，mode=“undirected”）

#去掉冗餘的邊（

multiple edges

、loop edges）；

g1<-< span=“”>simplify（g1）

#生成網路圖的結點標籤（OTU id）和degree屬性；

V（g1）$label <-< span=“”> V（g1）$name

V（g1）$degree <-< span=“”> degree（g1）

#檢視網路圖的物件結構；

print（g1）

第一行：‘U’ 表示 undirected ，相反， ‘D’ 表示 directed 網路圖；‘N’ 表示網路圖具有name屬性；‘W’表示weighted graphs；接下來的‘-’表示type屬性不存在；接下來是網路共有240個結點和1068條邊。

匯出graph物件

#將網路圖匯出為“graphml”、“gml”格式，方便匯入Gephi中使用；

write_graph（g1， “g1。graphml”， format=“graphml”）

write_graph（g1， “g1。gml”， format=“gml”）

#支援的格式有“edgelist”， “pajek”， “ncol”， “lgl”，“graphml”， “dimacs”， “gml”， “dot”， “leda”）， 。。。

匯出“graphml”格式的網路物件，之後可用Gephi和Cytoscape等網路圖視覺化軟體直接開啟，做進一步調整、美化。

使用igraph繪製網格圖

雖然直接使用igraph繪製個性化網路圖沒有Gephi簡單，igraph主要用於計算網路圖的拓撲性質（topological properties），我這裡仍然做了一些嘗試。

#計算

群體結構

（short random walks）；

c<-< span=“”> cluster_walktrap（g1）

#使用預設顏色列表；

V（g1）$color <-< span=“”>c$membership+1

#繪製網路圖；

layout<-< span=“”> layout。fruchterman。reingold

#也可以嘗試其他的Layout；

layout<-< span=“”>

layout_as_tree

layout<-< span=“”> layout_nicely

layout<-< span=“”> layout_on_sphere

#……

#繪製網路圖；

plot。igraph（g1，vertex。size=4，

vertex。label=NA，

edge。curved

=F，

edge。size=1。5，

layout= layout。fruchterman。reingold）

#預設的結點大小為15，這裡改為4，預設邊的粗細為1，這裡改為1。5；edge。curved若改為TURE，可實現類似文章開頭那樣的效果。

繪製的效果如下：

當然如果想做個性化調整，使用Gephi直接開啟igraph匯出的“g1。graphml”檔案對網路圖進行視覺化，如下，具體的使用方法可參考昨天發的《如何用Gephi繪製漂亮的網路圖？》一文。

使用Gephi匯出網路圖的效果如下，使用同樣的資料和佈局方式，我們可以看到使用igraph和Gephi得到網路圖雖有差異，但在模組上還是比較一致的。

這次就分享到這裡啦！還不會的同學可以收藏起來對照學習~有什麼問題或者想法可以留言哦~

參考資料：

https：//

igraph。org/r/