Interpretation

Tools that consume attribution scores and accessibility signal: GPU chromVAR (the AP-1 activity score behind Bolero-Score), motif scanning, TF-MoDISco seqlet discovery and embedding, finemo genome-wide hit calling, and seqlet→PWM matching via TomTom.

chromVAR

compute_deviations

compute_deviations(adata: AnnData, chunk_size: int = 10000, device: str = 'cuda') -> AnnData

Computes the deviation of motif matches from the background for each cell.

Parameters:

Name	Type	Description	Default
adata	AnnData	The input AnnData object containing the count matrix, background peaks, and motif match information.	required
chunk_size	int	The size of chunks to process the data in. Default is 10000. It is recommended to set this such that there's no GPU memory overflow. Although our implementation would tolerate that, it will be much slower if there's overflow.	10000
device	str	The device to use for computation. Can be either "cuda" for GPU or "cpu" for CPU. Default is "cuda".	'cuda'

Returns:

Name	Type	Description
dev	AnnData	The AnnData object containing the deviation values for each cell and motif match.

get_peak_bias

get_peak_bias(adata: AnnData, genome: Genome, base_order: str | None = None) -> None

Compute per-peak GC content and the overall background nucleotide frequency.

This is the bolero-native replacement for scPrinter's get_peak_bias: it fills adata.var["gc_content"] (used, together with per-peak read depth, to sample matched background peaks via :func:sample_bg_peaks) and adata.uns["bg_freq"] (the genome-wide-over-peaks A/C/G/T frequency used as the MOODS background for motif matching in :func:scan_peak_motifs).

Peak sequences are read directly from the genome FASTA (genome.fasta_path) with pyfaidx; only A/C/G/T are counted (N and other symbols are ignored).

Parameters:

Name	Type	Description	Default
adata	AnnData	Peak object; peaks are taken from adata.var (see :func:_parse_peak_regions).	required
genome	Genome	Genome providing fasta_path for the same assembly the peaks are on.	required
base_order	str	Base order for bg_freq. Defaults to bolero.pp.seq.DEFAULT_ONE_HOT_ORDER ("ACGT"), matching the order expected by MOODS.	None

Returns:

Type	Description
None	adata is modified in place (var["gc_content"], uns["bg_freq"]).

scan_peak_motifs

scan_peak_motifs(adata: AnnData, genome: Genome, motif_path: str | None = None, motif_db: str = 'JASPAR2024_CORE_vertebrates', bg: tuple[float, ...] | str | None = None, pvalue: float = 5e-05, pseudocount: float = 0.8, n_jobs: int = 16, mode: str = 'motifmatchr', verbose: bool = True) -> None

Match motifs against each peak and store the hit matrix on adata.

This is the bolero-native replacement for scPrinter's Motifs.chromvar_scan: it scans the peaks with :class:bolero.tl.motif.scan.Motifs (a MOODS motifmatchr port) and fills adata.varm["motif_match"] (a float32 peak-by-motif 0/1 matrix) and adata.uns["motif_name"] (the JASPAR matrix ids), the inputs :func:compute_deviations needs.

Parameters:

Name	Type	Description	Default
adata	AnnData	Peak object; peaks are taken from adata.var (see :func:_parse_peak_regions).	required
genome	Genome	Genome providing fasta_path for the same assembly the peaks are on.	required
motif_path	str or Path	JASPAR-format PFM file. If None, a cached copy of motif_db is fetched via :func:bolero.tl.motif.jaspar.get_jaspar_motif_file.	None
motif_db	str	JASPAR database key used when motif_path is None. Default "JASPAR2024_CORE_vertebrates".	'JASPAR2024_CORE_vertebrates'
bg	tuple of float or str	Background nucleotide frequency (A/C/G/T). If None, uses adata.uns["bg_freq"] when available (from :func:get_peak_bias), else "even".	None
pvalue	float	Motif match p-value threshold. Default 5e-5 (motifmatchr default).	5e-05
pseudocount	float	Pseudocount for the PFMs. Default 0.8 (motifmatchr default).	0.8
n_jobs	int	Number of processes for scanning. Default 16.	16
mode	(motifmatchr, moods)	Motif-matching mode. Default "motifmatchr".	"motifmatchr"
verbose	bool	Whether to show a scanning progress bar. Default True.	True

