---

title: "Microeco"

author: "LGP"

date: "2021/4/25"

output: html_document

---

```{r setup, include=FALSE}

knitr::opts_chunk$set(echo = TRUE)

```

##R安装

```{r}

# Github地址：https://github.com/ChiLiubio/microeco

#install.packages('microeco') # 安装 

# 或者通过Github安装最新版 因为Cran发布会有延迟

#devtools::install_github("ChiLiubio/microeco")

library(microeco)

#?microtable # 展示所有的功能和细节描述

```

##导入数据

```{r #加载该包示例数据}

data("sample_info_16S") # 分类表

head(sample_info_16S)

data("otu_table_16S") # OTU丰度表

head(otu_table_16S)

data("taxonomy_table_16S") # 物种注释表

head(taxonomy_table_16S)

data("phylo_tree_16S")

#进化树还有一个为树文件phylo_tree_16S，一般计算加权Alpha多样性时候会用到。

#class(phylo_tree_16S) 

# “phylo”

data("env_data_16S") #  环境因子表

head(env_data_16S)

```

```{r #加载一些基础包}

library(magrittr)

library(ggplot2)

theme_set(theme_bw())

```

```{r #创建R6类型数据}

#创建微生物数据类型，相比于phyloseq类型，microeco包的对象创建更简单。注意otu表的列名和样本表的行名要一致。

dataset <- microtable$new(sample_table = sample_info_16S, otu_table = otu_table_16S, tax_table = taxonomy_table_16S, phylo_tree = phylo_tree_16S)

class(dataset)

print(dataset)

#移除没有比对到"k__Archaea"和"k__Bacteria"域的OTU

dataset$tax_table %<>% subset(Kingdom == "k__Archaea" | Kingdom == "k__Bacteria")

print(dataset)

#移除比对到"mitochondria"和"chloroplast"相关的OTU

dataset$filter_pollution(taxa = c("mitochondria", "chloroplast"))

print(dataset)

#然后为使otu_table，tax_table和phylo_tree中的OTU相同，使用tidy_dataset()。

dataset$tidy_dataset()

print(dataset)

#我们使用sample_sums检查每个样品的测序数量

dataset$sample_sums()%>% range

#为了减少不同样品物种数目对多样性测量的影响，使用rarefy_samples进行稀释抽样。

## 比如，稀释到10000条reads

dataset$rarefy_samples(sample.size = 10000)

dataset$sample_sums() %>% range

```

##保存物种丰度表

```{r}

#使用cal_abund计算每个分类水平物种丰度

dataset$cal_abund()

# return dataset$taxa_abund

class(dataset$taxa_abund)

dir.create("taxa_abund")

dataset$save_abund(dirpath = "taxa_abund")

```

#计算Alpha多样性

```{r #使用cal_alphadiv函数}

# 如想增加系统发育多样性，使用PD =TRUE，速度略慢

dataset$cal_alphadiv(PD = FALSE)

# 返回物种多样性

class(dataset$alpha_diversity)

## 保存结果

dir.create("alpha_diversity")

dataset$save_alphadiv(dirpath = "alpha_diversity")

```

#计算Beta多样性

```{r #使用cal_betadiv函数，支持的距离类型：Bray-curtis, Jaccard, weighted Unifrac and unweighted unifrac}

dataset$cal_betadiv(unifrac = TRUE)

# 返回结果

class(dataset$beta_diversity)

# 保存

dir.create("beta_diversity")

dataset$save_betadiv(dirpath = "beta_diversity")

```

#trans_abun对象

#用于转换生物分类丰度数据，以便使用ggplot2软件包绘制生物分类丰度。我们首先将此类用于条形图。

