手工绘制PCA载荷向量与PCAtools双图不一致问题求助

使用PCAtools对宏基因组基因计数数据做PCA分析时，其自带的biplot()绘制的双图中，载荷向量长度能匹配数据离散度（向量长度代表对离散的影响）；但直接提取p$rotated（样本得分）和p$loadings（变量载荷）用ggplot2手动绘图时，载荷向量极小，比如K20444在PCAtools双图中处于-1到-2.5之间，ggplot2图中却在0到-1之间，二者尺度明显不符。

原始PCAtools绘图代码 #

library(PCAtools)
p <- pca(ex,removeVar = 0.1)
biplot(p,lab=NULL,hline = 0, vline = 0,
       showLoadings = T)
# 此代码生成含3个载荷的PCAtools双图

原始ggplot2手动绘图代码 #

# 选择匹配的载荷
loadplot <- c("K02469","K10118","K20444")
# 提取PC得分和载荷
PC <- p$rotated
load <-p$loadings
load <- na.omit(load[loadplot,])

# ggplot2绘图
p.pcaNORM <- ggplot(data = PC, aes(x = PC1, y = PC2)) +
  geom_point(aes(),shape = 21,colour = "black", size = 5) +
  geom_text(data=load, aes(x=PC1, y=PC2, label=row.names(load)),
            inherit.aes = FALSE, size=3) +
  geom_segment(data=load, aes(x=0, xend=load$PC1, y=0, yend=load$PC2), inherit.aes = FALSE,
               arrow = arrow(length = unit(0.25, "cm")), colour = "grey") +
  labs(x = "PC1 27.8%", y = "PC2 12.3%", title="PCA of Gene Count Differential Abundances") +
  geom_hline(yintercept = 0, linetype="dashed",color="black") +
  geom_vline(xintercept=0, linetype="dashed",color="black") +
  theme_linedraw() +
  theme(legend.position  = "bottom", 
        legend.box       = "horizontal",
        legend.direction = "horizontal",
        panel.background = element_blank(),
        axis.text.x      = element_text(size=8),
        axis.text.y      = element_text(size=8),
        plot.title       = element_text(size = 10))
print(p.pcaNORM)

输入数据 #

ex <- structure(list(Ga0598246 = c(6.56539918955656, 6.07685285958297, 
9.80857543461825, 7.45380707677968, 5.50253732907563, 0, 1.00547341301176, 
9.57886787328451, 5.91763530872237, 9.13761113968035), Ga0598239 = c(5.37714039215198, 
4.44627250098116, 9.12527549535166, 7.37818173071451, 4.22343869037075,  0, 0, 9.41186693570199, 2.04265725161515, 9.06211714113395), 
    Ga0598242 = c(7.23329507168685, 3.30011916965994, 9.14571484611697, 
    6.63420100118826, 5.33389579265697, 0, 0, 9.23317518972483, 
    4.83457149887647, 9.02564207931812), Ga0598241 = c(6.42805923331295, 
    4.41343389152249, 9.17585857319441, 6.53724813980237, 5.95254094365946, 
    0, 1.55299220383005, 9.34416529172944, 5.9067619585435, 9.35272687344133
    ), Ga0598244 = c(6.51891536576828, 4.56535598153963, 9.40792503580225, 
    6.60909974507472, 5.37234302403178, 0, 0.989803762135471, 
    9.61255610440809, 5.76126935753228, 9.26497019795582), Ga0598240 = c(6.59592709583215, 
    5.0116546308929, 9.10583727581318, 6.45845808171092, 6.32419402380456, 
    0, 1.97474341675587, 9.59380242030398, 6.25206560327315, 
    9.37141719913754), Ga0598243 = c(6.48685136970719, 4.54940912022517, 
    8.95163946868835, 6.2487649672149, 6.21130294335921, 0, 1.97225049345733, 
    9.55911910161097, 6.0719500604804, 9.61218376859803), Ga0598247 = c(3.42017839411825, 
    5.22468821931174, 9.43811837974567, 8.49700801179772, 5.02416088318049, 
    0, 0, 10.2938603383254, 0, 8.70421802755197), Ga0598259 = c(6.85309401728228, 
    5.08842002313584, 9.16500402830579, 6.20640673224555, 5.24441419037537, 
    0, 0, 9.46180201100927, 5.5134145812396, 9.53162068178508
    ), Ga0598260 = c(6.7006085700961, 4.80510960111303, 9.06424555064897, 
    7.05413314135196, 5.56105321976388, 0, 0, 9.94666671304116, 
    5.99044411371985, 9.14986983395412)), row.names = c("K20444",  "K01711", "K10119", "K10118", "K02469", "K03046", "K14358", "K00945",  "K12454", "K01129"), class = "data.frame")

核心原因：PCAtools对载荷向量做了自动缩放 #

PCAtools的biplot()函数为了让载荷向量的尺度和样本得分的尺度匹配，方便直观观察变量对样本分组的影响，默认会对原始载荷值进行缩放——具体是将载荷乘以对应主成分的标准差（p$sdev），确保向量长度能反映变量对数据离散的贡献。

而直接从p$loadings提取的是原始未缩放的载荷值，这些值本身的量级很小（通常在-1到1之间），所以用ggplot2绘制时向量会显得极小，和PCAtools的双图尺度不符。

解决方法：手动对载荷进行缩放 #

要让ggplot2绘制的载荷向量和PCAtools双图一致，只需将原始载荷乘以对应PC的标准差（p$sdev），即可得到和PCAtools中一致的缩放后载荷。

修改后的ggplot2绘图代码如下：

# 选择匹配的载荷
loadplot <- c("K02469","K10118","K20444")
# 提取PC得分和载荷
PC <- p$rotated
load <- p$loadings
load <- na.omit(load[loadplot,])

# 关键步骤：对载荷进行缩放，匹配PCAtools的尺度
load_scaled <- load %*% diag(p$sdev[1:2])  # 只缩放PC1和PC2

# ggplot2绘图，使用缩放后的载荷
p.pcaNORM <- ggplot(data = PC, aes(x = PC1, y = PC2)) +
  geom_point(shape = 21, colour = "black", size = 5) +
  geom_text(data = load_scaled, aes(x = PC1, y = PC2, label = row.names(load_scaled)),
            inherit.aes = FALSE, size = 3) +
  geom_segment(data = load_scaled, aes(x = 0, xend = PC1, y = 0, yend = PC2), 
               inherit.aes = FALSE, arrow = arrow(length = unit(0.25, "cm")), colour = "grey") +
  labs(x = paste0("PC1 ", round(p$variance[1]*100,1), "%"), 
       y = paste0("PC2 ", round(p$variance[2]*100,1), "%"), 
       title = "PCA of Gene Count Differential Abundances") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  theme_linedraw() +
  theme(legend.position  = "bottom", 
        legend.box       = "horizontal",
        legend.direction = "horizontal",
        panel.background = element_blank(),
        axis.text.x      = element_text(size = 8),
        axis.text.y      = element_text(size = 8),
        plot.title       = element_text(size = 10))
print(p.pcaNORM)

验证说明 #

运行上述代码后，载荷向量的尺度会和PCAtools的双图完全一致：比如K20444的PC1值会放大到和样本得分同量级的范围（-1到-2.5之间），向量长度也能正确反映该变量对数据离散的影响程度。

内容的提问来源于stack exchange，提问作者Geomicro