scGPT 分析基因调控网络 (Gene Regulatory Network, GRN) 的核心思路是：它不是在模拟基因间的线性相关，而是利用预训练模型从海量数据中学习到的“基因语义”来推断调控关系。其基本假设是，功能或调控行为相似的基因，在 scGPT 的向量空间中也拥有相近的“基因嵌入表示”。
- 无需微调：GRN 推断可直接在预训练模型上进行 (Zero-shot)，效率很高。下面这些预训练的网络拿过来可以直接使用：
- 方法的局限性：基于嵌入的 GRN 推断可能更偏向于捕捉“共表达”模式，而非“因果调控”的关系，类似于WGCNA而不是SCENIC。
基于以上两点可知，scGPT的GRN分析用来识别特异的调控模块是一个不错的选择，而不能确定驱动特定调控网络的TF。
import warnings
warnings.filterwarnings('ignore')
import json
from pathlib import Path
import torch
from anndata import AnnData
import scanpy as sc
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import networkx as nx
import pandas as pd
import tqdm
import gseapy as gp
import scgpt as scg
from scgpt.tasks import GeneEmbedding
from scgpt.tokenizer.gene_tokenizer import GeneVocab
from scgpt.model import TransformerModel
from scgpt.preprocess import Preprocessor
from scgpt.utils import set_seed
set_seed(42)
pad_token = "<pad>"
special_tokens = [pad_token, "<cls>", "<eoc>"]
n_hvg = 1200
n_bins = 51
mask_value = -1
pad_value = -2
n_input_bins = n_bins
从预训练的模型中选择合适的模型：
# Specify model path; here we load the pre-trained scGPT blood model
model_dir = Path("../pretrained_model/blood")
model_config_file = model_dir / "args.json"
model_file = model_dir / "best_model.pt"
vocab_file = model_dir / "vocab.json"
vocab = GeneVocab.from_file(vocab_file)
for s in special_tokens:
if s not in vocab:
vocab.append_token(s)
# Retrieve model parameters from config files
with open(model_config_file, "r") as f:
model_configs = json.load(f)
print(f"Resume model from {model_file}, the model args will override the config {model_config_file}.")
embsize = model_configs["embsize"]
nhead = model_configs["nheads"]
d_hid = model_configs["d_hid"]
nlayers = model_configs["nlayers"]
n_layers_cls = model_configs["n_layers_cls"]
gene2idx = vocab.get_stoi()
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
ntokens = len(vocab) # size of vocabulary
model = TransformerModel(
ntokens,
embsize,
nhead,
d_hid,
nlayers,
vocab=vocab,
pad_value=pad_value,
n_input_bins=n_input_bins,
)
try:
model.load_state_dict(torch.load(model_file))
print(f"Loading all model params from {model_file}")
except:
# only load params that are in the model and match the size
model_dict = model.state_dict()
pretrained_dict = torch.load(model_file)
pretrained_dict = {k: v for k, v in pretrained_dict.items() if k in model_dict and v.shape == model_dict[k].shape}
for k, v in pretrained_dict.items():
print(f"Loading params {k} with shape {v.shape}")
model_dict.update(pretrained_dict)
model.load_state_dict(model_dict)
model.to(device)
读取并预处理数据，包括过滤细胞和基因、标准化、归一化、分箱：
adata = sc.read("Immune_ALL_human.h5ad") # 33506 × 12303
ori_batch_col = "batch"
adata.obs["celltype"] = adata.obs["final_annotation"].astype(str)
data_is_raw = False
preprocessor = Preprocessor(
use_key="X", # the key in adata.layers to use as raw data
filter_gene_by_counts=3, # step 1
filter_cell_by_counts=False, # step 2
normalize_total=1e4, # 3. whether to normalize the raw data and to what sum
result_normed_key="X_normed", # the key in adata.layers to store the normalized data
log1p=data_is_raw, # 4. whether to log1p the normalized data
result_log1p_key="X_log1p",
subset_hvg=n_hvg, # 5. whether to subset the raw data to highly variable genes
hvg_flavor="seurat_v3" if data_is_raw else "cell_ranger",
binning=n_bins, # 6. whether to bin the raw data and to what number of bins
result_binned_key="X_binned", # the key in adata.layers to store the binned data
)
preprocessor(adata, batch_key="batch")
提取基因embeddings，并构建GeneEmbedding类的实例对象：
gene_ids = np.array([id for id in gene2idx.values()])
gene_embeddings = model.encoder(torch.tensor(gene_ids, dtype=torch.long).to(device))
gene_embeddings = gene_embeddings.detach().cpu().numpy()
gene_embeddings = {gene: gene_embeddings[i] for i, gene in enumerate(gene2idx.keys()) if gene in adata.var.index.tolist()}
print('Retrieved gene embeddings for {} genes.'.format(len(gene_embeddings)))
embed = GeneEmbedding(gene_embeddings)
提取基因程序，对基因程序向量构成的矩阵进行图聚类和降维，过程与处理表达谱类：
gdata = embed.get_adata(resolution=20)
metagenes = embed.get_metagenes(gdata)
mgs = dict()
for mg, genes in metagenes.items():
if len(genes) > 4:
mgs[mg] = genes
embed对象的score_metagenes方法用来计算基因程序在每个细胞中的活性，底层调用的是scanpy.tl.score_genes。其本身也有一个方法plot_metagenes_scores可以绘制热图来可视化不同细胞类型中基因程序的平均活性，但由于画出的热图有点奇怪，图例与图间隔大，还会多画两个空白页，看起来有点影响心情，这里修改一下画图函数：
embed.score_metagenes(adata, metagenes)
def plot_metagenes_scores(adata, metagenes, column, plot=None, figsize=(6, 13)):
matrix = []
meta_genes = []
cfnum = 1
cfams = dict()
for cluster, vector in metagenes.items():
row = []
cts = []
for ct in set(adata.obs[column]):
sub = adata[adata.obs[column] == ct]
val = np.mean(sub.obs[str(cluster) + "_SCORE"].tolist())
row.append(val)
cts.append(ct)
matrix.append(row)
label = str(cluster) + "_SCORE: " + ", ".join(vector[:5])
if len(set(vector)) > 10:
label += "*"
meta_genes.append(label)
cfams[cluster] = label
cfnum += 1
matrix = np.array(matrix)
df = pd.DataFrame(matrix, index=meta_genes, columns=cts)
fig = sns.clustermap(
df,
figsize=figsize,
cbar_pos=(0.01, 0.908, 0.02, 0.08), # x, y, width%, height%
dendrogram_ratio=0.11,
cmap="mako",
yticklabels=True,
standard_scale=0,
)
if plot:
plt.savefig(plot)
else:
return fig
plot_metagenes_scores(adata, dict(list(mgs.items())[:20]), "celltype", figsize=(7, 7))
热图可以同时展示很多基因程序在各种细胞类型的情况，而网络图可以进一步展示感兴趣的基因程序中基因间的相互关系：
CD_genes = mgs['5']
print(CD_genes)
# Compute cosine similarities among genes in this gene program
list_of_dfs = []
for i in tqdm.tqdm(CD_genes):
df = embed.compute_similarities(i, CD_genes)
df['Gene1'] = i
list_of_dfs.append(df)
df_CD = pd.concat(list_of_dfs, ignore_index=True)
df_CD_sub = df_CD[df_CD['Similarity'] < 1].sort_values(by='Gene') # Filter out edges from each gene to itself
# Creates a graph from the cosine similarity network
input_node_weights = [(row['Gene'], row['Gene1'], round(row['Similarity'], 2)) for i, row in df_CD_sub.iterrows()]
G = nx.Graph()
G.add_weighted_edges_from(input_node_weights)
thresh = 0.4
plt.figure(figsize=(8, 8))
widths = nx.get_edge_attributes(G, 'weight')
elarge = [(u, v) for (u, v, d) in G.edges(data=True) if d["weight"] > thresh]
esmall = [(u, v) for (u, v, d) in G.edges(data=True) if d["weight"] <= thresh]
width_large = {}
width_small = {}
for i, v in enumerate(list(widths.values())):
if v > thresh:
width_large[list(widths.keys())[i]] = v*10
else:
width_small[list(widths.keys())[i]] = max(v, 0)*10
pos = nx.spring_layout(G, k=0.4, iterations=15, seed=3)
nx.draw_networkx_edges(G, pos, edgelist = width_small.keys(), width=list(width_small.values()), edge_color='lightblue', alpha=0.8)
nx.draw_networkx_edges(G, pos, edgelist = width_large.keys(), width = list(width_large.values()), alpha = 0.5, edge_color = "blue")
nx.draw_networkx_labels(G, pos, font_size=10) # # node labels
d = nx.get_edge_attributes(G, "weight") # edge weight labels
edge_labels = {k: d[k] for k in elarge}
nx.draw_networkx_edge_labels(G, pos, edge_labels, font_size=10)
plt.axis("off")
plt.show()
最后，将得到的基因程序做个富集分析方便知道其中蕴含的生物学意义：
# Perform pathway enrichment analysis using the gseapy package in the Reactome database
enr_Reactome = gp.enrichr(gene_list=CD_genes, gene_sets=['Reactome_2022'], organism='human', outdir=None, cutoff=0.5) # cutoff参数主要影响画图
out = enr_Reactome.results
out.head()
Gene_set Term Overlap P-value Adjusted P-value Old P-value Old Adjusted P-value Odds Ratio Combined Score Genes
0 Reactome_2022 Immunoregulatory Interactions Between A Lympho... 5/123 6.201738e-09 3.969112e-07 0 0 120.278450 2273.074626 IFITM1;CD8B;CD8A;CD3E;CD3D
1 Reactome_2022 Translocation Of ZAP-70 To Immunological Synap... 3/17 1.116740e-07 2.707479e-06 0 0 475.571429 7612.796330 LCK;CD3E;CD3D
2 Reactome_2022 Immune System R-HSA-168256 9/1943 1.269131e-07 2.707479e-06 0 0 28.005171 444.715482 IL32;IFITM1;CD8B;CD8A;LCK;LTB;CD3E;IL7R;CD3D
3 Reactome_2022 Phosphorylation Of CD3 And TCR Zeta Chains R-H... 3/20 1.870300e-07 2.792044e-06 0 0 391.588235 6066.483731 LCK;CD3E;CD3D
4 Reactome_2022 PD-1 Signaling R-HSA-389948 3/21 2.181285e-07 2.792044e-06 0 0 369.814815 5672.286820 LCK;CD3E;CD3D
