VinceLNeuro/scFM_PerturbationPredictor_PosteriorEst

c2s perturb-seq finetuned model - posterior estimations

0

stars

3

commits

Jupyter Notebook

primary language

Jan 17, 2026

updated

README

C2S

Author: Tianze (Vincent) Luo
Updated: 2025-12-04

Table of Contents


0. preprocessing data

  • ~30k cells (2 donors), ~36k genes from the Domínguez Conde et al. study

  • 5'-capture (TCR/BCR)

  • gene names in adata.var_names --> will be used to create the cell sentences by the C2S code base functions later on

  • should start from raw counts (counts, not continuous (normalized)).


Step1. Data processing + Normalization: C2S only deviates from the standard preprocessing and normalization pipeline in that the log transformation is done with a base of 10 rather than natural logarithm

Step2. PCA > gKNN (nPC=50) > calcUMAP > plotUMAP


1. Cell Sentence Conversion & Reconstruction

1.1. Conversion Workflow

  1. AnnData object: containing our single-cell dataset
  2. Huggingface PyArrow dataset (arrow_ds: cell meta, cell sentence) && vocabulary (vocabulary/features/genes)
    • cs.CSData.adata_to_arrow
  3. CSData object: wraper of arrow dataset; data for inference or finetuning
    • cs.CSData.csdata_from_arrow
    • /ix/ccdg/storage3/til177/ParkLab/Project_C2S_scFM/code/tutorials/dominguez_immune_tissue_tutorial1

