Scoring a variant

Scoring a variant#

Scoring summarizes the difference between reference (REF) and alternate (ALT) allele predictions as one number per track (or, in the case of gene-specific scorers, one number per gene per track). This tutorial scores chr2:219243701:C>T near GLB1L in CD14-positive monocytes in blood (UBERON:0000178_CL:0001054) with recommended and custom variant scorers.

Key activities:

Set up#

# @title Install AlphaGenome

# @markdown Run this cell to install AlphaGenome.
from IPython.display import clear_output
! pip install alphagenome
clear_output()

# @title Imports
from alphagenome import colab_utils
from alphagenome.data import genome
from alphagenome.models.v2 import dna_client
from alphagenome.models.v2 import variant_scorers

Store your API key in Colab Secrets for persistent access across sessions (see Installation).

dna_model = dna_client.create(colab_utils.get_api_key())

Score a variant#

score_variant predicts the REF and ALT sequences and aggregates their difference with one or more variant scorers. It takes a condition, a genome.Variant, a genome.Interval that contains it, and a list of scorers.

Set condition#

Use a reference condition CURIE from Choosing a condition.

# Set condition.
target_condition = 'UBERON:0000178_CL:0001054'  # CD14-positive monocytes in blood.

Tip: For expression conditioning (cell_conditioning= alongside cell_condition_name=) with your own human scRNA-seq counts, use variant_scorers.RECOMMENDED_VARIANT_SCORERS['SC_GENE_EXPRESSION_DENOISED'] to score denoised single-cell gene expression; see Conditioning on your own single-cell data.

Define variant and interval#

Define a variant by creating a genome.Variant object (see Selecting a genomic region):

# Set variant.
variant = genome.Variant(
    chromosome='chr2',
    position=219243701,
    reference_bases='C',
    alternate_bases='T',
)
print(f'Variant: {variant}')

Variant: chr2:219243701:C>T

Center a 1 Mb interval on the variant with variant.reference_interval.resize(dna_client.SEQUENCE_LENGTH_1MB) (see Selecting a genomic region):

# Fetch a 1 Mb interval around the above variant.
interval = variant.reference_interval.resize(dna_client.SEQUENCE_LENGTH_1MB)

Choose scorers#

variant_scorers.RECOMMENDED_VARIANT_SCORERS maps each scored output type across bulk and single-cell assays to a recommended scorer, plus _ACTIVE scorers for quantifying absolute activity where applicable. If you omit variant_scorers, score_variant uses the recommended scorers for the organism and condition mode. See How variant scoring works for more details.

# List recommended scoring methods.
for name, scorer in variant_scorers.RECOMMENDED_VARIANT_SCORERS.items():
  print(f"{name:<20} -> {scorer}")

ATAC                 -> CenterMaskScorer(requested_output=ATAC, width=501, aggregation_type=DIFF_LOG2_SUM)
ATAC_ACTIVE          -> CenterMaskScorer(requested_output=ATAC, width=501, aggregation_type=ACTIVE_SUM)
CAGE                 -> CenterMaskScorer(requested_output=CAGE, width=501, aggregation_type=DIFF_LOG2_SUM)
CAGE_ACTIVE          -> CenterMaskScorer(requested_output=CAGE, width=501, aggregation_type=ACTIVE_SUM)
CHIP_HISTONE         -> CenterMaskScorer(requested_output=CHIP_HISTONE, width=2001, aggregation_type=DIFF_LOG2_SUM)
CHIP_HISTONE_ACTIVE  -> CenterMaskScorer(requested_output=CHIP_HISTONE, width=2001, aggregation_type=ACTIVE_SUM)
CHIP_TF              -> CenterMaskScorer(requested_output=CHIP_TF, width=501, aggregation_type=DIFF_LOG2_SUM)
CHIP_TF_ACTIVE       -> CenterMaskScorer(requested_output=CHIP_TF, width=501, aggregation_type=ACTIVE_SUM)
CONTACT_MAPS         -> ContactMapScorer()
DNASE                -> CenterMaskScorer(requested_output=DNASE, width=501, aggregation_type=DIFF_LOG2_SUM)
DNASE_ACTIVE         -> CenterMaskScorer(requested_output=DNASE, width=501, aggregation_type=ACTIVE_SUM)
METHYLATION          -> MethylationScorer(width=1201)
POLYADENYLATION      -> PolyadenylationScorer()
RNA_SEQ              -> GeneMaskLFCScorer(requested_output=RNA_SEQ, aggregation_type=DIFF_MEAN)
RNA_SEQ_ACTIVE       -> GeneMaskActiveScorer(requested_output=RNA_SEQ)
SC_ATAC              -> CenterMaskScorer(requested_output=SC_ATAC, width=1201, aggregation_type=DIFF_LOG2_SUM)
SC_ATAC_ACTIVE       -> CenterMaskScorer(requested_output=SC_ATAC, width=1201, aggregation_type=ACTIVE_SUM)
SC_GENE_EXPRESSION   -> GeneMaskLFCScorer(requested_output=SC_GENE_EXPRESSION, aggregation_type=DIFF_SUM)
SC_GENE_EXPRESSION_ACTIVE -> GeneMaskActiveScorer(requested_output=SC_GENE_EXPRESSION)
SC_GENE_EXPRESSION_DENOISED -> GeneMaskLFCScorer(requested_output=SC_GENE_EXPRESSION_DENOISED, aggregation_type=DIFF_SUM)
SPLICE_JUNCTIONS     -> SpliceJunctionScorer()
SPLICE_SITES         -> GeneMaskSplicingScorer(requested_output=SPLICE_SITES, width=None)
SPLICE_SITE_USAGE    -> GeneMaskSplicingScorer(requested_output=SPLICE_SITE_USAGE, width=None)

