STR/VNTR across populations

STR/VNTR across populations#

This is a demonstration of the downstream analysis after using vampire anno to annotate a TR locus across population.

[!NOTE] You can download this notebook by clicking the Download button in the upper-right corner of this page and selecting the .ipynb format.

import numpy as np
import vampire as vp
import logging
import plotly.io as pio
pio.renderers.default = "notebook_connected"

Set up logging for intermediate information.

vampire_logger = logging.getLogger("vampire")
vampire_logger.setLevel(logging.INFO) # set level to "DEBUG", "INFO", "WARNING"
if not vampire_logger.handlers:
    handler = logging.StreamHandler()
    handler.setFormatter(
        logging.Formatter(
            "[%(asctime)s] %(message)s",
            datefmt="%Y-%m-%d %H:%M:%S",
        )
    )
    vampire_logger.addHandler(handler)

Set the default plot style for all the plots in this notebook. You can set font_size, font_family, line_width, height, width, showgrid to set default parameters.

vp.anno.pl.set_default_plotstyle(font_size=14)

The dataset wdr7_hprc used in this tutorial was collected from the Yang et al., 2026. It contains annotations of a 69-bp motif VNTR locus (chr18:57,226,379–57,227,527 in T2T-CHM13v2.0) located within an intron of the WDR7 gene across 95 human haplotype assemblies, including 47 individuals from Phase 1 of the Human Pangenome Reference Consortium (HPRC) and the T2T-CHM13v2.0 reference genome.

Importing Data#

adata = vp.datasets.wdr7_hprc()
adata

AnnData object with n_obs × n_vars = 95 × 24
    obs: 'length', 'copy_number', 'score'
    var: 'motif', 'motif_length', 'copy_number', 'label'
    uns: 'motif_array', 'orientation_array', 'sequence'
    varp: 'motif_distance', 'rc_motif_distance'

Data Structure#

`X`#

Motif abundance / copy-number matrix with shape (n_obs × n_var).

X[i, j] = copy number of motif j in sample i

`obs`#

Sample-level metadata with shape (n_obs × metadata).

Column	Type	Description
`length`	`int`	Length of the TR sequence
`copy_number`	`float`	Total copy number of the locus
`score`	`int`	Annotation score

`var`#

Motif-level metadata with shape (n_var × metadata).

Column	Type	Description
`motif`	`str`	Motif sequence
`motif_length`	`int`	Length of the motif
`copy_number`	`float`	Total motif copy number across samples
`label`	`str`	Motif label

`varp`#

Motif-level pairwise relations with shape (n_var × n_var).

Key	Type	Description
`motif_distance`	`int`	Pairwise motif distance matrix
`rc_motif_distance`	`int`	Reverse-complement motif distance matrix

`uns`#

Unstructured genomic annotations not aligned to X.

Key	Type	Description
`sequence`	`Dict[str, str]`	Original TR sequences
`motif_array`	`Dict[str, List[str]]`	Motif ID arrays
`orientation_array`	`Dict[str, List[str]]`	Motif orientation arrays

Adding ancestry information#

adata.obs["sample"] = (
    adata.obs.index
    .astype(str)
    .str.replace(r"[\.][12]$", "", regex=True)
)

ancestry = vp.datasets.ancestry()
adata.obs["ancestry"] = adata.obs["sample"].map(ancestry)

adata.obs

	length	copy_number	score	sample	ancestry
sample
CHM13Y	1147	16.7	2276	CHM13Y	NA
HG002.1	666	9.7	1314	HG002	EUR
HG002.2	2183	31.7	4357	HG002	EUR
HG00438.1	1905	27.7	3783	HG00438	EAS
HG00438.2	3080	44.7	6151	HG00438	EAS
...	...	...	...	...	...
NA19240.2	732	10.7	1446	NA19240	AFR
NA20129.1	3977	57.7	7909	NA20129	AFR
NA20129.2	389	5.7	769	NA20129	AFR
NA21309.1	1147	16.7	2276	NA21309	AMR
NA21309.2	526	7.7	1025	NA21309	AMR

95 rows × 5 columns

Visualization across population#

You can set row_annotation="ancestry" to annotate samples according to the ancestry variable stored in .obs.

# sort samples by copy number
sample_order = list(
    adata.obs.sort_values(by="copy_number").index
)

