本文由阿童木根據(jù)實(shí)踐經(jīng)驗(yàn)而整理,，希望對大家有幫助。

原創(chuàng)微文,，歡迎轉(zhuǎn)發(fā)轉(zhuǎn)載,。

完成以上步驟將得到如下結(jié)果：

1 模擬豐度矩陣

set.seed(1995)

# 設(shè)置隨機(jī)種子：為了能多次獲取同一組隨機(jī)數(shù)

data=matrix(abs(round(rnorm(200, mean=1000, sd=500))), 20, 10)

# 獲取隨機(jī)正整數(shù)矩陣：20個樣本，10個菌種

colnames(data)=paste("Species", 1:10, sep=" ")

# 給細(xì)菌命名

rownames(data)=paste("Sample", 1:20, sep=" ")

# 給樣本命名,，結(jié)果如下：

2 模擬分組文件

group=c("A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "B", "B", "B", "B", "B", "B", "B", "B", "B", "B")

# 把樣本分為A,、B兩組

sample_id=rownames(data)

# 取樣本名

data_group=data.frame(sample_id, group)

# 做分組文件,，如下：

3 標(biāo)準(zhǔn)化豐度

data_norm=data

# 新建data_norm數(shù)據(jù)框，賦初值

for(i in 1:20){

    sample_sum=apply(data, 1, sum)

    # 計(jì)算每個樣本的菌種總豐度

    for(j in 1:10){

        data_norm[i,j]=data[i,j]/sample_sum[i]

        # [每個樣本的每個物種的豐度] / [該樣本物種總豐度] = 相對豐度

    }

}

apply(data_norm, 1, sum)

# 檢查每個樣本物種總豐度是否成功標(biāo)準(zhǔn)化到1,，如下：

4 計(jì)算Alpha多樣性

library(vegan)

# 加載VEGAN包

data_norm_shannon=diversity(data_norm, "shannon")

# 使用VEGAN包中的diversity函數(shù)計(jì)算每個樣品的Shannon Alpha多樣性指數(shù)

data_ggplot=data.frame(data_norm_shannon)

# 多樣性計(jì)算結(jié)果如下： 

data_ggplot=data.frame(data_ggplot, data_group["group"])

# 添加分組信息,，如下：

5 T檢驗(yàn)

成組設(shè)計(jì)T檢驗(yàn)需要兩個條件：1）個體間相互獨(dú)立；2）兩組樣本均來自正態(tài)分布總體,；3）方差齊,。我們可以利用直方圖和QQ圖觀察樣本是否呈正態(tài)分布，當(dāng)然也可以直接用夏皮羅-威爾克（Shapiro-Wilk）假設(shè)檢驗(yàn)的方法判斷樣本數(shù)據(jù)的分布,。如果樣本數(shù)據(jù)符合正態(tài)分布,，那么可以接著用Bartlett檢驗(yàn)不同組樣本的方差齊性。兩個條件都滿足的話就可以做T檢驗(yàn)了,。我這里使用的是隨機(jī)產(chǎn)生的數(shù)據(jù),，所以分析結(jié)果就不要太在意了。

group_A=data_ggplot$data_norm_shannon[data_ggplot$group=='A']

# 獲取A組的多樣性數(shù)據(jù)

group_B=data_ggplot$data_norm_shannon[data_ggplot$group=='B']

# 獲取B組的多樣性數(shù)據(jù)

hist(group_A)

hist(group_B)

# 繪制多樣性指數(shù)分布直方圖

# 觀察形狀若為倒鐘形那便是接近正態(tài)分布的

qqnorm(group_A)

qqnorm(group_B)

# 繪制多樣性指數(shù)QQ圖

# 觀察形狀是一條連接主對角線的線那便是接近正態(tài)分布

shapiro.test(data_ggplot$data_norm_shannon)

# 夏皮羅-威爾克（Shapiro-Wilk）檢驗(yàn)正態(tài)性,，p>0.05,，接受原假設(shè)，符合正態(tài)分布

bartlett.test(data_norm_shannon~group,data=data_ggplot)

# 巴特利特（Bartlett）檢驗(yàn)方差齊性,，p>0.05,，接受原假設(shè)，即兩樣本數(shù)據(jù)方差齊

with(data_ggplot,t.test(formula=data_norm_shannon~group,conf.level=0.95))

# T檢驗(yàn),，p>0.05,，沒法拒絕原假設(shè)，兩組Shannon多樣性指數(shù)差異不顯著

6 ggplot繪制箱形圖

library(ggplot2)

# 加載R包ggplot2

alpha_boxplot=ggplot(data_ggplot, aes(x=group, y=data_norm_shannon, fill=group))+

# 添加數(shù)據(jù),、xy值,、 顏色參數(shù)給畫圖函數(shù)ggplot

geom_boxplot()+

# 盒圖

labs(title="Alpha diversity", x="Group", y="Shannon index")+

# 標(biāo)題

theme(plot.title=element_text(hjust=0.5), legend.title=element_blank())

# 標(biāo)題居中

pdf('result.pdf')

alpha_boxplot

dev.off()

# 保存結(jié)果，打開result.pdf文件,，結(jié)果如下：

參考

1）R語言-t檢驗(yàn)

https://blog.csdn.net/weixin_38322363/article/details/89811964

2）Shapiro-Wilk檢驗(yàn)

https://www.jianshu.com/p/e202069489a6

3) 方差齊性檢驗(yàn)的原理

https://zhuanlan.zhihu.com/p/26943946

感謝您的閱讀~

你可能還喜歡