跟着Nature学数据分析：R语言vegan包做主坐标分析（PCoA）并使用ggplot2可视化结果

2022-05-23 15:45:48 浏览数 (1)

论文是

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这个数据集是表型数据，行是样本，列是变量

读取数据集

代码语言：javascript复制
inMBss<-read.table("20220420/nature/inMBss.tsv",
                   row.names = 1,
                   header=T)
dim(inMBss)

使用vegan包做PCoA分析

代码语言：javascript复制
library(vegan)
distMatrix <- vegdist(inMBsss,method = "bray")
pCoa <- cmdscale(distMatrix, eig = T,k = 2 )

这里计算距离的时候会遇到警告信息

他这个是模拟数据，我们可以忽略这个警告，但是真实数据如果遇到这个情况应该怎么处理呢？欢迎大家留言讨论

计算每个坐标的贡献

代码语言：javascript复制
varExp <- (eigenvals(pCoa)/sum(eigenvals(pCoa)))[1:2]
xVar <- as.numeric(varExp[1]*100)
yVar <- as.numeric(varExp[2]*100)

整理作图数据

代码语言：javascript复制
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绘图

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

下面是用来制作封面图，与本文内容无关

代码语言：javascript复制
library(patchwork)
g g plot_layout(guide="collect")

image.png

下期推文介绍如何将主坐标分析的结果关联表型数据

https github typescript 网络安全 git

0 人点赞