Returns:

Type	Description
None	adata is modified in place (varm["motif_match"], uns["motif_name"]).

Motif scanning

Motifs

A class for motif matching based on MOODS

Parameters:

Name	Type	Description	Default
ref_path_motif	str | Path	Path to the motif file, in JASPAR format	required
ref_path_fa	str | Path	Path to the reference genome fasta file. An eazy way would be passing, the genome.fetch_fa() function	required
bg	Literal['even', 'genome'] | tuple	Background nucleotide frequency, by default 'even' stands for equal frequency of A, C, G, T 'genome' stands for the frequency of A, C, G, T in the reference genome When passing a tuple of size 4, they would be used as background frequency	'even'
pseudocount	float	Pseudocount for the motif, by default 0.8 (same as motifmatchR)	0.8
pvalue	float	P-value threshold for motif matching, by default 5e-5	5e-05
n_jobs	int	Number of cores to use, by default 32	32
split_tf	bool	Whether to split the TFs by comma (for cistarget motif collection), by default True	True
mode	Literal['motifmatchr', 'moods']	The mode of the motif matching, by default 'motifmatch	'motifmatchr'

prep_scanner

prep_scanner(tf_genes: list[str] | None = None, pseudocount: float = 0.8, pvalue: float = 5e-05, window: int = 7)

Prepare the MOODS scanner for motif matching

Parameters:

Name	Type	Description	Default
tf_genes	list[str] | None	List of TFs to be used for motif matching, by default None, which means all TFs. Note that if you pass motif_name_func and/or split_tf, you should specify TFs by their names after being processed by those (so likely the true TF names).	None
pseudocount	float	Pseudocount for the motif, by default 0.8 (same as motifmatchR)	0.8
pvalue	float	P-value threshold for motif matching, by default 5e-5	5e-05
window	int	Window size for motif matching, by default 7 (same as motifmatchR), passed to MOODS	7

Returns:

Name	Type	Description
scanner	Scanner

scan

scan(region_bed: DataFrame, verbose: bool = True, mutation_table: DataFrame = None, use_genotype: Literal['ref', 'alt'] = 'ref')

Perform motif scanning on the given peaks and store the results in the AnnData object.

Parameters:

Name	Type	Description	Default
region_bed	DataFrame	The region bed file containing the region information.	required
verbose	bool	Whether to display a progress bar during the motif scanning process. Default is True.	True
mutation_table	DataFrame	Table containing mutation information. Position is 1-based. Columns: ["Chromosome", "Position", "Ref", "Alt"] Currently only supports SNV mutations. Mutation table index is the same as the region bed index. If provided, the DNA sequence will be mutated in the SNVs before motif matching.	None
use_genotype	Literal['ref', 'alt']	Whether to use the reference or alternative genotype for the mutation. Default is "ref".	'ref'

Returns:

Type	Description
DataFrame	A dataframe containing the motif match information for each region. The index is the region names, the columns are the motif names.

collect_child_process

collect_child_process(p_list, verbose, bar, motifs_length, name2id, maps, break_on_min_jobs)

Collect completed processes from the process list and process their motif mapping results.

Parameters:

Name	Type	Description	Default
p_list	list	A list of Process objects representing the child processes.	required
verbose	bool	A flag indicating whether to display a progress bar.	required
bar	tqdm	A progress bar object for displaying the progress.	required
motifs_length	list	A list containing the lengths of the motifs.	required
name2id	dict	A dictionary mapping motif names to their indices.	required
maps	dict	A dict to store the results of the motif scanning process.	required
break_on_min_jobs	bool	A flag indicating whether to break the loop when the number of remaining processes is less than or equal to the number of workers.	required

Returns:

Type	Description
None	The function modifies the maps list in-place by adding the results of the completed processes.

scan_motif

scan_motif(peaks_iter, clean: bool = False, verbose: bool = False, split_tfs: bool = True, strand: bool = True, mutation_table: DataFrame = None, use_genotype: Literal['ref', 'alt'] = 'ref')

Scan motifs in the peaks iterator.

Parameters:

