R统计绘图-使用rgl或pca3D包绘制3DPCA图

虽然PCA和RDA分析及绘图都写过教程，但是对于结果的解释都没有写的很详细，刚好最近有人询问怎样使用FactoMineR factoextra包进行PCA分析。所以使用R统计绘图-环境因子相关性热图中的不同土壤环境因子数据进行PCA绘图和结果解读推文。

一、 PCA分析过程

1.1 数据准备

# 1.1.1 设置工作路径
setwd("D:\\EnvStat\\PCA\\3DPCA")
getwd()#获取工作路径

# 1.1.2 调用R包
library(factoextra)
library(scales)
library(ggsci)
library(pca3d)
library(rgl)

# 1.1.3 导入数据
data = read.csv("env.csv",
                header = TRUE,row.names=1,
                stringsAsFactors = FALSE,
                check.names = FALSE) 
##设置check.names=FALSE，读入数据是不会更改行名和列名中的特殊符号(-、空格等符号)
head(data)
dim(data) # 1-3列为分组信息，4-14列为分析变量信息。

图1|环境因子及分组信息表，env.csv。行为样品名称，列为环境因子名称和分组信息，共有11个环境变量，3个分组信息。

1.2 主成分分析(PCA)

使用R中能进行PCA分析的R函数(vegan::rda(),prcomp()或princomp)就行，pca3D需要prcomp()输出结果，所以这里使用prcomp()进行PCA分析。

## 1.2.1 使用prcomp()进行PCA分析
#library(vegan)
#pca=rda(data[4:14],scale=TRUE)
#pca = prcomp(data[4:14], scale.= TRUE)
pca <- princomp(data[4:14],cor=TRUE,scores=TRUE) # cor = TRUE使用相关性矩阵进行相关性分析。

## 1.2.2 提取样本PC1-3坐标数据
ind = data.frame(get_pca_ind(pca)$coord[,1:3])
colnames(ind)= c("PC1","PC2","PC3")
ind = cbind(ind,data[2:3])
ind
####设置因子水平---基于tillage和depth分组信息绘图
ind$tillage = factor(ind$tillage,levels = c("CK","NT","PT","WL"))
ind$depth =  factor(ind$depth,levels = c("0-10 cm","10-20 cm","20-30 cm"))

## 1.2.3 提取变量PC1-PC3数据
var = data.frame(get_pca_var(pca)$coord[,1:3])
colnames(var)= c("PC1","PC2","PC3")
var

注：这里只展示绘图过程，PCA分析的详细分析及结果解释流程都写在R统计绘图-PCA分析绘图及结果解读(误差线，多边形，双Y轴图、球形检验、KMO和变量筛选等)中。

图2|样本坐标数据及分组信息表,ind。

图3|变量坐标数据表,var。

二、绘制3D PCA图

3D PCA样本图其实就是3D散点图，如果要绘制3D PCA双序图，可以以箭头等其他形式展示变量。绘图原理都是差不多的，本质上都需要提供坐标信息，向图中添加各种需要的元素并设置图中的各种参数美化图形。可以绘制3D图形的R包挺多的，我觉得rgl是最强大的，其他很多可以绘制3D图的函数，很多都是调用的rgl包。另外还有一个pca3d专门用于绘制3D PCA图。这里主要介绍使用rgl或pca3d绘制3D PCA图以及输出pdf格式图。

2.1 pca3d绘制3D PCA图

样本直接用散点显示(样本名称太长，重叠情况严重，所以不添加样本标签)，使用分组信息分别着色。变量使用红色箭头显示。

## 2.1.1 准备绘图颜色
###数据分析时，最好统一绘图结果的色系，所以一般同一篇文章的数据分析结果，我都先根据需要的颜色数量，选择ggsci中的期刊颜色。之后都统一使用这一颜色集。这里选择d3颜色集
cols = pal_d3("category10")(10)
show_col(cols) # 查看颜色的序号、十六进制名称和颜色。可以自己手动选择自己想要的颜色。
### 如果有想要使用的RGB颜色，也可以使用rgb()函数转为十六进制颜色值
### col2rgb()可以将颜色名称或十六进制颜色转换为RGB值。
### grDevices::convertColor() # 可以更灵活的在各个颜色空间进行转换。
### 例：RGB(80,196,143);RGB(54.133,254);SD:RGB(245,97,111)
#col2rgb('blue', alpha=T) # 颜色名称转为RGB值。
#col2tgb("hexadecimal",alpha = TRUE) # 16进制转为RGB值。
#cols = c(rgb(red=80, green=196, blue=143, alpha=255, max=255),
         #rgb(red=54, green=133, blue=254, alpha=255, max=255),
         #rgb(red=245, green=97, blue=111, alpha=255, max=255))
         
