lyao222lll的个人博客分享 http://blog.sciencenet.cn/u/lyao222lll

博文

R语言绘制PCA、RDA等排序图的一些示例

已有 7517 次阅读 2019-12-8 13:41 |个人分类:R作图|系统分类:科研笔记| _R语言, _统计图表, _排序分析

R语言绘制PCA、RDA等排序图的一些示例

 

前文已分享了PCACADCAPCoANMDSRDAdb-RDACCA等排序方法在R语言中的计算过程。对于排序图的绘制,用于计算的R包中一般会提供可视化函数,如vegan包的ordiplot()plot.cca()等,ade4scatter()等,便于我们在计算后快速观测数据特征。如果有必要,之后可再细调一些参数,让排序图更美观些。

排序图的本质是散点图,所以更普遍的做图方法,就是将排序坐标导出,然后用专业的作图包(如ggplot2)绘制,实现更好的可视化效果。本篇展示一些ggplot2绘制二维排序图的示例,以及再列举一些可绘制三维排序图的R包。

 

示例文件、R脚本等的百度盘链接:

https://pan.baidu.com/s/1KUA5owf4y13y1RJbff-l7g

pca_site.txt”是某PCA分析的样方(对象)排序坐标,“pca_var.txt”是变量的坐标。

cap_site.txt”是某db-RDACAP)分析的样方(对象)排序坐标,“cap_sp.txt”是物种(响应变量)排序坐标,“cap_env.txt”是环境变量(解释变量)的排序坐标。

 


二维散点图,ggplot2的一些作图示例


二维散点图只展示了ggplot2的方法,我觉得下面的这些示例,帮助掌握作图方法应该足够了。文献中大多也是ggplot2的风格,其它情况下,参考文献中的样式仿着来就可以了。

 

PCA等排序图的绘制


先以某PCA排序图为例吧。

library(ggplot2)
 
#数据
pca_site <- read.delim('pca_site.txt', sep = '\t', stringsAsFactors = FALSE)
pca_site$group_all <- paste(pca_site$group1, pca_site$group2, sep = '_')
 
pca_var <- read.delim('pca_var.txt', sep = '\t', stringsAsFactors = FALSE)

 

常规的只展示对象的排序图,PCA前两轴,颜色区分分组,这里假设原分析中前两轴的贡献率分别为30%20%

#常规仅样方的点图,单组样式
p <- ggplot(data = pca_site, aes(x = PC1, y = PC2)) +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.key = element_rect(fill = 'transparent')) +
labs(x = 'PCA1: 30%', y = 'PCA2: 20%')
 
p + geom_point(aes(color = group1)) +
scale_color_manual(values = c('red2', 'purple2', 'green3'))

1.png

 

如果存在多分组,可以分别以更多的颜色区分,或者点的形状等。

#多组可以颜色表示
p + geom_point(aes(color = group_all)) +
scale_color_manual(values = c('red', 'purple', 'green', 'red4', 'purple4', 'green4'), limits = c('A_l', 'B_l', 'C_l', 'A_h', 'B_h', 'C_h'))
 
#或者使用点的颜色+形状区分
p + geom_point(aes(color = group1, shape = group2)) +
scale_color_manual(values = c('red2', 'purple2', 'green3'), limits = c('A', 'B', 'C')) +
scale_shape_manual(values = c(17, 15), limits = c('l', 'h'))
 
#渐变色表示连续数值
p + geom_point(aes(color = group3, shape = group1)) +
scale_shape_manual(values = c(19, 17, 15), limits = c('A', 'B', 'C')) +
scale_color_gradientn(colors = colorRampPalette(c('#C0D857', '#3EB897', '#663490'))(7))

2.png

 

一些背景网格线的常选样式。

#以上述某样式继续为例展示
p <- p + geom_point(aes(color = group1, shape = group2)) +
scale_color_manual(values = c('red2', 'purple2', 'green3'), limits = c('A', 'B', 'C')) +
scale_shape_manual(values = c(17, 15), limits = c('l', 'h'))
 
