PCA降维原理以及用R绘制PCA图

PCA图的原理

基本概念 - PCA（Principal Component Analysis）即主成分分析，是一种常用的数据分析方法。它的主要目的是通过线性变换将原始数据转换为一组新的变量，这些新变量被称为主成分。主成分是原始变量的线性组合，并且它们是按照方差大小依次排列的。 - 比如，假设有一组二维数据点$(x_1,y_1),(x_2,y_2),\cdots,(x_n,y_n)$，PCA的目标就是找到两个新的变量（主成分），使得这些数据点在新的变量表示下能够最大程度地保留数据的方差信息。
原理推导
- 协方差矩阵
  - 首先考虑数据的协方差矩阵。对于有 $m$ 个特征（变量）的数据集$X$（$n$个样本），其协方差矩阵$\Sigma$的元素$\sigma_{ij}$表示第$i$个特征和第$j$个特征之间的协方差，计算公式为$\sigma_{ij}=\frac{1}{n - 1}\sum_{k = 1}^{n}(x_{ki}-\overline{x}i)(x{kj}-\overline{x}_j) $，其中$$\overline{x}_i$$是第$i$个特征的均值。
  - 协方差矩阵是一个对称矩阵，它衡量了不同特征之间的相关性。例如，如果两个特征的协方差为正，表示它们正相关；协方差为负，表示负相关；协方差为0，表示不相关。
- 特征值和特征向量
  - PCA的核心是求解协方差矩阵的特征值和特征向量。对于一个$m\times m$的协方差矩阵$\Sigma$，如果存在一个非零向量$v$和一个标量$\lambda$，使得$\Sigma v=\lambda v$，那么$\lambda$就是协方差矩阵的特征值，$v$就是对应的特征向量。
  - 特征向量代表了数据的方向，而特征值代表了数据在该方向上的方差大小。例如，较大的特征值对应的特征向量所指向的方向是数据方差较大的方向。
- 主成分的选取
  - 在求出协方差矩阵的所有特征值和特征向量后，将特征值从大到小排序。选择前$k$个最大特征值对应的特征向量，这些特征向量构成了一个$m\times k$ 的矩阵$W$。
  - 原始数据$X$（$n\times m$）通过线性变换$Y = XW$得到新的数据$Y$（$n\times k$），这里的$Y$就是数据在主成分空间中的表示。新的数据$Y$的每个维度（主成分）是相互正交（不相关）的，并且按照方差大小依次排列，前几个主成分包含了原始数据的大部分方差信息。
降维应用原理
- 假设原始数据有10个特征，但通过PCA分析发现，前3个主成分已经能够解释原始数据80%以上的方差。那么就可以将原始数据从10维降到3维。
- 这是因为在实际应用中，高维数据可能存在很多冗余信息，通过PCA降维可以去除这些冗余，同时保留数据的主要结构信息。例如，在图像识别中，图像数据可能有很高的维度（很多像素点），PCA可以帮助提取图像的主要特征，减少后续处理的数据量，提高计算效率。
可视化应用原理
- 当数据维度大于3时，很难直接可视化。但如果通过PCA将数据降到2维或3维，就可以将数据点绘制在平面或三维空间中。
- 例如，在基因表达数据分析中，有成千上万个基因（特征），通过PCA将其降到2维或3维后，可以直观地观察不同样本（数据点）之间的关系，如不同疾病组和健康组样本是否能够在PCA图中明显区分开来，从而帮助分析疾病相关的基因表达模式。

PCA图在区分转录组样本时的简单判断依据

PCA图主要用来降维和分类，PC1 和 PC2 实际上是基于这些计算得到的特征向量所确定的新的坐标轴方向，它们分别代表了数据在不同维度上的最大方差方向，也就是能最大程度地区分样本的方向。通俗一点的理解就是，PCA 的作用就是找到一些新的坐标轴，将这些原本在高维空间的样本点投影到这个由 PC1 和 PC2 构成的二维平面上，这样我们就可以通过观察这些样本点在这个二维平面上的分布情况，比如它们之间的距离远近等，来判断样本之间的相似程度或者是否能分成不同的类别等。

所以问题来了：怎么利用PCA去除离群样本？PCA不好要不要重测？

首先观察PC1和PC2占比多少，PC1+PC2之和最好大于80%比较有说服力。
- 情景1：PC1（大于90%）远大于PC2。这时候主要观察样本在PC1（横坐标的距离），PC1距离接近，但是PC2距离较远，样本之间的相似度还是比较接近的。反之，PC1距离较远，不管PC2多接近，样本之间的相似度可能是比较远的。
- 情景2：PC1和PC2接近。通过两点之间的最短距离来判断相似度。

注意注意！PCA只能用来相对地看样本之间的差异，如果没有对照，只有重复的三个处理样本，PCA可能会放大不重要的差异！

在R中简单绘制PCA（基于prcomp）

# 清除当前环境中的所有对象，确保没有残留变量干扰
rm(list = ls())

# 加载必要的包
library(ggplot2)

# 读取 data.csv 文件，添加错误处理
data <- read.csv("data.csv", row.names = 1)
head(data)
#               wt1   wt2   wt3   ko1   ko2   ko3   kd1   kd2   kd3
# MSTRG.10004 2.639 1.849 2.087 0.149 0.506 0.816 0.382 0.571 0.736
# MSTRG.10008 1.420 1.341 1.350 0.564 0.390 0.682 0.248 0.540 0.637
# MSTRG.10158 1.470 1.655 1.081 0.000 0.771 0.528 0.136 0.000 0.787
# MSTRG.10160 0.699 0.611 0.483 1.014 0.724 0.478 0.219 0.307 0.580
# MSTRG.10171 0.558 1.025 0.221 0.896 0.790 0.649 1.983 2.475 2.239
# MSTRG.10179 0.565 0.545 0.573 1.392 1.717 1.368 1.334 0.831 1.215
data_t = t(data)
dim(data_t)

# 过滤掉全为0的行和列
data_t <- data_t[,colSums(data_t)>0]
data_t <- data_t[,apply(data_t, 1, var)!=0]

# 执行主成分分析
pca_result <- prcomp(data_t, scale. = TRUE)

# 计算主成分的百分比解释方差
pca_percentage <- (pca_result$sdev^2 / sum(pca_result$sdev^2)) * 100

# 可视化
library(factoextra)
fviz_eig(pca_result, addlabels = TRUE)
fviz_pca_ind(pca_result,
  col.ind = rownames(data_t), mean.point = F,
  addEllipses = T, legend.title = "Groups"
)

# 或者ggplot2
p <- ggplot(
  pca_result$x,
  aes(x = PC1, y = PC2, color = rownames(data_t))
) +
  geom_point() +
  labs(
    title = "PCA Plot",
    x = "Principal Component 1",
    y = "Principal Component 2",
    color = "Group"
  ) +
  # 添加坐标轴
  xlab(paste0("PC1 (", round(pca_percentage[1], 2), "%)")) +
  ylab(paste0("PC2 (", round(pca_percentage[2], 2), "%)"))
p

NEXTUpset图的原理以及在富集分析中简单的应用