## 2.1.2 pca3d绘制3D PCA图
#browseVignettes("pca3d")# 查看帮助手册
#??pca3d

### 默认参数绘图
pca3d(pca) # 图5

### 添加分组信息和自定义形状和颜色
listShapes() # 可以查看pca2d()和pca3d可以使用的形状
pca3d(pca, # 图6
      components = 1:3, # 选择PC1-PC3用于绘图
      group=ind[,4],# 选择的分组信息,不能使用ind[4]会报“Error in matchShapes(shape) : Incorrect shapes: NA”错误
      palette=cols,# 设置色板
      shape ="s", # 需要同时设置形状，不然默认根据分组信息分配不同的形状(tetrahedrons, cubes, spheres和octahedrons)。
      radius = 2 # 类似于cex参数，设置形状大小。
      ) ## 不使用分组信息，向参数palette和shape提供一个与样本数等长的颜色和性状矢量也可指定颜色和形状。

###  添加分组椭圆
pca3d(pca, # 图7
      components = 1:3, 
      group=ind[,4],
      palette=cols, 
      shape ="s",
      show.plane = TRUE, # 在y=0处添加水平平面  
      show.shadows = TRUE, # 绘制棒棒图，点到y=0平面的垂直线和阴影。
      show.ellipses = TRUE, # 添加分组椭圆
      ellipse.ci = 0.95,# 椭圆集合区间水平，默认0.95。
      show.axes = TRUE,
      show.scale =TRUE, # 显示轴刻度
      show.centroids=TRUE, # 显示聚类质心
      show.group.labels=TRUE,# 显示分组标签，可以自定义字符，显示在质心处。
      axes.color= "black",#设置轴线颜色
      radius = 2,
      bg = "white",# 设置图背景颜色
      axe.titles = c("PC1","PC2","PC3") # 可用于更改轴标题。
      )  # 图6

###  绘制3D PCA双序图---样本以形状显示
pca3d(pca, # 图8
      components = 1:3, 
      group = ind[,4],
      biplot = TRUE, # 绘制双序图 ，变量以红色箭头显示。
      biplot.vars = 5, # 需要展示的变量数量，默认是5。
      palette = cols, # 设置色板
      shape ="s",
      new= TRUE, # 在新的窗口展示图，不覆盖之前的绘图窗口。
      #fancy = TRUE,# 相当于设置同时展示样本标签，质心、样本阴影和分组标签
      bg= "white",
      show.axe.titles=TRUE,
      legend="right",
      show.axes = TRUE,
      show.scale = TRUE,
      axes.color= "black")# 设置图例位置，

###  绘制3D PCA双序图---样本以样本名展示
pca3d(pca, # 图9
      components = 1:3, 
      biplot=TRUE, 
      biplot.vars=6,
      group=ind[,4],
      palette=cols[1:4], # 会按照分组因子水平按序分配颜色。
      show.ellipses = TRUE,
      show.shapes = FALSE, 
      show.labels=TRUE, # 显示样本标签
      labels.col = c(rep(cols[1],9),
                     rep(cols[2],9),
                     rep(cols[3],9),
                     rep(cols[4],9)), # 设置样本标签颜色
      #shape ="s",
      bg= "white",
      legend="right",
      show.axes = TRUE,
      axes.color= "black", 
      show.scale = FALSE,
      new= TRUE
      )

###绘图后可用鼠标将图调整到需要的角度，然后保存该角度到本地
snapshotPCA3d(file="pca3d.png")

图4|ggsci中d3的category10颜色集。

图5|pca3d默认参数绘图结果。用鼠标可以调整图形展示角度。

图6|添加分组信息和自定义形状和颜色。