```{r #绘制群落组成图}

# 创造trans_abun对象

# 使用排名前10 的门

t1 <- trans_abund$new(dataset = dataset, taxrank = "Phylum", ntaxa = 10)

# 在x轴上删除了样本名称，并添加了分面以根据组显示丰度。

t1$plot_bar(others_color = "grey70", facet = "Group", xtext_keep = FALSE, legend_text_italic = FALSE)

#以组分开绘制

# 以组平均值绘制群落柱状图

t1 <- trans_abund$new(dataset = dataset, taxrank = "Phylum", ntaxa = 10, groupmean = "Group")

t1$plot_bar(others_color = "grey70", legend_text_italic = FALSE)

#绘制冲击图

t1 <- trans_abund$new(dataset = dataset, taxrank = "Phylum", ntaxa = 10)

# use_alluvium = TRUE 绘制冲击图, clustering =TRUE 重新排序样本

t1$plot_bar(use_alluvium = TRUE, clustering = TRUE, xtext_type_hor = FALSE, xtext_size = 6)

#箱线图

# 展示纲水平15个物种

t1 <- trans_abund$new(dataset = dataset, taxrank = "Class", ntaxa = 15)

t1$plot_box(group = "Group")

# 展示属水平前40物种

t1 <- trans_abund$new(dataset = dataset, taxrank = "Genus", ntaxa = 40)

t1$plot_heatmap(facet = "Group", xtext_keep = FALSE, withmargin = FALSE)

#饼图

t1 <- trans_abund$new(dataset = dataset, taxrank = "Phylum", ntaxa = 6, groupmean = "Group")

# all pie chart in one row

t1$plot_pie(facet_nrow = 1)

```

#trans_venn对象

#Trans_venn对象用于venn分析，为了分析组间独有和共有OTUs，我们首先按组合并样品。### Venn图

```{r}

# merge samples as one community for each group

dataset1 <- dataset$merge_samples(use_group = "Group")

# 创建trans_venn对象

t1 <- trans_venn$new(dataset1, ratio = "seqratio")

t1$plot_venn()

#整数是OTU数量，括号中数据为测序数/总测序量

#花瓣图

#当分组太多时候，可以使用花瓣图

# use "Type" column in sample_table

dataset1 <- dataset$merge_samples(use_group = "Type")

t1 <- trans_venn$new(dataset1)

t1$plot_venn(petal_plot = TRUE)

#有时，我对独有或共有物种感兴趣，换句话说，独特或共有物种的组成可能解释了群体之间生态特征的不同和相似部分，可以绘制柱状图展示这些部分的组成。

dataset1 <- dataset$merge_samples(use_group = "Group")

t1 <- trans_venn$new(dataset1)

# 将venn结果转换为sample-species表，此处不考虑丰度，仅考虑存在/不存在信息。

t2 <- t1$trans_venn_com(use_OTUs_frequency = TRUE)

# t2是一个新对象

class(t2)

# 计算丰度比例

t2$cal_abund()

# 转变绘图

t3 <- trans_abund$new(dataset = t2, taxrank = "Genus", ntaxa = 10)

t3$plot_bar(bar_type = "part", legend_text_italic = T, ylab_title = "Frequency (%)", xtext_type_hor = FALSE)

#使用饼图展示每一部分门水平成分

t3 <- trans_abund$new(dataset = t2, taxrank = "Phylum", ntaxa = 8)

t3$plot_pie(facet_nrow = 3, use_colors = rev(c(RColorBrewer::brewer.pal(8, "Dark2"), "grey50")))

```

#trans_alpha对象

#可以使用trans_alpha类转换和绘制Alpha多样性结果。创建trans_alpha对象可以返回两个数据帧：alpha_data和alpha_stat。

```{r #stat}

t1 <- trans_alpha$new(dataset = dataset, group = "Group")

# return t1$alpha_stat

t1$alpha_stat[1:5, ]

#多组差异检验

#然后，我们使用KW秩和检验和多重矫正的Anova进行组间差异检验

t1$cal_diff(method = "KW")

# return t1$res_alpha_diff

t1$res_alpha_diff[1:5, ]

#使用Anova检验

t1$cal_diff(method = "anova")

# return t1$res_alpha_diff

t1$res_alpha_diff

#误差图

#绘制误差图展示每组Alpha多样性，然后增添ducan检验

t1$plot_alpha(add_letter = TRUE, measure = "Chao1")

#箱线图

#箱线图展示各组差异比较结果

t1$plot_alpha(pair_compare = TRUE, measure = "Chao1")

```

```{r,#trans_beta对象}

#β多样性的距离矩阵可以使用trans_beta类进行变换和绘制。此类中提到的beta多样性的分析主要包括排序，组距离，聚类和manova。我们首先展示使用PCoA的排序。

# 创建对象选择PCoA对象

t1 <- trans_beta$new(dataset = dataset, group = "Group", measure = "bray", ordination = "PCoA")

# t1$res_ordination is the ordination result list

class(t1$res_ordination)

## 绘制PCoA

t1$plot_ordination(plot_color = "Group", plot_shape = "Group", plot_group_ellipse = TRUE)

#然后我们绘制箱线图比较组间距离Beta多样性

# 计算绘制组内距离

t1$cal_group_distance()

# 计算绘制组间距离

t1$cal_group_distance(within_group = FALSE)

# 返回

t1$plot_group_distance(distance_pair_stat = TRUE)

#层级聚类

# 使用replace_name设置标签名称，使用组参数设置颜色

t1$plot_clustering(group = "Group", replace_name = c("Saline", "Type"))

#perMANOVA分析

# 对所有样本进行manova 分析

t1$cal_manova(cal_manova_all = TRUE)

t1$res_manova$aov.tab

# 对每对样品分组间进行manova分析

t1$cal_manova(cal_manova_paired = TRUE)

t1$res_manova

# 对特定的分组进行manova分析: 

t1$cal_manova(cal_manova_set = "Group + Type")

t1$res_manova$aov.tab

```

