跟着NatureCommunication学数据分析:R语言相对丰度数据主坐标分析(PcoA)
2023-02-19 12:27:38 时间
论文
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
读取数据集
读取相对丰度数据
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')
论文中提供的代码出图效果如下
image.png
这个图和最终论文中的图还是有些差别的,主要是图例的位置和边框,如何用代码把图例的位置调整到右下角并添加一个边框,这个另起推文来介绍吧
相关文章
- 学生数据库管理系统
- SpringDataJpa 用MySQL语句怎么分页,spring全家桶SpringDataJpa 用MySQL语句怎么分页
- Docker创建MySQL容器模板命令
- Elasticsearch对应MySQL的对应关系
- 使用SpringDataJpa保存(save)报错误:SQL Error: 1062, SQLState: 23000 控制台会报:Duplicate entry ‘数‘ for key ‘PRIMA
- Navicat Premium 连接sqlserver数据库时提示安装Client失败,解决方案
- Mysql查询当前用户所有数据库语句(SHOW DATABASES)
- MySQL语句-查看当前数据库有哪些表(SHOW TABLES)
- MySQL5.0版本以上新增的 information_schema 数据库是什么?
- MariaDB数据库备份之逻辑备份
- MariaDB数据库创建用户
- MariaDB数据库给用户授权
- MariaDB数据库刷新权限表命令
- MariaDB数据库删除用户命令
- PhpStudy 2016搭建-sqli-libs靶场
- MySQL手动注入步骤
- Pikachu靶场-SQL注入-数字型注入(post)过关步骤
- Pikachu靶场-SQL注入-字符型注入(get)过关步骤
- 利用SQL注入漏洞实现MySQL数据库读写文件
- Kali-工具-sqlmap常见用法