图7|设置置信椭圆、质心、轴刻度标签等信息。

图8|3D PCA双序图。变量以红色箭头显示。

图9|3D PCA双序图。样本以字符展示。样本太多，就不适合用字符标签显示，可用于探索数据。

2.2 rgl包绘制3D PCA图

rgl是一个很强大的绘制3D图的R包。许多绘制3D图的函数都需要调用rgl包。比如car::scatter3d();compositions::plot3D等函数。使用rgl可以绘制多种类型的3D图形(hist3D(),points3d(),lines3d(),segments3d(),triangles3d()和quads3d()等),可以更方便的添加绘图元素和调整图形参数。使用说明：http://127.0.0.1:16323/library/rgl/doc/rgl.html

## 2.2.1 rgl包绘制3D PCA双序图---单个分组因子
###browseVignettes("rgl") # 查看帮助信息
## 设置颜色
cols = pal_d3("category10")(10)
col = cols[c(7,1,5,2)] # 选择四个颜色，为每个分组指定特定颜色。
get_colors <- function(groups, cols = palette()){
  groups <- as.factor(groups)
  ngrps <- length(levels(groups))
  if(ngrps > length(cols)) 
  cols <- rep(cols, ngrps)
  color <- cols[as.numeric(groups)]
  names(color) <- as.vector(groups)
  return(color)
}
labels.col = get_colors(ind[,4],col)
labels.col
### 或者
ind[,4] = factor(ind[,4],levels = c("CK","NT","PT","WL")) # 设置分组的因子水平
labels.col = col[as.numeric(ind[,4])]
names(labels.col) <- as.vector(ind[,4])
labels.col
### 或者，三种设置颜色方法，效果是一致的
labels.col =  c(rep(cols[7],9),
                     rep(cols[2],9),
                     rep(cols[1],9),
                     rep(cols[5],9)) # 根据tillage分组为样本赋予颜色，这个要根据样本的排布顺序设置颜色。
open3d() # 打开一个绘图窗口；图10.
par3d(family="serif",
      cex=1.2,
      font=1) # 设置图形参数，应用于后面的所有绘图对象,需要与绘图参数一起运行。

### 绘制样本散点3D图
plot3d(ind[,1:3],
       type="p",# "p"表示点，"s"表示球，"l"表示线,"h"表示线段
       col=labels.col,
       lwd = 1,
       size=12,
       ticktype = "detailed",
       #theta = 20,
       #phi = 10, # theta和phi是极坐标
       #fov = 60, # 以度为单位的视场角
       #zoom = 1 # zoom表示缩放系数,这些都可以使用绘图后，使用鼠标调整。
       )

## 3D图在屏幕的视图也可以用view3d()设置
#view3d(theta = 20,
       #phi = 10, # theta和phi是极坐标
       #fov = 90, # 以度为单位的视场角
       #zoom = 1 )
### 默认图形以立方体形式展示，不设置网格线
axes3d() # 只显示x y和z 3个轴
grid3d(c("x","y","z")) #x y和z 3个平面添加网格线

### 使用arrow3d()添加变量箭头
####需要准备箭头的起始点(0,0,0)数据表
origin = data.frame(PC1= rep(0,11),
                    PC2= rep(0,11),
                    PC3= rep(0,11),
                    row.names = rownames(var))
#### 绘制变量箭头
####11个环境因子数量有点多，这里只在图中展示5个环境因子可以使用vegan包先筛选环境因子。
####为了图美观，环境因子坐标都*3.
for (i in 1:nrow(var[1:5,])) {
  arrow3d(p0 = origin[i,], 
        p1 = var[i,]*3, 
        type = "lines", #箭头类型可选 "extrusion"，"lines","rotation","flat"
        width = 3,# 箭头末端倒钩的宽度
        thickness = 0.618 * width, # 箭头倒钩的厚度
        theta = 0.45, # 箭头末端倒钩的角度
        s = 0.2,# 箭头末端倒钩的长度
        n = 2, # 箭头末端倒钩的个数，默认为3
        col = "red",
        plot = TRUE) # plot = FALSE则返回绘图参数，但不绘图
}