#一些常见背景样式
p + theme(axis.line = element_line(colour = 'black'), panel.background = element_blank())
 
p + geom_vline(xintercept = 0, color = 'gray', size = 0.5) +
geom_hline(yintercept = 0, color = 'gray', size = 0.5)
 
p + theme(panel.grid.major = element_line(color = 'gray', size = 0.2))

3.png

 

有时能够看到直接在原图中标记对象名称的。

#添加标签,常规方法
p + geom_text(aes(label = name, color = group1), size = 3, vjust = 2, show.legend = FALSE)
 
#防止重叠的添加标签方法
library(ggrepel)
 
p + geom_text_repel(aes(label = name, color = group1), size = 3, box.padding = unit(0.5, 'lines'), show.legend = FALSE)
 
p + geom_label_repel(aes(label = name, color = group1), size = 3, box.padding = unit(0.5, 'lines'), show.legend = FALSE)

4.png

 

这种在图中,根据样本添加分组椭圆或者多边形的样式,经常能看到。

#添加置信椭圆,注意它不是聚类
p + stat_ellipse(aes(color = group1), level = 0.95, linetype = 2, show.legend = FALSE)
 
p + stat_ellipse(aes(fill = group1), geom = 'polygon', level = 0.95, alpha = 0.1, show.legend = FALSE) +
scale_fill_manual(values = c('red', 'purple', 'green'))
 
#“伪聚类”,我也不知道怎么称呼,既不是真正的聚类也不是置信区间,仅用于标记分组而已
#这种样式也经常看到吧
#参考自:https://stackoverflow.com/questions/42575769/how-to-modify-the-backgroup-color-of-label-in-the-multiple-ggproto-using-ggplot2
source('geom_enterotype.r')  #见网盘附件
 
p + geom_enterotype(aes(fill = group1, color = group1, label = group1), show.legend = FALSE) +
scale_fill_manual(values = c('#ffa6a9', '#e8b4fd', '#c7ffc4'))
 
#还有分组多边形这种
library(plyr)
 
group1_border <- ddply(pca_site, 'group1', function(df) df[chull(df[[2]], df[[3]]), ])
 
p + geom_polygon(data = group1_border, aes(fill = group1, color = group1), alpha = 0.1, show.legend = FALSE) +
scale_fill_manual(values = c('red', 'purple', 'green'))

5.png

 

图中额外添加一些标识符、标签或者特殊形状之类的,一般推荐后续用AIPSP图上去更为方便;如果用代码添加的话,以一些文字为例。

#额外添加一些标签等,例如
p + stat_ellipse(aes(fill = group1), geom = 'polygon', level = 0.95, alpha = 0.1, show.legend = FALSE) +
scale_fill_manual(values = c('red', 'purple', 'green')) +
annotate('text', label = 'A', x = -0.15, y = -0.10, size = 5, color = 'red2') +
annotate('text', label = 'B', x = 0.20, y = 0.09, size = 5, color = 'purple2') +
annotate('text', label = 'C', x = -0.20, y = -0.05, size = 5, color = 'green3') +
annotate('text', label = 'PERMANOVA: P<0.001', x = 0.2, y = 0.14, size = 3)

6.png

 

以上均展示了对象的位置,平常见到的大多数PCA排序图也只展示了对象;如果期望将变量也呈现出来。

#添加变量,PCA 这种线性模型中,变量常以向量形式展示
p + geom_segment(data = pca_var, aes(x = 0, y = 0, xend = PC1, yend = PC2), arrow = arrow(length = unit(0.1, 'cm')), color = 'blue', size = 0.3) +
geom_text(data = pca_var, aes(x = PC1 * 1.1, y = PC2 * 1.1, label = name), color = 'blue', size = 3)
 
