R 语言|绘制主坐标图以分析 PCoA 图
最编程
2024-03-31 17:55:48
...
主坐标分析(Principal Coordinates Analysis,PCoA)是一种研究数据相似性或差异性的可视化方法,把复杂的数据用一系列的特征值和特征向量进行排序后,选择排在前面的几位特征值,用来表示样品之间的关系并以坐标的形式展现,结果是数据矩阵的旋转,通过PCA 可以找到距离矩阵中最主要的坐标,观察个体或者群体之间的差异。
PCoA分析比较适用样本数目比较少,而物种数目比较多的情况。
小编今天给大家分享R包vegan绘制主坐标分析PCoA图。
具体步骤如下:
1.加载vegan和ggplot2包;
#导入作图包
library(vegan)
library(ggplot2)
#颜色选择
color=c( "#3C5488B2","#00A087B2",
"#F39B7FB2","#91D1C2B2",
"#8491B4B2", "#DC0000B2",
"#7E6148B2","yellow",
"darkolivegreen1", "lightskyblue",
"darkgreen", "deeppink", "khaki2",
"firebrick", "brown1", "darkorange1",
"cyan1", "royalblue4", "darksalmon",
"darkgoldenrod1", "darkseagreen", "darkorchid")
2.读取otu数据文件,计算样方间距离,并读入此距离测度作为PCoA的输入。
(小编这里没有现成的距离文件,所以需要先根据OTU的丰度组成来计算样方间距离)
#读取otu数据文件
otu <- read.delim('otu_table.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
otu <- data.frame(t(otu))
#根据物种组成计算样方距离,结果以 dist 数据类型存储
bray_dis <- vegdist(otu, method = 'bray')
#样方排序坐标
pcoa <- cmdscale(bray_dis, k = (nrow(otu) - 1), eig = TRUE)
site <- data.frame(pcoa$point)[1:2]
site$name <- rownames(site)
3.读取分组文件,并将分组文件和距离文件合并;
#读取分组文件
map<-read.table('group.txt',header=T,sep="\t",row.names=1)
#site$group <- c(rep('A', 12), rep('B', 12), rep('C', 12))
#将分组文件和数据文件以行名合并
merged=merge(site,map,by="row.names",all.x=TRUE)
4.使用 wascores() 被动添加物种得分(坐标),丰度加权平均方法;
#使用 wascores() 被动添加物种得分(坐标),丰度加权平均方法
species <- wascores(pcoa$points[,1:4], otu)
5.这里小编主要展示top10 丰度物种;
#所有物种太多了,挑主要的展示,如 top10 丰度物种
#计算 top10 丰度物种
abundance <- apply(otu, 2, sum)
abundance_top10 <- names(abundance[order(abundance, decreasing = TRUE)][1:10])
species_top10 <- data.frame(species[abundance_top10,1:2])
species_top10$name <- rownames(species_top10)
pcoa_exp <- pcoa$eig/sum(pcoa$eig)
pcoa1 <- paste('PCoA axis1 :', round(100*pcoa_exp[1], 2), '%')
pcoa2 <- paste('PCoA axis2 :', round(100*pcoa_exp[2], 2), '%')
6.使用ggplot2包作图;
#ggplot2 作图
library(ggplot2)
p <- ggplot(data = merged, aes(X1, X2)) +
geom_point(aes(color = group)) +
stat_ellipse(aes(fill = group), geom = 'polygon', level = 0.95, alpha = 0.1, show.legend = FALSE) + #添加置信椭圆,注意不是聚???
scale_color_manual(values =color[1:length(unique(map$group))]) +
scale_fill_manual(values =color[1:length(unique(map$group))]) +
theme(panel.grid.major = element_line(color = 'gray', size = 0.2), panel.background = element_rect(color = 'black', fill = 'transparent'),
plot.title = element_text(hjust = 0.5)) +
geom_vline(xintercept = 0, color = 'gray', size = 0.5) +
geom_hline(yintercept = 0, color = 'gray', size = 0.5) +
geom_text(data = species_top10, aes(label = name), color = "lightskyblue", size = 3) + #??? top10 丰度物种标签
labs(x = pcoa1, y = pcoa2, title = 'PCoA top10丰度物种')
p
“作图帮”公众号,免费分享绘图代码及数据,快来关注吧!
上一篇: 基于图形的植物泛基因组
下一篇: C 语言的二叉树遍历
推荐阅读
-
YUV420数据格式详解,以图文形式揭示 - 提示:加强阅读理解,配合图片观看更加易懂
-
AVCodecContext和AVCodec的关系分析
-
分析Opencv中BGR、YUV、YUV_I420和NV12颜色空间
-
深入解析YUV420、YUV420(YUY2)、YUV422(YVYU):以图解方式详解
-
深入探究音视频(1):对YUV颜色空间的初步分析
-
学会分析YUV数据:视频采集与处理入门教程
-
分析YUV数据解码方法
-
C语言内存管理的重要性
-
视频格式的YUV和RGB分析及常见形式
-
10bit YUV(P010)的存储结构和处理-随着计算机处理信息的能力越来越厉害,这种能表现更高动态范围的图像存储格式将会逐渐成为主流,但是现在很多算法都不能直接处理 10bit 的 YUV ,都是先将其转换为 8bit YUV ,然后再进行处理,这实际上是丢弃了 10bit YUV 的图像高动态范围优势。 令人遗憾的是在渲染图像时,目前 OpenGL 也无法直接对 10bit YUV 进行渲染,也是需要先转换为 8bit YUV 。 接下来以一种常见的 10bit YUV (P010) 格式为例,介绍一下 10bit YUV 到 8bit YUV 的转换过程。 P010 最早是微软定义的格式,表示的是 YUV 4:2:0 的采样方式,也就是说 P010 表示的是一类 YUV 格式,它的内存排布方式可能是 NVNVYUYV12 。