#### 添加变量标签
text3d(var[1:5,]*3, 
       texts=rownames(var[1:5,]),
       #adj = c(1,1,1), # 调整字符标签相对于坐标的位置.(0,0)将标签放在左下角落，(0.5,0.5) #中心，(1,1)右上
       offset = 0.5, # 与pos联用,设置标签与坐标点的距离。
       pos = 3, # 设置标签位置。0=坐标所在位置, 1=below,2= left,3=above,4=right,5=in front of,6=behind
       col="black",
       #family = "serif",
       cex=1.5)

#### 添加图例
legend3d("right",
         levels(ind$tillage),
         pch=16, # 设置形状
         col=cols[1:4] # 设置颜色
         )

#### 从open3d()开始运行完命令后，调整为自己想要的角度和缩放程度后，保存pdf格式图片。图例没有输出
rgl.postscript("rgl.pca1.pdf", fmt = "pdf", drawText = T)

close3d() # 关闭3D图的屏幕显示

图10|rgl包绘制3D PCA双序图---单个分组因子，rgl.pca1.pdf。

上面只使用颜色显示了tillage的分组信息，还要一个分组信息depth可以使用形状显示。下面介绍使用pch3d()根据tillage设置样本点颜色，根据depth设置样本点形状。

## 2.2.2 rgl包绘制3D PCA双序图---双分组因子
### 样本根据tillage和depth分组信息设置不同的颜色和形状
#### 设置形状：?points，可以查看可以使用的平面点形状(原点，三角等)。
#### 立体图形有：cube3d();tetrahedron3d();octahedron3d();icosahedron3d();dodecahedron3d();cuboctahedron3d()。
## 设置颜色：颜色上面已经设置好，就不再重复了。
open3d() # 图11
pch3d(ind[,1:3],
       pch = rep(c(15,16,17),# ?points,可以查看可使用的shape
                 time = 12),# 根据depth设置不同的形状
       col=labels.col, # 根据tillage设置形状颜色
       cex = 0.5, # 形状相对于图形的大小
       radius = 0.25 # 坐标中绘制图形大小
      )
bg3d("lightgrey") # 设置图形背景颜色为灰色
axes3d(edges = "bbox",
       box = TRUE, # box = FALSE，则只显示立方体三面的框线。
       col = "black",# 设置轴为黑色
       lwd = 2,# 设置轴线宽度
       family = "serif",
       font = 1,
       marklen = 0.3,
       marklen.rel = FALSE) # 设置刻度线宽度,marklen.rel = TREU,则刻度线宽度 = 1/marklen * axis length。
aspect3d(1, 1, 1) # 设置x,y和z 轴的长度比率为1:1:1，图形近似一个正方体。

title3d(xlab = "PC1",
        ylab = "PC2",
        zlab = "PC3",
        #main = "3D PCA biplot", # 设置图标题
        col = "black")
for (i in 1:nrow(var[1:5,])) {
  arrow3d(p0 = origin[i,], 
        p1 = var[i,]*3, 
        type = "lines",
        width = 3,
        thickness = 0.618 * width, 
        theta = 0.45, 
        s = 0.2,
        n = 2, 
        lwd = 2,
        col = "red",
        plot = TRUE
        ) 
}
text3d(var[1:5,]*3, 
       texts=rownames(var[1:5,]),
       offset = 0.5, 
       pos = 3,
       col="black",
       family = "serif",
       cex=1.2)

## 调整好3D图形展示角度之后，输出pdf图。设置的灰色图形背景没有输出到图中。
rgl.postscript("rgl.pca2.pdf", fmt = "pdf", drawText = T)

图11|rgl包绘制3D PCA双序图---双分组因子，rgl.pca2.pdf。样本点用形状和颜色显示不同分组信息。设置刻度线长度，轴标签等图形参数。

## 2.2.3 rgl包绘制3D PCA双序图---设置xyz坐标轴在(0,0,0)相交。
### rgl实现类似pca3d()的绘图结果
open3d() # 图12
#par3d() # 3d绘图参数设置
#layout3d(matrix(1:4, 2,2)) # 3D图排版设置
###heights表示图布局中行的高度，widths表示图布局中列的宽度，matrix()设置多张图依次在屏幕上的显示位置。
#mfrow3d(2,2,byrow = TRUE,sharedMouse = TRUE) # 设置绘图参数，画布被按行分为4份。

