解決疑難：排序分析三個重複的置信橢圓

寫在前麪

排序圖添加一個置信橢圓可以很好的對觀察不同分組微生物數據的變化槼律。

實戰

對ps對象取子集，每個分組衹畱下三個樣本，用於測試三個重複添加置信橢圓的函數是否可靠

library(vegan)
library(tidyverse)
library(ellipse)
library(ggClusterNet)
library(phyloseq)

ps =  ps %>% 
subset_samples.wt("ID",c("OE1","OE2","OE3","WT1","WT2","WT3","KO1","KO2","KO3"))

這裡對數據進行一個pcoa排序

otu = ps %>%
  vegan_otu() %>% t() %>%
as.data.frame()
map = sample_data(ps) 
head(map)

bray_curtis =vegan:: vegdist(t(otu),method ="bray", na.rm=TRUE)
pcoa <- otu %>% 
  t() %>%
  vegan:: vegdist(method ="bray", na.rm=TRUE) %>%
cmdscale(k=2,eig=T)

points <- pcoa$points %>% as.data.frame() %>%
  # as.tibble() %>%
  dplyr::rename(.,x ="V1",y ="V2" )
head(points)

eig = pcoa$eig
eig
points = cbind(points, map[match(rownames(points), map$ID), ])
ado = adonis(bray_curtis~ map$Group,permutations = 999,method="bray")
a = round(as.data.frame(ado$aov.tab[5])[1,1],3)
R2 <- paste("adonis:R",a, sep ="")
b = as.data.frame(ado$aov.tab[6])[1,1]
p_v = paste("p:",b, sep ="")
title = paste(R2,"",p_v, sep ="")
title

設定置信橢圓範圍爲0.95

level = 0.95
head(points)
centroids <- aggregate(cbind(x,y)~ Group,points,mean)
head(centroids)
row.names(centroids) = centroids$Group

conf.rgn  <- do.call(rbind,lapply(unique(points$Group),function(t)
data.frame(Group=as.character(t),
ellipse::ellipse(cov(points[points$Group==t,1:2]),
centre=as.matrix(centroids[t,2:3]),
level=level),
stringsAsFactors=FALSE)))

p4 <- ggplot(points, aes(x=x, y=y, fill = Group))  
  geom_point(alpha=.7, size=5, pch = 21)  
  labs(x=paste("PCoA 1 (", format(100 * eig[1] / sum(eig), digits=4),"%)", sep=""),
       y=paste("PCoA 2 (", format(100 * eig[2] / sum(eig), digits=4),"%)", sep=""),
title=title)
  geom_path(data=conf.rgn,linetype = 2,aes(color = Group))  
theme_classic()
p4

設定置信區間爲0.7；

level = 0.7
head(points)
centroids <- aggregate(cbind(x,y)~ Group,points,mean)
head(centroids)
row.names(centroids) = centroids$Group

conf.rgn  <- do.call(rbind,lapply(unique(points$Group),function(t)
data.frame(Group=as.character(t),
ellipse::ellipse(cov(points[points$Group==t,1:2]),
centre=as.matrix(centroids[t,2:3]),
level=level),
stringsAsFactors=FALSE)))

p4 <- ggplot(points, aes(x=x, y=y, fill = Group))  
  geom_point(alpha=.7, size=5, pch = 21)  
  labs(x=paste("PCoA 1 (", format(100 * eig[1] / sum(eig), digits=4),"%)", sep=""),
       y=paste("PCoA 2 (", format(100 * eig[2] / sum(eig), digits=4),"%)", sep=""),
title=title)
  geom_path(data=conf.rgn,linetype = 2,aes(color = Group))  
theme_classic()
p4
根際互作生物學研究室 簡介