畫一個帶統(tǒng)計檢驗的PCoA分析結(jié)果
前情回顧
PERMANOVA原理解釋：這個統(tǒng)計檢驗可用于判斷PCA/PCoA等的分群效果是否顯著！
經(jīng)過前面的鋪墊，下面來實戰(zhàn)一下，理論應(yīng)用于實際看看會出現(xiàn)什么問題？
PERMANOVA 實戰(zhàn) (一)
采用vegan包自帶的一套數(shù)據(jù)（也解釋了如何自己準(zhǔn)備數(shù)據(jù)）看下PERMANOVA的具體代碼和應(yīng)用。
dune數(shù)據(jù)集描述
dune是一套包含了20個樣品和30個物種豐度數(shù)據(jù)的統(tǒng)計表。其格式是常見OTU表轉(zhuǎn)置后的格式，每一行是一個樣品，每一列是一個物種 (檢測指標(biāo))。
library(vegan)
data(dune)
dim(dune)
## [1] 20 30
head(dune)
## Achimill Agrostol Airaprae Alopgeni Anthodor Bellpere Bromhord Chenalbu Cirsarve Comapalu Eleopalu Elymrepe Empenigr Hyporadi
## 1 1 0 0 0 0 0 0 0 0 0 0 4 0 0
## 2 3 0 0 2 0 3 4 0 0 0 0 4 0 0
## 3 0 4 0 7 0 2 0 0 0 0 0 4 0 0
## 4 0 8 0 2 0 2 3 0 2 0 0 4 0 0
## 5 2 0 0 0 4 2 2 0 0 0 0 4 0 0
## 6 2 0 0 0 3 0 0 0 0 0 0 0 0 0如果我們有一個OTU豐度表，怎么轉(zhuǎn)成這個格式呢？
otu_table <- read.table("otutable_rare",sep="\t", row.names=1, header=T)
as.data.frame(t(otu_table))
## OTU1 OTU2 OTU3
## Samp1 2 12 22
## Samp2 13 13 10
## Samp3 14 8 14
## Samp4 15 10 11dune.env是元數(shù)據(jù)信息，包含數(shù)據(jù)的分子信息、生存環(huán)境信息等，記錄了5個因素 (同時包含連續(xù)變量信息和分組變量信息)：
A1: 土壤厚度信息 a numeric vector of thickness of soil A1 horizon.Moisture: 濕度等級信息，分4個等級，1 < 2 < 4 < 5.Management: 分組信息，不同的管理方式 a factor with levels: BF (Biological farming), HF (Hobby farming), NM (Nature Conservation Management), and SF (Standard Farming).Use: 一個分組信息 an ordered factor of land-use with levels: Hayfield < Haypastu < Pasture.Manure: 一個分組信息，0 < 1 < 2 < 3 < 4.
data("dune.env")
head(dune.env)
## A1 Moisture Management Use Manure
## 1 2.8 1 SF Haypastu 4
## 2 3.5 1 BF Haypastu 2
## 3 4.3 2 SF Haypastu 4
## 4 4.2 2 SF Haypastu 4
## 5 6.3 1 HF Hayfield 2
## 6 4.3 1 HF Haypastu 2
summary(dune.env)
## A1 Moisture Management Use Manure
## Min. : 2.800 1:7 BF:3 Hayfield:7 0:6
## 1st Qu.: 3.500 2:4 HF:5 Haypastu:8 1:3
## Median : 4.200 4:2 NM:6 Pasture :5 2:4
## Mean : 4.850 5:7 SF:6 3:4
## 3rd Qu.: 5.725 4:3
## Max. :11.500這個文件就是我們常用的metadata文件，組織格式也一致，每一行是一個樣品，每一列對應(yīng)樣品的不同屬性。
繪制一個PcOA的圖看一下
# 計算加權(quán)bray-curtis距離
dune_dist <- vegdist(dune, method="bray", binary=F)
dune_pcoa <- cmdscale(dune_dist, k=3, eig=T)
dune_pcoa_points <- as.data.frame(dune_pcoa$points)
sum_eig <- sum(dune_pcoa$eig)
eig_percent <- round(dune_pcoa$eig/sum_eig*100,1)
colnames(dune_pcoa_points) <- paste0("PCoA", 1:3)
dune_pcoa_result <- cbind(dune_pcoa_points, dune.env)
head(dune_pcoa_result)
## PCoA1 PCoA2 PCoA3 A1 Moisture Management Use Manure
## 1 -0.35473182 -0.25667235 0.31129225 2.8 1 SF Haypastu 4
## 2 -0.29462318 -0.18609437 0.03355954 3.5 1 BF Haypastu 2
## 3 -0.07276681 -0.29087086 -0.01169171 4.3 2 SF Haypastu 4
## 4 -0.06925423 -0.26419764 -0.01634735 4.2 2 SF Haypastu 4
## 5 -0.30706200 0.03031589 -0.09124310 6.3 1 HF Hayfield 2
## 6 -0.25302974 0.09420852 0.02814297 4.3 1 HF Haypastu 2
library(ggplot2)
ggplot(dune_pcoa_result, aes(x=PCoA1, y=PCoA2, color=Management)) +
labs(x=paste("PCoA 1 (", eig_percent[1], "%)", sep=""),
y=paste("PCoA 2 (", eig_percent[2], "%)", sep="")) +
geom_point(size=4
) + stat_ellipse(level=0.6) +
theme_classic()
## Too few points to calculate an ellipse
## Warning: Removed 1 row(s) containing missing values (geom_path).
樣品中重復(fù)太少了，做不出置信橢圓。換個方式，用ggalt包中的geom_encircle把樣品包起來。
# install.packages("ggalt")
library(ggalt)
ggplot(dune_pcoa_result, aes(x=PCoA1, y=PCoA2, color=Management, group = Management)) +
labs(x=paste("PCoA 1 (", eig_percent[1], "%)", sep=""),
y=paste("PCoA 2 (", eig_percent[2], "%)", sep="")) +
geom_point(size=5) +
geom_encircle(aes(fill=Management), alpha = 0.1, show.legend = F) +
theme_classic() + coord_fixed(1)
那么不同管理風(fēng)格對物種組成是否有顯著影響呢？
關(guān)注不同管理風(fēng)格對物種組成是否有顯著影響
假如關(guān)注的問題是：不同的管理風(fēng)格對物種組成是否有顯著影響？這就是一個典型的單因素非參多元方差分析。因素就是Management。
基于bray-curtis距離進行PERMANOVA分析，代碼和結(jié)果如下：
dune是轉(zhuǎn)置后的物種豐度表 (抽平或相對比例都行)Management是dune.env中的列名字，代表一列信息，可以是任意樣品屬性信息或分組信息permutations設(shè)置置換次數(shù)method指定距離計算方法R2值顯示Management可以解釋總體差異的34.2%，且P<0.05，表示不同的管理風(fēng)格下的物種組成差異顯著。當(dāng)然還有
65.8%的差異是其它因素造成的。這通常是我們對
PcOA等降維圖標(biāo)記統(tǒng)計檢驗P值的常用方式。
注意：因為是隨機置換，在未指定隨機數(shù)種子時，每次執(zhí)行的結(jié)果都會略有不同，但通常對結(jié)論沒有影響。
# 基于bray-curtis距離進行計算
dune.div <- adonis2(dune ~ Management, data = dune.env, permutations = 999, method="bray")
dune.div
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dune ~ Management, data = dune.env, permutations = 999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Management 3 1.4686 0.34161 2.7672 0.004 **
## Residual 16 2.8304 0.65839
## Total 19 4.2990 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1注意：因為是隨機置換，在未指定隨機數(shù)種子時，每次執(zhí)行的結(jié)果都會略有不同，但通常對結(jié)論沒有影響。也可以如下設(shè)置隨機數(shù)種子，則結(jié)果穩(wěn)定。
# 基于bray-curtis距離進行計算
set.seed(1)
dune.div <- adonis2(dune ~ Management, data = dune.env, permutations = 999, method="bray")
dune.div
## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = dune ~ Management, data = dune.env, permutations = 999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## Management 3 1.4686 0.34161 2.7672 0.002 **
## Residual 16 2.8304 0.65839
## Total 19 4.2990 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1把統(tǒng)計檢驗結(jié)果加到PcOA的圖上。
dune_adonis <- paste0("adonis R2: ",round(dune.div$R2,2), "; P-value: ", dune.div$`Pr(>F)`)
# install.packages("ggalt")
library(ggalt)
ggplot(dune_pcoa_result, aes(x=PCoA1, y=PCoA2, color=Management, group = Management)) +
labs(x=paste("PCoA 1 (", eig_percent[1], "%)", sep=""),
y=paste("PCoA 2 (", eig_percent[2], "%)", sep=""),
title=dune_adonis) +
geom_point(size=5) +
geom_encircle(aes(fill=Management), alpha = 0.1, show.legend = F) +
theme_classic() + coord_fixed(1)
整體有差異了，后面就看看具體那兩組之間有差異，哪兩組之間無差異~~~
往期精品(點擊圖片直達文字對應(yīng)教程)
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
后臺回復(fù)“生信寶典福利第一波”或點擊閱讀原文獲取教程合集