par3d(windowRect = 50 + c( 0, 0, 640, 640)) # 设置绘图窗口大小
pch3d(ind[,1:3],
       pch = rep(c(15,16,17),
                 time = 12),
       col=labels.col, 
       cex = 0.5, 
       radius = 0.25 
      )
### 设置xyz轴线在(0,0,0)相交，轴线颜色设置为灰色
axis3d("x",pos = c(NA,0,0),labels = TRUE,color = "lightgrey") # pos决定了轴线在另两个轴的位置，不与level同时使用。
axis3d("y",pos = c(0,NA,0),labels = TRUE,color = "lightgrey")
axis3d("z",pos = c(0,0,NA),labels = TRUE,color = "lightgrey")

### 添加轴标签
#### 与轴线相同方向的平行线，刻度线从line = 0到line=1，刻度标签放在line=2。level表示与轴线正交的线，
title3d(xlab = "PC1",
        ylab = "PC2",
        zlab = "PC3",
        col = "black")

### 绘制变量箭头与名称
for (i in 1:nrow(var[1:5,])) {
  arrow3d(p0 = origin[i,], 
        p1 = var[i,]*3, 
        type = "lines",
        width = 3,
        thickness = 0.618 * width, 
        theta = 0.45, 
        s = 0.2,
        n = 2, 
        lwd = 2,
        col = "black", # 变量箭头设置为黑色
        plot = TRUE
        ) 
}
text3d(var[1:5,]*3, 
       texts=rownames(var[1:5,]),
       offset = 0.5, 
       pos = 3,
       col="black",
       family = "serif",
       cex=1.2)

## 调整好3D图形展示角度之后，输出pdf图。
#rgl.postscript("rgl.pca3.pdf", fmt = "pdf", drawText = T)


### 添加框面
bbox3d(draw_front = FALSE, # 只绘制立方体的三个边界
       marklen = 0.3,
       marklen.rel = FALSE, 
       color = "#333377", # 设置边界框颜色
        emission="#333377", #emission, specular, shininess: properties for lighting calculation
         specular="#3333FF", shininess=5, alpha=0.8
       ) # bbox3d()与aexs3d()功能有重叠，先运行的会被后运行的结果覆盖。

#highlevel(integer()) # 以rglwidget显示，运行命令前后，显示的图没什么差别。

## 调整好3D图形展示角度之后，输出pdf图。设置的狂免颜色亮度没有输出到图中。
rgl.postscript("rgl.pca3.pdf", fmt = "pdf",drawText = TRUE)

图12|rgl模仿pca3d()绘图结果，rgl.pca3.pdf。

## 2.2.4 rgl包绘制3D PCA双序图---根据分组信息添加集合椭圆/球
names(col) = c("CK","NT","PT","WL") # 因子水平与样本排列顺序不一致，在设置颜色和形状时要特别注意。样本点的颜色与集合椭圆/球的颜色保持一致。
open3d() #图13
pch3d(ind[,1:3],
       pch = rep(c(15,16,17),# ?points,可以查看可使用的shape
                 time = 12),# 根据depth设置不同的形状
       col=labels.col, # 根据tillage设置形状颜色
       cex = 0.5, # 形状相对于图形的大小
       radius = 0.25 # 坐标中绘制图形大小
      )
axes3d(edges = "bbox",labels = TRUE,
       box = FALSE, # box = TRUE，则只显示立方体六面的框线。
       col = "black",# 设置轴为黑色
       lwd = 2,# 设置轴线宽度
       family = "serif",
       font = 1,
       marklen = 0.3,
       marklen.rel = FALSE) # 设置刻度线宽度,marklen.rel = TREU,则刻度线宽度 = 1/marklen * axis length。
aspect3d(1, 1, 1)

### 添加轴标签
#### 与轴线相同方向的平行线，刻度线从line = 0到line=1，刻度标签放在line=2。level表示与轴线正交的线，
title3d(xlab = "PC1",
        ylab = "PC2",
        zlab = "PC3",
        col = "black")

