前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >跟着NatureCommunication学数据分析:R语言相对丰度数据主坐标分析(PcoA)

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

作者头像
用户7010445
发布2023-01-06 18:35:10
7730
发布2023-01-06 18:35:10
举报

论文

Microbiomes in the Challenger Deep slope and bottom-axis sediments

https://www.nature.com/articles/s41467-022-29144-4#code-availability

对应代码链接

https://github.com/ucassee/Challenger-Deep-Microbes

论文里提供了大部分图的数据和代码,很好的学习材料,感兴趣的同学可以找来参考,今天的推文重复一下论文中的Figure2b

image.png

部分数据集截图如下

相对丰度数据

image.png

分组数据

image.png

读取数据集

读取相对丰度数据

代码语言:javascript
复制
otu<-read.delim("data/20220602/dat01.txt",
                sep="\t",
                header = TRUE,
                row.names = 1,
                check.names = FALSE,
                stringsAsFactors = FALSE)
dim(otu)
head(otu)

对数据转置

代码语言:javascript
复制
otu <- data.frame(t(otu))
otu[1:6,1:6]

读取分组数据

代码语言:javascript
复制
group<-read.delim("data/20220602/dat02.txt",
                  header=TRUE,
                  sep="\t",
                  stringsAsFactors = FALSE)
head(group)

这个分组数据和论文中提供的代码的分组信息还少一些内容,我们再给它增加几列

代码语言:javascript
复制
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

这个分组信息可能和原文中有差别

主坐标分析代码

代码语言:javascript
复制
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)

构造作图数据

代码语言:javascript
复制
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()

给准备好的数据赋予因子水平

代码语言:javascript
复制
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'))

构造画边界的数据

代码语言:javascript
复制
group_border <- plyr::ddply(sample_site, 'Site', 
                      function(df) df[chull(df[[2]], df[[3]]), ]) 

画图代码

代码语言:javascript
复制
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')

论文中提供的代码出图效果如下

image.png

这个图和最终论文中的图还是有些差别的,主要是图例的位置和边框,如何用代码把图例的位置调整到右下角并添加一个边框,这个另起推文来介绍吧

本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2022-06-02,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 小明的数据分析笔记本 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 论文
  • 对应代码链接
  • 部分数据集截图如下
  • 读取数据集
  • 主坐标分析代码
  • 构造作图数据
  • 给准备好的数据赋予因子水平
  • 构造画边界的数据
  • 画图代码
相关产品与服务
图数据库 KonisGraph
图数据库 KonisGraph(TencentDB for KonisGraph)是一种云端图数据库服务,基于腾讯在海量图数据上的实践经验,提供一站式海量图数据存储、管理、实时查询、计算、可视化分析能力;KonisGraph 支持属性图模型和 TinkerPop Gremlin 查询语言,能够帮助用户快速完成对图数据的建模、查询和可视化分析。
领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档