「微生物组」基于microeco包进行共现性网络分析！！！

为了避免各位错过最新的推文教程，强烈建议大家将“科研后花园”设置为“星标”！

##加载R包
library(microeco) # Microbial Community Ecology Data Analysis
library(WGCNA) # Weighted Correlation Network Analysis
library(magrittr) # A Forward-Pipe Operator for R

####加载数据，包括OTU特征表、样本信息以及物种分类表
#这里我们使用microeco包自带数据集
data(sample_info_16S)
data(otu_table_16S)
data(taxonomy_table_16S)
data(phylo_tree_16S)
data(env_data_16S)
#转换为数据框格式
sample_info_16S <- as.data.frame(sample_info_16S)#确保行名是样本名称，重要
otu_table_16S <- as.data.frame(otu_table_16S)
taxonomy_table_16S <- as.data.frame(taxonomy_table_16S)
phylo_tree_16S <- as.data.frame(phylo_tree_16S)
env_data_16S <- as.data.frame(env_data_16S)

#统一分类信息
taxonomy_table_16S %<>% tidy_taxonomy

#确保样本信息表行名是样本名称
sample_info_16S[1:5, ]
# row.names(sample_info_16S) <- sample_info_16S$SampleID

## 构造microtable
df <- microtable$new(sample_table = sample_info_16S,
                     otu_table = otu_table_16S,
                     tax_table = taxonomy_table_16S,
                     phylo_tree = phylo_tree_16S,
                     auto_tidy = F)#是否使所有文件信息统一

##去除不属于非古菌和细菌的OTU
# df$tax_table %<>% subset(Kingdom == "k__Archaea" | Kingdom == "k__Bacteria")
df$tax_table %<>% .[grepl("Bacteria|Archaea", .$Kingdom), ]

##去除“线粒体”和“叶绿体”污染
df$filter_pollution(taxa = c("mitochondria", "chloroplast"))

##统一各数据的样本和OTU信息
df$tidy_dataset()
df

network <- trans_network$new(dataset = df, 
                              cor_method = "spearman",
                              use_WGCNA_pearson_spearman = TRUE,
                              filter_thres = 0.0001)

network$cal_network(COR_p_thres = 0.01, #设置p阈值
                     COR_cut = 0.6)#设置R值阈值
head(network$res_network)

6 x 1378 sparse Matrix of class "dgCMatrix"
  [[ suppressing 35 column names ‘OTU_1556’, ‘OTU_32’, ‘OTU_3692’ ... ]]
                                                                                
OTU_1556 . . . . . . . . . . . . . . . . . .         . . . . . . . .         . .
OTU_32   . . . . . . . . . . . . . . . . . .         . . . . . . . .         . .
OTU_3692 . . . . . . . . . . . . . . . . . .         . . . . . . . .         . .
OTU_50   . . . . . . . . . . . . . . . . . 0.7030676 . . . . . . . 0.6314738 . .
OTU_999  . . . . . . . . . . . . . . . . . .         . . . . . . . .         . .
OTU_305  . . . . . . . . . . . . . . . . . .         . . . . . . . .         . .
                             
OTU_1556 . . . . . . . ......
OTU_32   . . . . . . . ......
OTU_3692 . . . . . . . ......
OTU_50   . . . . . . . ......
OTU_999  . . . . . . . ......
OTU_305  . . . . . . . ......

 .....suppressing 1343 columns in show(); maybe adjust 'options(max.print= *, width = *)'
 ..............................

###将网络划分模块
network$cal_module(method = "cluster_fast_greedy")

###计算模块的特征基因，即基于PCA分析的第一主成分，以及方差百分比
network$cal_eigen()

###计算网络属性，即拓扑机构参数
network$cal_network_attr()

###获取节点属性
network$get_node_table(node_roles = TRUE)

###获取边属性
network$get_edge_table() 
network$get_adjacency_matrix()

#install.packages("rgexf")
network$save_network(filepath = "network.gexf")

##看枢纽节点数据
network$res_node_table
##可视化
network$plot_taxa_roles(use_type = 1)