Title:

• 16s多樣性分析和可視化（16s diversity analysis and visualization ）

Author:

• GCC

Date:

• 2019年4月

1. clean數(shù)據(jù)

加載第三方包

#setup
knitr::opts_chunk$set(warning = FALSE, message = FALSE,tidy=TRUE, tidy.opts = list(blank = TRUE, width.cutoff=23),dpi = 600)
options(stringsAsFactors = FALSE)

#加載第三方包
library(dplyr)
library(ggplot2)
library(reshape)
library(vegan)
library(tidyr)
library(formatR)

導(dǎo)入備注信息和種豐度

#導(dǎo)入數(shù)據(jù)
setwd('C:\\Users\\GCC\\FMT')
annotations <-
        readxl::read_xlsx('./annotation.xlsx') %>% select(., oldID, ID)
species <- read.delim('./種豐度.txt', sep = '\t')

重新標(biāo)注樣本信息

#重新對樣本命名
newNames <- annotations$I
names(newNames) <- annotations$oldID
newNames <- c(FeatureID = 'FeatureID', newNames)
speciesNew <- plyr::rename(species, newNames)
speciesNew <-
        speciesNew[, stringr::str_detect(colnames(speciesNew), pattern = '(negativeControl)|(FeatureID)|(treat)|(wt)|(model).*')]

2. alpha多樣性分析

Simpson指數(shù)：Simpson反映的是優(yōu)勢種在群落中的地位和作用，若一個群落中優(yōu)勢種占的多，其他非優(yōu)勢物種所占的比例則會減少，那么Simpson 指數(shù)值較大，這說明群落多樣性較低。
Shannon指數(shù)：Shannon值越大，說明群落多樣性越高。

#多樣性分析
set.seed(0408)
H <- diversity(t(speciesNew[, -1]))
#計算Shannon-Wiener指數(shù)
shannon = diversity(t(speciesNew[,-1]), "shannon")
#計算Simpson指數(shù)
simpson = diversity(t(speciesNew[,-1]), "simpson")
#計算Inverse Simpson指數(shù)
insimpson = diversity(t(speciesNew[,-1]), "inv")

pairs(cbind(H, simpson, insimpson), pch = "+", col = "blue")
# Species richness (S) and Pielou's evenness (J):
S <- specnumber(t(speciesNew[, -1])) 
J <- H / log(S)
diversitys <-
        data.frame(shannon = shannon,
                        simpson = simpson,
                        insimpson = insimpson)
diversitys <-
        diversitys[stringr::str_detect(rownames(diversitys), '[^(_)].*'), ]
diversitys <-
        within(diversitys,
        {
        time <-
        stringr::str_split(rownames(diversitys), '_', simplify = T)[, 2] %>% factor(., levels = c(1:10))
        group <-
        stringr::str_extract(rownames(diversitys), '[^(_)]*')
        }) %>% gather(.,
        key = 'index',
        value = 'values',
        c(shannon, simpson, insimpson))

3. 作圖

#index-plotting, fig.align='center', fig.cap = 'the diversity index in different group'
p1 <-
        diversitys %>% ggplot(., aes(time, values, fill = group)) + geom_boxplot(alpha =
        0.8, show.legend = FALSE) + geom_point(
        position = position_jitter(),
        size = 1.5,
        alpha = 0.8,
        aes(color = group),
        show.legend = FALSE
        ) + facet_grid(index ~ group, scales = 'free', labeller = as_labeller(
        c(
        model = 'model',
        negativeControl = 'negative',
        treat = 'treat',
        wt = 'wild type',
        shannon = 'Shannon',
        simpson = 'Simpson',
        insimpson = 'Insimpson'
        )
        )) + labs(x = '', y = 'Alpha Diversity Measure') + theme(
        panel.grid = element_blank(),
        panel.background = element_blank(),
        panel.border = element_rect(color = 'black', fill = NA)
        ) #+ ggpubr::stat_compare_means(aes(group=interaction(time,group)))

print(p1)