GeneKI：原文说“吊打”scTenifoldKnk｜仅需WT scRNA数据直接进行目标基因虚拟敲除

#创建环境
#安装GenKI
#GenKI requires Python ≥ 3.10. 
conda create -n geneki python=3.10
conda activate geneki
#安装
pip install scanpy -i https://pypi.tuna.tsinghua.edu.cn/simple
pip install GenKI -i https://pypi.tuna.tsinghua.edu.cn/simple

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scanpy as sc

import GenKI as gk
from GenKI.preprocessing import build_adata
from GenKI.dataLoader import DataLoader
from GenKI.train import VGAE_trainer
from GenKI import utils

#读入数据
adata = build_adata("./microglial_seurat_WT.h5ad")
adata = adata[:, :].copy()
adata

adata.obs['celltype'] = 'microglial'

#后续需要使用多线程，需要先安装：
pip install "GenKI[ray]"

# load data — build the GRN and prepare WT / virtual-KO graph data
data_wrapper = DataLoader(
    adata,                  # AnnData object,前面处理好的表达矩阵anndata
    target_gene=["TREM2"],  # KO gene name(s)，需要敲除的任意基因名称
    target_cell=None,       # obs label for cell-type subset; None = all cells
                            #也即是说，你的anndata是所有celltype的时候，可以设置需要分析的celltype
                            #无需提取数据，这里因为是提取的数据，都是一种celltype，所以设置None

    obs_label="ident",      # adata.obs column for cell-type labels，细胞注释列
    GRN_file_dir="GRNs",    # folder for cached GRN .npz，储存GRN .npz的文件夹名称
    rebuild_GRN=True,       # build GRN if not cached
    n_cpus=8,               # Ray sharding; needs `pip install "GenKI[ray]"`
)
data_wt = data_wrapper.load_data()
data_ko = data_wrapper.load_kodata()

#初始化并配置一个VGAE模型的训练器，用于在WT数据上训练模型
hyperparams = {"epochs": 100, 
               "lr": 7e-4, 
               "beta": 1e-4, 
               "seed": 8096}
log_dir=None 
sensei = VGAE_trainer(data_wt,
                     epochs=hyperparams["epochs"], 
                     lr=hyperparams["lr"], 
                     log_dir=log_dir, 
                     beta=hyperparams["beta"], 
                     seed=hyperparams["seed"],
                     verbose=True,
                     )

#模型训练
sensei.train()

#保存训练模型
#sensei.save_model('model_example')

# get distance between wt and ko
z_mu_wt, z_std_wt = sensei.get_latent_vars(data_wt)
z_mu_ko, z_std_ko = sensei.get_latent_vars(data_ko)
dis = gk.utils.get_distance(z_mu_ko, z_std_ko, z_mu_wt, z_std_wt, by="KL")
print(dis.shape)

#根据扰动得分对基因进行排序，排名较前的说明敲除目标基因后收到的扰动越大
res_raw = utils.get_generank(data_wt, dis, rank=True)
res_raw.head()

#permutation test
null = sensei.pmt(data_ko, n=100, by="KL")
res = utils.get_generank(data_wt, dis, null,
                        save_significant_as = 'MC_KOgenes'#将扰动基因保存到文件
                        )
res

import pandas as pd
import gseapy as gp
import matplotlib.pyplot as plt

gene_list = res.index.to_list()
enr = gp.enrichr(gene_list=gene_list, # gene形式需要与数据库保持一致
                 gene_sets=['GO_Biological_Process_2025','KEGG_2021_Human'],#需要分析的基因集,可以同时分析GO、KEGG两个gene set
                 organism='human', #物种
                 outdir=None)

from gseapy import barplot, dotplot

# categorical scatterplot
ax = dotplot(enr.results,#分析结果
              column="Adjusted P-value",#按照adj pvalue展示结果
              x='Gene_set', # set x axis, so you could do a multi-sample/library comparsion
              size=10,#点大小
              top_term=10,#展示的top5terms
              figsize=(3,5),#图长宽
              title = "Enrichment",#标题
              xticklabels_rot=45, # rotate xtick labels
              show_ring=False, # set to False to revmove outer ring
              marker='o',#plot图形形式，参照https://matplotlib.org/stable/api/markers_api.html
             )

# categorical barplot
ax = barplot(enr.results,
              column="Adjusted P-value",
              group='Gene_set', # set group, so you could do a multi-sample/library comparsion
              size=10,
              top_term=10,
              figsize=(3,5),
              color = {'GO_Biological_Process_2025': 'salmon','KEGG_2021_Human':'darkblue'}
             )