R语言如何绘制三维PCA图

前言

本文是R包scatterplot3d绘制三维PCA的教程。

什么是PCA？

人眼一般能感知的空间为二维和三维。高维数据可视化的重要目标就是将高维数据呈现于二维或三维空间中。高维数据变换就是使用降维度的方法，使用线性或非线性变换把高维数据投影到低维空间，去掉冗余属性，但同时尽可能地保留高维空间的重要信息和特征。

PCA, 主成分分析法，也被称为主分量分析法，是很常用的一种数据降维方法。主成分分析法采用一个线性变换将数据变换到一个新的坐标系统，使得任何数据点投影到第一个坐标（第一主成分）的方差最大，在第二个坐标（第二主成分）的方差为第二大，以此类推。因此，主成分分析可以减少数据的维数，并保留对方差贡献最大的特征。

三维PCA和普通PCA的区别就是多了一个维度。

是主成分分析的PC1，PC2和PC3的结果，xyz坐标分别为前三个主成分，括号内的百分比为该主成分能解释的变量的百分比。PCA得分图能将对照组和实验组样本区分开。在PCA图中，如果样本之间聚集在一起，说明这些样本差异性小；反之样本之间距离越远，说明样本之间差异性越大。不同颜色的散点表示不同实验分组的样本。

绘图前的数据准备

绘图数据

数据来源一般是搜库结果定量表。包含2个维度的数据，一般情况下，每一行是一个基因，每一列是一个样本。

demo数据可以在https://www.r2omics.cn/res/demodata/PCA/data.txt下载。

分组数据

demo数据可以在https://www.r2omics.cn/res/demodata/PCA/sample.class.txt下载。

R语言如何绘制三维PCA图

# 代码来源：https://www.r2omics.cn/
# 加载R包，没有安装请先安装  install.packages("包名") 
library(scatterplot3d)

# 读取PCA数据文件
df = read.delim("https://www.r2omics.cn/res/demodata/PCA/data.txt",# 这里读取了网络上的demo数据，将此处换成你自己电脑里的文件
                header = T,    # 指定第一行是列名
                row.names = 1  # 指定第一列是行名
)
df=t(df) # 对数据进行转置，如果想对基因分组则不用转置

# 读取样本分组数据文件
dfGroup = read.delim("https://www.r2omics.cn/res/demodata/PCA/sample.class.txt",
                     header = T,
                     row.names = 1
)

# PCA计算
pca_result <- prcomp(df,
                     scale=T  # 一个逻辑值，指示在进行分析之前是否应该将变量缩放到具有单位方差
)
pca_result$x<-data.frame(pca_result$x)
# 设置颜色，有几个分组就写几个颜色
colors <- c("#ff007a","#ffc700","#1dd66d")
colors <- colors[as.numeric(as.factor(dfGroup[,1]))]
# 设置点形状，仅供参考
# shape<-16:18
# shape<-shape[as.numeric(as.factor(dfGroup[,1]))]

# 计算PC值，并替换列名，用来替换坐标轴上的标签
pVar <- pca_result$sdev^2/sum(pca_result$sdev^2)
pVar = round(pVar,digits = 3)
colnames(pca_result$x)[1:3] = c(
  paste0("PC1 (",as.character(pVar[1] * 100 ),"%)"),
  paste0("PC2 (",as.character(pVar[2] * 100 ),"%)"),
  paste0("PC3 (",as.character(pVar[3] * 100 ),"%)")
)

# 绘图
s3d <- scatterplot3d(pca_result$x[,1:3],
                     pch = 16,       # 点形状
                     color=colors,   # 点颜色
                     cex.symbols = 2 # 点大小
)
# 设置图例
legend("top",
       legend = unique(dfGroup[,1]),
       col =  c("#ff007a","#ffc700","#1dd66d"),
       pch = 16,
       inset = -0.1,
       xpd = TRUE,
       horiz = TRUE)
# 设置文字标注
text(s3d$xyz.convert(pca_result$x[,c(1,2,3)] + 2),
     labels = row.names(pca_result$x),
     cex = 0.8,col = "black")