论文是
Environmental factors shaping the gut microbiome in a Dutch population
数据和代码的github主页链接
https://github.com/GRONINGEN-MICROBIOME-CENTRE/DMP
这个也是数据代码的下载链接，可以看目录结构
https://zenodo.org/record/5910709#.YmAcp4VBzic
论文没有提供真实数据，而是提供了一个模拟数据，可能是考虑到数据量比较大，后面的分析只用了模拟数据的前两千行。
论文里也提供了数据预处理的代码，这里我就不展示了，我直接用处理好的数据
两个数据集
inMBss.tsv这个数据行是样本，列是物种，但是数值代表的是什么意思我没有看懂phe.tsv这个数据集是表型数据，行是样本，列是变量
读取数据集
inMBss<-read.table("20220420/nature/inMBss.tsv",
row.names = 1,
header=T)
dim(inMBss)
使用vegan包做PCoA分析
library(vegan)
distMatrix <- vegdist(inMBsss,method = "bray")
pCoa <- cmdscale(distMatrix, eig = T,k = 2 )
这里计算距离的时候会遇到警告信息
image.png
他这个是模拟数据，我们可以忽略这个警告，但是真实数据如果遇到这个情况应该怎么处理呢？欢迎大家留言讨论
计算每个坐标的贡献
varExp <- (eigenvals(pCoa)/sum(eigenvals(pCoa)))[1:2]
xVar <- as.numeric(varExp[1]*100)
yVar <- as.numeric(varExp[2]*100)
整理作图数据
pCoaVecs <- as.data.frame(pCoa$points)
colnames(pCoaVecs) <- paste0("PCo",c(1:2))
pCoaVecs$ID <- row.names(pCoaVecs)
inMBssID <- inMBsss
inMBssID$ID <- row.names(inMBssID)
inMBssPrev <- inMBssID[,c("ID","k__Bacteria.p__Bacteroidetes.c__Bacteroidia.o__Bacteroidales.f__Prevotellaceae.g__Prevotella.s__Prevotella_copri")]
colnames(inMBssPrev) <- c("ID","P.copri")
pCoaVecsTaxa <- merge(pCoaVecs,inMBssPrev,by="ID")
这里的代码他为什么挑P.copri这个种来映射颜色我没有看懂，还需要仔细看原文
ggplot2绘图
library(ggplot2)
g <- ggplot(pCoaVecsTaxa,aes(x=PCo1,y=PCo2,color=P.copri)) +
geom_point(size=1.25) +
theme_classic() +
scale_color_gradientn(colours = rainbow(5)) +
xlab(paste0('PCo1 ',round(xVar,2),' %')) +
ylab(paste0('PCo2 ',round(yVar,2),' %')) +
theme(text = element_text(size = 16))
print(g)
image.png
下面是用来制作封面图，与本文内容无关
library(patchwork)
g+g+plot_layout(guide="collect")
image.png
下期推文介绍如何将主坐标分析的结果关联表型数据
本期推文的示例数据和代码可以先
点赞，然后点击在看，最后在公众号后台留言20220421获取
欢迎大家关注我的公众号
小明的数据分析笔记本
小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！
