PCA双图绘制问题：基因箭头与样本表达关联不符

开发者阿泽

2026-6-12

PCA双图箭头指向错误的修正方案

需求说明

我要复现一张PCA双图，要求如下：

散点代表样本，按pb.obs中的sample列着色
箭头代表基因，指向该基因高表达样本所在的区域：若基因与特定样本簇关联，箭头指向该簇；若在两类样本中高表达，箭头指向两类簇之间，依此类推

问题描述

已生成带箭头的PCA图（样本聚类正常），但箭头无实际意义——部分基因箭头指向未表达该基因的样本簇，需要修正箭头使其符合表达关联。

原代码（PCA聚类正常，但箭头无意义）

import numpy as np
import matplotlib.pyplot as plt
from adjustText import adjust_text
import scanpy as sc

sc.pp.normalize_total(pb, target_sum=1e6)  # counts per million
sc.pp.log1p(pb)

sc.pp.highly_variable_genes(pb, n_top_genes=1000)
pb = pb[:, pb.var['highly_variable']].copy()

sc.tl.pca(pb)
sc.pp.neighbors(pb, n_neighbors=5, n_pcs=30)
sc.tl.umap(pb)

PCS        = (1, 2)          # PCs
N_PER_SIDE = 15          
COLOR_BY   = 'sample'         # column in pb.obs for dot colours
GENE_CLR   = 'firebrick'
ARROW_KW   = dict(ls='-', lw=1.4, alpha=.9)
LABEL_PAD  = 1.08            # push gene names just beyond arrow tips
LEN_FRAC   = .85             # longest arrow reaches 85 % of sample radius
# ─────────────────────────────────────────────────────────────────────

# 1 pull scores & loadings
c1, c2    = (i - 1 for i in PCS)
S         = pb.obsm['X_pca'][:, [c1, c2]]            
L_raw     = pb.varm['PCs'][:, [c1, c2]]             
eigvals   = pb.uns['pca']['variance'][[c1, c2]]
L_corr    = L_raw * np.sqrt(eigvals)                  # correlations

# 2 pick genes: top N each side of each PC
idxs = []
for dim in (0, 1):                      # PC-1 then PC-2
    for sign in (+1, -1):
        # sort by signed loading
        order = np.argsort(sign * L_corr[:, dim])
        idxs.extend(order[-N_PER_SIDE:])              # grab N from the end
idxs = list(dict.fromkeys(idxs))                      # keep order unique

# 3 rescale arrows so they fill the plot
max_r_gene  = np.sqrt((L_corr[idxs]**2).sum(1)).max()
max_r_score = np.abs(S).max()
scale       = (LEN_FRAC * max_r_score) / max_r_gene
V           = L_corr[idxs] * scale

# 4 plot
fig, ax = plt.subplots(figsize=(15, 13))

# dots
codes = pb.obs[COLOR_BY].astype('category').cat.codes
scat  = ax.scatter(S[:, 0], S[:, 1], c=codes, cmap='tab20',
                   s=200, edgecolor='k', alpha=.85)
# legend to the right
handles, _ = scat.legend_elements(prop='colors')
labels     = pb.obs[COLOR_BY].astype('category').cat.categories
ax.legend(handles, labels, title=COLOR_BY, loc='center left',
          bbox_to_anchor=(1.05, .5), frameon=False)

# arrows and text 
texts = []
for (dx, dy), i in zip(V, idxs):
    ax.arrow(0, 0, dx, dy, color=GENE_CLR, **ARROW_KW,
             head_width=.03*scale/max_r_score,
             head_length=.04*scale/max_r_score,
             length_includes_head=True, zorder=3)
    texts.append(ax.text(dx*LABEL_PAD, dy*LABEL_PAD,
                         pb.var_names[i], color=GENE_CLR,
                         fontsize=12, ha='center', va='center', zorder=4))


adjust_text(texts, arrowprops=dict(arrowstyle='-', color=GENE_CLR))


vr = pb.uns['pca']['variance_ratio'][[c1, c2]] * 100
ax.set_xlabel(f'PC{PCS[0]}  ({vr[0]:.1f} % var)', fontsize=13)
ax.set_ylabel(f'PC{PCS[1]}  ({vr[1]:.1f} % var)', fontsize=13)
ax.axhline(0, lw=.5, color='grey'); ax.axvline(0, lw=.5, color='grey')
ax.set_aspect('equal')
ax.set_title('PCA biplot – top gene drivers per axis')
plt.tight_layout(); plt.show()

模拟数据生成代码

import numpy as np
import pandas as pd
import scanpy as sc

np.random.seed(42)

