【实战指南】如何用Python绘制带置信椭圆的PCA双标图？附完整代码

1. PCA双标图与置信椭圆的核心价值主成分分析PCA是数据降维的经典方法但在生物信息学和多组学数据分析中单纯展示样本点在二维空间的分布往往不够。我第一次处理基因表达数据时就踩过坑——虽然能看到样本聚类但完全不明白哪些基因在驱动这种分布模式。这就是为什么我们需要双标图Biplot和置信椭圆的组合可视化。双标图的精妙之处在于它同时呈现了两个维度的信息样本分布展示降维后的数据点聚类情况变量载荷用箭头表示原始高维特征对主成分的贡献度而置信椭圆则像给每个类别画了个势力范围椭圆大小反映类内离散程度。实测下来这种组合能快速回答三个关键问题不同组别的样本是否显著分离哪些原始变量如基因、代谢物主导了这种分离组内变异程度如何2. 环境准备与数据模拟2.1 安装必要库建议使用conda创建虚拟环境conda create -n pca_plot python3.8 conda activate pca_plot pip install numpy matplotlib scikit-learn pandas2.2 模拟多组学数据我们模拟三类代谢组学样本每类包含20个代谢物特征import numpy as np from sklearn.decomposition import PCA # 设置随机种子保证可复现 np.random.seed(42) # 定义三类样本的均值和协方差矩阵 class_params { A: {mean: [5, 10], cov: [[2, 1.5], [1.5, 2]]}, B: {mean: [10, 5], cov: [[1, -0.8], [-0.8, 1]]}, C: {mean: [15, 15], cov: [[3, 0], [0, 1]]} } # 生成模拟数据每类50个样本 samples [] for class_name, params in class_params.items(): samples.append( np.random.multivariate_normal( meanparams[mean], covparams[cov], size50 ) ) X np.vstack(samples) y np.array([k for k in class_params.keys() for _ in range(50)]) # 添加20个模拟代谢物特征 n_features 20 feature_loadings np.random.randn(n_features, 2) * 2 # 随机生成特征载荷3. PCA计算与双标图绘制3.1 执行PCA降维pca PCA(n_components2) X_pca pca.fit_transform(X) # 获取主成分解释方差比 explained_var pca.explained_variance_ratio_ print(f主成分解释方差: PC1{explained_var[0]:.1%}, PC2{explained_var[1]:.1%})3.2 绘制基础双标图import matplotlib.pyplot as plt plt.figure(figsize(10, 8), dpi120) # 绘制样本散点 colors {A: #FF6B6B, B: #4ECDC4, C: #45B7D1} for class_name, color in colors.items(): mask y class_name plt.scatter( X_pca[mask, 0], X_pca[mask, 1], colorcolor, labelclass_name, alpha0.7 ) # 绘制特征载荷箭头 for i in range(n_features): plt.arrow( 0, 0, # 从原点出发 feature_loadings[i, 0] * 8, feature_loadings[i, 1] * 8, # 放大显示 colorgray, alpha0.5, head_width0.2 ) plt.text( feature_loadings[i, 0] * 8.5, feature_loadings[i, 1] * 8.5, fF{i1}, fontsize9, colorblack ) plt.xlabel(fPC1 ({explained_var[0]:.1%})) plt.ylabel(fPC2 ({explained_var[1]:.1%})) plt.title(PCA Biplot with Feature Loadings) plt.legend() plt.grid(alpha0.3) plt.show()4. 置信椭圆的计算与绘制4.1 椭圆参数计算原理置信椭圆基于多元正态分布假设关键参数包括中心点类内样本均值长/短轴由协方差矩阵的特征值决定旋转角度由特征向量方向决定数学表达式 $$(x-\mu)^T \Sigma^{-1} (x-\mu) \chi^2_{crit}(p)$$ 其中$\chi^2_{crit}$是卡方分布的临界值对应置信水平。4.2 Python实现代码from matplotlib.patches import Ellipse from scipy.stats import chi2 def plot_confidence_ellipse(ax, points, color, alpha0.3, level0.95): # 计算协方差矩阵 cov np.cov(points.T) # 计算椭圆参数 lambda_, v np.linalg.eig(cov) angle np.degrees(np.arctan2(v[1,0], v[0,0])) width, height 2 * np.sqrt(lambda_ * chi2.ppf(level, df2)) # 绘制椭圆 ellipse Ellipse( xynp.mean(points, axis0), widthwidth, heightheight, angleangle, colorcolor, alphaalpha, lw2, fillFalse, zorder0 ) ax.add_patch(ellipse) # 在之前图形上添加椭圆 fig, ax plt.subplots(figsize(10, 8), dpi120) for class_name, color in colors.items(): mask y class_name class_points X_pca[mask] ax.scatter(class_points[:,0], class_points[:,1], colorcolor, labelclass_name) plot_confidence_ellipse(ax, class_points, colorcolor) # 添加特征箭头略 plt.title(PCA with 95% Confidence Ellipses) plt.show()5. 高级定制技巧5.1 双标图箭头优化原始箭头可能重叠严重可以只显示重要特征载荷绝对值前10%添加透明度避免遮挡# 计算特征重要性 loading_norms np.linalg.norm(feature_loadings, axis1) threshold np.percentile(loading_norms, 90) # 取前10% for i in range(n_features): if loading_norms[i] threshold: plt.arrow(0, 0, feature_loadings[i,0]*7, feature_loadings[i,1]*7, colorpurple, alpha0.3, width0.05) plt.text(feature_loadings[i,0]*7.2, feature_loadings[i,1]*7.2, fF{i1}, fontsize8, colordarkred)5.2 多椭圆样式定制ellipse_styles { A: {color: #FF6B6B, linestyle: --, linewidth: 1.5}, B: {color: #4ECDC4, linestyle: -, linewidth: 2}, C: {color: #45B7D1, linestyle: :, linewidth: 2} } for class_name, style in ellipse_styles.items(): mask y class_name points X_pca[mask] # 计算椭圆参数同上 ... ellipse Ellipse(..., **style) ax.add_patch(ellipse)5.3 3D PCA可视化扩展from mpl_toolkits.mplot3d import Axes3D pca_3d PCA(n_components3) X_3d pca_3d.fit_transform(X) fig plt.figure(figsize(10, 8)) ax fig.add_subplot(111, projection3d) for class_name, color in colors.items(): mask y class_name ax.scatter(X_3d[mask,0], X_3d[mask,1], X_3d[mask,2], colorcolor, labelclass_name, alpha0.7) ax.set_xlabel(PC1) ax.set_ylabel(PC2) ax.set_zlabel(PC3) plt.legend() plt.show()6. 实战注意事项数据标准化至关重要from sklearn.preprocessing import StandardScaler X_scaled StandardScaler().fit_transform(X)当特征量远大于样本量时考虑稀疏PCAfrom sklearn.decomposition import SparsePCA spca SparsePCA(n_components2, alpha0.1) X_spca spca.fit_transform(X_scaled)置信水平选择建议探索性分析95%默认严格验证99%调整方法修改chi2.ppf()中的level参数图形标注最佳实践主成分轴标注解释方差比特征箭头用不同颜色标识正负贡献添加图例说明椭圆置信水平我在实际项目中遇到过椭圆形状异常的情况通常是因为样本量过少每类建议至少30个样本数据存在极端异常值协方差矩阵计算不稳定解决方法包括增加样本量、进行数据清洗或使用稳健协方差估计方法。