跟着NatureCommunication学数据分析：R语言相对丰度数据主坐标分析（PcoA)

otu<-read.delim("data/20220602/dat01.txt",
                sep="\t",
                header = TRUE,
                row.names = 1,
                check.names = FALSE,
                stringsAsFactors = FALSE)
dim(otu)
head(otu)

otu <- data.frame(t(otu))
otu[1:6,1:6]

group<-read.delim("data/20220602/dat02.txt",
                  header=TRUE,
                  sep="\t",
                  stringsAsFactors = FALSE)
head(group)

library(tidyverse)

group %>% 
  mutate(Site=case_when(
    Group == "Slope" ~ "None-CD",
    TRUE ~ "CD"
  ),
  high=case_when(
    `Depth (m)`< 6000 ~ '5k',
    #`Depth (m)`>=6000 & `Depth (m)` < 7000 ~ '6k',
    `Depth (m)`>=6000 & `Depth (m)` < 8000 ~ '7k',
    `Depth (m)`>=8000 & `Depth (m)` < 9000 ~ '8k',
    `Depth (m)`>=9000 & `Depth (m)` < 10000 ~ '9k',
    `Depth (m)`>=10000 ~ '10k'
  ),
  position=case_when(
    `Depth (m)`< 8000 ~ 'North',
    `Depth (m)`>= 8000 & `Depth (m)`< 10000 ~ 'South',
    TRUE ~ 'axis'
  )) -> new.group

library(vegan)
distance <- vegdist(otu, method = 'bray')
pcoa <- cmdscale(distance, k = (nrow(otu) - 1), eig = T)


ordiplot(scores(pcoa)[ ,c(1, 2)], type = 't')

summary(pcoa)

point <- data.frame(pcoa$point)
point %>% head()

species <- wascores(pcoa$points[,1:2], otu)
species %>% head()
pcoa_eig <- (pcoa$eig)[1:2] / sum(pcoa$eig)
pcoa_eig

sample_site <- data.frame({pcoa$point})[1:2]
sample_site$Sample <- rownames(sample_site)
names(sample_site)[1:2] <- c('PCoA1', 'PCoA2')
sample_site <- merge(sample_site, new.group, by = 'Sample', all.x = T)
sample_site %>% head()

sample_site$Site <- factor(sample_site$Site, 
                           levels = c('None-CD', 'CD'))
sample_site$high <- factor(sample_site$high, 
                           levels = c('5k', '7k', '8k', '9k', '10k'))
sample_site$Position <- factor(sample_site$position, 
                               levels = c('North', 'South', 'axis'))

group_border <- plyr::ddply(sample_site, 'Site', 
                      function(df) df[chull(df[[2]], df[[3]]), ]) 

ggplot(sample_site, aes(PCoA1, PCoA2, group = Site)) +
  theme(panel.grid = element_line(color = 'gray', 
                                  linetype = 2, size = 0.1), 
        panel.background = element_rect(color = 'black', 
                                        fill = 'transparent'), 
        legend.key = element_rect(fill = 'transparent')) + #去掉背景框
  geom_vline(xintercept = 0, color = 'gray', size = 0.4) + 
  geom_hline(yintercept = 0, color = 'gray', size = 0.4) +
  geom_polygon(data = group_border, aes(fill = Site),alpha=0.1) +
  guides(fill = guide_legend(order = 1), 
         shape = guide_legend(order = 2), 
         color = guide_legend(order = 3)) +
  scale_shape_manual(values = c(17, 16,15,12,10)) + 
  geom_point(aes(color = position, shape = high), 
             size = 2.5, alpha = 0.8) + 
  labs(x = paste('PCoA axis1: ', 
                 round(100 * pcoa_eig[1], 2), '%'),
       y = paste('PCoA axis2: ', 
                 round(100 * pcoa_eig[2], 2), '%')) +
  annotate('text', label = 'Slope', 
           x = -0.31, y = -0.15, 
           size = 5, colour = '#C673FF') +
  annotate('text', label = 'Bottom-axis', 
           x = 0.32, y = 0.05, 
           size = 5, colour = '#FF985C')