#!/usr/bin/env R

#定义函数

library(picante)       #picante 包加载时默认同时加载 vegan

alpha <- function(x, tree = NULL, base = exp(1)) {

        est <- estimateR(x)

        Richness <- est[1, ]

        Chao1 <- est[2, ]

        ACE <- est[4, ]

        Shannon <- diversity(x, index = 'shannon', base = base)

        Simpson <- diversity(x, index = 'simpson')    #Gini-Simpson 指数

        Pielou <- Shannon / log(Richness, base)

        goods_coverage <- 1 - rowSums(x == 1) / rowSums(x)

        result <- data.frame(Richness, Shannon, Simpson, Pielou, Chao1, ACE, goods_coverage)

        if (!is.null(tree)) {

                PD_whole_tree <- pd(x, tree, include.root = FALSE)[1]

                names(PD_whole_tree) <- 'PD_whole_tree'

                result <- cbind(result, PD_whole_tree)

        }

        result

}

#现在直接使用定义好的命令 alpha ()，一步得到多种 Alpha 多样性指数

#加载 OTU 丰度表和进化树文件

otu <- read.delim('feature-table.tsv', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)

otu <- t(otu)

tree <- read.tree('tree.nwk')

#不包含谱系多样性，无需指定进化树；Shannon 公式的 log 底数我们使用 2

##alpha_all <- alpha(otu, base = 2)

#包含谱系多样性时，指定进化树文件；Shannon 公式的 log 底数我们使用 2

alpha_all <- alpha(otu, tree, base = 2)

#输出保存在本地

write.table(alpha_all, 'alpha.txt', sep="\t", quote = FALSE)

#!/usr/bin/env R

# 基于 github 安装包，需要 devtools，检测是否存在，不存在则安装

if (!requireNamespace("devtools", quietly = TRUE))

    install.packages("devtools")

# 加载 github 包安装工具

library(devtools)

# 检测 amplicon 包是否安装，没有从源码安装

if (!requireNamespace("amplicon", quietly = TRUE))

    install_github("liaochenlanruo/amplicon")

# 提示升级，选择 3 None 不升级；升级会容易出现报错

# library 加载包，suppress 不显示消息和警告信息

suppressWarnings(suppressMessages(library(amplicon)))

# 设置工作目录，根据实际情况设定

setwd("用户自己的工作目录")

# 读取元数据，参数指定包括标题行 (TRUE)，列名为 1 列，制表符分隔，无注释行，不转换为因子类型

metadata = read.table("sample-metadata.tsv", header=T, row.names=1, sep="\t", comment.char="", stringsAsFactors = F)

# 读取第一步中 vegan 计算的 6 种 alpha 多样性指数

#alpha_div = read.table("alpha.txt", header=T, row.names=1, sep="\t", comment.char="")

alpha_div <- read.delim("alpha.txt", quote="", stringsAsFactors=FALSE)

# 绘制各组香农指数分布，外层 () 可对保存的图形同时预览

(p1 = alpha_boxplot(alpha_div, index = "Shannon", metadata, groupID = "Depth2"))

## 物种丰富度 Richness 指数

(p2 = alpha_boxplot(alpha_div, index = "Richness", metadata, groupID = "Depth2"))

#Gini-Simpson 指数（我们平时常用的 Simpson 指数即为 Gini-Simpson 指数）

(p3 = alpha_boxplot(alpha_div, index = "Simpson", metadata, groupID = "Depth2"))

#Chao1 指数

(p4 = alpha_boxplot(alpha_div, index = "Chao1", metadata, groupID = "Depth2"))

#ACE 指数

(p5 = alpha_boxplot(alpha_div, index = "ACE", metadata, groupID = "Depth2"))

#Shannon 均匀度（Pielou 均匀度）

(p6 = alpha_boxplot(alpha_div, index = "Pielou", metadata, groupID = "Depth2"))

## 谱系多样性

(p7 = alpha_boxplot(alpha_div, index = "PD_whole_tree", metadata, groupID = "Depth2"))

# 保存图片，指定图片为 pdf 格式方便后期修改，图片宽 89 毫米，高 75 毫米

ggsave(paste0("alpha_boxplot_shannon.pdf"), p1, width=89, height=75, units="mm")

ggsave(paste0("alpha_boxplot_Richness.pdf"), p2, width=89, height=75, units="mm")

ggsave(paste0("alpha_boxplot_Simpson.pdf"), p3, width=89, height=75, units="mm")

ggsave(paste0("alpha_boxplot_Chao1.pdf"), p4, width=89, height=75, units="mm")

ggsave(paste0("alpha_boxplot_ACE.pdf"), p5, width=89, height=75, units="mm")

ggsave(paste0("alpha_boxplot_Pielou.pdf"), p6, width=89, height=75, units="mm")

ggsave(paste0("alpha_boxplot_PD_whole_tree.pdf"), p7, width=89, height=75, units="mm")

setwd("用户自己的路径")

if (!requireNamespace("devtools", quietly=TRUE))

    install.packages("devtools")

library(devtools)

