xQTL Model Training Tutorial

xQTL Model Training Tutorial#

This notebook implements the sc-EMS (single-cell Expression Modifier Scores) methodology from our previous research (FIXME). The original research demonstrated that machine learning can predict which genetic variants affect gene expression in specific brain cell types, particularly for Alzheimer’s disease research.

Motivation#

Traditional genetic studies only explain 5% of Alzheimer’s disease heritability. The remaining 95% “missing heritability” occurs because bulk tissue approaches cannot capture regulatory effects that occur in specific cell types. Our machine learning approach identifies genetic variants that affect gene expression specifically in brain cell types like microglia, potentially explaining previously undetected disease mechanisms.

  • Our Motivation is to create a reproducible pipeline that accurately predicts functional genetic variants using cell-type specific genomic annotations, transforming complex research findings into a standardized tool for the broader scientific community.

Methods Overview#

Feature-Weighted CatBoost Algorithm#

We use a single CatBoost (Categorical Boosting) gradient boosting model optimized for genomic data. CatBoost builds an ensemble of decision trees sequentially, where each tree learns from previous errors, making it particularly effective for high-dimensional biological datasets with mixed data types.

Our Training Strategy: Single feature-weighted model (Model 5) that emphasizes biology-informed features while maintaining comprehensive genomic context.

Training Data Construction#

Data Source: Chromosome 2 variants with non-overlapping splits

  • Training set: 3,056 variants (80% of available data)

  • Testing set: 761 variants (20% of available data)

  • Critical validation: Zero variant overlap between training and testing sets

Label Definition (Y):

  • Y = 1: Functional eQTL (PIP > 0.1) - variant significantly affects gene expression

  • Y = 0: Non-functional variant (PIP < 0.01) - no detectable expression effect

  • Class distribution: 9% positive rate reflects biological reality that most variants are non-functional

Input Data: Feature Categories from Manuscript Table C#

Based on our previous research (FIXME) (Figure 1, Table C), we use these broad feature categories. Each category represents a branch of features that can be applied to many different contexts and datasets.

# Load gene constraint data
import pandas as pd
gene_data = pd.read_excel('data/41588_2024_1820_MOESM4_ESM.xlsx')
print(gene_data.head(3))

Sample output showing GeneBayes constraint scores:

ensg            hgnc        chrom  obs_lof  exp_lof   post_mean    post_lower_95  post_upper_95
ENSG00000198488 HGNC:24141  chr11    12      8.9777   6.46629E-05  7.16256E-06   0.00017805
ENSG00000164363 HGNC:26441  chr5     31      28.55    0.00016062   2.59918E-05   0.00044175
ENSG00000159337 HGNC:30038  chr15    28      41.84    0.00016978   0.000018674   0.00053317

Feature Categories#

1. Distance Features#

Biological Rationale: Physical proximity determines regulatory potential

  • abs_distance_TSS_log: Log-transformed distance to transcription start site

  • Example: Value 8.2 = ~160,000 base pairs from gene start

  • Impact: Closest predictor of regulatory effects (importance: 23.58)

2. Cell-Type Regulatory Features#

Biological Rationale: Microglia-specific regulatory mechanisms

  • ABC_score_microglia: Activity-by-Contact regulatory activity score

  • Example: Score 0.73 indicates high regulatory potential in microglia

  • diff_32_ENCFF140MBA: Cell-type chromatin accessibility change

  • Example: Value 1.83 indicates strong accessibility difference

  • Impact: Multiple diff_32_* features contribute consistently (0.5-1.8 importance)

3. Population Genetics Features#

Biological Rationale: Evolutionary constraint and population frequency

  • gnomad_MAF: Minor allele frequency from population database

  • Example: 0.15 = variant found in 15% of human population

  • Impact: Secondary predictor (importance: 0.88)

4. Conservation Features#

Biological Rationale: Evolutionary preservation indicates functional importance

  • Cross-species conservation scores and constraint metrics

  • Impact: Moderate contribution to overall prediction

5. Deep Learning Predictions#

Biological Rationale: Sequence-to-function model predictions

  • Enformer and BPNet chromatin accessibility predictions

  • Impact: Complementary information beyond experimental features

Output Data Generated#

Output Type

Description

Research Use

Trained models

Single feature-weighted CatBoost classifier

Variant effect prediction

Performance metrics

AP/AUC scores

Model comparison

Feature importance

Genomic feature rankings

Biological interpretation

What We’re Trying to Predict (Y)#

Binary Classification:

  • Y = 1: Functional variant (affects gene expression)

  • Y = 0: Non-functional variant (no effect on genes)

Training Data Breakdown:

  • Total: 3,056 variants

  • Functional (Y=1): 278 variants (9%)

  • Non-functional (Y=0): 2,778 variants (91%)

  • This 9% rate reflects biological reality - most genetic variants don’t affect genes

Description of Toy Dataset#

Dataset Overview#

This tutorial demonstrates the sc-EMS (single-cell Expression Modifier Scores) methodology from our previous research (FIXME), which predicts brain cell type eQTLs for Alzheimer’s disease research. Our toy dataset represents a subset of the full study data:

Dataset specifications:

  • Training variants: 3,056 genetic variants from chromosome 2

  • Testing variants: 761 genetic variants from chromosome 2 (non-overlapping)

  • Feature dimensions: 4,839 genomic annotations per variant

  • Target prediction: Binary classification (functional vs. non-functional eQTLs)

What makes this a “toy” dataset:

  • We use a smaller, simpler version of the full research study

  • The full study used all 22 chromosomes with thousands of people

  • Our tutorial uses only chromosome 2 data from a smaller sample