n_cells   = 180
n_genes   = 500
groups    = np.repeat(['A', 'B', 'C'], repeats=n_cells//3)     # 60 each

# counts
X = np.random.poisson(lam=1.5, size=(n_cells, n_genes)).astype(float)

# marker genes: first 50 for A, next 50 for B, next 50 for C
X[groups == 'A',  :50] += np.random.poisson(4,  ( (groups == 'A').sum(), 50))
X[groups == 'B', 50:100] += np.random.poisson(4, ((groups == 'B').sum(), 50))
X[groups == 'C',100:150] += np.random.poisson(4, ((groups == 'C').sum(), 50))

pb = sc.AnnData(X,
                obs  = pd.DataFrame({'sample': groups},
                                    index=[f'cell{i}' for i in range(n_cells)]),
                var  = pd.DataFrame(index=[f'gene{j}' for j in range(n_genes)]))

sc.pp.normalize_total(pb, target_sum=1e6, inplace=True)
sc.pp.log1p(pb)
sc.pp.pca(pb, n_comps=20)

问题原因与修正方案

问题核心

原代码按单个PC轴的loading极端值筛选基因，导致选中的基因不一定是样本簇的特异性标记基因；即使L_corr的计算逻辑正确（代表基因表达与PCA得分的相关性），非标记基因的箭头也无法对应到高表达样本区域。

修正步骤

筛选样本簇标记基因：通过差异分析获取每个样本簇的特异性标记基因，确保箭头对应的是与簇高度关联的基因
保留箭头方向正确性：标记基因的箭头方向由基因表达与PCA得分的相关性决定，自然会指向该基因高表达的样本簇区域

修正后代码

import numpy as np
import matplotlib.pyplot as plt
from adjustText import adjust_text
import scanpy as sc

# --------------- 数据预处理（保留原逻辑）---------------
sc.pp.normalize_total(pb, target_sum=1e6)  # counts per million
sc.pp.log1p(pb)

sc.pp.highly_variable_genes(pb, n_top_genes=1000)
pb = pb[:, pb.var['highly_variable']].copy()

sc.tl.pca(pb)
sc.pp.neighbors(pb, n_neighbors=5, n_pcs=30)
sc.tl.umap(pb)

# --------------- 参数设置 ---------------
PCS        = (1, 2)          # 要展示的PC轴
COLOR_BY   = 'sample'         # 样本着色列
GENE_CLR   = 'firebrick'
ARROW_KW   = dict(ls='-', lw=1.4, alpha=.9)
LABEL_PAD  = 1.08            # 基因标签偏移量
LEN_FRAC   = .85             # 最长箭头占样本分布范围的比例
TOP_MARKERS_PER_GROUP = 5    # 每个样本簇选Top5标记基因

# --------------- 1. 计算样本簇标记基因 ---------------
sc.tl.rank_genes_groups(pb, groupby=COLOR_BY, method='wilcoxon')
markers = sc.get.rank_genes_groups_df(pb, group=None)
# 每个簇选Top5显著标记基因
top_markers = markers.groupby('group').head(TOP_MARKERS_PER_GROUP)['names'].tolist()
# 获取标记基因在var中的索引
idxs = [pb.var_names.get_loc(gene) for gene in top_markers if gene in pb.var_names]

# --------------- 2. 提取PCA得分与基因相关性向量 ---------------
c1, c2    = (i - 1 for i in PCS)
S         = pb.obsm['X_pca'][:, [c1, c2]]            
L_raw     = pb.varm['PCs'][:, [c1, c2]]             
eigvals   = pb.uns['pca']['variance'][[c1, c2]]
L_corr    = L_raw * np.sqrt(eigvals)                  # 基因表达与PCA得分的相关性

# --------------- 3. 缩放箭头到合适范围 ---------------
max_r_gene  = np.sqrt((L_corr[idxs]**2).sum(1)).max()
max_r_score = np.abs(S).max()
scale       = (LEN_FRAC * max_r_score) / max_r_gene
V           = L_corr[idxs] * scale

# --------------- 4. 绘图（保留原绘图逻辑）---------------
fig, ax = plt.subplots(figsize=(15, 13))

# 绘制样本散点
codes = pb.obs[COLOR_BY].astype('category').cat.codes
scat  = ax.scatter(S[:, 0], S[:, 1], c=codes, cmap='tab20',
                   s=200, edgecolor='k', alpha=.85)
# 图例设置
handles, _ = scat.legend_elements(prop='colors')
labels     = pb.obs[COLOR_BY].astype('category').cat.categories
ax.legend(handles, labels, title=COLOR_BY, loc='center left',
          bbox_to_anchor=(1.05, .5), frameon=False)

# 绘制基因箭头与标签
texts = []
for (dx, dy), gene_name in zip(V, top_markers):
    ax.arrow(0, 0, dx, dy, color=GENE_CLR, **ARROW_KW,
             head_width=.03*scale/max_r_score,
             head_length=.04*scale/max_r_score,
             length_includes_head=True, zorder=3)
    texts.append(ax.text(dx*LABEL_PAD, dy*LABEL_PAD,
                         gene_name, color=GENE_CLR,
                         fontsize=12, ha='center', va='center', zorder=4))

adjust_text(texts, arrowprops=dict(arrowstyle='-', color=GENE_CLR))

# 坐标轴与标题设置
vr = pb.uns['pca']['variance_ratio'][[c1, c2]] * 100
ax.set_xlabel(f'PC{PCS[0]}  ({vr[0]:.1f} % 方差)', fontsize=13)
ax.set_ylabel(f'PC{PCS[1]}  ({vr[1]:.1f} % 方差)', fontsize=13)
ax.axhline(0, lw=.5, color='grey'); ax.axvline(0, lw=.5, color='grey')
ax.set_aspect('equal')
ax.set_title('PCA双图 - 样本簇标记基因')
plt.tight_layout(); plt.show()

效果说明

修正后的代码会：

自动筛选每个样本簇的Top标记基因，确保箭头对应的是与簇高度关联的基因
箭头方向严格对应基因高表达样本所在的PCA区域（比如A簇的标记基因箭头会指向A簇样本的聚集区域）
保留原代码的绘图美观性，同时解决箭头指向错误的问题

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

PCA双图绘制问题：基因箭头与样本表达关联不符

需求说明 #

问题描述 #

原代码（PCA聚类正常，但箭头无意义） #

模拟数据生成代码 #

问题原因与修正方案 #

问题核心 #

修正步骤 #

修正后代码 #

效果说明 #