# visualize the TR locus across population
# motifs exceeding the color limit will be displayed in black
vp.anno.pl.waterfall(
    adata,
    colormap = "rainbow", # choose from "rainbow", "glasbey" and "sequential", or you can specify colormap yourself
    sample_order = sample_order, # apply custom sample ordering
    row_annotation = "ancestry",
    track_name_dx = -0.01, # adjust horizontal position of sample labels
    # keyword arguments passed to plotly.graph_objects.Figure.update_layout(), including margin, font, etc.
    font=dict(size=8), # adjust to decrease the figure size
)

[2026-06-11 21:19:21] Number of unique labels (24) exceeds colormap size (11), the overflow labels will be colored in #0c0c0c

Set deduplicate=True to collapse identical tandem repeat structures and display each structure only once.

vp.anno.pl.waterfall(
    adata,
    colormap = "rainbow", # choose from "rainbow", "glasbey" and "sequential", or you can specify colormap yourself
    sample_order = sample_order, # apply custom sample ordering
    row_annotation = "ancestry",
    track_name_dx = -0.01, # adjust horizontal position of sample labels
    deduplicate = True, 
    # keyword arguments passed to plotly.graph_objects.Figure.update_layout(), including margin, font, etc.
    font=dict(size=8), # adjust to decrease the figure size
)

[2026-06-11 21:09:16] Number of unique labels (24) exceeds colormap size (11), the overflow labels will be colored in #0c0c0c
[2026-06-11 21:09:16] unique_group not found in adata.obs. vp.anno.pp.markdup() has not been run. Running it automatically.
[2026-06-11 21:09:16] markdup completed. 88 unique groups identified from 95 samples.

Generate a legend for the waterfall plot. Use waterfall_legend() with the same parameters as waterfall() except the color parameter. To display motif sequences, set color as "motif". To display motif IDs, set color as "id".

vp.anno.pl.waterfall_legend(
    adata,
    color="motif", 
    colormap = "rainbow",
    track_name_dx = -0.01,
    font=dict(size=12)
)

[2026-06-11 21:09:44] Number of unique labels (24) exceeds colormap size (11), the overflow labels will be colored in #0c0c0c

Motif Variation Pattern#

vp.anno.tl.motif_msa(adata)

[2026-06-11 21:09:46] Performing MSA of 24 motifs...
[2026-06-11 21:09:56] Refinement converged after 1 iteration(s).
[2026-06-11 21:09:56] MSA completed. Aligned length: 69. Found 55 variants.

AnnData object with n_obs × n_vars = 95 × 24
    obs: 'length', 'copy_number', 'score', 'sample', 'ancestry', 'unique_group'
    var: 'motif', 'motif_length', 'copy_number', 'label'
    uns: 'motif_array', 'orientation_array', 'sequence', 'motif_msa'
    varp: 'motif_distance', 'rc_motif_distance'

By default, in vp.anno.pl.logo(), conserved positions are shown in gray. Set conserved_color to None to disable this feature or specify a color to highlight conserved positions.

vp.anno.pl.logo(
    adata,
    feature = "information", # choose from "count", "probability" and "information"
    #drop_gap=True
)

vp.anno.pl.logo(
    adata,
    feature = "information", # choose from "count", "probability" and "information"
    conserved_color = None,
)

Use vp.anno.pl.motif_msa() to visualize sequence variants for each motif. The first upper sequence is the consensus motif in msa mode or the reference motif in pairwise mode. Bases same as the reference will be ignore and substitution/insertion/deletion will be shown. By default, the inserted bases will not show on the plot, set show_ins_bases to True to show the inserted bases.

Use vp.anno.pl.motif_msa() to visualize sequence variants for each motif. The top sequence represents the consensus motif in msa mode and the reference motif in pairwise mode. Bases identical to the reference are omitted from the display, while substitutions, insertions, and deletions are highlighted. By default, inserted bases are hidden to reduce visual clutter; set show_ins_bases=True to display the inserted sequences.

vp.anno.pl.motif_msa(
    adata,
    show_ins_bases=True,
)

[2026-06-11 21:10:05] You are using alignment data generated in msa mode.

Alternatively, you can specify a reference to perform pairwise alignment. Set reference=<motif_id> (an integer motif ID) to use the corresponding motif sequence as the reference, or provide a sequence string directly as a custom reference.

vp.anno.tl.motif_msa(
    adata,
    reference = 4
)

[2026-06-11 21:10:07] Reference specified; performing pairwise alignments of 24 motifs against 'AGGGATATGTGACTATGATAGTTATGTCCCTTACAGGATTTCACATATCCCTATTCTCAGACAGGAAT'
[2026-06-11 21:10:07] Aligned 24 motifs against reference 'AGGGATATGTGACTATGATAGTTATGTCCCTTACAGGATTTCACATATCCCTATTCTCAGACAGGAAT'. Found 74 variants.