Why we can still prevent data leakage:

  • We split chromosome 2 data into two completely separate groups

  • Training data: 3,056 variants with looser selection criteria

  • Test data: 761 variants with stricter selection criteria

  • No variant appears in both groups

Data Split Strategy#

Critical Design: Our dataset is smaller and simpler than the real research, but it teaches the same machine learning method and maintains scientific rigor by keeping training and testing data completely separate

Key Parameters#

Model 5 Hyperparameters#

  • Tree depth: 6 (balances complexity vs. overfitting)

  • L2 regularization: 3.0 (controls genomic noise)

  • Learning rate: 0.03 (conservative for stability)

  • Iterations: 1000 (sufficient for convergence)

Feature Weighting Strategy:#

  • Standard features: Weight = 1.0 (distance, conservation, population genetics)

  • High-priority regulatory features: Weight = 10.0

    • chrombpnet_positive: Chromatin accessibility predictions

    • tf_positive: Transcription factor binding sites

    • diff_*: Cell-type specific differential signals

Biological Rationale: Experimentally validated regulatory features receive higher weight than purely computational predictions, reflecting greater confidence in biology-based annotations.

Cross-Validation Configuration#

Leave-One-Chromosome-Out (LOCO) Procedure#

The sc-EMS methodology uses chromosome-level cross-validation to prevent data leakage:

  • Training approach: For chromosome i, train on all other chromosomes j where j ≠ i

  • Testing approach: Evaluate performance on held-out chromosome i

  • Implementation: 22 separate models trained (one per chromosome)

Mathematical Framework#

For chromosome validation:

Detailed Workflow#

This section explains how to use our xQTL training pipeline

Step 1: Running the Pipeline#

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

```python
# file check

import os
required_files = [
    'model_training_simplified.py',
    'data_config.yaml',
    'model_config.yaml'
]
print("Configuration File Check:")
for file in required_files:
    status = "✅ Found" if os.path.exists(file) else "❌ Missing"
    print(f"  {file}: {status}")

Step 2: Understanding Training Process Overview#

  1. Load leak-free training data (3,056 variants)

  2. Train single feature-weighted CatBoost model

  3. Evaluate on held-out test set (761 variants)

  4. Generate feature importance rankings and predictions

Step 3: Results & Results Analysis#

# Performance results 
print("Model Performance Results:")
print("=" * 50)

# Single model results (corrected from previous 4-model approach)
corrected_results = {
    'Feature-Weighted Model (Model 5)': {'AP': 0.5050, 'AUC': 0.8978}
}

for model, metrics in corrected_results.items():
    print(f"{model}:")
    print(f"  Average Precision: {metrics['AP']:.4f}")
    print(f"  AUC Score: {metrics['AUC']:.4f}")
    print()

print("Feature Importance (Top 10):")
print("-" * 40)
top_features = [
    ("abs_distance_TSS_log", 23.58),
    ("diff_32_ENCFF140MBA", 1.83),
    ("diff_32_ENCFF457MJJ", 1.12),
    ("gnomad_MAF", 0.88),
    ("diff_32_ENCFF455KGB", 0.84),
    ("diff_32_ENCFF579IKK", 0.80),
    ("ABC_score_microglia", 0.54),
    ("microglia_enhancer_promoter_union_atac", 0.58),
    ("diff_32_ENCFF098QMC", 0.66),
    ("diff_32_ENCFF757EYT", 0.66)
]

for i, (feature, importance) in enumerate(top_features, 1):
    print(f"{i:2d}. {feature:<35} {importance:>6.2f}")

Results Analysis#

🎯 Strong Performance Indicators (Our Results):#

  • AUC: 0.8978 - Excellent discrimination (89.78% > target of 75%)

  • Average Precision: 0.5050 - Strong precision (50.5% > target of 40%)

  • Distance dominance: Clear biological pattern (23.58 importance)

  • Feature diversity: Multiple regulatory signals contribute (0.5-1.8 importance)

✅ What This Means:#

  • Realistic genomic performance - Not artificially perfect, indicates genuine learning

  • Biology makes sense - Distance matters most, regulatory features add value

  • Model works correctly - Can distinguish functional from non-functional variants

🔬 Key Biological Findings:#

  • Distance rules everything: 20x more important than other features

  • Cell-type features help: Microglia-specific signals provide additional information

  • Feature weighting worked: 10x weighting successfully elevated regulatory signals

  • Population genetics secondary: MAF contributes but less than regulatory data

Step 4: Generated Output Files#

Model File:

  • model_feature_weighted_chr_chr2_NPR_10.joblib - Trained CatBoost classifier

Analysis Results:

  • summary_dict_catboost_1model_chr_chr2_NPR_10.pkl - Performance metrics and validation statistics

  • features_importance_1model_chr_chr2_NPR_10.csv - Complete feature importance rankings

  • predictions_1model_chr2.tsv - Per-variant prediction probabilities

Model Application#

The trained model predicts functional genetic variants with 89.78% discriminative accuracy, enabling:

  • Identification of disease-relevant variants in cell-type specific contexts

  • Prioritization of variants for experimental validation

  • Understanding of regulatory mechanisms in neurodegenerative diseases

Step 5: Using Your Trained Models#

Once training is complete, load the trained models for predictions. Please refer to EMS Predictions (https://statfungen.github.io/xqtl-protocol/code/xqtl_modifier_score/ems_prediction.html) for detailed prediction workflows and variant scoring.

Contents