#CA 这种单峰模型中,变量则常直接展示为点,这里顺便以该示例数据作个展示,假设它是个 CA 的数据
p + geom_text(data = pca_var, aes(x = PC1, y = PC2, label = name), color = 'blue', size = 3) +
labs(x = 'CA1', y = 'CA2')

7.png

 

组合一下上述提到的方法,随便来一个。

#象征性来个组合样式的
p <- ggplot(data = pca_site, aes(x = PC1, y = PC2)) +
geom_point(aes(color = group1, shape = group2)) +
geom_enterotype(aes(fill = group1, color = group1, label = group1), show.legend = FALSE) +
scale_color_manual(values = c('red2', 'purple2', 'green3'), limits = c('A', 'B', 'C')) +
scale_shape_manual(values = c(17, 15), limits = c('l', 'h')) +
scale_fill_manual(values = c('#ffa6a9', '#e8b4fd', '#c7ffc4')) +
guides(color = FALSE) +
theme(panel.grid.major = element_line(color = 'gray', size = 0.2), panel.background = element_rect(color = 'black', fill = 'transparent'), 
        legend.title = element_blank(), legend.key = element_rect(fill = 'transparent'), legend.position = c(0.92, 0.87), legend.background= element_blank()) +
labs(x = 'PCA1: 30%', y = 'PCA2: 20%') +
annotate('text', label = 'PERMANOVA: P<0.001', x = 0.2, y = 0.14, size = 3)
 
p

8.png

 

以及偶尔也能看到和其它图形组合在一起的,例如结合箱线图表示对象的分布更为直观。

#在图的两侧添加箱线图展示点的分布
pc1_boxplot <- ggplot(data = pca_site, aes(x = group1, y = PC1)) +
geom_boxplot(aes(color = group1)) +
scale_color_manual(values = c('red2', 'purple2', 'green3'), limits = c('A', 'B', 'C')) +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), 
        legend.position= 'none', axis.title = element_text(color = 'transparent')) +
coord_flip()
 
pc2_boxplot <- ggplot(data = pca_site, aes(x = group1, y = PC2)) +
geom_boxplot(aes(color = group1)) +
scale_color_manual(values = c('red2', 'purple2', 'green3'), limits = c('A', 'B', 'C')) +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), 
        legend.position= 'none', axis.title = element_text(color = 'transparent'))
 
library(grid)
 
grid.newpage()
pushViewport(viewport(layout = grid.layout(3, 3)))
print(p, vp = viewport(layout.pos.row = 1:2, layout.pos.col = 1:2))
print(pc1_boxplot, vp = viewport(layout.pos.row = 3, layout.pos.col = 1:2))
print(pc2_boxplot, vp = viewport(layout.pos.row = 1:2, layout.pos.col = 3))
 
 #这里用了 grid 包实现组合,还有很多 R 包可用于组合,还请自行了解了
#如果觉得对齐效果不太好,在 print() 函数里调试参数对齐
#或者后面用 AI、PS 调整下更方便

9.png

 

RDA等排序图的绘制


这里以某db-RDA(CAP)的计算坐标为例吧。

基本的作图方法上文已作了描述,这里直接先展示一例只包含样方和环境变量的双序图,假设原分析中前两个约束轴的解释率分别为25%15%

library(ggplot2)
 
#数据
cap_site <- read.delim('cap_site.txt', sep = '\t', stringsAsFactors = FALSE)
cap_sp <- read.delim('cap_sp.txt', sep = '\t', stringsAsFactors = FALSE)
cap_env <- read.delim('cap_env.txt', sep = '\t', stringsAsFactors = FALSE)
 
#作图,双序图
p <- ggplot(data = cap_site, aes(x = CAP1, y = CAP2)) +
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 = c('red2', 'purple2', 'green3'), limits = c('A', 'B', 'C')) +
scale_shape_manual(values = c(17, 15), limits = c('l', 'h')) +
scale_fill_manual(values = c('red', 'purple', 'green')) +
theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), 
        legend.title = element_blank(), legend.key = element_rect(fill = 'transparent')) +