AnnData object with n_obs × n_vars = 95 × 24
    obs: 'length', 'copy_number', 'score', 'sample', 'ancestry', 'unique_group'
    var: 'motif', 'motif_length', 'copy_number', 'label'
    uns: 'motif_array', 'orientation_array', 'sequence', 'motif_msa'
    varp: 'motif_distance', 'rc_motif_distance'

vp.anno.pl.motif_msa(
    adata,
    show_ins_bases=True,
)

[2026-06-11 21:10:07] You are using alignment data generated in pairwise mode.

Motif Abundance#

Tandem repeats often exhibit substantial sequence variation at the motif level across populations. Beyond differences in total copy number, individual alleles may vary in motif composition/abundance, and organization, reflecting diverse mutational and evolutionary processes such as motif duplication, substitution, and local rearrangement. Population-level motif abundance analysis aims to characterize these patterns by quantifying the relative abundance of distinct motifs across samples and identifying shared or divergent repeat architectures among individuals.

vp.anno.pl.motif_abundance_heatmap(
    adata,
    cluster_rows = True,
    cluster_cols = True,
    standard_scale = "obs", # use Max-Min normalization on "obs" row or "var" column， or set to None to disable normalization
    row_annotation = "ancestry",
    font = dict(size=10), # adjust to decrease the figure size
)

Except for heatmap visualization, motif composition patterns can also be quantitatively summarized using principal component analysis (PCA). By representing each tandem repeat allele as a motif composition vector, PCA projects population-level variation into a low-dimensional space, allowing major axes of motif composition diversity to be identified and visualized. Samples with similar motif usage patterns tend to cluster together in the embedding space, while distinct haplotypes or repeat architectures may separate along principal components. This provides a compact and interpretable representation of motif composition variation across populations.

vp.anno.tl.motif_abundance_pca(
    adata,
    n_components = 10,
)

[2026-06-11 21:11:02] PCA on motif abundance: 10 components. Explained variance: 36.0%, 18.0%, 13.9%, 8.6%, 6.1%, 4.9%, 2.8%, 2.5%, 1.9%, 1.2%

AnnData object with n_obs × n_vars = 95 × 24
    obs: 'length', 'copy_number', 'score', 'sample', 'ancestry', 'unique_group'
    var: 'motif', 'motif_length', 'copy_number', 'label'
    uns: 'motif_array', 'orientation_array', 'sequence', 'motif_msa', 'motif_abundance_pca'
    obsm: 'X_motif_abundance_pca'
    varm: 'motif_abundance_PCs'
    varp: 'motif_distance', 'rc_motif_distance'

vp.anno.pl.motif_abundance_pca_variance(adata)

vp.anno.pl.motif_abundance_pca(
    adata,
    shape_by = "ancestry",
    color_by = "copy_number",
    components = (1, 2), # choose which principal components to plot
)

Haplotype Clustering#

vp.anno.tl.sample_msa() performs a multiple sequence alignment (MSA) of motif arrays across samples. It first builds a substitution matrix from pairwise motif distances and lengths, where similar motifs receive higher alignment scores.

By default, the aligned motif arrays, orientation arrays, and consensus sequence are stored in .uns["aligned_motif_array"], .uns["aligned_orientation_array"] and .uns["aligned_consensus"].

vp.anno.tl.sample_msa(
    adata,
    match_score = 2,
    mismatch_penalty = -3,
    gap_open_penalty = -5,
    gap_extend_penalty = -1,
)

[2026-06-11 21:11:03] Aligning 95 samples: 88 unique sequences (7 duplicates removed).
[2026-06-11 21:11:07] Refinement converged after 1 iteration(s).
[2026-06-11 21:11:07] Aligned 95 samples (88 unique). Original lengths: [17, 10, 32, 28, 45, 27, 27, 10, 26, 24, 29, 18, 10, 7, 2, 18, 17, 31, 33, 2, 19, 30, 10, 7, 34, 24, 16, 13, 5, 16, 44, 24, 20, 2, 43, 4, 10, 18, 29, 27, 27, 39, 38, 8, 44, 25, 9, 26, 51, 54, 7, 26, 35, 5, 34, 6, 3, 20, 6, 4, 4, 13, 55, 54, 19, 6, 7, 9, 17, 45, 31, 11, 11, 7, 16, 7, 5, 16, 3, 8, 5, 12, 17, 9, 15, 5, 3, 4, 3, 9, 11, 58, 6, 17, 8]. Aligned length: 58.