Name	Type	Description	Default
peaks_iter	iterable	An iterable yielding a list [chrom, start, end], e.g., a pandas dataframe with three columns. This function currently only supports a specific format of regions.	required
clean	bool	Whether to clean the output. If True, overlapping motif hits of the same TF will be merged into one hit. Default is False.	False
verbose	bool	Whether to display a progress bar. Default is False.	False
split_tfs	bool	Whether to split the output by transcription factors. Default is True.	True
strand	bool	Whether to consider the strand information. Default is True.	True
mutation_table	DataFrame	Table containing mutation information. Position is 1-based. Columns: ["Chromosome", "Position", "Ref", "Alt"] Currently only supports SNV mutations. If provided, the DNA sequence will be mutated in the SNVs before motif matching.	None
use_genotype	Literal['ref', 'alt']	Whether to use the reference or alternative genotype for the mutation. Default is "ref".	'ref'

Returns:

Type	Description
list or ndarray	The output of the motif scanning process. The format depends on the values of the clean, split_tfs, and strand parameters.

TF-MoDISco & hit calling

ModiscoResults

Represents the results of a Modisco analysis.

Attributes:

Name	Type	Description
output_dir (str)	The output directory of the analysis.	_pos_patterns (list[ModiscoPattern]): The positive patterns identified in the analysis. _neg_patterns (list[ModiscoPattern]): The negative patterns identified in the analysis. _dna_one_hot (np.ndarray): The one-hot encoded DNA sequences. _attr (np.ndarray): The attribute scores. _attr_1d (np.ndarray): The 1D attribute scores. _tfbs (np.ndarray): The TFBS scores. _region (pd.DataFrame): The region information. _hits (pd.DataFrame): The motif hits.

neg_patterns property

neg_patterns

Modisco patterns with negative attribution scores.

neg_seqlets property

neg_seqlets

Get the representative seqlet for each negative pattern.

pos_patterns property

pos_patterns

Modisco patterns with positive attribution scores.

pos_seqlets property

pos_seqlets

Get the representative seqlet for each positive pattern.

dna_one_hot property

dna_one_hot: ndarray

Get the one-hot encoded DNA sequences.

Returns:

Type	Description
np.ndarray: The one-hot encoded DNA sequences.

attr property

attr: ndarray

Get the attribute scores.

Returns:

Type	Description
np.ndarray: The attribute scores.

attr_1d property

attr_1d: ndarray

Get the 1D attribute scores.

Returns:

Type	Description
np.ndarray: The 1D attribute scores.

tfbs property

tfbs: ndarray

Get the TFBS scores.

Returns:

Type	Description
np.ndarray: The TFBS scores.

region property

region: DataFrame

Get the region information.

Returns:

Type	Description
pd.DataFrame: The region information.

get_motif_hits

get_motif_hits(tfbs=False, attr_1d=False, slop=30) -> pd.DataFrame

Get the motif hits.

Returns:

Type	Description
pd.DataFrame: The motif hits.

modisco_seqlets_embedding

modisco_seqlets_embedding(seqlets, nearest_neighbors_to_compute=500, sign='pos', min_overlap_while_sliding=0.7, affmat_correlation_threshold=0.15, tsne_perplexity=30.0, corr_filter=False, n_jobs=1, tsne=True)

Calculate affinity matrix between seqlets and embed them in 2D space using t-SNE.

Adapted from modisco.tfmodisco.seqlets_to_patterns function. https://github.com/jmschrei/tfmodisco-lite/blob/main/modiscolite/tfmodisco.py

Parameters:

Name	Type	Description	Default
seqlets	list of modisco.core.Seqlet	List of seqlets to embed.	required
nearest_neighbors_to_compute	int	Number of nearest neighbors to compute for cosine similarity as an approximate start, by default 500.	500
sign	str	Sign of the cosine similarity, by default "pos".	'pos'
min_overlap_while_sliding	float	Minimum overlap while sliding, by default 0.7.	0.7
affmat_correlation_threshold	float	Threshold for filtering by correlation, by default 0.15.	0.15
tsne_perplexity	float	Perplexity for t-SNE, by default 10.0.	30.0
corr_filter	bool	Whether to filter seqlet by correlation, by default False.	False
n_jobs	int	Number of parallel jobs, by default 1.	1
tsne	bool	Whether to compute the 2D t-SNE embedding, by default True.	True