Tip

For an example of scoring a splicing variant, see Analyzing a splicing variant.

To customize scoring, you can construct a new scorer manually. The one below averages the ALT minus REF DNase-seq difference over a 501 bp window centered on the variant:

# Create a variant scoring method.
dnase_variant_scorer = variant_scorers.CenterMaskScorer(
    requested_output=dna_client.OutputType.DNASE,
    width=501,
    aggregation_type=variant_scorers.AggregationType.DIFF_MEAN,
)

We will also include the recommended scorer for single-cell gene expression:

# Fetch a recommended scoring method for single-cell gene expression.
sc_gex_scorer = variant_scorers.RECOMMENDED_VARIANT_SCORERS[
    'SC_GENE_EXPRESSION'
]

Run score_variant#

Score variant effects using score_variant:

# Score variant with both the recommended scRNA-seq scorer and custom DNase scorer.
variant_scores = dna_model.score_variant(
    interval=interval,
    variant=variant,
    variant_scorers=[sc_gex_scorer, dnase_variant_scorer],
    condition=target_condition,
)

The returned variant_scores list contains one AnnData object per requested scorer:

len(variant_scores)

2

Each score object is an AnnData instance (see Key objects), storing numerical scores in .X, track metadata in .var, and gene metadata (for gene-centric scorers) in .obs:

sc_gex_scores = variant_scores[0]
sc_gex_scores

AnnData object with n_obs × n_vars = 73 × 7
    obs: 'gene_id', 'strand', 'gene_name', 'gene_type'
    var: 'name', 'strand', 'biosample_curie', 'imputation_mode', 'biosample_type', 'biosample_name', 'Assay title', 'suspension_type', 'data_pipeline', 'cellxgene_dataset_id', 'paper_doi'
    uns: 'interval', 'variant', 'variant_scorer'
    layers: None (.X)

For improved readability, we recommend flattening AnnData score objects using variant_scorers.tidy_scores, which returns a DataFrame with one row per combination of (variant, gene, scorer, track). By default (match_gene_strand=True) it excludes stranded tracks on the opposite strand to the gene for gene-centric scorers, and filter_protein_coding=True keeps only protein-coding genes:

