1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
| # =============================================================================
# 矫正sample_id + 考虑sample type（更精细）
# =============================================================================
import scvi
import scanpy as sc
adata = all_obj[all_obj.obs["cancer_type"] == "BRCA"].copy()
print(adata)
# ============ 1. 检查数据 ============
print("=== 数据概览 ===")
print("Dataset分布:\n", adata.obs['dataset'].value_counts())
print("\nSample分布:\n", adata.obs['sample_id'].value_counts())
print("\nTissue分布:\n", adata.obs['sample_type'].value_counts())
# ============ 2. 数据准备 ============
# 确保X是原始counts
adata.layers["counts"] = adata.X.copy()
# 预处理：筛选高变基因
sc.pp.highly_variable_genes(
adata,
n_top_genes=3000,
subset=True,
layer="counts",
flavor="seurat_v3",
batch_key="sample_id"
)
print(f"\n保留了 {adata.n_vars} 个高变基因")
# ============ 3. 设置 scVI 模型 ============
scvi.model.SCVI.setup_anndata(
adata,
layer="counts",
batch_key="sample_id", #主批次效应
categorical_covariate_keys=["sample_type"] #协变量
)
# ============ 4. 建立模型 ============
model = scvi.model.SCVI(
adata,
n_layers=2,
n_latent=30,
gene_likelihood="nb"
)
print("\n=== 模型信息 ===")
print(model)
# ============ 5. 训练模型 ============
print("\n=== 开始训练 ===")
model.train(
max_epochs=400,
early_stopping=True,
early_stopping_patience=20,
train_size=0.9,
batch_size=128
)
# 查看训练历史
import matplotlib.pyplot as plt
train_history = model.history
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.plot(train_history['elbo_train'], label='Train ELBO')
plt.plot(train_history['elbo_validation'], label='Validation ELBO')
plt.xlabel('Epoch')
plt.ylabel('ELBO')
plt.legend()
plt.title('Training History')
plt.subplot(1, 2, 2)
plt.plot(train_history['reconstruction_loss_train'], label='Train Recon Loss')
plt.plot(train_history['reconstruction_loss_validation'], label='Val Recon Loss')
plt.xlabel('Epoch')
plt.ylabel('Reconstruction Loss')
plt.legend()
plt.title('Reconstruction Loss')
plt.tight_layout()
plt.savefig("F:/Work/Age_LN/AgeLN_BRCA/scVI_training_history.pdf", dpi=300)
plt.show()
# ============ 6. 获取结果 ============
#adata.layers["scvi_normalized"] = model.get_normalized_expression(library_size=1e4)
# 6. 获取潜在表示（去批次后的embedding）
adata.obsm["X_scVI"] = model.get_latent_representation()
print(adata.obsm["X_scVI"].shape) # 应该是 (细胞数, 30)
# 7. 保存模型和数据
model.save("F:/Work/Age_LN/AgeLN_BRCA/scVI_model_BATCHid_COvatype", overwrite=True)
adata.raw = None
adata.write_h5ad("F:/Work/Age_LN/AgeLN_BRCA/data/AgeLN_qc_BRCA_BATCHsampleid_COvasampletype.h5ad")
# ============ 读取完整基因矩阵作为 adata_raw (因为adata做scvi只取了3000个高变基因)============
#adata = sc.read_h5ad("F:/Work/Age_LN/AgeLN_BRCA/data/AgeLN_qc_BRCA_scvi.h5ad")
#adata_ba = sc.read_h5ad("F:/Work/Age_LN/AgeLN_BRCA/data/AgeLN_qc_BRCA_nobatch.h5ad")
adata_ba = adata_raw
print("adata_ba 基因名:", adata_ba.var_names[:10])
assert (adata.obs_names == adata_ba.obs_names).all(), "细胞顺序不一致！"
print("细胞顺序验证通过，共", len(adata.obs_names), "个细胞")
adata_ba.obsm["X_scVI"] = adata.obsm["X_scVI"].copy()
for col in adata.obs.columns:
if col not in adata_ba.obs.columns:
adata_ba.obs[col] = adata.obs[col].values
# ============ 7. 下游分析 ============
sc.pp.neighbors(adata_ba, use_rep='X_scVI', n_neighbors=15)
sc.tl.umap(adata_ba, min_dist=0.3)
sc.tl.leiden(adata_ba, resolution=0.8)
# ============ 8. 可视化检查 ============
fig, axes = plt.subplots(2, 2, figsize=(18, 12))
sc.pl.umap(adata_ba, color='dataset', ax=axes[0, 0], show=False,
title='Dataset (should be mixed)')
sc.pl.umap(adata_ba, color='sample_id', ax=axes[0, 1], show=False,
title='Sample ID', legend_loc=None)
sc.pl.umap(adata_ba, color='sample_type', ax=axes[1, 0], show=False,
title='Tissue type (should be separated)')
sc.pl.umap(adata_ba, color='leiden', ax=axes[1, 1], show=False,
title='Clusters', legend_loc='on data')
plt.tight_layout()
plt.savefig("F:/Work/Age_LN/AgeLN_BRCA/scVI_integration_QC.png", dpi=300)
plt.show()
sc.pl.umap(adata_ba, color='sample_id', show=False,
title='Sample ID')
plt.savefig("F:/Work/Age_LN/AgeLN_BRCA/scVI_integration_QC_sampleid.png", dpi=300, bbox_inches="tight")
|