arrow_ds

  • 29,773 cells have now been converted into rows of a Dataset object with an additional cell_sentence column
    • The cell sentence contains a sentence of gene names ordered by descending expression level, giving a rank-based gene name representation of the cell.
    • Each row is a dict for a cell
    • e.g.
      {'cell_name': 'Pan_T7935490_AAACCTGCAAATTGCC',
       'cell_sentence': 'RPLP1 ACTB EEF1A1 HSP90AA1 TMSB4X B2M FTH1 KLF6 HSPA1B MALAT1 RPS12 HSPA8 RPL13 MT-CO1 ATF3 MT-CO2 RPL41 TPT1 MT-CO3 ..., ...}
      

vocabulary

  • vocabulary is an OrderedDict of gene features, corresponding to the original 23944 genes in our adata object. The OrderedDict denotes the gene features present in our single-cell dataset, and also stores the number of cells that gene was expressed in.
  • e.g.
    [('RP11-34P13', 38),
     ('RP11-34P13-3', 106),
      ...
    ]
    

1.2. Cell Sentence Conversion Benchmarking

Aim: Know how well the conversion did, and how much expression information was lost when we switched to a rank ordering of genes rather than exact expression values.

  • Paper Fig 10: linear relationship was found between the Log-Rank of a gene and its Log-Norm expression value

  • benchmark_expression_conversion()

    • Fit a linear model on the ranks and expression of the original data, which can be used to reconstruct expression from rank
    • Save plots of (1) log rank vs log expression and (2) log expression vs reconstructed expression from rank
    logNorm vs logRank Reconstructed Expression (from cs) vs Original Expression

1.3. Reconstruct Cell Expression Matrix From Cell Sentences

reconstruct_expression_from_cell_sentence()

  • Need

    • cell_sentences_list (from csdata)
    • vocab_list
    • benchmarking slope & intercept
  • Predict logNorm expression vector && Convert back to Anndata

    • predicted_expression = intercept + (slope * log(rank_of_gene))
  • Very successful inverse transformtion, in terms of rebuilding expression_vector -> rebuilding anndata -> UMAP comparison

umap_comparison.png


2. Cell Embedding with C2S Foundation Models without finetuning

Rationale

  • By loading/defining a CSModel object (pretrained model), the model parameters are completely frozen during embedding extraction. The model simply performs a full forward pass through the pretrained transformer to obtain the last-layer learned latent hidden states , which are then average pooled the latents into a single embedding vector per cell.

    • manuscript BioRxiv L448
    • 1024-d
  • Note: given the pretraining process is using large-scale data across batches/datasets with many tissues, donors, cell types --> the forward pass will show biological structure and implicitly reduce batch effect.

  • By converting cells into learned embeddings, we create a compact representation that captures the essential information from the cell sentences. Cell embeddings are crucial for downstream tasks such as clustering, visualization, and classification

Workflow

  1. Reload preprocessed immune tissue single-cell dataset (preprocessed in tutorial notebook 0, two sample donors) -> create a CSData() wrapper around it

  2. Load a pretrained C2S model (preferably, C2S models which have been trained to do cell type or tissue prediction) to create a CSModel object.

    • Check here for models to use (not include the 2-27B Gemma model): https://github.com/vandijklab/cell2sentence?tab=readme-ov-file#model-zoo
      • vandijklab/C2S-Pythia-410m-diverse-single-and-multi-cell-tasks
      • downloaded @ /ihome/hpark/til177/.cache/huggingface/hub/models--vandijklab--C2S-Pythia-410m-diverse-single-and-multi-cell-tasks/snapshots/51f7c9d46776273ea4732ddaf494d1db733ca5d6/README.md
  3. Embed the cells using the specific C2S model

    • embed_cells()
    • Input: CSData, CSModel,number of genes to use per cell sentence
    • How it worked (details here):
      • Load the C2S model and data
      • Format the cell sentences into prompts for task="cell_type_prediction" (same as a later task)
      • Run the prompts through the C2S model
      • Uses the internal hidden states (via csmodel.embed_cells_batched) as cell embeddings, instead of sampling generated text (normal task="cell_type_prediction")
  4. Visualize the cell embeddings to gain insights into the data

    • Steps:

      • SKIP PCA (given that this latent space is already a kind of nonlinear dimension reduction -- low‑dimensional, model-learned representation)

      • Do gKNN construction -> UMAP

    • Results:

      • Retain distinct clusters separated by cell type and tissue (i.e., cells with similar expression programs and similar cell types should end up close to each other in this refined-learned-representation space) --> all over the place is a warning message
      • Similar cell-type clustering patterns in cell embedding UMAP, when comparing with original UMAP
      • Batch being corrected

      Cell embedding colored by batch Cell embedding colored by tissue

      Original UMAP

      Cell Embedding UMAP


3. [key] Finetuning on a New Single-Cell Dataset

Question: Why do we do cell embedding on a Pre-trained model (given that we are fine-tuning the model, why not just do cell embedding on the fine-tuned model)?

Use Pre-trained model (FM) as "Initialization" / "feature extractor" (better than PCA)

  • M_pre has already learned a good representation of immune cells (basic cell types, marker relationships, etc.).

  • Those embeddings are richer and more biologically meaningful than raw gene counts or PCA, especially if M_pre saw a huge amount of data.

  • So, embeddings from M_pre may be used as a fixed reference space:

    • For QC (outlier detection, batch inspection)
    • For mapping new datasets into the same space later.

Fine‑tuning

  • Next-token prediction objective (as pretraining) with prompts formatted to match each task --> so, even if the cells are from the same dataset, the loss function and goal are different.

  • Code-wise: same underlying loss function for all tasks (Hugging Face causal LM cross‑entropy on tokens).

    • What changes by task is:
      • What you ask the model to do in the prompt (model_input).
      • What you treat as the “ground‑truth” response (response).
      • Whether you take loss on response only vs prompt + response.

Workflow

  1. Load an preprocessed immune tissue single-cell dataset (two sample donors)

  2. (Optional) Custom Prompt Formatting: Format the dataset using a CustomPromptFormatter object, which prepares the data for the fine-tuning process.

  3. (get CSData & CSModel) Load a pretrained C2S model.

  4. Fine-tune the C2S model to improve its performance on cell type prediction.

Output

Training time (GPU): 4296 seconds

  • Model_step3600: {'eval_loss': 1.392207384109497}

  • Model_step3700: {'eval_loss': 1.3921808004379272}

  • *Model_step3725: {'eval_loss': 1.3921139240264893}

    see ../csmodel_tutorial_3/2025-11-24-18_42_05_finetune_cell_type_prediction/checkpoint-3725/loss_curves.png

4. Cell Type Prediction/Annotation (Using Fine-Tuned Model)

Workflow

  1. Load an preprocessed immune tissue single-cell dataset (two sample donors)

  2. Load fine-tuned CSModel (training_task = "cell_type_prediction") CHECKPOINT

    • Load the data_split_indices_dict.pkl associated with the fine-tuned model
      • Containing indices
      • train/val/test split = 80/10/10
  3. C2S conversion (AnnData [full] -> Arrow [full] -> CSData [test only])

  4. INFERENCE: Predict Cell Types

    • Input: finetuned cell type prediction model CSModel and have our test set CSData
    • predict_cell_types_of_data()
  5. Calculate accuracy

    • accuracy = 0.81

7. [key] Custom prompt templates (for PromptEngineering)

  1. Load an preprocessed immune tissue single-cell dataset (two sample donors)

  2. AnnData -> Arrow

  3. Custom Prompt Formatting

    • Create a subclass of the PromptFormatter class, which has to define a format_hf_ds method that takes in a cell sentence arrow dataset and returns a formatted dataset

    • prompt_formatter = CustomPromptFormatter()

    • It can be beneficial to provide variations of prompt templates to provide some diversity - simply create several templates and choose one when formatting each sample in the formatting function!

    • Note: csmodel.fine_tune(prompt_formatter=prompt_formatter) function will do the formatting on the full dataset for us --> no need to do formatting here

      # Example Arrow
      Dataset({
          features: ['cell_name', 'cell_sentence', 'cell_type', 'tissue', 'batch_condition', 'organism', 'sex'],
          num_rows: 10
      })
      
      # Example Formatted Arrow
      Dataset({
          features: ['sample_type', 'model_input', 'response'],
          num_rows: 10
      })
      
  4. Arrow -> CSData

  5. Load C2S FM (https://huggingface.co/collections/vandijklab/cell2sentence-models)

  6. Fine-tune on new task!


10. [key] Finetuning for Perturbation Response Prediction

Rationale

Perturbation Response:

how a cell's gene expression profile changes in response to a specific perturbation (e.g., a genetic knockout or a drug treatment)

We will treat this as a "translation" task in natural language: translating a cell (in cell sentence format) from its basal (control) state to its perturbed state, conditioned on the perturbation applied.

Workflow

At a high level, we will:

  1. Load a public single-cell perturbation dataset.

    • Data requirement

      • AnnData object
      • .obs dataframe must contain:
        • A column that distinguishes control cells from perturbed cells, e.g., adata.obs['condition']
    • Data used in this analysis

      • Original Paper: https://www-nature-com.pitt.idm.oclc.org/articles/s41588-025-02169-3#Sec10

        • "A second CRISPRi Jurkat cell line expressing the optimized UCOE-EF1α-Zim3-dCas9-P2A-mCherry CRISPRi construct was generated as previously described and was used for Perturb-seq."

          • Same cell line
        • "For the Jurkat Perturb-seq experiment, Jurkat cells expressing Zim3-dCas9-P2A-mCherry were transduced with dJR092 library lentivirus by spinfection (1,000g) with polybrene (8 µg ml−1; Sigma-Aldrich) with a targeted low infection rate of ~10%. This low rate was chosen to reduce the chances of a single cell being infected by several viruses."

          • Controlled the doses intentionally
      • Perturb-seq data/experiments in Jurkat cells: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE264667

        • Different cells get different sgRNAs (targeting different genes or control sgRNAs) --> Pooled all the edited cells together to do scRNA-seq

        • The downloaded data (/ix/ccdg/storage3/til177/ParkLab/Project_C2S_scFM/code/tutorials/data_PerturbSeq/GSE264667_jurkat.h5ad) already passed filtering used in original paper

          • Filtering in original paper -> 262,956 cells retained
            1. min_umi >1,750
            2. max_mito = 14
    • Our basic filtering did not remove cells/features; main concern is sequencing depth -> used advanced filtering (obs: 262956 -> 257412)

      # For Jurkat Perturb-seq with median UMI ~10k
      min_umi = 1000    # ~10% of median
      max_umi = 40000   # ~4x median (doublet filter)
      min_genes = 500
      max_genes = 6000
      max_mito = 15
      
      # Apply filters
      adata = adata[
          (adata.obs['UMI_count'] > min_umi) &
          (adata.obs['UMI_count'] < max_umi) &
          (adata.obs['n_genes'] > min_genes) &
          (adata.obs['n_genes'] < max_genes) &
          (adata.obs['mitopercent'] < max_mito)
      ].copy()
      
      • [Added 12/11/2025] Filtering functioning here are:
        • max_umi = 40000
        • max_genes = 6000
    • Normalization

  2. Write a custom prompt template for perturbation prediction.

    • Subclass the PromptFormatter class (ABC) to create pairs of control and perturbed cells.
      • format_hf_ds method (will be applied automatically in csmodel.fine_tune())
      • Output: formatted HF Dataset
      • Note: Using top 200 genes for this example. For real applications, ideal to use all nonzero expressed genes if possible.
  3. Load and Finetune a pretrained C2S-Scale model on this new task.

    • Note: For this tutorial, we'll run for a small number of steps (max_steps=500). For a full finetuning run, you would typically train for several epochs.

    • loss_on_response_only=True: use input 'control cell' & 'purturbation gene name' as condition (p( perturbed_cell | control_cell, perturbation ))

      • We only want to compute loss on the predicted perturbed cell sentence
      • Do not need to waste resources on the control cell-sentence prediction
  4. [c2s_tvl_11] Generate a prediction with our new finetuned model to see it in action.


11. [key] Posterior Over Responses Instead of Point Estimation

See c2s_tvl_11_perturbation_ftEval_PosteriorEst.ipynb for details

Created two helper functions for empirical statistics calculation

  • sample_perturbation_posterior_sentences()
  • posterior_sentences_to_expression()

11v2/v3. [key] Sampling 100 times for 1 inference prompt (comparing with two temperatures)

Data stored at /ix/ccdg/storage3/til177/ParkLab/Project_C2S_scFM/code/tutorials/data_PerturbSeq/perturbation_predictor_finetuned_final_benchmarking/posterior_samples_meta_script11v2.pkl

  • (see 'c2s_tvl_11v2_perturbation_PosteriorEst.py')

Caveat

If using top_k_genes=200 but not the whole list of genes (e.g., in v2 PerturbationPromptFormatter), then reconstruct_expression_from_cell_sentence() will create zeros for those genes NOT among the top_k_genes.

Use top_k_genes=200 model for visualization to quality check if temperature difference during generation will create a different distribution of the samples (i.e., higher temperature would have a more spread distribution on the shared UMAP)

see perturbation_predictor_finetuned_final_benchmarking/PCA__top_k_genes200_PosteriorSamplesByTemp.png see perturbation_predictor_finetuned_final_benchmarking/umap__top_k_genes200_PosteriorSamplesByTemp.png

In v3, we switched to (full-length) top 2048 genes (Pythia-1B with 8192-token context limit -> 2048 genes)

  • best_model_checkpoint: /ix/ccdg/storage3/til177/ParkLab/Project_C2S_scFM/code/tutorials/data_PerturbSeq/finetunedModel_2025-12-26-04_18_42_FullLengthFinetune_perturbation_prediction/checkpoint-36500

    • No overfitting
    • Seems that the loss curve is too noisy --> might need larger eff. batch size
    • The loss curve decreases slower after 15k steps --> may consider early stop see loss_curves_FullLength.png
  • Problematic in the generation step

    • The top_k_genes=2048 model did not exceed the 8192 context length. However, the generation length is limited given the context length cap.

    • prompt_len_tokens: 7015 model_max: 8192 max_new_tokens: 1113

      • The current generation script is beyond the model_max (total=7015+8192) --> The posterior_samples_meta_script11v3Run3.pkl is low quality from the PCA/UMAP (expected degraded beyond model limit)

      • Truncation (to within the model limit) improves the result but still invalid... -> should regenerate if want to use the top_k_genes=2048 model

13. Head & Tail Gene Generation QC (top-200-genes model)

See code/tutorials/c2s_tvl_pipeline/c2s_tvl_13_diffPerturbGene_QC.ipynb

Check all unique perturbations (genes) in the test set

  • In total 2339 unique perturbations/genes in the test set (which makes sense when there are many duplicates before train-eval-test split)

Check if 'head' and 'tail' genes (HJ informed) are in the test samples && OOD-OOV issues

HJ: Check the genes in "well-studied “head” genes (e.g., TP53, EGFR, MKI67, ribosomal RPL/RPS; interferon ISG15/IFI6) acquire rich, stable embeddings, while rare or lineage-restricted “tail” genes (e.g., tuft-cell POU2F3, mTEC AIRE, hair-cell ATOH1".

  • A lot of ribosomal RPL/RPS genes in the test samples (prompts)

  • If we directly swap the 'model_input' in the formatted_test_ds as 'TP53' ==> this is OOD/zero-shot perturbation label, leading to weaker/inconsistent conditioning on the perturbation label, or outputs resembling an “average perturbation” or nearest seen perturbations.

  • If we manually use OOV (wrt this HF dataset) genes (e.g., provided 'tail' genes), the model will still accept the prompt (tokenizer can encode it).

    • But it’s OOD wrt fine-tuning (the model likely didn’t see that token in gene lists).
    • It can confuse generation slightly.
    • And importantly: the reconstruction/post-processing process will treat it as NOT_A_GENE / removed, so it won’t be represented in expression space.

*Newly defined 'head' and 'tail' genes

Definition

  • 'Head' genes are those have broad and essential biological and cellular functions and has a ubiquitous expression pattens (such as genes/proteins important to transcription and translation).
  • 'Tail' genes are lineage/tissue restricted, usually developmentally related TFs (I double checked using GTEx)
  • All these genes were selected from the test set perturbations/genes

Potential 'head' genes

  • TARDBP (ALS-related, broadly studied RNA-binding protein/gene vital for RNA processing, splicing, stability, and gene regulation)
  • EIF4B (translation initiation factor)
  • RPL/RPS pair
    • RPL13A('housekeeping gene')
    • RPS13('stable reference gene')

Potential 'tail' genes:

  • FOXL2 (TF in ovary - granulosa cells)
  • GATA1 (crucial TF, essentially the master regulator, of erythroid (red blood cell), megakaryocytic (platelet precursor) development)
  • KRT10 (GTEx; HPA)(later added for its presence in the perturb-seq data): skin keratinocyte (very non–T-cell-like).
  • HMX3 (GTEx; HPA)(later added for its presence in the perturb-seq data): Developmental TF in Tuft cells.

In total, 8 genes for generation

=== Generation Step ===

  • See directory: code/tutorials/c2s_tvl_pipeline/c2s_tvl_13v2_diffPerturbGenes_Gen_QC for python script and slurm array submission code

    • Basically, we locate the sample indices for those eight genes in formatted_test_ds (by searching across inference prompts -> many replications of single unique perturbs ->) and further random.choice one inference prompt for each target gene.

    • The index dictionary is saved at code/tutorials/data_PerturbSeq/finetunedModel_2025-12-09-16_44_46_testFinetune_perturbation_prediction/idx__formatted_test_ds__selected_head-tail-genes_SEED1234.pkl with README inside.

      defaultdict(<class 'list'>, {'TARDBP': 2660, 'EIF4B': 455, 'RPL13A': 19826, 'RPS13': 13270, 'FOXL2': 23190, 'GATA1': 21352, 'KRT10': 18395, 'HMX3': 15316})
      
  • See directory code/tutorials/data_PerturbSeq/perturbation_predictor_finetuned_final_benchmarking/c2s_tvl_13v2_diffPerturbGenes_Gen for stored generation results

=== Plot Step ===


Other Optional Tasks

5. Cell Generation

Workflow

  • Similar pipeline
  • Need a model finetuned to do cell generation (from tutorial 3)
  • generate_cells_conditioned_on_cell_type()
  • Post-processing and Reconstruction
    • Remove words that are not gene names in the dataset (vocabulary)
    • Deal with duplicated genes
    • post_processed_sentence, num_genes_replaced = post_process_generated_cell_sentences()
    • Reconstruction (see tutorial 1)
  • Visualization
    • Compare the UMAP seperately
    • Compare heatmap
    • Plot in shared space --> should be very similar in the shared embedding space

6. Cell Type Annotation with C2S Foundation Model

Workflow very similar to tutorial 4

8 & 9 [not yet implemented] Multi-cell

Contributors

VinceLNeuro

3 commits

VinceLNeuro/scFM_PerturbationPredictor_PosteriorEst

c2s perturb-seq finetuned model - posterior estimations

0

stars

3

commits

Jupyter Notebook

primary language

Jan 17, 2026

updated

README

C2S

Author: Tianze (Vincent) Luo
Updated: 2025-12-04

Table of Contents


0. preprocessing data

  • ~30k cells (2 donors), ~36k genes from the Domínguez Conde et al. study

  • 5'-capture (TCR/BCR)

  • gene names in adata.var_names --> will be used to create the cell sentences by the C2S code base functions later on

  • should start from raw counts (counts, not continuous (normalized)).


Step1. Data processing + Normalization: C2S only deviates from the standard preprocessing and normalization pipeline in that the log transformation is done with a base of 10 rather than natural logarithm

Step2. PCA > gKNN (nPC=50) > calcUMAP > plotUMAP


1. Cell Sentence Conversion & Reconstruction

1.1. Conversion Workflow

  1. AnnData object: containing our single-cell dataset
  2. Huggingface PyArrow dataset (arrow_ds: cell meta, cell sentence) && vocabulary (vocabulary/features/genes)
    • cs.CSData.adata_to_arrow
  3. CSData object: wraper of arrow dataset; data for inference or finetuning
    • cs.CSData.csdata_from_arrow
    • /ix/ccdg/storage3/til177/ParkLab/Project_C2S_scFM/code/tutorials/dominguez_immune_tissue_tutorial1

arrow_ds

  • 29,773 cells have now been converted into rows of a Dataset object with an additional cell_sentence column
    • The cell sentence contains a sentence of gene names ordered by descending expression level, giving a rank-based gene name representation of the cell.
    • Each row is a dict for a cell
    • e.g.
      {'cell_name': 'Pan_T7935490_AAACCTGCAAATTGCC',
       'cell_sentence': 'RPLP1 ACTB EEF1A1 HSP90AA1 TMSB4X B2M FTH1 KLF6 HSPA1B MALAT1 RPS12 HSPA8 RPL13 MT-CO1 ATF3 MT-CO2 RPL41 TPT1 MT-CO3 ..., ...}
      

vocabulary

  • vocabulary is an OrderedDict of gene features, corresponding to the original 23944 genes in our adata object. The OrderedDict denotes the gene features present in our single-cell dataset, and also stores the number of cells that gene was expressed in.
  • e.g.
    [('RP11-34P13', 38),
     ('RP11-34P13-3', 106),
      ...
    ]
    

1.2. Cell Sentence Conversion Benchmarking

Aim: Know how well the conversion did, and how much expression information was lost when we switched to a rank ordering of genes rather than exact expression values.

  • Paper Fig 10: linear relationship was found between the Log-Rank of a gene and its Log-Norm expression value

  • benchmark_expression_conversion()

    • Fit a linear model on the ranks and expression of the original data, which can be used to reconstruct expression from rank
    • Save plots of (1) log rank vs log expression and (2) log expression vs reconstructed expression from rank
    logNorm vs logRank Reconstructed Expression (from cs) vs Original Expression

1.3. Reconstruct Cell Expression Matrix From Cell Sentences

reconstruct_expression_from_cell_sentence()

  • Need

    • cell_sentences_list (from csdata)
    • vocab_list
    • benchmarking slope & intercept
  • Predict logNorm expression vector && Convert back to Anndata

    • predicted_expression = intercept + (slope * log(rank_of_gene))
  • Very successful inverse transformtion, in terms of rebuilding expression_vector -> rebuilding anndata -> UMAP comparison

umap_comparison.png


2. Cell Embedding with C2S Foundation Models without finetuning

Rationale

  • By loading/defining a CSModel object (pretrained model), the model parameters are completely frozen during embedding extraction. The model simply performs a full forward pass through the pretrained transformer to obtain the last-layer learned latent hidden states , which are then average pooled the latents into a single embedding vector per cell.

    • manuscript BioRxiv L448
    • 1024-d
  • Note: given the pretraining process is using large-scale data across batches/datasets with many tissues, donors, cell types --> the forward pass will show biological structure and implicitly reduce batch effect.

  • By converting cells into learned embeddings, we create a compact representation that captures the essential information from the cell sentences. Cell embeddings are crucial for downstream tasks such as clustering, visualization, and classification

Workflow

  1. Reload preprocessed immune tissue single-cell dataset (preprocessed in tutorial notebook 0, two sample donors) -> create a CSData() wrapper around it

  2. Load a pretrained C2S model (preferably, C2S models which have been trained to do cell type or tissue prediction) to create a CSModel object.

    • Check here for models to use (not include the 2-27B Gemma model): https://github.com/vandijklab/cell2sentence?tab=readme-ov-file#model-zoo
      • vandijklab/C2S-Pythia-410m-diverse-single-and-multi-cell-tasks
      • downloaded @ /ihome/hpark/til177/.cache/huggingface/hub/models--vandijklab--C2S-Pythia-410m-diverse-single-and-multi-cell-tasks/snapshots/51f7c9d46776273ea4732ddaf494d1db733ca5d6/README.md
  3. Embed the cells using the specific C2S model

    • embed_cells()
    • Input: CSData, CSModel,number of genes to use per cell sentence
    • How it worked (details here):
      • Load the C2S model and data
      • Format the cell sentences into prompts for task="cell_type_prediction" (same as a later task)
      • Run the prompts through the C2S model
      • Uses the internal hidden states (via csmodel.embed_cells_batched) as cell embeddings, instead of sampling generated text (normal task="cell_type_prediction")
  4. Visualize the cell embeddings to gain insights into the data

    • Steps:

      • SKIP PCA (given that this latent space is already a kind of nonlinear dimension reduction -- low‑dimensional, model-learned representation)

      • Do gKNN construction -> UMAP

    • Results:

      • Retain distinct clusters separated by cell type and tissue (i.e., cells with similar expression programs and similar cell types should end up close to each other in this refined-learned-representation space) --> all over the place is a warning message
      • Similar cell-type clustering patterns in cell embedding UMAP, when comparing with original UMAP
      • Batch being corrected

      Cell embedding colored by batch Cell embedding colored by tissue

      Original UMAP

      Cell Embedding UMAP


3. [key] Finetuning on a New Single-Cell Dataset

Question: Why do we do cell embedding on a Pre-trained model (given that we are fine-tuning the model, why not just do cell embedding on the fine-tuned model)?

Use Pre-trained model (FM) as "Initialization" / "feature extractor" (better than PCA)

  • M_pre has already learned a good representation of immune cells (basic cell types, marker relationships, etc.).

  • Those embeddings are richer and more biologically meaningful than raw gene counts or PCA, especially if M_pre saw a huge amount of data.

  • So, embeddings from M_pre may be used as a fixed reference space:

    • For QC (outlier detection, batch inspection)
    • For mapping new datasets into the same space later.

Fine‑tuning

  • Next-token prediction objective (as pretraining) with prompts formatted to match each task --> so, even if the cells are from the same dataset, the loss function and goal are different.

  • Code-wise: same underlying loss function for all tasks (Hugging Face causal LM cross‑entropy on tokens).

    • What changes by task is:
      • What you ask the model to do in the prompt (model_input).
      • What you treat as the “ground‑truth” response (response).
      • Whether you take loss on response only vs prompt + response.

Workflow

  1. Load an preprocessed immune tissue single-cell dataset (two sample donors)

  2. (Optional) Custom Prompt Formatting: Format the dataset using a CustomPromptFormatter object, which prepares the data for the fine-tuning process.

  3. (get CSData & CSModel) Load a pretrained C2S model.

  4. Fine-tune the C2S model to improve its performance on cell type prediction.

Output

Training time (GPU): 4296 seconds

  • Model_step3600: {'eval_loss': 1.392207384109497}

  • Model_step3700: {'eval_loss': 1.3921808004379272}

  • *Model_step3725: {'eval_loss': 1.3921139240264893}

    see ../csmodel_tutorial_3/2025-11-24-18_42_05_finetune_cell_type_prediction/checkpoint-3725/loss_curves.png

4. Cell Type Prediction/Annotation (Using Fine-Tuned Model)

Workflow

  1. Load an preprocessed immune tissue single-cell dataset (two sample donors)

  2. Load fine-tuned CSModel (training_task = "cell_type_prediction") CHECKPOINT

    • Load the data_split_indices_dict.pkl associated with the fine-tuned model
      • Containing indices
      • train/val/test split = 80/10/10
  3. C2S conversion (AnnData [full] -> Arrow [full] -> CSData [test only])

  4. INFERENCE: Predict Cell Types

    • Input: finetuned cell type prediction model CSModel and have our test set CSData
    • predict_cell_types_of_data()
  5. Calculate accuracy

    • accuracy = 0.81

7. [key] Custom prompt templates (for PromptEngineering)

  1. Load an preprocessed immune tissue single-cell dataset (two sample donors)

  2. AnnData -> Arrow

  3. Custom Prompt Formatting

    • Create a subclass of the PromptFormatter class, which has to define a format_hf_ds method that takes in a cell sentence arrow dataset and returns a formatted dataset

    • prompt_formatter = CustomPromptFormatter()

    • It can be beneficial to provide variations of prompt templates to provide some diversity - simply create several templates and choose one when formatting each sample in the formatting function!

    • Note: csmodel.fine_tune(prompt_formatter=prompt_formatter) function will do the formatting on the full dataset for us --> no need to do formatting here

      # Example Arrow
      Dataset({
          features: ['cell_name', 'cell_sentence', 'cell_type', 'tissue', 'batch_condition', 'organism', 'sex'],
          num_rows: 10
      })
      
      # Example Formatted Arrow
      Dataset({
          features: ['sample_type', 'model_input', 'response'],
          num_rows: 10
      })
      
  4. Arrow -> CSData

  5. Load C2S FM (https://huggingface.co/collections/vandijklab/cell2sentence-models)

  6. Fine-tune on new task!


10. [key] Finetuning for Perturbation Response Prediction

Rationale

Perturbation Response:

how a cell's gene expression profile changes in response to a specific perturbation (e.g., a genetic knockout or a drug treatment)

We will treat this as a "translation" task in natural language: translating a cell (in cell sentence format) from its basal (control) state to its perturbed state, conditioned on the perturbation applied.

Workflow

At a high level, we will:

  1. Load a public single-cell perturbation dataset.

    • Data requirement

      • AnnData object
      • .obs dataframe must contain:
        • A column that distinguishes control cells from perturbed cells, e.g., adata.obs['condition']
    • Data used in this analysis

      • Original Paper: https://www-nature-com.pitt.idm.oclc.org/articles/s41588-025-02169-3#Sec10

        • "A second CRISPRi Jurkat cell line expressing the optimized UCOE-EF1α-Zim3-dCas9-P2A-mCherry CRISPRi construct was generated as previously described and was used for Perturb-seq."

          • Same cell line
        • "For the Jurkat Perturb-seq experiment, Jurkat cells expressing Zim3-dCas9-P2A-mCherry were transduced with dJR092 library lentivirus by spinfection (1,000g) with polybrene (8 µg ml−1; Sigma-Aldrich) with a targeted low infection rate of ~10%. This low rate was chosen to reduce the chances of a single cell being infected by several viruses."

          • Controlled the doses intentionally
      • Perturb-seq data/experiments in Jurkat cells: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE264667

        • Different cells get different sgRNAs (targeting different genes or control sgRNAs) --> Pooled all the edited cells together to do scRNA-seq

        • The downloaded data (/ix/ccdg/storage3/til177/ParkLab/Project_C2S_scFM/code/tutorials/data_PerturbSeq/GSE264667_jurkat.h5ad) already passed filtering used in original paper

          • Filtering in original paper -> 262,956 cells retained
            1. min_umi >1,750
            2. max_mito = 14
    • Our basic filtering did not remove cells/features; main concern is sequencing depth -> used advanced filtering (obs: 262956 -> 257412)

      # For Jurkat Perturb-seq with median UMI ~10k
      min_umi = 1000    # ~10% of median
      max_umi = 40000   # ~4x median (doublet filter)
      min_genes = 500
      max_genes = 6000
      max_mito = 15
      
      # Apply filters
      adata = adata[
          (adata.obs['UMI_count'] > min_umi) &
          (adata.obs['UMI_count'] < max_umi) &
          (adata.obs['n_genes'] > min_genes) &
          (adata.obs['n_genes'] < max_genes) &
          (adata.obs['mitopercent'] < max_mito)
      ].copy()
      
      • [Added 12/11/2025] Filtering functioning here are:
        • max_umi = 40000
        • max_genes = 6000
    • Normalization

  2. Write a custom prompt template for perturbation prediction.

    • Subclass the PromptFormatter class (ABC) to create pairs of control and perturbed cells.
      • format_hf_ds method (will be applied automatically in csmodel.fine_tune())
      • Output: formatted HF Dataset
      • Note: Using top 200 genes for this example. For real applications, ideal to use all nonzero expressed genes if possible.
  3. Load and Finetune a pretrained C2S-Scale model on this new task.

    • Note: For this tutorial, we'll run for a small number of steps (max_steps=500). For a full finetuning run, you would typically train for several epochs.

    • loss_on_response_only=True: use input 'control cell' & 'purturbation gene name' as condition (p( perturbed_cell | control_cell, perturbation ))

      • We only want to compute loss on the predicted perturbed cell sentence
      • Do not need to waste resources on the control cell-sentence prediction
  4. [c2s_tvl_11] Generate a prediction with our new finetuned model to see it in action.


11. [key] Posterior Over Responses Instead of Point Estimation

See c2s_tvl_11_perturbation_ftEval_PosteriorEst.ipynb for details

Created two helper functions for empirical statistics calculation

  • sample_perturbation_posterior_sentences()
  • posterior_sentences_to_expression()

11v2/v3. [key] Sampling 100 times for 1 inference prompt (comparing with two temperatures)

Data stored at /ix/ccdg/storage3/til177/ParkLab/Project_C2S_scFM/code/tutorials/data_PerturbSeq/perturbation_predictor_finetuned_final_benchmarking/posterior_samples_meta_script11v2.pkl

  • (see 'c2s_tvl_11v2_perturbation_PosteriorEst.py')

Caveat

If using top_k_genes=200 but not the whole list of genes (e.g., in v2 PerturbationPromptFormatter), then reconstruct_expression_from_cell_sentence() will create zeros for those genes NOT among the top_k_genes.

Use top_k_genes=200 model for visualization to quality check if temperature difference during generation will create a different distribution of the samples (i.e., higher temperature would have a more spread distribution on the shared UMAP)

see perturbation_predictor_finetuned_final_benchmarking/PCA__top_k_genes200_PosteriorSamplesByTemp.png see perturbation_predictor_finetuned_final_benchmarking/umap__top_k_genes200_PosteriorSamplesByTemp.png

In v3, we switched to (full-length) top 2048 genes (Pythia-1B with 8192-token context limit -> 2048 genes)

  • best_model_checkpoint: /ix/ccdg/storage3/til177/ParkLab/Project_C2S_scFM/code/tutorials/data_PerturbSeq/finetunedModel_2025-12-26-04_18_42_FullLengthFinetune_perturbation_prediction/checkpoint-36500

    • No overfitting
    • Seems that the loss curve is too noisy --> might need larger eff. batch size
    • The loss curve decreases slower after 15k steps --> may consider early stop see loss_curves_FullLength.png
  • Problematic in the generation step

    • The top_k_genes=2048 model did not exceed the 8192 context length. However, the generation length is limited given the context length cap.

    • prompt_len_tokens: 7015 model_max: 8192 max_new_tokens: 1113

      • The current generation script is beyond the model_max (total=7015+8192) --> The posterior_samples_meta_script11v3Run3.pkl is low quality from the PCA/UMAP (expected degraded beyond model limit)

      • Truncation (to within the model limit) improves the result but still invalid... -> should regenerate if want to use the top_k_genes=2048 model

13. Head & Tail Gene Generation QC (top-200-genes model)

See code/tutorials/c2s_tvl_pipeline/c2s_tvl_13_diffPerturbGene_QC.ipynb

Check all unique perturbations (genes) in the test set

  • In total 2339 unique perturbations/genes in the test set (which makes sense when there are many duplicates before train-eval-test split)

Check if 'head' and 'tail' genes (HJ informed) are in the test samples && OOD-OOV issues

HJ: Check the genes in "well-studied “head” genes (e.g., TP53, EGFR, MKI67, ribosomal RPL/RPS; interferon ISG15/IFI6) acquire rich, stable embeddings, while rare or lineage-restricted “tail” genes (e.g., tuft-cell POU2F3, mTEC AIRE, hair-cell ATOH1".

  • A lot of ribosomal RPL/RPS genes in the test samples (prompts)

  • If we directly swap the 'model_input' in the formatted_test_ds as 'TP53' ==> this is OOD/zero-shot perturbation label, leading to weaker/inconsistent conditioning on the perturbation label, or outputs resembling an “average perturbation” or nearest seen perturbations.

  • If we manually use OOV (wrt this HF dataset) genes (e.g., provided 'tail' genes), the model will still accept the prompt (tokenizer can encode it).

    • But it’s OOD wrt fine-tuning (the model likely didn’t see that token in gene lists).
    • It can confuse generation slightly.
    • And importantly: the reconstruction/post-processing process will treat it as NOT_A_GENE / removed, so it won’t be represented in expression space.

*Newly defined 'head' and 'tail' genes

Definition

  • 'Head' genes are those have broad and essential biological and cellular functions and has a ubiquitous expression pattens (such as genes/proteins important to transcription and translation).
  • 'Tail' genes are lineage/tissue restricted, usually developmentally related TFs (I double checked using GTEx)
  • All these genes were selected from the test set perturbations/genes

Potential 'head' genes

  • TARDBP (ALS-related, broadly studied RNA-binding protein/gene vital for RNA processing, splicing, stability, and gene regulation)
  • EIF4B (translation initiation factor)
  • RPL/RPS pair
    • RPL13A('housekeeping gene')
    • RPS13('stable reference gene')

Potential 'tail' genes:

  • FOXL2 (TF in ovary - granulosa cells)
  • GATA1 (crucial TF, essentially the master regulator, of erythroid (red blood cell), megakaryocytic (platelet precursor) development)
  • KRT10 (GTEx; HPA)(later added for its presence in the perturb-seq data): skin keratinocyte (very non–T-cell-like).
  • HMX3 (GTEx; HPA)(later added for its presence in the perturb-seq data): Developmental TF in Tuft cells.

In total, 8 genes for generation

=== Generation Step ===

  • See directory: code/tutorials/c2s_tvl_pipeline/c2s_tvl_13v2_diffPerturbGenes_Gen_QC for python script and slurm array submission code

    • Basically, we locate the sample indices for those eight genes in formatted_test_ds (by searching across inference prompts -> many replications of single unique perturbs ->) and further random.choice one inference prompt for each target gene.

    • The index dictionary is saved at code/tutorials/data_PerturbSeq/finetunedModel_2025-12-09-16_44_46_testFinetune_perturbation_prediction/idx__formatted_test_ds__selected_head-tail-genes_SEED1234.pkl with README inside.

      defaultdict(<class 'list'>, {'TARDBP': 2660, 'EIF4B': 455, 'RPL13A': 19826, 'RPS13': 13270, 'FOXL2': 23190, 'GATA1': 21352, 'KRT10': 18395, 'HMX3': 15316})
      
  • See directory code/tutorials/data_PerturbSeq/perturbation_predictor_finetuned_final_benchmarking/c2s_tvl_13v2_diffPerturbGenes_Gen for stored generation results

=== Plot Step ===


Other Optional Tasks

5. Cell Generation

Workflow

  • Similar pipeline
  • Need a model finetuned to do cell generation (from tutorial 3)
  • generate_cells_conditioned_on_cell_type()
  • Post-processing and Reconstruction
    • Remove words that are not gene names in the dataset (vocabulary)
    • Deal with duplicated genes
    • post_processed_sentence, num_genes_replaced = post_process_generated_cell_sentences()
    • Reconstruction (see tutorial 1)
  • Visualization
    • Compare the UMAP seperately
    • Compare heatmap
    • Plot in shared space --> should be very similar in the shared embedding space

6. Cell Type Annotation with C2S Foundation Model

Workflow very similar to tutorial 4

8 & 9 [not yet implemented] Multi-cell

Contributors

VinceLNeuro

3 commits

Languages

Jupyter Notebook

63.4%

HTML

36.5%