if (!requireNamespace("amplicon", quietly=TRUE))

    install_github("microbiota/amplicon")

suppressWarnings(suppressMessages(library(amplicon)))

# 从文件读取元数据，特征表和物种注释

metadata=read.table("sample-metadata.tsv", header=T, row.names=1, sep="\t", comment.char="", stringsAsFactors=F)

taxonomy=read.table("taxonomy.tsv", header=T, row.names=1, sep="\t", comment.char="", stringsAsFactors=F)

otutab2=read.table("level-2.txt", header=T, row.names=1, sep="\t", comment.char="", stringsAsFactors=F)

#堆叠柱状图。以分组均值绘制，展示丰度最高的 19 个门，其余归类为其他 (Other)

# 门水平物种组成表和元数据作为输入，分组列名为 Group，默认显示前 19 个分类，按丰度排序

otutab2=t(otutab2)

(ps=tax_stackplot(otutab2, metadata, groupID="Site", topN=20, style="sample", sorted="abundance"))

ps2 = ps + scale_fill_manual(values = rev(c("#1f77b4", "#aec7e8", "#ff7f0e", "#ffbb78", "#2ca02c", "#98df8a", "#d62728", "#ff9896", "#9467bd", "#c5b0d5", "#8c564b", "#c49c94", "#e377c2", "#f7b6d2", "#7f7f7f", "#c7c7c7", "#bcbd22", "#dbdb8d", "#17becf", "#9edae5"))) + theme(axis.text.x = element_text(angle = 90, hjust = 0.5, vjust = 0.5))

ps2

setwd("用户的工作目录")

library("vegan")

# 读取数据，一份 otu.table 文件和一份分组信息文件

otu <- read.table("feature-table.tsv",row.names=1, header=T,sep="\t",check.names=F)

design <- read.table("sample-metadata.tsv",header=T,sep="\t",row.names=1,check.names=F)

# 数据调整为列名是 OTU，行名是样本名

otu=t(otu)

# 标准化，当然 method 有很多，可以通过？decostand 查看其它 method

vare.hel<-decostand(otu,method="hellinger")

# 计算 Bray-curtis 距离矩阵

vare.dis <- vegdist(vare.hel,method="bray")

# 使用 NMDS 的方法

vare.mds <- metaMDS(vare.hel,distance = "bray")

# 首先提取前两轴坐标

point = scores(vare.mds)

# 将分组文件和得分文件合并

index = merge(design, point,by="row.names",all=F)

# 查看 Stress 值

# Stress 值是反映模型合适程度的指标，NMDS 会多次打乱数据计算 Stress 值，直到找到最合适的模型，也就是最低的 Stress 值；理想状况下，Stress 值为 0，一般 Stress 值低于 0.1 较为合理。

vare.mds

# Stress:     0.07498046 #记住这个数字，后面会打印在图中

# 显着性检验；anosim 本质是基于排名的算法，更加适合 NMDS，这里基于元数据中的深度（Depth）计算

anosim.result<-anosim(vare.dis, design$Depth,permutations =999)

summary(anosim.result)

# ANOSIM statistic R: 0.8036 #记住这个数字，后面会打印在图中

#       Significance: 0.001 #记住这个数字，后面会打印在图中

# tiff 输出图形，适合大部分出版刊物，入门级别分辩率 300,18*14 的长宽；

#tiff(file="beta_bray_NMDS.tiff", res = 300, compression ="none", width=180,height=140,units= "mm")

# 在这里我选择输出为 PDF，后续可以进一步的编辑

pdf(file="beta_bray_NMDS.pdf", width=180,height=140)

#开始出图，代码如下：

library("ggplot2")

p = ggplot(index, aes(x=NMDS1, y=NMDS2, color=as.factor(Depth))) +

      geom_point(size=0) +

     #scale_colour_manual (values = c ("red","blue", "green")) +# 这行代码是用于自定义颜色的，有几个组就用几个颜色，少于 3 个组的话不支持自定义

     labs(x=paste("NMDS1"), y=paste("NMDS2"), title="")

#置信区间当然要加上，有三种方式，线条类型也可以更改。将上面得到的三个指标在图中更换 Stress，R，p。

p+stat_ellipse(type = "t", linetype = 2, level = 0.95, show.legend = TRUE)+

    annotate("text",x=-2.5,y=-1.25,parse=TRUE,size=4,label="'R: '*0.8036",family="serif",fontface="italic",colour="darkred",hjust = 0)+

    annotate("text",x=-2.5,y=-1.35,parse=TRUE,size=4,label="'p: '*0.001",family="serif",fontface="italic",colour="darkred",hjust = 0)+

    annotate("text",x=-2.5,y=-1.45,parse=TRUE,size=4,label="'Stress: '*0.075",family="serif",fontface="italic",colour="darkred",hjust = 0) + geom_point(aes(x=NMDS1, y=NMDS2, shape=Site), size=2)+

	#scale_shape_manual (values = c (15, 19, 17))+# 这行代码是用于自定义形状的，有几个组就用几个形状，少于 3 个组的话不支持自定义

    theme_classic()

dev.off()