scEEMS Prediction - FunGen-xQTL Consortium

Under active development. This module is a work in progress: its interface, parameters and outputs may still change, and it is not yet covered by the automated test suite.

Description¶

Pipeline Execution Workflow¶

Step 1: Input Data Preparation¶

Create a tab-separated file with variant identifiers:

variant_id
2:12345:A:T
2:67890:G:C
2:11111:T:A

Format requirements:

Chromosome notation: Numeric (e.g., 2) or with prefix (e.g., chr2)
Position: 1-based genomic coordinate (GRCh38/hg38 reference)
Alleles: Reference and alternate alleles (ACGT notation)

Step 2: Execute Prediction Pipeline¶

cd ~/xqtl-protocol/code/SoS/xqtl_modifier_score/
python model_training_model5_only.py Mic_mega_eQTL 2 \
  --data-config data_config.yaml \
  --model-config model_config.yaml

What happens during execution:

Variant parsing and coordinate validation
Feature annotation (distance, regulatory, population genetics, conservation, deep learning predictions)
Gene constraint and MAF integration with imputation using training statistics
Feature subsetting and absolute value transformations
Model inference with 10x weighting for regulatory features
Probability calibration and output generation

Step 3: Output Files Generated¶

File	Description	Use Case
`model_standard_subset_weighted_chr_chr2_NPR_10.cbm`	Serialized CatBoost model	Future predictions, model inspection
`predictions_weighted_model_chr2.tsv`	Per-variant predictions with scores	Primary analysis, variant prioritization
`summary_dict_catboost_weighted_model_chr_chr2_NPR_10.pkl`	AP/AUC metrics, class distributions	Model validation, quality control
`features_importance_model5_chr_chr2_NPR_10.csv`	Feature importance rankings	Biological interpretation

Load and Inspect Results¶

# Inspect the EMS training predictions (functional-score table written by [train]).
R:
    suppressPackageStartupMessages(library(vroom))
    results_dir <- "output/xqtl_modifier_score/protocol_example/model_results"
    results <- vroom(file.path(results_dir, "predictions_parquet_catboost",
                               "predictions_weighted_model_chr2.tsv"), show_col_types = FALSE)
    cat("Loaded predictions for", nrow(results), "variants\n")
    cat("Columns: variant_id, standard_subset_weighted_pred_prob (EMS score 0-1),\n")
    cat("         standard_subset_weighted_pred_label (0/1), actual_label\n")

Example output:

Model Configuration:
  Algorithm: CatBoostClassifier
  Features: 4839
  Classes: [0, 1]
  Tree depth: 6
  Iterations: 1000

Loaded predictions for 761 variants

Available columns:
  - variant_id: Genomic coordinates
  - standard_subset_weighted_pred_prob: EMS functional score (0-1)
  - standard_subset_weighted_pred_label: Binary classification (0/1)
  - actual_label: True label from test set (if available)

Statistical Summary of Predictions¶

The standard_subset_weighted_pred_prob column contains continuous EMS scores representing the probability that each variant has functional regulatory impact on gene expression.

R:
    suppressPackageStartupMessages(library(vroom))
    results <- vroom("tests/fixtures/ems_training/expected/predictions_weighted_model_chr2.tsv", show_col_types = FALSE)
    s <- results$standard_subset_weighted_pred_prob
    cat("SCORE DISTRIBUTION ANALYSIS\n"); cat(strrep("=", 50), "\n", sep = "")
    cat("Total variants analyzed:", length(s), "\n\n")
    cat("Priority Classification:\n")
    cat(sprintf("  High priority (>0.8):      %5d (%5.1f%%)\n", sum(s > 0.8), 100*mean(s > 0.8)))
    cat(sprintf("  Medium priority (0.5-0.8): %5d (%5.1f%%)\n", sum(s > 0.5 & s <= 0.8), 100*mean(s > 0.5 & s <= 0.8)))
    cat(sprintf("  Low priority (<=0.5):      %5d (%5.1f%%)\n", sum(s <= 0.5), 100*mean(s <= 0.5)))
    cat(sprintf("\nScore Statistics: mean=%.4f median=%.4f sd=%.4f min=%.4f max=%.4f\n",
                mean(s), median(s), sd(s), min(s), max(s)))
    for (q in c(90, 95, 99)) cat(sprintf("  %dth percentile: %.4f\n", q, quantile(s, q/100)))

