Selecting a genomic region

Selecting a genomic region#

AlphaGenome 2 makes predictions for a DNA sequence of fixed length, centered on the region you care about. After choosing a condition, the second step of an analysis is to describe that input: a genome.Interval on the reference genome, an interval around a specific genome.Variant, or a custom DNA string.

Here we set up the genomic region around the GLB1L gene and a nearby variant, chr2:219243701:C>T.

Key activities:

  1. Create a genome.Interval and resize it to the recommended 1 Mb context length (dna_client.SEQUENCE_LENGTH_1MB).

  2. Look up the interval of a gene from GENCODE annotations.

  3. Define a genome.Variant and center the 1 Mb interval on the variant.

  4. Prepare a custom DNA sequence string for predict_sequence.

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.data import gene_annotation
from alphagenome.data import genome
from alphagenome.models.v2 import dna_client
import pandas as pd

This tutorial does not call the model, so it does not need an API key.

Supported sequence lengths#

AlphaGenome 2 accepts four input sequence lengths: 16 kb (16,384 bp, dna_client.SEQUENCE_LENGTH_16KB), 131,072 bp (dna_client.SEQUENCE_LENGTH_100KB), 524,288 bp (dna_client.SEQUENCE_LENGTH_500KB), and 1 Mb (1,048,576 bp, dna_client.SEQUENCE_LENGTH_1MB):

list(dna_client.SUPPORTED_SEQUENCE_LENGTHS.keys())

['SEQUENCE_LENGTH_16KB',
 'SEQUENCE_LENGTH_100KB',
 'SEQUENCE_LENGTH_500KB',
 'SEQUENCE_LENGTH_1MB']

Important

Use a 1 Mb interval (dna_client.SEQUENCE_LENGTH_1MB) whenever possible for best accuracy, and center variants inside the interval (variant.reference_interval.resize(dna_client.SEQUENCE_LENGTH_1MB)) so the model sees sequence context on both sides (see Start from a variant).

Define an interval#

A genome.Interval is a genomic region with a 0-based start and an exclusive end position:

interval = genome.Interval(chromosome='chr2', start=219_200_000, end=219_300_000)

print(interval)
print(f'Width: {interval.width:,} bp')

chr2:219200000-219300000:.
Width: 100,000 bp

resize expands (or shrinks) the interval around its center to the requested width, adding reference genome sequence on both sides. Use dna_client.SEQUENCE_LENGTH_1MB for full context:

interval_1mb = interval.resize(dna_client.SEQUENCE_LENGTH_1MB)

print(interval_1mb)
print(f'Width: {interval_1mb.width:,} bp')

chr2:218725712-219774288:.
Width: 1,048,576 bp

Start from a gene#

Often, you will want to start from a gene rather than from coordinates. GENCODE annotations (in GTF format) give the location of every gene:

# GENCODE v46 gene annotations for the hg38 human genome assembly.
gtf = pd.read_feather(
    'https://storage.googleapis.com/alphagenome/reference/gencode/'
    'hg38/gencode.v46.annotation.gtf.gz.feather'
)

gene_interval = gene_annotation.get_gene_interval(gtf, gene_symbol='GLB1L')
gene_interval

Interval(chromosome='chr2', start=219236597, end=219245478, strand='-', name='GLB1L')

get_gene_interval also records whether the gene lies on the + or - strand. Remember to resize the interval to dna_client.SEQUENCE_LENGTH_1MB before making predictions:

gene_interval.resize(dna_client.SEQUENCE_LENGTH_1MB)

Interval(chromosome='chr2', start=218716749, end=219765325, strand='-', name='GLB1L')

Start from a variant#

A genome.Variant is specified by its chromosome, 1-based position (as in VCF files), and reference and alternate alleles:

variant = genome.Variant(
    chromosome='chr2',
    position=219243701,
    reference_bases='C',
    alternate_bases='T',
)
print(variant)

chr2:219243701:C>T

reference_interval covers the reference allele. Resizing reference_interval to dna_client.SEQUENCE_LENGTH_1MB centers the 1 Mb prediction window on the variant, which is recommended for variant effect predictions and scoring:

interval = variant.reference_interval.resize(dna_client.SEQUENCE_LENGTH_1MB)
interval

Interval(chromosome='chr2', start=218719413, end=219767989, strand='.', name='')

The centered 1 Mb window around the variant contains the whole GLB1L gene:

interval.contains(gene_interval)

True

Define a DNA sequence#

For dna_model.predict_sequence, provide a string of A, C, G, T, and N whose length matches one of dna_client.SUPPORTED_SEQUENCE_LENGTHS. Pad shorter motifs or synthetic constructs with 'N', centering your sequence of interest (16 kb via dna_client.SEQUENCE_LENGTH_16KB or 1 Mb via dna_client.SEQUENCE_LENGTH_1MB). Note that extensively 'N'-padded sequences have not been benchmarked and may affect prediction performance:

sequence = 'GATTACA'.center(dna_client.SEQUENCE_LENGTH_16KB, 'N')
print(f'Sequence length: {len(sequence):,} bp')

Sequence length: 16,384 bp

Next steps#

You now have the inputs for the next step: a 1 Mb interval (optionally, with a variant centered inside it) or a custom DNA sequence.

For more operations on intervals and variants, such as indels and overlap checks, see Key objects.

Next in the core workflow, we will predict tracks for an interval or sequence, predict a variant, or score a variant.