{ "cells": [ { "cell_type": "markdown", "id": "4885915e", "metadata": {}, "source": [ "# In-Silico Perturbation" ] }, { "cell_type": "code", "execution_count": null, "id": "e54fe7f9", "metadata": {}, "outputs": [], "source": [ "import CancerStFormer as cstformer\n", "from cstformer.perturbation.perturb import Perturber\n", "from cstformer.perturbation.perturb_stats import PerturberStats\n", "import pandas as pd \n", "import pickle\n", "import os\n", "import re" ] }, { "cell_type": "markdown", "id": "3ad23a7c", "metadata": {}, "source": [ "### 0.1 Parameter Descriptions" ] }, { "cell_type": "raw", "id": "6aea863b", "metadata": { "vscode": { "languageId": "raw" } }, "source": [ "**Parameters:**\n", " mode: {'spot','neighborhood'} \n", " | Type of dataset tokenization option \n", " | 'spot': if dataset is tokenized genes from only spot \n", " | 'neighborhood': if dataset is tokenized from spot and surrounding neighbors\n", " perturb_type : {\"delete\", \"overexpress\", \"inhibit\", \"activate\"}\n", " | Type of perturbation.\n", " | \"delete\": delete gene from rank value encoding\n", " | \"overexpress\": move gene to front of rank value encoding\n", " combos : {0,1}\n", " | Whether to perturb genes individually (0) or in pairs (1).\n", " anchor_gene : None, str\n", " | ENSEMBL ID of gene to use as anchor in combination perturbations.\n", " | For example, if combos=1 and anchor_gene=\"ENSG00000148400\":\n", " | anchor gene will be perturbed in combination with each other gene.\n", " model_type : {\"Pretrained\", \"GeneClassifier\", \"CellClassifier\"}\n", " | Whether model is the pretrained Geneformer or a fine-tuned gene or cell classifier.\n", " num_classes : int\n", " | If model is a gene or cell classifier, specify number of classes it was trained to classify.\n", " | For the pretrained Geneformer model, number of classes is 0 as it is not a classifier.\n", " emb_mode : {\"cell\", \"cell_and_gene\"}\n", " | Whether to output impact of perturbation on cell and/or gene embeddings.\n", " | Gene embedding shifts only available as compared to original cell, not comparing to goal state.\n", " cell_emb_style : \"mean_pool\"\n", " | Method for summarizing cell embeddings.\n", " | Currently only option is mean pooling of gene embeddings for given cell.\n", " filter_data : None, dict\n", " | Default is to use all input data for in silico perturbation study.\n", " | Otherwise, dictionary specifying .dataset column name and list of values to filter by.\n", " cell_states_to_model : None, dict\n", " | Cell states to model if testing perturbations that achieve goal state change.\n", " | Four-item dictionary with keys: state_key, start_state, goal_state, and alt_states\n", " | state_key: key specifying name of column in .dataset that defines the start/goal states\n", " | start_state: value in the state_key column that specifies the start state\n", " | goal_state: value in the state_key column taht specifies the goal end state\n", " | alt_states: list of values in the state_key column that specify the alternate end states\n", " | For example: {\"state_key\": \"disease\",\n", " | \"start_state\": \"dcm\",\n", " | \"goal_state\": \"nf\",\n", " | \"alt_states\": [\"hcm\", \"other1\", \"other2\"]}\n", " state_embs_dict : None, dict\n", " | Embedding positions of each cell state to model shifts from/towards (e.g. mean or median).\n", " | Dictionary with keys specifying each possible cell state to model.\n", " | Values are target embedding positions as torch.tensor.\n", " | For example: {\"nf\": emb_nf,\n", " | \"hcm\": emb_hcm,\n", " | \"dcm\": emb_dcm,\n", " | \"other1\": emb_other1,\n", " | \"other2\": emb_other2}\n", " max_ncells : None, int\n", " | Maximum number of cells to test.\n", " | If None, will test all cells.\n", " cell_inds_to_perturb : \"all\", list\n", " | Default is perturbing each cell in the dataset.\n", " | Otherwise, may provide a dict of indices of cells to perturb with keys start_ind and end_ind.\n", " | start_ind: the first index to perturb.\n", " | end_ind: the last index to perturb (exclusive).\n", " | Indices will be selected *after* the filter_data criteria and sorting.\n", " | Useful for splitting extremely large datasets across separate GPUs.\n", " emb_layer : {-1, 0}\n", " | Embedding layer to use for quantification.\n", " | 0: last layer (recommended for questions closely tied to model's training objective)\n", " | -1: 2nd to last layer (recommended for questions requiring more general representations)\n", " forward_batch_size : int\n", " | Batch size for forward pass.\n", " nproc : int\n", " | Number of CPU processes to use.\n", " token_dictionary_file : Path\n", " | Path to pickle file containing token dictionary (Ensembl ID:token)" ] }, { "cell_type": "markdown", "id": "0539f30a", "metadata": {}, "source": [ "### 1.1 Option A: Delete Single Gene and View both Gene and Cell Embedding Shifts" ] }, { "cell_type": "markdown", "id": "441eaa32", "metadata": {}, "source": [ "Here we look at the following option:\n", "1. Delete Gene: PDCD1\n", "2. Load in pretrained model on masked learning objective\n", "3. want to look at the cosine similarity to null distribution when perturbing PDCD1 in both per gene and per cell embedding\n", "4. batch our perturbations at each model forward pass\n" ] }, { "cell_type": "code", "execution_count": null, "id": "aa0ecf37", "metadata": {}, "outputs": [], "source": [ "isp = Perturber(\n", " mode = 'extended', #perturbing either spot or extended model (based on tokenization)\n", " perturb_type=\"delete\", # in-silico deletion\n", " genes_to_perturb=['PDCD1'], #immune checkpiont\n", " perturb_rank_shift = None,\n", " model_type=\"GeneClassifier\", # GeneClassifier Model (can use CellClassifier, GeneClassifier, Pretrained (ML))\n", " num_classes=2, # classes used in classifier\n", " emb_mode=\"cell_and_gene\", # cell and gene shifts \n", " cell_emb_style='mean_pool',\n", " max_ncells=1000, #num cells use in perturbation\n", " emb_layer=-1, \n", " forward_batch_size=80, #minibatch size\n", " nproc=1, # num threads\n", " token_dictionary_file='new_token_dictionary.pickle', # gene -> token mapping\n", " )\n", " \n", "isp.perturb_dataset(\n", " model_directory='run-8eb93bdf/checkpoint-1000', #which final model to choose\n", " input_data_file='STFormer_TNBC_neighbor.dataset', #tokenized dataset\n", " output_directory='output/perturb/cell_shift', \n", " output_prefix='perturb_extended')" ] }, { "cell_type": "markdown", "id": "06ab3e15", "metadata": {}, "source": [ "We want to create stats to evaluate the effect of PDCD1 gene expression on:\n", "1. Cell embeddings similarity to unperturbed cells\n", "2. Gene embedding similarity to unperturbed genes, where we aggregate the gene shifts for our cells for each other token in our dataset\n", " - ex: how does deletion of PDCD1 effect the expression of another gene like CCR5AS\n", " - These results are ranked by mean cosine similarity (lowest -> highest) where low similarity signifies a greater impact on gene expression and a high dependency on the gene expression of our perturbed gene\n", " - provides standeard deviation of cosine similarity between gene embeddings and the number of detections" ] }, { "cell_type": "code", "execution_count": null, "id": "1a9bc21e", "metadata": {}, "outputs": [], "source": [ "from cstformer.perturbation.perturb_stats import PerturberStats\n", "\n", "ispstats = PerturberStats(\n", " mode='aggregate_gene_shifts', # view gene shifts due to perturbation\n", " gene_id_name_dict='gene_id_dictionary.pickle', # maps gene -> ENS (prints in gene_name)\n", " token_dictionary_file = 'new_token_dictionary.pickle' # gene -> token dict\n", ")\n", "ispstats.compute_stats(\n", " input_data_dir = 'output/perturb', # where embeddings are located\n", " null_data_dir=None, # compare to null distribution\n", " output_directory = 'output/perturb/perturb_stats',\n", " output_prefix= 'perturb_extended'\n", ")" ] }, { "cell_type": "code", "execution_count": null, "id": "70c8f02e", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
PerturbedGene_nameEnsembl_IDAffectedAffected_gene_nameAffected_Ensembl_IDCosine_sim_meanCosine_sim_stdevN_Detections
0PDCD1PDCD1PDCD1NEK11NEK11NEK110.5671510.02
1PDCD1PDCD1PDCD1SYT17SYT17SYT170.6329330.02
2PDCD1PDCD1PDCD1CASS4CASS4CASS40.6451200.02
3PDCD1PDCD1PDCD1CAPN10-DTCAPN10-DTCAPN10-DT0.6796280.02
4PDCD1PDCD1PDCD1PDE1BPDE1BPDE1B0.7060050.02
..............................
15370PDCD1PDCD1PDCD1AL133453.1AL133453.1AL133453.10.9996780.02
15371PDCD1PDCD1PDCD1ACP7ACP7ACP70.9997010.02
15372PDCD1PDCD1PDCD1AL589182.1AL589182.1AL589182.10.9997100.02
15373PDCD1PDCD1PDCD1PIH1D2PIH1D2PIH1D20.9997180.02
15374PDCD1PDCD1PDCD1FAM218AFAM218AFAM218A0.9997740.02
\n", "

15375 rows × 9 columns

\n", "
" ], "text/plain": [ " Perturbed Gene_name Ensembl_ID Affected Affected_gene_name \\\n", "0 PDCD1 PDCD1 PDCD1 NEK11 NEK11 \n", "1 PDCD1 PDCD1 PDCD1 SYT17 SYT17 \n", "2 PDCD1 PDCD1 PDCD1 CASS4 CASS4 \n", "3 PDCD1 PDCD1 PDCD1 CAPN10-DT CAPN10-DT \n", "4 PDCD1 PDCD1 PDCD1 PDE1B PDE1B \n", "... ... ... ... ... ... \n", "15370 PDCD1 PDCD1 PDCD1 AL133453.1 AL133453.1 \n", "15371 PDCD1 PDCD1 PDCD1 ACP7 ACP7 \n", "15372 PDCD1 PDCD1 PDCD1 AL589182.1 AL589182.1 \n", "15373 PDCD1 PDCD1 PDCD1 PIH1D2 PIH1D2 \n", "15374 PDCD1 PDCD1 PDCD1 FAM218A FAM218A \n", "\n", " Affected_Ensembl_ID Cosine_sim_mean Cosine_sim_stdev N_Detections \n", "0 NEK11 0.567151 0.0 2 \n", "1 SYT17 0.632933 0.0 2 \n", "2 CASS4 0.645120 0.0 2 \n", "3 CAPN10-DT 0.679628 0.0 2 \n", "4 PDE1B 0.706005 0.0 2 \n", "... ... ... ... ... \n", "15370 AL133453.1 0.999678 0.0 2 \n", "15371 ACP7 0.999701 0.0 2 \n", "15372 AL589182.1 0.999710 0.0 2 \n", "15373 PIH1D2 0.999718 0.0 2 \n", "15374 FAM218A 0.999774 0.0 2 \n", "\n", "[15375 rows x 9 columns]" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import pandas as pd\n", "pd.read_csv('perturb_extended_emb_aggregate_gene_shifts.csv')" ] }, { "cell_type": "markdown", "id": "4735250d", "metadata": {}, "source": [ "### 1.2 Option B: Delete all Genes and compare Cell State Shifts (Group A -> Group B)" ] }, { "cell_type": "markdown", "id": "48b214d8", "metadata": {}, "source": [ "\n", "1. **Data Filtering & Batching** \n", " - `perturb_dataset()` loads and filters your input dataset. \n", " - In `isp_perturb_set()` or `isp_perturb_all()`, cells are grouped into batches of size `forward_batch_size`.\n", "\n", "2. **Define Transition** \n", " - `get_state_embs` \n", "\n", "3. **Embedding Extraction** \n", " - **Original** (`full_orig`) and **perturbed** (`full_pert`) token sequences are passed through the model (via `get_embs`) to collect hidden‐state tensors. \n", " - Because overexpressed tokens were prepended, the **first N** positions of `full_pert` are your newly-inserted genes; the rest align with the original token order.\n", "\n", "4. **Cosine‐Similarity Quantification** \n", " - **Gene-level**: `pu.quant_cos_sims(pert_emb, original_emb, …, emb_mode=\"gene\")` computes how each remaining gene embedding shifts when your target gene is overexpressed. \n", " - **Cell-level** (if `cell_states_to_model` is set): average the non-padding embeddings before & after overexpression and quantify that shift against your state-embedding targets.\n", "\n", "5. **Aggregation & Output** \n", " - For each cell or state, cosine-similarities are averaged (or bucketed) and stored in a dictionary keyed by the perturbed gene(s). \n", " - Final results are written out as “_cell_embs_dict_…” and—if `emb_mode=\"cell_and_gene\"`—also “_gene_embs_dict_…” files you can analyze with `in_silico_perturber_stats`." ] }, { "cell_type": "code", "execution_count": null, "id": "25423ce4", "metadata": {}, "outputs": [], "source": [ "cell_states_to_model = {'state_key':'classification',\n", " 'start_state': 'Invasive cancer',\n", " 'goal_state': 'Invasive cancer + lymphocytes',\n", " 'alt_states': ['Invasive cancer + stroma','Invasive cancer + stroma + lymphocytes']}" ] }, { "cell_type": "code", "execution_count": null, "id": "5804f4f5", "metadata": {}, "outputs": [], "source": [ "from cstformer.tokenization.embedding_extractor import EmbExtractor\n", "embex = EmbExtractor(\n", " model_type = 'Pretrained',\n", " num_classes=2,\n", " filter_data = None,\n", " max_ncells = 1000,\n", " emb_layer=0,\n", " summary_stat='exact_mean',\n", " forward_batch_size=100,\n", " token_dictionary_file='output/token_dictionary.pickle',\n", " nproc = 16\n", ")\n", "\n", "state_embs_dict = embex.get_state_embs(\n", " cell_states_to_model,\n", " model_directory='output/spot/models/250422_102707_cstformer_L6_E3/final',\n", " input_data_file='output/spot/visium_spot.dataset',\n", " output_dir='output/perturb/cell_shift',\n", " output_prefix='spot_state_shift'\n", ")" ] }, { "cell_type": "markdown", "id": "9840fb81", "metadata": {}, "source": [ "visualize state embeddings" ] }, { "cell_type": "code", "execution_count": 11, "id": "655a0287", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Invasive cancer\n", "Invasive cancer + lymphocytes\n", "Invasive cancer + stroma\n", "Invasive cancer + stroma + lymphocytes\n" ] } ], "source": [ "for key,value in state_embs_dict.items():\n", " print(key)" ] }, { "cell_type": "markdown", "id": "523120d6", "metadata": {}, "source": [ "Perform Perturbtion of all genes and compare cosine similarity from start state to goal state and rank according to shift to goal state" ] }, { "cell_type": "code", "execution_count": null, "id": "ae1fd7be", "metadata": {}, "outputs": [], "source": [ "from stFormer.perturbation.stFormer_perturb import Perturber\n", "from stFormer.perturbation.perturb_stats import PerturberStats" ] }, { "cell_type": "code", "execution_count": null, "id": "e224b78e", "metadata": {}, "outputs": [], "source": [ "isp = Perturber(\n", " perturb_type=\"delete\",\n", " perturb_rank_shift=None,\n", " combos=0,\n", " anchor_gene=None,\n", " genes_to_perturb='all',\n", " cell_states_to_model=cell_states_to_model,\n", " state_embs_dict = state_embs_dict,\n", " model_type=\"Pretrained\",\n", " num_classes=0,\n", " emb_mode=\"cell\",\n", " cell_emb_style='mean_pool',\n", " max_ncells=None,\n", " emb_layer=0,\n", " forward_batch_size=100,\n", " nproc=12,\n", " token_dictionary_file='output/token_dictionary.pickle',\n", " )\n", " \n", "isp.perturb_dataset(\n", " model_directory='output/spot/models/250422_102707_cstformer_L6_E3/final',\n", " input_data_file='output/spot/visium_spot.dataset',\n", " output_directory='output/perturb/cell_shift',\n", " output_prefix='perturb_spot')" ] }, { "cell_type": "code", "execution_count": null, "id": "cca50e99", "metadata": {}, "outputs": [], "source": [ "ispstats = PerturberStats(mode=\"goal_state_shift\",\n", " genes_perturbed=\"all\",\n", " combos=0,\n", " anchor_gene=None,\n", " cell_states_to_model=cell_states_to_model,\n", " token_dictionary_file='output/token_dictionary.pickle',\n", " gene_id_name_dict='output/ensembl_mapping_dict.pickle')\n", "ispstats.compute_stats(\n", " input_data_dir=\"output/perturb/cell_shift\", # this should be the directory \n", " null_data_dir=None,\n", " output_dir=\"output/perturb/cell_shift/perturb_stats\",\n", " output_prefix=\"visium_spot\")\n" ] }, { "cell_type": "markdown", "id": "9fcd71e3", "metadata": {}, "source": [ "Structure of this dataset follows:\n", "1. `Gene`: Token ID \n", "2. `Gene Name`: hgnc name for gene\n", "3. `Ensembl_ID`: matching ensembl gene id \n", "4. `Shift to goal end`: cosine shift from start state towards goal end state in response to given perturbation\n", "5. `Goal end vs random pvalue`: pvalue of cosine shift from start state towards goal end state by Wilcoxon to random distribution of max 10,000 cells\n", "6. `N Detections`: Number of cells where the perturbed gene was detected/perturbed\n", "7. `Shift to alternate end state`: Cosine shift from start state towards alternate end state in response to perturbation\n", "8. `Goal end FDR`: Benjamini Hochberg correction of Goal State vs Null Pvalue\n", "9. `Alt End FDR`: Multiple Hypothesis Test Correction of Alternate State End vs Random Pvalue \n", "10. `Sig`: Binarized False Discovery Rate below significant (FDR < 0.05)" ] }, { "cell_type": "code", "execution_count": null, "id": "304e769d", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
GeneGene_nameEnsembl_IDShift_to_goal_endGoal_end_vs_random_pvalN_DetectionsShift_to_alt_end_Invasive cancer + stromaAlt_end_vs_random_pval_Invasive cancer + stromaShift_to_alt_end_Invasive cancer + stroma + lymphocytesAlt_end_vs_random_pval_Invasive cancer + stroma + lymphocytesGoal_end_FDRAlt_end_FDR_Invasive cancer + stromaAlt_end_FDR_Invasive cancer + stroma + lymphocytesSig
29793907ANKRD37ENSG000001863520.0010140.02774150.0007120.1743320.0007500.0901400.9246710.9135860.8733150
1164617551SPINK4ENSG000001227110.0009100.17544260.0011640.1327970.0011390.1335000.9535330.8864580.8886110
788410285ARID4AENSG000000322190.0008240.023335100.0003480.6019560.0004030.4247660.9246710.9821600.9598370
949712278HOXB-AS1ENSG000002301480.0008180.02490960.0007500.0442140.0007700.0269930.9246710.8235520.7902710
1064413690PBX4ENSG000001057170.0007780.00364370.0012210.0054390.0011540.0050400.9246710.7346570.7159310
.............................................
5364TTC34ENSG00000215912-0.0007630.01367390.0000240.509069-0.0000800.6304170.9246710.9747610.9746860
50666672IDO1ENSG00000131203-0.0007820.0554777-0.0004790.125197-0.0005220.0745210.9246710.8864580.8538430
58567701CAVIN3ENSG00000170955-0.0008790.0131256-0.0004910.008692-0.0005450.0054540.9246710.7346570.7159310
479604FYB2ENSG00000187889-0.0009250.00585711-0.0001710.891572-0.0002690.6942670.9246710.9961870.9852340
46096032SHROOM2ENSG00000146950-0.0011320.0157316-0.0006310.063054-0.0007290.0407950.9246710.8480500.8282480
\n", "

11664 rows × 14 columns

\n", "
" ], "text/plain": [ " Gene Gene_name Ensembl_ID Shift_to_goal_end \\\n", "2979 3907 ANKRD37 ENSG00000186352 0.001014 \n", "11646 17551 SPINK4 ENSG00000122711 0.000910 \n", "7884 10285 ARID4A ENSG00000032219 0.000824 \n", "9497 12278 HOXB-AS1 ENSG00000230148 0.000818 \n", "10644 13690 PBX4 ENSG00000105717 0.000778 \n", "... ... ... ... ... \n", "53 64 TTC34 ENSG00000215912 -0.000763 \n", "5066 6672 IDO1 ENSG00000131203 -0.000782 \n", "5856 7701 CAVIN3 ENSG00000170955 -0.000879 \n", "479 604 FYB2 ENSG00000187889 -0.000925 \n", "4609 6032 SHROOM2 ENSG00000146950 -0.001132 \n", "\n", " Goal_end_vs_random_pval N_Detections \\\n", "2979 0.027741 5 \n", "11646 0.175442 6 \n", "7884 0.023335 10 \n", "9497 0.024909 6 \n", "10644 0.003643 7 \n", "... ... ... \n", "53 0.013673 9 \n", "5066 0.055477 7 \n", "5856 0.013125 6 \n", "479 0.005857 11 \n", "4609 0.015731 6 \n", "\n", " Shift_to_alt_end_Invasive cancer + stroma \\\n", "2979 0.000712 \n", "11646 0.001164 \n", "7884 0.000348 \n", "9497 0.000750 \n", "10644 0.001221 \n", "... ... \n", "53 0.000024 \n", "5066 -0.000479 \n", "5856 -0.000491 \n", "479 -0.000171 \n", "4609 -0.000631 \n", "\n", " Alt_end_vs_random_pval_Invasive cancer + stroma \\\n", "2979 0.174332 \n", "11646 0.132797 \n", "7884 0.601956 \n", "9497 0.044214 \n", "10644 0.005439 \n", "... ... \n", "53 0.509069 \n", "5066 0.125197 \n", "5856 0.008692 \n", "479 0.891572 \n", "4609 0.063054 \n", "\n", " Shift_to_alt_end_Invasive cancer + stroma + lymphocytes \\\n", "2979 0.000750 \n", "11646 0.001139 \n", "7884 0.000403 \n", "9497 0.000770 \n", "10644 0.001154 \n", "... ... \n", "53 -0.000080 \n", "5066 -0.000522 \n", "5856 -0.000545 \n", "479 -0.000269 \n", "4609 -0.000729 \n", "\n", " Alt_end_vs_random_pval_Invasive cancer + stroma + lymphocytes \\\n", "2979 0.090140 \n", "11646 0.133500 \n", "7884 0.424766 \n", "9497 0.026993 \n", "10644 0.005040 \n", "... ... \n", "53 0.630417 \n", "5066 0.074521 \n", "5856 0.005454 \n", "479 0.694267 \n", "4609 0.040795 \n", "\n", " Goal_end_FDR Alt_end_FDR_Invasive cancer + stroma \\\n", "2979 0.924671 0.913586 \n", "11646 0.953533 0.886458 \n", "7884 0.924671 0.982160 \n", "9497 0.924671 0.823552 \n", "10644 0.924671 0.734657 \n", "... ... ... \n", "53 0.924671 0.974761 \n", "5066 0.924671 0.886458 \n", "5856 0.924671 0.734657 \n", "479 0.924671 0.996187 \n", "4609 0.924671 0.848050 \n", "\n", " Alt_end_FDR_Invasive cancer + stroma + lymphocytes Sig \n", "2979 0.873315 0 \n", "11646 0.888611 0 \n", "7884 0.959837 0 \n", "9497 0.790271 0 \n", "10644 0.715931 0 \n", "... ... ... \n", "53 0.974686 0 \n", "5066 0.853843 0 \n", "5856 0.715931 0 \n", "479 0.985234 0 \n", "4609 0.828248 0 \n", "\n", "[11664 rows x 14 columns]" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import pandas as pd\n", "cell_shift_stats = pd.read_csv('output/perturb/cell_shift/perturb_stats/visium_spot.csv',index_col=0)\n", "cell_shift_stats" ] }, { "cell_type": "markdown", "id": "1bbee5f1", "metadata": {}, "source": [ "### 1.3 Option C: Overexpress Gene" ] }, { "cell_type": "markdown", "id": "309994e2", "metadata": {}, "source": [ "\n", "1. **Data Filtering & Batching** \n", " - `perturb_dataset()` loads and filters your input dataset. \n", " - In `isp_perturb_set()` or `isp_perturb_all()`, cells are grouped into batches of size `forward_batch_size`.\n", "\n", "2. **Applying Overexpression** \n", " - For each example, `pu.overexpress_tokens(...)` takes the token indices of your gene(s) and **inserts** them at the **front** of the token list. \n", "\n", "3. **Embedding Extraction** \n", " - **Original** (`full_orig`) and **perturbed** (`full_pert`) token sequences are passed through the model (via `get_embs`) to collect hidden‐state tensors. \n", " - Because overexpressed tokens were prepended, the **first N** positions of `full_pert` are your newly-inserted genes; the rest align with the original token order.\n", "\n", "4. **Cosine‐Similarity Quantification** \n", " - **Gene-level**: `pu.quant_cos_sims(pert_emb, original_emb, …, emb_mode=\"gene\")` computes how each remaining gene embedding shifts when your target gene is overexpressed. \n", " - **Cell-level** (if `cell_states_to_model` is set): average the non-padding embeddings before & after overexpression and quantify that shift against your state-embedding targets.\n", "\n", "5. **Aggregation & Output** \n", " - For each cell or state, cosine-similarities are averaged (or bucketed) and stored in a dictionary keyed by the perturbed gene(s). \n", " - Final results are written out as “_cell_embs_dict_…” and—if `emb_mode=\"cell_and_gene\"`—also “_gene_embs_dict_…” files you can analyze with `in_silico_perturber_stats`." ] }, { "cell_type": "code", "execution_count": null, "id": "c9e6a363", "metadata": {}, "outputs": [], "source": [ "from cstformer.perturbation.cstformer_perturb import Perturber\n", "from cstformer.perturbation.perturb_stats import PerturberStats" ] }, { "cell_type": "code", "execution_count": null, "id": "7d148937", "metadata": {}, "outputs": [], "source": [ "isp = Perturber(\n", " perturb_type=\"overexpress\", # simulate gene overexpression\n", " mode = 'spot',\n", " genes_to_perturb=['ENSG00000091831'], #ESR1,\n", " #perturb_rank_shift = None,\n", " model_type=\"Pretrained\", # ML Trained Model\n", " filter_data={'subtype':'TNBC'}, #include filtering to TNBC to see overexpression of ESR1 on these cells\n", " num_classes=0,\n", " emb_mode=\"cell_and_gene\", # cell and gene embeddings \n", " cell_emb_style='mean_pool',\n", " max_ncells=1000,\n", " emb_layer=0,\n", " forward_batch_size=50,\n", " nproc=12,\n", " token_dictionary_file='output/token_dictionary.pickle',\n", " )\n", " \n", "isp.perturb_dataset( \n", " model_directory='output/spot/models/250422_102707_cstformer_L6_E3/final',\n", " input_data_file='output/spot/visium_spot.dataset',\n", " output_directory='output/perturb/overexpress',\n", " output_prefix='perturb_spot')" ] }, { "cell_type": "code", "execution_count": null, "id": "3f712f5a", "metadata": {}, "outputs": [], "source": [ "ispstats = PerturberStats(mode='aggregate_gene_shifts',\n", " token_dictionary_file='output/token_dictionary.pickle',\n", " gene_name_id_dictionary_file='output/ensembl_mapping_dict.pickle'\n", " )\n", "ispstats.compute_stats(\n", " input_data_dir='output/perturb/overexpress',\n", " null_data_dir=None,\n", " output_dir='output/perturb/overexpress/stats',\n", " output_prefix='visium_spot'\n", ")" ] } ], "metadata": { "kernelspec": { "display_name": "stFormer", "language": "python", "name": "stformer" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.12" } }, "nbformat": 4, "nbformat_minor": 5 }