#trans_diff对象

#主要用于差异分析，设置了三个比较流行的差异分析方法：metastat、Lefse、random forest分析

```{r #metastat分析}

#Metastat取决于排列和t检验，并且在稀疏数据上表现良好。用于两组之间的比较。

# metastat analysis 属水平

t1 <- trans_diff$new(dataset = dataset, method = "metastat", group = "Group", metastat_taxa_level = "Genus")

# t1$res_metastat is the result

# t1$res_metastat_group_matrix is the group comparisons order for plotting

# plot the first paired groups, choose_group = 1

t1$plot_metastat(use_number = 1:10, qvalue = 0.05, choose_group = 1)

#Lefse分析

#Lefse分析结合了非参数检验和线性判别分析，我们在这个包含了这个途径代替python

t1 <- trans_diff$new(dataset = dataset, method = "lefse", group = "Group", alpha = 0.01, lefse_subgroup = NULL)

# t1$res_lefse is the LEfSe result

# t1$res_abund is the abundance information

t1$plot_lefse_bar(LDA_score = 4) # 设置LDA分值

# 表格展示lefse结果

t1$res_lefse[1:5, ]

# 我们也可以基于Lefse分析的结果绘制差异柱状图

t1$plot_diff_abund(use_number = 1:30)

#也可以展示差异分析树状图，

# clade_label_level 5代表门水平

# 加载ggtree包

t1$plot_lefse_cladogram(use_taxa_num = 200, use_feature_num = 50, clade_label_level = 5)

#随机森林分析

#当前方法可以像LEfSe中的方法那样通过自举计算随机森林，并且仅使用重要功能。在分析中选择MeanDecreaseGini作为指标值。

# 使用 rf_taxa_level设置属水平，如果想要使用全部的分类，改变成all

# nresam = 1 and boots = 1 代表无自举和直接使用全部样品

t1 <- trans_diff$new(dataset = dataset, method = "rf", group = "Group", rf_taxa_level = "Genus")

# t1$res_rf 是这个对象的分类结果

# plot the result

t2 <- t1$plot_diff_abund(use_number = 1:20, only_abund_plot = FALSE)

gridExtra::grid.arrange(t2$p1, t2$p2, ncol=2, nrow = 1, widths = c(2,2))

# 中间的星号代表显著性

```

#trans_env对象

```{r }

#RDA分析

#环境变量对于分析微生物群落结构和组装机制非常有用。首先显示RDA分析（db-RDA和RDA）。

# add_data 用于增加环境数据

t1 <- trans_env$new(dataset = dataset, add_data = env_data_16S[, 4:11])

# use bray-curtis 距离进行dbrda

t1$cal_rda(use_dbrda = TRUE, use_measure = "bray")

# t1$res_rda 提取结果

t1$trans_rda(adjust_arrow_length = TRUE, max_perc_env = 10)

# t1$res_rda_trans 转变结果绘图

t1$plot_rda(plot_color = "Group")

#使用属水平

t1$cal_rda(use_dbrda = FALSE, taxa_level = "Genus")

#由于RDA的主要结果与不同箭头之间的投影和角度有关，因此我们使用几个参数来调整箭头的长度以清楚地显示它们。

t1$trans_rda(show_taxa = 10, adjust_arrow_length = TRUE, max_perc_env = 1500, max_perc_tax = 3000, min_perc_env = 200, min_perc_tax = 300)

# t1$res_rda_trans 绘图

t1$plot_rda(plot_color = "Group")

#Mantel 检验

#Mantel测试可用于检查环境变量与距离矩阵之间是否存在显着相关性。

t1$cal_mantel(use_measure = "bray")

# return t1$res_mantel

t1$res_mantel

#环境变量与分类单元之间的相关性对于分析和推断影响社区结构的因素很重要。在此示例中，我们首先执行差异丰度测试和随机森林分析以获得重要属。然后，我们使用这些分类单元进行相关分析。

# 先创造trans_diff对象

t2 <- trans_diff$new(dataset = dataset, method = "rf", group = "Group", rf_taxa_level = "Genus")

# 然后创建trans_env对象

t1 <- trans_env$new(dataset = dataset, add_data = env_data_16S[, 4:11])

# 计算相关性

t1$cal_cor(use_data = "other", p_adjust_method = "fdr", other_taxa = t2$res_rf$Taxa[1:40])

# return t1$res_cor 

# 使用热图绘制相关性代码

t1$plot_corr()

# 添加聚类

t1$plot_corr(pheatmap = TRUE)

# 分组绘制热图展示相关性

t1$cal_cor(by_group = "Group", use_data = "other", p_adjust_method = "fdr", other_taxa = t2$res_rf$Taxa[1:40])

# return t1$res_cor

t1$plot_corr()

```