Returns:

Name	Type	Description
dist_mat	ndarray	Distance matrix between seqlets.
embeddings	DataFrame	Embeddings of seqlets in 2D space.
tsne	TSNE	t-SNE object.

Finemo

plot_pattern

plot_pattern(name, **kwargs)

Plot a modisco pattern by name.

run

run(all_sequences, all_contribs, all_attr_region, max_round=3)

Run the finemo pipeline.

SeqletTomtom

Class to perform motif comparison using Tomtom.

adjust_seqlets_idx_

adjust_seqlets_idx_(seqlets: list[Seqlet], chunk_start: int)

Adjust seqlet indices to account for the chunk start. This is necessary to maintain global indexing across chunks.

tomtom

tomtom(seqlets: list[Seqlet], blacklist: list[str] = MOTIF_BLACKLIST) -> xr.DataArray

Perform Tomtom motif comparison on seqlets.

simple_tomtom

simple_tomtom(attr_scores: ndarray, score_names: list[str] = None, blacklist: list[str] = MOTIF_BLACKLIST) -> xr.DataArray

Perform Tomtom motif comparison on attribute scores.

Parameters:

Name	Type	Description	Default
attr_scores	ndarray	Attribute scores. Shape (n, 4, seq_len).	required
score_names	list[str]	Score names. Default is None.	None

Returns:

Type	Description
DataArray	Tomtom results. Shape (n, motif_rank, value_type). value_type: ["p_values", "scores", "offsets", "overlaps", "strands", "idxs"].

annotate_seqlet_ds

annotate_seqlet_ds(seqlet_ds: Dataset, seqlet_info: DataFrame, flank_seqlet_size: int = 2) -> pd.DataFrame

Annotate a tengermeme extracted seqlet dataset with motif information.

Parameters:

Name	Type	Description	Default
seqlet_ds	Dataset	Seqlet dataset.	required
seqlet_info	DataFrame	Seqlet information.	required
flank_seqlet_size	int	Flank size of the seqlet cluster. Default is 2.	2

run

run(one_hot: str | ndarray, hypothetical_contribs: str | ndarray, output_dir: str, regions: DataFrame | None = None, pseudobulk_ids: Series | ndarray | None = None)

Run the entire Modisco Tomtom pipeline.

First using modisco extract seqlets from the dna and attribution scores, then perform Tomtom motif comparison, and finally save the results to a Zarr file.

Zarr file structure: Dimensions: - seqlet, motif_id, base, position, value_type, motif_rank Coordinates: - seqlet: Seqlet name. If regions are provided, the seqlet name will be the genomic region. - motif_id: Motif ID. - base: DNA base A, C, G, T. - value_type: Type of value from Tomtom ["p_values", "scores", "offsets", "overlaps", "strands", "idxs"]. Data variables: - attr_region: relative attr region coords of each seqlets: {input_region_idx}{rel_start}{rel_end} (seqlet,) - motif_name: motif_name (motif_id,) - seqlets_score: attr score for each seqlet (seqlet, base, position) - seqlets_seq: DNA one hot for each seqlet (seqlet, base, position) - seqlets_sign: sign of the seqlet (seqlet,) - seqlets_tomtom: Tomtom results for each seqlet (seqlet, motif_rank, value_type)

Parameters:

Name	Type	Description	Default
one_hot	str | ndarray	Path to .npz file or numpy array of one-hot encoded DNA sequences. Shape (n, seq_len, 4).	required
hypothetical_contribs	str | ndarray	Path to .npz file or numpy array of hypothetical contribution scores. Shape (n, seq_len, 4).	required
output_dir	str	Path to the output directory.	required
regions	DataFrame | None	DataFrame of regions. If provided, the seqlet index will be annotated with the regions. First three columns should be "chrom", "start", "end".	None

Returns:

Type	Description
None

AggregateModiscoResults

annotate_patterns

annotate_patterns()

Annotate patterns using TOMTOM.

add_pattern_corr

add_pattern_corr(pattern_annot)

For each TOMTOM hits, further calculate correlation with motif PWM. tomtom implementation doesn't support corr, so calculate it here as a post filtering step.

run

run()

Run the aggregation and save the results to a file.