geneformer | 基因扰动分析

如果某个基因被删除或者过表达，细胞内会发生什么变化？哪些通路会被抑制或者激活？哪些基因表达会随之上升或下降？这类问题以前只能靠实验来验证，成本高、周期长。现在随着工具的发展，各种虚拟扰动的工具被开发出来，我们可以在计算机中对任意基因进行“虚拟扰动”，无需基因编辑、不用培养细胞，只需scRNA-seq数据作为输入，依靠软件在电脑里模拟基因的删除或者过表达，来预测下游调控、细胞命运与功能通路。

有什么好的想法，通过虚拟扰动来做预实验去验证可行性，相比于湿实验来说，大大地降低了成本和周期。今天咱们就来体验一下用geneformer做基因扰动。

geneformer做扰动分析大体分成两种情况：有参考状态，即扰动后细胞是否有初始状态偏向目标状态；无参考状态，即扰动后细胞嵌入状态的变化。底层原理的本质就是通过改变基因token的排序来计算扰动后细胞状态的改变。

参数emb_mode初始化时有三个选项：cls、cell、gene，模型版本为V1时选择cell，V2时可以选择cls。

from geneformer import InSilicoPerturber
from geneformer import InSilicoPerturberStats
from geneformer import EmbExtractor

cell_states_to_model = {"state_key": "disease", "start_state": "dcm", "goal_state": "nf", "alt_states": ["hcm"]}
filter_data_dict = {"cell_type": ["Cardiomyocyte1","Cardiomyocyte2","Cardiomyocyte3"]}

# OF NOTE: model_version should match version of model to be used (V1 or V2) to use the correct token dictionary
embex = EmbExtractor(model_type="CellClassifier",
                     num_classes=3,
                     emb_mode = 'cell',
                     filter_data=filter_data_dict,
                     max_ncells=1000,
                     emb_layer=0,
                     summary_stat="exact_mean",
                     forward_batch_size=256,
                     model_version="V1",
                     nproc=16)

state_embs_dict = embex.get_state_embs(cell_states_to_model,
                                       "../Geneformer-V1-10M",
                                       "../genecorpus-30m/example_input_files/cell_classification/disease_classification/human_dcm_hcm_nf.dataset",
                                       "output_perturbation",
                                       "my_perturbation")

上面的代码用来获取细胞嵌入，仅在有参考状态时需要，如果没有参考状态可以省略，并且无参考也需要省略cell_states_to_model、state_embs_dict两个参数。

gene_list = ['ENSG00000155657', 'ENSG00000198626', 'ENSG00000154556', 'ENSG00000129116', 'ENSG00000163110', 'ENSG00000115641']
isp = InSilicoPerturber(perturb_type="delete",
                        perturb_rank_shift=None,
                        genes_to_perturb=gene_list,
                        model_type="CellClassifier",
                        num_classes=3,
                        emb_mode="cell",
                        cell_emb_style="mean_pool",
                        filter_data=filter_data_dict,
                        cell_states_to_model=cell_states_to_model,
                        state_embs_dict=state_embs_dict,
                        max_ncells=2000,
                        emb_layer=0,
                        forward_batch_size=400,
                        model_version="V1",
                        nproc=16)

默认会扰动所有基因，即基因token列表里所有的基因。如果没有GPU，扰动时会非常耗时。可以提供一个需要扰动的基因列表，但需要注意的是，基因列表会作为一个整体，如果想单独扰动一个基因可以用循环。基因名只能用ensembl ID，这个由基因token列表决定的。

isp.perturb_data("../Geneformer-V1-10M",
                 "../genecorpus-30m/example_input_files/cell_classification/disease_classification/human_dcm_hcm_nf.dataset",
                 "output_perturbation",
                 "my_perturbation")
Some weights of BertForSequenceClassification were not initialized from the model checkpoint at ../Geneformer-V1-10M and are newly initialized: ['bert.pooler.dense.bias', 'bert.pooler.dense.weight', 'classifier.bias', 'classifier.weight']
100%|████████████████████████████████████████████████| 5/5 [04:42<00:00, 44.52s/it]

上面这一步是核心步骤，也是耗时最长的步骤。如果是CPU环境，设置16个核，扰动一个基因需要5分钟左右。如果需要全局扰动来筛选的话，还是找个GPU环境吧。

ispstats = InSilicoPerturberStats(mode="goal_state_shift", genes_perturbed=gene_list, cell_states_to_model=cell_states_to_model, model_version="V1")
ispstats.get_stats("output_perturbation", None, "output_perturbation", "my_perturbation")

接着，对扰动结果进行统计，mode参数决定了如何统计扰动效应，根据需要选择即可。

pst = pu.pd.read_csv('output_perturbation/my_perturbation.csv', index_col=0)
pst
   Gene Gene_name       Ensembl_ID  Shift_to_goal_end  Shift_to_alt_end_hcm  Goal_end_vs_random_pval  Alt_end_vs_random_pval_hcm  Goal_end_FDR  Alt_end_FDR_hcm  N_Detections  Sig
0  9837       TTN  ENSG00000155657          -0.000469             -0.000307                      1.0                         1.0           1.0              1.0          2000    0

扰动结果是一张表格，可以查看扰动是否引起细胞状态的显著性变化。从上面的这个结果可以看出，删除这个基因并没有引起显著性变化。

有无参考状态，以及扰动模式都会影响输出的结果，以下是CSV输出文件中可能列的定义。

注意：并非所有列都会出现在所有输出文件中。某些列特定于特定的扰动模式。