AnnData object with n_obs × n_vars = 95 × 24
    obs: 'length', 'copy_number', 'score', 'sample', 'ancestry', 'unique_group'
    var: 'motif', 'motif_length', 'copy_number', 'label'
    uns: 'motif_array', 'orientation_array', 'sequence', 'motif_msa', 'motif_abundance_pca', 'aligned_motif_array', 'aligned_orientation_array', 'aligned_consensus', 'aligned_consensus_orientation'
    obsm: 'X_motif_abundance_pca'
    varm: 'motif_abundance_PCs'
    varp: 'motif_distance', 'rc_motif_distance'

Show aligned motif arrays by setting feature = "aligned_motif" in waterfall() plot.

vp.anno.pl.waterfall(
    adata,
    feature = "aligned_motif",
    colormap = "rainbow",
    sample_order = sample_order,
    margin = dict(l=120),
    track_name_dx = -0.01,
    font = dict(size=8),
)

[2026-06-11 21:11:07] Number of unique labels (24) exceeds colormap size (11), the overflow labels will be colored in #0c0c0c

vp.anno.tl.haplotype_neighbor() builds a fused k-nearest-neighbour (KNN) graph from aligned motif arrays for downstream haplotype clustering. It computes one or more pairwise distance matrices:

structural distance: distance of motif organization, including higher-order repeats and other shared substructures, between two samples
composition distance: the motif composition difference between two samples
CNV distance: the copy number difference between two samples

And then it converts each into an adaptive kNN similarity graph using a Gaussian kernel with a per-node local scale. The individual graphs are symmetrised and then fused by element-wise aggregation (max (default), mean, or min); the resulting connectivities are stored in .obsp["haplotype_connectivities"].

vp.anno.tl.haplotype_leiden() clusters samples into haplotypes by running the Leiden community-detection algorithm on the fused graph produced by vp.anno.tl.haplotype_neighbor(). After obtaining cluster labels, it constructs a consensus motif sequence for each haplotype by taking the most frequent non-gap motif at every alignment position among cluster members. The haplotype labels (e.g., H1, H2) are stored in .obs["haplotype"], and the consensus sequences are saved in .uns["haplotype_consensus"].

# select from "structural", "composition", "cnv", or you can pass self-defined functions
vp.anno.tl.haplotype_neighbor(adata, metrics=["structural", "composition"])

[2026-06-11 21:11:10] 95 samples collapsed to 88 unique sequences. (7 duplicates removed)
[2026-06-11 21:11:10] Computed built-in metric 'structural'
[2026-06-11 21:11:10] Computed built-in metric 'composition'
[2026-06-11 21:11:10] Stored structural distance matrix in obsp['haplotype_structural_distance']
[2026-06-11 21:11:10] Stored composition distance matrix in obsp['haplotype_composition_distance']
[2026-06-11 21:11:10] Building adaptive kNN graphs (k=9, n=88)
[2026-06-11 21:11:10] Built kNN graph for 'structural': 525 undirected edges
[2026-06-11 21:11:10] Built kNN graph for 'composition': 558 undirected edges
[2026-06-11 21:11:10] Fusing 2 graphs using 'max'
[2026-06-11 21:11:10] Fused graph: 828 undirected edges
[2026-06-11 21:11:10] Haplotype neighbour graph ready: 88 unique nodes, 828 edges, 1 component(s)

AnnData object with n_obs × n_vars = 95 × 24
    obs: 'length', 'copy_number', 'score', 'sample', 'ancestry', 'unique_group'
    var: 'motif', 'motif_length', 'copy_number', 'label'
    uns: 'motif_array', 'orientation_array', 'sequence', 'motif_msa', 'motif_abundance_pca', 'aligned_motif_array', 'aligned_orientation_array', 'aligned_consensus', 'aligned_consensus_orientation', 'haplotype'
    obsm: 'X_motif_abundance_pca'
    varm: 'motif_abundance_PCs'
    obsp: 'haplotype_structural_distance', 'haplotype_composition_distance', 'haplotype_connectivities'
    varp: 'motif_distance', 'rc_motif_distance'

You can use vp.anno.tl.haplotype_leiden_res_scan to identify the optimal resolution parameter that yields the highest modularity score.

best_res = vp.anno.tl.haplotype_leiden_res_scan(adata)
vp.anno.pl.haplotype_leiden_res_scan(adata)

[2026-06-11 21:11:19] Resolution scan complete: 20 resolutions → best modularity 0.437 at res=1.00 (k=4)