Example output:

SCORE DISTRIBUTION ANALYSIS
==================================================
Total variants analyzed: 761

Priority Classification:
  High priority (>0.8):        12 (  1.6%)
  Medium priority (0.5-0.8):   45 (  5.9%)
  Low priority (<0.5):        704 ( 92.5%)

Score Statistics:
  Mean:   0.1245
  Median: 0.0823
  Std:    0.1567
  Min:    0.0012
  Max:    0.9234

Percentile Distribution:
  90th percentile: 0.3421 (76 variants)
  95th percentile: 0.5012 (38 variants)
  99th percentile: 0.7834 (8 variants)

Variant Prioritization and Export¶

Extract and rank high-confidence functional predictions for downstream experimental validation.

R:
    suppressPackageStartupMessages(library(vroom))
    results <- vroom("tests/fixtures/ems_training/expected/predictions_weighted_model_chr2.tsv", show_col_types = FALSE)
    hp <- results[results$standard_subset_weighted_pred_prob > 0.8, ]
    hp <- hp[order(-hp$standard_subset_weighted_pred_prob), ]
    cat("HIGH-PRIORITY VARIANTS (Score >0.8)\n"); cat(strrep("=", 50), "\n", sep = "")
    cat("Total:", nrow(hp), "variants\n")
    if (nrow(hp) > 0) {
        cols <- intersect(c("variant_id", "standard_subset_weighted_pred_prob", "actual_label"), names(hp))
        print(utils::head(hp[, cols], 10))
        vroom::vroom_write(hp, "high_priority_variants_validation.tsv", delim = "\t")
        cat("\nExported", nrow(hp), "high-priority variants to high_priority_variants_validation.tsv\n")
    } else cat("\nNo variants exceed 0.8 threshold\n")

Example output:

HIGH-PRIORITY VARIANTS (Score >0.8)
==================================================
Total: 12 variants

Top 10 by EMS score:
      variant_id  standard_subset_weighted_pred_prob  actual_label
   2:54321:A:G                              0.9234             1
   2:12345:T:C                              0.9102             1
   2:98765:C:A                              0.8956             1
   2:44444:G:T                              0.8723             0
   2:77777:A:C                              0.8612             1
   2:33333:T:G                              0.8501             1
   2:66666:C:T                              0.8398             1
   2:11111:G:A                              0.8267             0
   2:55555:A:T                              0.8145             1
   2:99999:T:A                              0.8034             1

  Exported to: high_priority_variants_validation.tsv
   Contains all 12 high-priority variants
   Ready for: CRISPR screens, luciferase assays, functional validation

Verify Model Performance¶

Review metrics from the held-out test set to understand model reliability.

# Model performance on the held-out test set (written by [train] as summary JSON).
R:
    suppressPackageStartupMessages(library(jsonlite))
    summ <- fromJSON("tests/fixtures/ems_training/expected/model_5_summary_chr_chr2_NPR_1.json")
    m <- summ$CatBoost$standard_subset_weighted
    cat("MODEL PERFORMANCE ON TEST SET\n"); cat(strrep("=", 50), "\n", sep = "")
    cat(sprintf("Average Precision (AP): %.4f\n", m$AP_test))
    cat(sprintf("AUC-ROC:                %.4f\n", m$AUC_test))
    npos <- summ$CatBoost$test_num_positive_labels; nneg <- summ$CatBoost$test_num_negative_labels
    cat(sprintf("\nTest set: %d positive, %d negative (%.1f%% positive)\n",
                npos, nneg, 100*npos/(npos+nneg)))

