寫在前面今年的畢業(yè)生應(yīng)該是6月底離校,，但是這位朋友到了九月份終于完善了手上的工作,，對(duì)這幾年來了一個(gè)完善的結(jié)束，感謝送給我的一方印章，也祝鎧鳴大展宏圖。再有就是ggclusterNet功能也完善了90%,，剩下的慢慢補(bǔ)充上來，尤其是代碼的優(yōu)化要抽出些時(shí)間多和大佬交流,，優(yōu)化,。要知道這部分實(shí)在是也很頭疼。原計(jì)劃今年將微生物組分析的四個(gè)R包夠?qū)懞?，想來想去似乎目前只有easystat和ggClusterNet可以順利完成,，我也會(huì)抓緊時(shí)間在今年結(jié)束的時(shí)候至少完成一個(gè)包（phyloBiota）的功能。 導(dǎo)入所需要的R包#--導(dǎo)入所需R包#-------
library(phyloseq)
library(igraph)
library(network)
library(sna)
library(ggplot2)
library(ggClusterNet)數(shù)據(jù)導(dǎo)入我在ggClusterNet包中加入了一個(gè)phyloseq封裝的微生物組數(shù)據(jù),，這份數(shù)據(jù)來自劉永鑫老師的示例文件,。這里直接使用data()提取即可 data(ps)網(wǎng)絡(luò)分析一網(wǎng)打盡為什么說是一網(wǎng)打盡呢，我使用了network函數(shù)打包了ggclusterNet封裝的各種功能,，一般情況下網(wǎng)絡(luò)所需要的內(nèi)容都基本上囊括了,，一條代碼即可分析完成網(wǎng)絡(luò)，并且可以按照分組同時(shí)分析完成多組網(wǎng)絡(luò),。完成多組網(wǎng)絡(luò)的共同比對(duì)和輸出,。 network 函數(shù)計(jì)算相關(guān)矩陣生成邊和節(jié)點(diǎn)文件保存為Gephi格式； 使用sne中十多種layout方法并用ggplotv出圖,，組圖 計(jì)算節(jié)點(diǎn)性質(zhì) 計(jì)算網(wǎng)絡(luò)性質(zhì) 計(jì)算隨機(jī)網(wǎng)絡(luò)和樣本網(wǎng)絡(luò)比對(duì)，檢測(cè) 計(jì)算zipi 分組網(wǎng)絡(luò)計(jì)算自動(dòng)保存
先看看出來的結(jié)果 
函數(shù)測(cè)試這個(gè)例子表示的是提取微生物相對(duì)豐度在千分之一以上的OTU,，,，然后對(duì)其做 皮爾森相關(guān)（推薦使用斯皮爾曼相關(guān)），相關(guān)閾值選擇大于0.6的,，并且顯著性閾值低于0.05的相關(guān)關(guān)系,，并且使用ggrapg可視化網(wǎng)絡(luò)，通過igraph計(jì)算網(wǎng)絡(luò)節(jié)點(diǎn)和整體屬性,，后計(jì)算zipi模塊化信息,，將全部結(jié)果保存在path中。 ```{R}
path = “./result”
dir.create(path)
result = network(ps = ps,N = 0.001,r.threshold=0.6,p.threshold=0.05,label = FALSE,path = path ,zipi = TRUE)
除了保存的結(jié)果,，這里還有全部分組的網(wǎng)絡(luò)可視化輸出,，共比較,，和全部網(wǎng)絡(luò)全局參數(shù)一起。
```{R}
# 全部樣本的網(wǎng)絡(luò)比對(duì)
p = result[[1]]
p
# 全部樣本網(wǎng)絡(luò)參數(shù)比對(duì)
data = result[[2]]
plotname1 = paste(path,"/network_all.jpg",sep = "")
ggsave(plotname1, p,width = 22,height = 22)
plotname1 = paste(path,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 22,height = 22)
tablename <- paste(path,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)結(jié)果展示-網(wǎng)絡(luò)可視化
結(jié)果展示-Gephi輸入的兩個(gè)文件上面的可視化如果你不喜歡,，就可以在gephi中出圖,，使用下面的這兩個(gè)文件。 KO_Gephi_edge.csv KO_Gephi_node.csv ··· 結(jié)果展示-網(wǎng)絡(luò)性質(zhì)同隨機(jī)網(wǎng)絡(luò)比對(duì)KO_net_VS_erdos_properties.csv WT_net_VS_erdos_properties.csv ··· 結(jié)果展示-節(jié)點(diǎn)屬性KO_node_properties.csv WT_node_properties.csv ··· 結(jié)果展示-zipi計(jì)算結(jié)果和圖形KOZiPi.csv KO_ZiPi.pdf 
結(jié)果展示-冪律分布這里看到不同網(wǎng)絡(luò)的冪律分布是不同的,，和隨機(jī)網(wǎng)絡(luò)之間的差距也不是很相同,，這表明了不同分組的網(wǎng)絡(luò)不一定都是無標(biāo)度網(wǎng)絡(luò)。  
二分網(wǎng)絡(luò)一網(wǎng)打盡二分網(wǎng)絡(luò)的介紹在之前我已經(jīng)講過了,，參見推送： 點(diǎn)擊跳轉(zhuǎn),。 或者在專輯中查看。 微生物結(jié)合環(huán)境因子等其他指標(biāo)網(wǎng)絡(luò)這里是我構(gòu)造的一組模擬環(huán)境因子數(shù)據(jù)
環(huán)境因子的格式和vegan中的一樣,，也就是樣本作為列,，指標(biāo)作為行。
—導(dǎo)入多個(gè)環(huán)境變量及其分組data1 <- read.delim(“..//ori_data/bios_network_normal//env.txt”) head(data1) row.names(data1) <- data1SampleID = NULL data1 <- as.data.frame(t(data1)) head(data1) data1$ID = row.names(data1) data1 <- data1 %>% select(ID,everything())
Gru <- read.delim(“..//ori_data/bios_network_normal//env_group.txt”)細(xì)菌真菌phyloseq對(duì)象我們已經(jīng)準(zhǔn)備好了,，這里可以參照學(xué)習(xí),， 如果不會(huì)構(gòu)建phyloseq對(duì)象，可以參見之前的推文,。```{R}ps16 = readRDS("..//ori_data/bios-network_micro/16S/ps_16S.rds")psIT = readRDS("..//ori_data/bios-network_micro/ITS/ps_ITS.rds")ps.merge <- merge16S_ITS (ps16,psIT,N16s = 0.005,NITS = 0.005)ps.merge
同樣是一條函數(shù)即可計(jì)算得到類似于上面的全部內(nèi)容,，由于zipi的計(jì)算不夠快,，所以我這里設(shè)置：zipi = FALSE。下面是結(jié)果文件的總覽,，少了zipi的部分： library(ggpubr)
path = "./result_bio/"
dir.create(path)
result <- corBionetwork(ps = ps.merge,
N = 0.001,
r.threshold = 0.9, # 相關(guān)閾值
p.threshold = 0.05,
label = FALSE,
group = "Group",
env = data1, # 環(huán)境指標(biāo)表格
envGroup = Gru,# 環(huán)境因子分組文件表格
# layout = "fruchtermanreingold",
path = path,# 結(jié)果文件存儲(chǔ)路徑
fill = "group", # 出圖點(diǎn)填充顏色用什么值
size = "igraph.degree", # 出圖點(diǎn)大小用什么數(shù)據(jù)
scale = TRUE, # 是否要進(jìn)行相對(duì)豐度標(biāo)準(zhǔn)化
bio = TRUE, # 是否做二分網(wǎng)絡(luò)
zipi = FALSE, # 是否計(jì)算ZIPI
step = 100, # 隨機(jī)網(wǎng)絡(luò)抽樣的次數(shù)
width = 12,
height = 10
)p = result[[1]]
p
# 全部樣本網(wǎng)絡(luò)參數(shù)比對(duì)
data = result[[2]]
plotname1 = paste(path,"/network_all.jpg",sep = "")
ggsave(plotname1, p,width = 22,height = 22)
plotname1 = paste(path,"/network_all.pdf",sep = "")
ggsave(plotname1, p,width = 22,height = 22)
tablename <- paste(path,"/co-occurrence_Grobel_net",".csv",sep = "")
write.csv(data,tablename)結(jié)果展示這里我就簡(jiǎn)單展示一下出圖,，使用的ggclusterNet的函數(shù)出圖，展示不同模塊之間的關(guān)系：
|