vp.anno.tl.haplotype_leiden(
    adata,
    resolution = best_res
)

[2026-06-11 21:11:19] Running Leiden clustering (resolution=1.00, random_state=0)
[2026-06-11 21:11:19] Leiden found 4 cluster(s)
[2026-06-11 21:11:19] Sorted haplotypes by sample_size
[2026-06-11 21:11:19] Built consensus for 4 haplotype(s)
[2026-06-11 21:11:19] Haplotype assignment complete: 95 sample(s) → 4 cluster(s). Sizes: {'H1': 33, 'H2': 26, 'H3': 22, 'H4': 14}

AnnData object with n_obs × n_vars = 95 × 24
    obs: 'length', 'copy_number', 'score', 'sample', 'ancestry', 'unique_group', 'haplotype'
    var: 'motif', 'motif_length', 'copy_number', 'label'
    uns: 'motif_array', 'orientation_array', 'sequence', 'motif_msa', 'motif_abundance_pca', 'aligned_motif_array', 'aligned_orientation_array', 'aligned_consensus', 'aligned_consensus_orientation', 'haplotype', 'haplotype_evaluation', 'haplotype_consensus'
    obsm: 'X_motif_abundance_pca'
    varm: 'motif_abundance_PCs'
    obsp: 'haplotype_structural_distance', 'haplotype_composition_distance', 'haplotype_connectivities'
    varp: 'motif_distance', 'rc_motif_distance'

vp.anno.pl.haplotype_distance_heatmap() visualises a pairwise distance matrix produced by vp.anno.tl.haplotype_neighbor() as an interactive heatmap. It retrieves the requested (structural, cnv, or composition) distance matrix from .obsp, optionally reorders rows and columns by haplotype label so that members of the same haplotype appear adjacent.

vp.anno.pl.haplotype_distance_heatmap(
    adata,
    metric = "structural", # choose from "structural", "composition" or "cnv"
    reorder = True,
    vmax = 0.5,
    font = dict(size=4),
)

vp.anno.pl.haplotype_distance_heatmap(
    adata,
    metric = "structural", # choose from "structural", "composition" or "cnv"
    cluster = True,
    vmax = 0.5,
    font = dict(size=4),
)

Then let’s group the samples by haplotype and mark the haplotype groups on waterfall plot.

adata_plot = adata.copy() #[adata.obs["ancestry"].isin(["AFR", "AMR", "EAS"])].copy()

# sort haplotype
adata_plot.obs["_hap_num"] = (
    adata_plot.obs["haplotype"]
    .astype(str)
    .str.extract(r"(\d+)")
    .astype(int)
)
# sort samples by haplotype and copy number
sample_order = (
    adata_plot.obs
    .sort_values(["_hap_num", "copy_number"])
    .index
    .tolist()
)

# plot
vp.anno.pl.waterfall(
    adata_plot,
    feature = "motif", # "aligned_motif" or "motif"
    colormap = "rainbow",
    sample_order = sample_order,
    row_annotation = ["ancestry", "haplotype"],
    font = dict(size=8),
    track_name_dx = -0.01,
)

[2026-06-11 21:11:19] Number of unique labels (24) exceeds colormap size (11), the overflow labels will be colored in #0c0c0c

Copy Number Variation#

By default, vp.anno.pl.copy_number_violin() displays the distribution of total copy numbers across different haplotypes, ancestries, or other categorical groups.

vp.anno.pl.copy_number_violin(
    adata,
    group_by = "haplotype",
    show_box = True,
    show_points = True,
)

vp.anno.pl.copy_number_violin(
    adata,
    group_by = "ancestry",
    show_box = True,
    show_points = True,
)

Set motif={motif_id} to display the copy number distribution of a specific motif across haplotypes, ancestries, or other categorical groups.

vp.anno.pl.copy_number_violin(
    adata,
    group_by = "haplotype",
    motif = 2,
    show_box = True,
    show_points = True,
)

vp.anno.pl.copy_number_violin(
    adata,
    group_by = "ancestry",
    motif = 2,
    show_box = True,
    show_points = True,
    show_counts = True,
)

You can also use vp.anno.pl.copy_number_stacked_violin() to visualize the copy number distribution of all motifs across different haplotypes, ancestries, or other categorical groups. By default, group_by=None, meaning that samples are displayed without grouping.

vp.anno.pl.copy_number_stacked_violin(
    adata,
    group_by = None,
    motifs = None,
    show_counts = True,
)

vp.anno.pl.copy_number_stacked_violin(
    adata,
    group_by = "haplotype",
    motifs = None,
)