Example output:

MODEL PERFORMANCE ON TEST SET
==================================================
Average Precision (AP): 0.5050
AUC-ROC:                0.8978

Test Set Composition:
  Positive labels (functional eQTLs): 68
  Negative labels (non-functional):   693
  Positive rate:                       8.9%

   These metrics reflect performance on 761 held-out chromosome 2 variants
   with stricter selection criteria (PIP >0.9 for positives) than training data.

Application to New Variant Lists¶

Workflow for Novel Variants¶

To score additional variants not included in the original training/test sets:

1. Prepare input file matching the format above (chr:pos:ref:alt)

2. Update configuration (data_config.yaml):

Point training_data.base_dir to directory containing your variant annotations
Ensure gene constraint file (GeneBayes scores) is accessible
Verify MAF file matches your chromosome

3. Execute pipeline (same command, different input data):

python model_training_model5_only.py [cohort] [chromosome] \
  --data-config data_config.yaml \
  --model-config model_config.yaml

The pipeline will:

Generate all 4,839 genomic features for your variants
Apply the same preprocessing (subsetting, absolute values, imputation)
Use the trained model for inference
Output predictions in identical format

4. Analyze results using the code blocks above

Steps¶

Step 1. Score variants with a trained feature-weighted CatBoost scEEMS model for one cohort / chromosome. The toy command below scores the microglia cohort (protocol_example) on chromosome 2 using a model produced by the training workflow.

Timing: ~varies on typical compute infrastructure.

sos run pipeline/ems_prediction.ipynb predict \
    --cohort protocol_example \
    --chromosome 2 \
    --model_path output/xqtl_modifier_score/protocol_example/model_results/model_standard_subset_weighted_chr_chr2_NPR_1.cbm \
    --data_config code/SoS/xqtl_modifier_score/data_config.yaml \
    --cwd output/ems_prediction

Command Interface¶

sos run pipeline/ems_prediction.ipynb -h

Workflow implementation¶

The predict step wraps the gems_pipeline.R predict engine. The pipeline loads the trained model, annotates the input variants with the full feature set, applies the 10x deep-learning feature weighting used at training time, and writes per-variant EMS scores. New datasets / cell types are scored by editing the data_config YAML only - no code changes are required.

[global]
import os
# Work directory & output directory
parameter: cwd = path('output/ems_prediction')
# Cohort / cell type to score (must match an entry in the data config)
parameter: cohort = 'protocol_example'
# Chromosome to score
parameter: chromosome = '2'
# Trained CatBoost model (.cbm) from the EMS Training workflow
parameter: model_path = path('output/xqtl_modifier_score/protocol_example/model_results/model_standard_subset_weighted_chr_chr2_NPR_1.cbm')
# Data configuration YAML (cohort, variant list, feature set)
parameter: data_config = path('code/SoS/xqtl_modifier_score/data_config.yaml')
# Directory whose relative config paths resolve against (holds the config YAMLs)
parameter: pipeline_dir = path('code/SoS/xqtl_modifier_score')
# Directory holding the modular analysis scripts (code/script)
parameter: modular_script_dir = path('code/script')
parameter: job_size = 1
parameter: mem = '60G'
parameter: walltime = '24h'

[predict]
bash: expand = "$[ ]", workdir = pipeline_dir
    Rscript $[modular_script_dir:a]/xqtl_modifier_score/gems_pipeline.R \
        --step predict \
        --cohort $[cohort] \
        --chromosome $[chromosome] \
        --model-path $[model_path:a] \
        --data-config $[data_config:a]

Anticipated Results¶

Running the predict step produces per-variant scEEMS prediction scores for the supplied variant list. See the Pipeline Execution Workflow walkthrough above for the statistical summary, variant prioritization, and model-performance verification of these outputs.