PCA со стрелками направления, показывающими гены, связанные с образцами ⇐ Python

[Anonymous]

Я пытаюсь воспроизвести биплат PCA, где
dots = образцы, окрашенные образец столбца в pb.obs и
стрелки = гены, указывая на область (s) на графике, где каждый ген наиболее высоко экспрессирован. Pd.x < /code> < /p>
Итак - если ген связан с определенной выборкой, на графике PCA направление стрелки будет направлено на эту выборку.
Если ген высокий в двух образцах, я ожидаю, что стрелка окажется между двумя областями образцов; Если он высокий в трех моделях, стрелка должна измениться соответственно .. и т. Д. Пример фотографии, который я сделал из Google, прекрасно показывает, что мне нужно:

Код, который я пробовал, но результаты в стрелках без значения во всех. Сам самостоятельно в порядке: < /p>
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 pickk 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, 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()

Dumy Data, которые вы можете использовать, pb выглядит нечто похожее на это:
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)


Подробнее здесь: https://stackoverflow.com/questions/796 ... th-samples