geom_vline(xintercept = 0, color = 'gray', size = 0.5) +
geom_hline(yintercept = 0, color = 'gray', size = 0.5) +
geom_segment(data = cap_env, aes(x = 0, y = 0, xend = CAP1, yend = CAP2), arrow = arrow(length = unit(0.2, 'cm')), color = 'blue', size = 0.3) +
geom_text(data = cap_env, aes(x = CAP1 * 1.1, y = CAP2 * 1.1, label = name), color = 'blue', size = 3) +
labs(x = 'CAP1: 25%', y = 'CAP2: 15%')
 
p

10.png

 

带样方、物种和环境变量的三序图。

#考虑将物种变量加入,展示为三序图
p + geom_segment(data = cap_sp, aes(x = 0, y = 0, xend = CAP1, yend = CAP2), arrow = arrow(length = unit(0.1, 'cm')), color = 'orange', size = 0.2) +
geom_text(data = cap_sp, aes(x = CAP1 * 1.1, y = CAP2 * 1.1, label = name), color = 'orange', size = 2)
 
#物种也用颜色区分下分类
#物种数量太多了,这里选部分物种演示下好了
cap_sp_select <- cap_sp[1:10, ]
 
p + geom_segment(data = cap_sp_select, aes(x = 0, y = 0, xend = CAP1, yend = CAP2, color = sp_class), 
        arrow = arrow(length = unit(0.1, 'cm')), size = 0.2, show.legend = FALSE) +
geom_text(data = cap_sp_select, aes(x = CAP1 * 1.1, y = CAP2 * 1.1, label = name, color = sp_class), size = 2) +
scale_color_manual(values = c('red2', 'purple2', 'green3', '#8DD3C7', '#BEBADA', '#FB8072', '#80B1D3'), 
        limits = c('A', 'B', 'C', 'pa', 'pb', 'pr', 'a'))
 
#物种直接以点展示,而非向量
p + geom_point(data = cap_sp_select, aes(x = CAP1, y = CAP2), color = 'orange', size = 1)
 
p + geom_text(data = cap_sp_select, aes(x = CAP1, y = CAP2, label = name), color = 'orange', size = 2)

11.png

 


三维散点图,一些可选R


二维点图简单,使用ggplot2作图挺方便的。

三维点图表示难画3.gif,就不整这么细了,下面列几种常用的R包,实际使用时看看帮助文档调调参数啥的,或再配合AIPS等看着修一修。

以下还是以某PCA结果为例吧。

 

scatterplot3d包绘制三维散点图。

pca_site.txtgroup1区分颜色,按group2区分点的形状。

#数据
pca_site <- read.delim('pca_site.txt', sep = '\t', stringsAsFactors = FALSE)
pca_var <- read.delim('pca_var.txt', sep = '\t', stringsAsFactors = FALSE)
 
#scatterplot3d 包,详情 ?scatterplot3d
library(scatterplot3d)
 
scatterplot3d(pca_site$PC1, pca_site$PC2, pca_site$PC3, color = rep(c('red2', 'purple2', 'green3'), c(8, 8, 8)), pch = rep(rep(c(17, 15), c(4, 4)), 3), 
        xlab = paste('PCA1: 20%'), ylab = paste('PCA2: 12%'), zlab = paste('PCA3: 8%'))
 
scatterplot3d(pca_site[2:4], color = rep(c('red2', 'purple2', 'green3'), c(8, 8, 8)), lwd = 3, type = 'h',
        xlab = paste('PCA1: 20%'), ylab = paste('PCA2: 12%'), zlab = paste('PCA3: 8%'))

12.png

 

rgl包绘制三维散点图,交互式风格便于观测数据。

pca_site.txtgroup1区分颜色

#rgl 包,交互式,详情 ?plot3d
library(rgl)
 