#环境因子与Alpha多样性之间的相关性热图

```{r }

t1 <- trans_env$new(dataset = dataset, add_data = env_data_16S[, 4:11])

# use add_abund_table parameter to add the extra data table

t1$cal_cor(add_abund_table = dataset$alpha_diversity)

t1$plot_corr()

```

#trans_network对象

#网络是研究微生物生态学中共生模式的常用方法。在这一部分中，我们描述了trans_network类中的所有核心内容。网络构建方法可以分为两种类型：基于相关性和非基于相关性。可以使用几种方法来计算相关性和重要性。首先介绍基于相关性的网络。trans_network中的参数cal_cor用于选择相关性计算方法。

```{r }

# 使用R基础函数cor.test, 比较慢

t1 <- trans_network$new(dataset = dataset, cal_cor = "base", taxa_level = "OTU", filter_thres = 0.0001, cor_method = "spearman")

# return t1$res_cor_p list; one table: correlation; another: p value

## SparCC方法，需要加载SpiecBasi包

#install.packages("SpiecBasi")

#library(SpiecBasi)

## 该方法非常慢，所以要考虑过滤低丰度序列

#t1 <- trans_network$new(dataset = dataset, cal_cor = "SparCC", #taxa_level = "OTU", filter_thres = 0.001, SparCC_simu_num = 100)

## 当OTU数目比较多时候，使用R包WGCNA提代R中基础函数计算相关性

# 加载WGCNA包

#install.packages("WGCNA")

#library(WGCNA)

#t1 <- trans_network$new(dataset = dataset, cal_cor = "WGCNA", #taxa_level = "OTU", filter_thres = 0.0001, cor_method = "spearman")

#参数COR_cut可以被用来选择相关性阈值，此外，COR_optimization = TRUE代表使用RMT理论去寻找最右的相关性阈值代替COR_cut自定义阈值。

# 构建网络，需要igraph包

#t1$cal_network(p_thres = 0.01, COR_optimization = TRUE)

# 返回t1$res_network

## 使用任意系数阈值构建网络

t1$cal_network(p_thres = 0.01, COR_cut = 0.7)

## 保存网络为gexf格式。

## 可用Gephi(https://gephi.org/)打开该文件

## 加载rgexf包

library(rgexf)

t1$save_network(filepath = "network.gexf")

#绘制网络并根据Gephi中计算的模块显示节点颜色。

#现在，我们用门水平信息显示节点颜色，并用正相关和负相关显示边缘颜色。所有使用的数据已存储在network.gexf文件中，包括模块分类，Phylum信息和边分类。

```

#计算网络的拓扑属性

```{r}

t1$cal_network_attr()

# return t1$res_network_attr

#节点分类

t1$cal_node_type()

# return t1$res_node_type

# we retain the file for the following example in trans_func part

network_node_type <- t1$res_node_type

#Zi-Pi识别网络或模块hub

t1$plot_taxa_roles(use_type = 1)

#根据门信息绘制节点的类别

t1$plot_taxa_roles(use_type = 2)

#显示模块的特征基因模块的特征基因，即PCA的第一个主要成分，

#代表了模块种类中丰度的主要变化。

t1$cal_eigen()

# return t1$res_eigen

#进行相关热图，以显示特征基因与环境因素之间的关系。

# 创造trans_env 对象

t2 <- trans_env$new(dataset = dataset, add_data = env_data_16S[, 4:11])

# 计算相关性

t2$cal_cor(add_abund_table = t1$res_eigen)

#绘制相关性热图

t2$plot_corr()

```