# Flatten scores from both scorers into a tidy DataFrame.
tidy_scores = variant_scorers.tidy_scores(
    variant_scores, match_gene_strand=True, filter_protein_coding=True
)
tidy_scores
variant_id scored_interval gene_id gene_name gene_type gene_strand junction_Start junction_End output_type variant_scorer track_name track_strand imputation_mode biosample_curie biosample_type biosample_name Assay title suspension_type data_pipeline cellxgene_dataset_id paper_doi data_source genetically_modified endedness raw_score
0 chr2:219243701:C>T chr2:218719413-219767989:. ENSG00000054356 PTPRN protein_coding - None None SC_GENE_EXPRESSION GeneMaskLFCScorer(requested_output=SC_GENE_EXP... UBERON:0000178_CL:0001054_EasySci-easy_sci-nuc... . imputed_supported UBERON:0000178_CL:0001054 primary_cell blood CD14-positive monocyte EasySci nuclear easy_sci NaN NaN NaN NaN NaN 0.013146
1 chr2:219243701:C>T chr2:218719413-219767989:. ENSG00000054356 PTPRN protein_coding - None None SC_GENE_EXPRESSION GeneMaskLFCScorer(requested_output=SC_GENE_EXP... UBERON:0000178_CL:0001054_SHARE-seq-share_seq-... . imputed_supported UBERON:0000178_CL:0001054 primary_cell blood CD14-positive monocyte SHARE-seq nuclear share_seq NaN NaN NaN NaN NaN 0.008913
2 chr2:219243701:C>T chr2:218719413-219767989:. ENSG00000054356 PTPRN protein_coding - None None SC_GENE_EXPRESSION GeneMaskLFCScorer(requested_output=SC_GENE_EXP... UBERON:0000178_CL:0001054_SPLIT-seq-star_solo-... . imputed_supported UBERON:0000178_CL:0001054 primary_cell blood CD14-positive monocyte SPLIT-seq nuclear encode NaN NaN NaN NaN NaN 0.011879
3 chr2:219243701:C>T chr2:218719413-219767989:. ENSG00000054356 PTPRN protein_coding - None None SC_GENE_EXPRESSION GeneMaskLFCScorer(requested_output=SC_GENE_EXP... UBERON:0000178_CL:0001054_scRNA-seq 10X 3-cell... . imputed_supported UBERON:0000178_CL:0001054 primary_cell blood CD14-positive monocyte scRNA-seq 10X 3 cell cellranger NaN NaN NaN NaN NaN 0.011667
4 chr2:219243701:C>T chr2:218719413-219767989:. ENSG00000054356 PTPRN protein_coding - None None SC_GENE_EXPRESSION GeneMaskLFCScorer(requested_output=SC_GENE_EXP... UBERON:0000178_CL:0001054_scRNA-seq 10X 3-cell... . imputed_supported UBERON:0000178_CL:0001054 primary_cell blood CD14-positive monocyte scRNA-seq 10X 3 nuclear cellranger NaN NaN NaN NaN NaN 0.009927
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
275 chr2:219243701:C>T chr2:218719413-219767989:. ENSG00000286143 ENSG00000286143 protein_coding + None None SC_GENE_EXPRESSION GeneMaskLFCScorer(requested_output=SC_GENE_EXP... UBERON:0000178_CL:0001054_SPLIT-seq-star_solo-... . imputed_supported UBERON:0000178_CL:0001054 primary_cell blood CD14-positive monocyte SPLIT-seq nuclear encode NaN NaN NaN NaN NaN 0.002575
276 chr2:219243701:C>T chr2:218719413-219767989:. ENSG00000286143 ENSG00000286143 protein_coding + None None SC_GENE_EXPRESSION GeneMaskLFCScorer(requested_output=SC_GENE_EXP... UBERON:0000178_CL:0001054_scRNA-seq 10X 3-cell... . imputed_supported UBERON:0000178_CL:0001054 primary_cell blood CD14-positive monocyte scRNA-seq 10X 3 cell cellranger NaN NaN NaN NaN NaN 0.001796
277 chr2:219243701:C>T chr2:218719413-219767989:. ENSG00000286143 ENSG00000286143 protein_coding + None None SC_GENE_EXPRESSION GeneMaskLFCScorer(requested_output=SC_GENE_EXP... UBERON:0000178_CL:0001054_scRNA-seq 10X 3-cell... . imputed_supported UBERON:0000178_CL:0001054 primary_cell blood CD14-positive monocyte scRNA-seq 10X 3 nuclear cellranger NaN NaN NaN NaN NaN 0.004684
278 chr2:219243701:C>T chr2:218719413-219767989:. ENSG00000286143 ENSG00000286143 protein_coding + None None SC_GENE_EXPRESSION GeneMaskLFCScorer(requested_output=SC_GENE_EXP... UBERON:0000178_CL:0001054_scRNA-seq 10X 3-star... . imputed_supported UBERON:0000178_CL:0001054 primary_cell blood CD14-positive monocyte scRNA-seq 10X 3 nuclear encode NaN NaN NaN NaN NaN 0.001952
279 chr2:219243701:C>T chr2:218719413-219767989:. ENSG00000286143 ENSG00000286143 protein_coding + None None SC_GENE_EXPRESSION GeneMaskLFCScorer(requested_output=SC_GENE_EXP... UBERON:0000178_CL:0001054_scRNA-seq 10X 5-cell... . observed UBERON:0000178_CL:0001054 primary_cell blood CD14-positive monocyte scRNA-seq 10X 5 cell cellranger 3a8a77a6-f069-479c-a2c1-252feafd0106 10.1016/j.cell.2025.02.017 NaN NaN NaN 0.000587

280 rows × 25 columns

Filter to observed tracks for GLB1L and any interval-level scorers (where gene_name is NaN):

tidy_scores.query(
    '(gene_name == "GLB1L" or gene_name.isna()) and imputation_mode == "observed"'
)
variant_id scored_interval gene_id gene_name gene_type gene_strand junction_Start junction_End output_type variant_scorer track_name track_strand imputation_mode biosample_curie biosample_type biosample_name Assay title suspension_type data_pipeline cellxgene_dataset_id paper_doi data_source genetically_modified endedness raw_score
188 chr2:219243701:C>T chr2:218719413-219767989:. ENSG00000163521 GLB1L protein_coding - None None SC_GENE_EXPRESSION GeneMaskLFCScorer(requested_output=SC_GENE_EXP... UBERON:0000178_CL:0001054_scRNA-seq 10X 5-cell... . observed UBERON:0000178_CL:0001054 primary_cell blood CD14-positive monocyte scRNA-seq 10X 5 cell cellranger 3a8a77a6-f069-479c-a2c1-252feafd0106 10.1016/j.cell.2025.02.017 NaN NaN NaN -1.835707

Score multiple variants#

To score multiple variants across intervals and conditions in parallel, pass matching lists of intervals, variants, and conditions to dna_model.score_variants, then flatten the returned nested list of AnnData objects with variant_scorers.tidy_scores:

variants = [
    genome.Variant(
        chromosome='chr2',
        position=219243701,
        reference_bases='C',
        alternate_bases='T',
    ),
    genome.Variant(
        chromosome='chr2',
        position=219243701,
        reference_bases='C',
        alternate_bases='G',
    ),
]
intervals = [
    v.reference_interval.resize(dna_client.SEQUENCE_LENGTH_1MB)
    for v in variants
]

batch_scores = dna_model.score_variants(
    intervals=intervals,
    variants=variants,
    variant_scorers=[sc_gex_scorer],
    conditions=[target_condition] * len(variants),
)
batch_tidy = variant_scorers.tidy_scores(
    batch_scores, match_gene_strand=True, filter_protein_coding=True
)
batch_tidy.query('gene_name == "GLB1L" and imputation_mode == "observed"')
variant_id scored_interval gene_id gene_name gene_type gene_strand junction_Start junction_End output_type variant_scorer track_name track_strand imputation_mode biosample_curie biosample_type biosample_name Assay title suspension_type data_pipeline cellxgene_dataset_id paper_doi raw_score
266 chr2:219243701:C>T chr2:218719413-219767989:. ENSG00000163521 GLB1L protein_coding - None None SC_GENE_EXPRESSION GeneMaskLFCScorer(requested_output=SC_GENE_EXP... UBERON:0000178_CL:0001054_scRNA-seq 10X 5-cell... . observed UBERON:0000178_CL:0001054 primary_cell blood CD14-positive monocyte scRNA-seq 10X 5 cell cellranger 3a8a77a6-f069-479c-a2c1-252feafd0106 10.1016/j.cell.2025.02.017 -1.835707
546 chr2:219243701:C>G chr2:218719413-219767989:. ENSG00000163521 GLB1L protein_coding - None None SC_GENE_EXPRESSION GeneMaskLFCScorer(requested_output=SC_GENE_EXP... UBERON:0000178_CL:0001054_scRNA-seq 10X 5-cell... . observed UBERON:0000178_CL:0001054 primary_cell blood CD14-positive monocyte scRNA-seq 10X 5 cell cellranger 3a8a77a6-f069-479c-a2c1-252feafd0106 10.1016/j.cell.2025.02.017 -1.867211

Next in the core workflow, we will explore how to plot predictions.