plot3d(pca_site[2:4], col = rep(c('red2', 'purple2', 'green3'), c(8, 8, 8)), size = 5)
 
#参考自:https://stackoverflow.com/questions/24282143/pca-multiplot-in-r
#有其它需要(如颜色选择等),自定义修改即可
library(rgl)
 
plotPCA <- function(x, nGroup) {
    n <- ncol(x) 
    if(!(n %in% c(2,3))) { # check if 2d or 3d
        stop("x must have either 2 or 3 columns")
    }
 
    fit <- hclust(dist(x), method="complete") # cluster
    groups <- cutree(fit, k=nGroup)
 
    if(n == 3) { # 3d plot
        plot3d(x, col=groups, type="s", size=1, axes=F)
        axes3d(edges=c("x--", "y--", "z"), lwd=3, axes.len=2, labels=FALSE)
        grid3d("x")
        grid3d("y")
        grid3d("z")
    } else { # 2d plot
        maxes <- apply(abs(x), 2, max)
        rangeX <- c(-maxes[1], maxes[1])
        rangeY <- c(-maxes[2], maxes[2])
        plot(x, col=groups, pch=19, xlab=colnames(x)[1], ylab=colnames(x)[2], xlim=rangeX, ylim=rangeY)
        lines(c(0,0), rangeX*2)
        lines(rangeY*2, c(0,0))
    }
}
 
plotPCA(pca_site[2:4], nGroup = 3)

13.png

 

car包结合rgl包绘制三维散点图,交互式风格便于观测数据。

pca_site.txtgroup1区分颜色,同时按group1添加分组椭圆(注意这个椭圆不是置信区间,仅为按分组拟合的曲面)

#car+rgl 包,交互式,详情 ?scatter3d
library(car)
library(rgl)
 
scatter3d(pca_site$PC1, pca_site$PC2, pca_site$PC3, surface = FALSE, point.col = rep(c('red2', 'purple2', 'green3'), c(8, 8, 8)), 
        xlab = paste('PCA1: 20%'), ylab = paste('PCA2: 12%'), zlab = paste('PCA3: 8%'))
 
scatter3d(pca_site$PC1, pca_site$PC2, pca_site$PC3, groups = as.factor(pca_site$group1), surface = FALSE, ellipsoid = TRUE, 
        xlab = paste('PCA1: 20%'), ylab = paste('PCA2: 12%'), zlab = paste('PCA3: 8%'))

14.png

 

还有广泛熟知的plotly包绘制三维散点图,交互式网页的风格便于观测数据。

pca_site.txtgroup1区分颜色,以及添加变量向量等。

#plotly 包,网页交互式,详情 ?plot_ly
library(plotly)
 
pca_site$color <- rep(c('red2', 'purple2', 'green3'), c(8, 8, 8))
plot_ly(pca_site, labels = ~name, x = ~PC1, y = ~PC2, z = ~PC3, color = ~color)
 
#添加变量向量
#参考自:https://stackoverflow.com/questions/44393823/3d-biplot-in-plotly-r
library(plotly)
 
p <- plot_ly() %>%
  add_trace(x = pca_site$PC1, y = pca_site$PC2, z = pca_site$PC3,
            type="scatter3d", mode="markers",
            marker = list(color=rep(c(1,2,3), each = 8), opacity = 0.7)) 
 
for (k in 1:nrow(pca_var)) {
   p <- p %>% add_trace(x=c(0, pca_var[k,'PC1']), y=c(0, pca_var[k,'PC2']), z=c(0, pca_var[k,'PC3']),
      type="scatter3d", mode="lines",
      line = list(width=3),
      opacity = 1) 
}
 
p

15.png

 



http://blog.sciencenet.cn/blog-3406804-1209233.html

上一篇:R包vegan的典范对应分析(CCA)
下一篇:R包vegan的Mantel tests

1 许亚东

该博文允许注册用户评论 请点击登录 评论 (2 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2020-8-8 20:42

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部