### 绘制变量箭头与名称
for (i in 1:nrow(var[1:5,])) {
  arrow3d(p0 = origin[i,], 
        p1 = var[i,]*3, 
        type = "lines",
        width = 3,
        thickness = 0.618 * width, 
        theta = 0.45, 
        s = 0.2,
        n = 2, 
        lwd = 2,
        col = "black", # 变量箭头设置为黑色
        plot = TRUE
        ) 
}
text3d(var[1:5,]*3, 
       texts=rownames(var[1:5,]),
       offset = 0.5, 
       pos = 3,
       col="black",
       family = "serif",
       cex=1.2)

### 根据tillage绘制分组集合区域
for (i in 1:length(levels(ind[,4]))){
    group <- levels(ind[,4])[i]
    xyz <- ind[which(ind[,4] == group),] 
    ellips <- ellipse3d(cov(cbind(xyz[,1:3])), 
              centre=c(mean(xyz[,1]), xyz[,2], mean(xyz[,3])), level = 0.95) 
    #shade3d(ellips, col = col[i], alpha = 0.2, lit = FALSE) # 阴影椭圆形式展示集合区间
    wire3d(ellips, col = col[i],  lit = FALSE) # 网格椭球形式
    # show group labels
    texts3d(apply(xyz[1:3], 2, mean),
            text = group, # 在椭圆中心添加分组标签
            col= col[i], cex = 2)
}
## 调整好3D图形展示角度之后，输出pdf图。设置的框面颜色亮度没有输出到图中。
rgl.postscript("rgl.pca4.pdf", fmt = "pdf",drawText = TRUE)


## 添加单个集合区间
### 获得一个mesh3d对象，用于后续绘图。
ellips <- ellipse3d(cov(ind[,1:3]), # 相关性或协方差矩阵
            centre=c(mean(ind[,1]), mean(ind[,2]), mean(ind[,3])),# 椭圆中心 
            level = 0.95) # level设置集合区间水平，用于控制椭球的大小。
### 绘制单个集合区间，阴影或网格形式
shade3d(ellips, col = "#D95F02", alpha = 0.1, lit = FALSE) #阴影平面
wire3d(ellips, col = "#D95F02",  lit = FALSE) # 网格椭球
#### shade3d()与wire3d()同时运行，会形成网格椭圆集合。

图13|添加集合区间3D PCA图，rgl.pca4.pdf。

数据表和代码可以QQ交流群文件夹中下载，或微信公众号后台发送“PCA3D”获取。

参考资料

http://www.sthda.com/english/wiki/a-complete-guide-to-3d-visualization-device-system-in-r-r-software-and-data-visualization#d-scatter-plothttps://dmurdoch.github.io/rgl/dev/reference/material.htmlhttps://dmurdoch.github.io/rgl/dev/reference/polygon3d.htmlhttps://dmurdoch.github.io/rgl/dev/reference/rglwidget.html

推荐阅读

R绘图-RDA排序分析

R统计绘图-VPA(方差分解分析)

R统计绘图-RDA分析、Mantel检验及绘图

R统计-PCA/PCoA/db-RDA/NMDS/CA/CCA/DCA等排序分析教程R统计绘图-PCA分析绘图及结果解读(误差线，多边形，双Y轴图、球形检验、KMO和变量筛选等)

R统计-微生物群落结构差异分析及结果解读

R统计绘图-PCA分析及绘制双坐标轴双序图

R统计绘图-分子生态相关性网络分析

R中进行单因素方差分析并绘图R统计-多变量单因素参数、非参数检验及多重比较

R绘图-相关性分析及绘图

R绘图-相关性系数图

R统计绘图-环境因子相关性热图

R统计绘图-corrplot绘制热图及颜色、字体等细节修改

R统计绘图-corrplot热图绘制细节调整2(更改变量可视化顺序、非相关性热图绘制、添加矩形框等)

R数据可视化之美-节点链接图R统计绘图-rgbif包下载GBIF数据及绘制分布图

R统计-正态性分布检验[Translation]

R统计-数据正态分布转换[Translation]

R统计-方差齐性检验[Translation]

R统计-Mauchly球形检验[Translation]R统计绘图-单、双、三因素重复测量方差分析[Translation]

R统计绘图-混合方差分析[Translation]

R统计绘图-协方差分析[Translation]

R统计绘图-One-Way MANOVA