Comparison and clustering¶

Comparing motifs with one another, merging similar ones, and building trees from them.

compare_motifs ¶

compare_motifs(
    motifs: Iterable[Motif],
    compare_to: int
    | str
    | Sequence[int | str]
    | None = None,
    qvalue: float = 0.1,
    min_overlap: int = 5,
    rc: bool = True,
    null: str = "empirical",
    bkg: Sequence[float]
    | Mapping[str, float]
    | None = None,
    matrix_out: str = "score",
    nthreads: int = 1,
    complete: bool = False,
)

Compare motifs.

Motifs of any alphabet can be compared. Each target is slid along the query, and with rc along the reverse complement of a DNA or RNA query too, and the alignment with the highest sum of column correlations is kept, among those overlapping at least min_overlap columns, or all of the shorter motif's columns when it is narrower than that. Column correlations are Pearson correlations of the probabilities. Each P-value comes from the exact distribution of that sum under a null model of columns, Bonferroni-corrected for the alignments tried. The "empirical" null uses the columns of all the motifs given, while the "parametric" null, for DNA and RNA only, uses a Dirichlet grid around bkg. Q-values are Benjamini-Hochberg adjusted within each query. P-values are worked out for each query's own columns and alignments, so a P-value matrix need not be symmetric even where the scores are.

By default only the overlap counts, so a short target that matches part of a long query can score as well as one that covers all of it. With complete, whole queries are scored instead, as TOMTOM does by default: each query column's correlations are measured against the median of its correlations with the null columns, so that a typical column scores zero and query columns hanging over the ends of the target count as typical, and the P-value of the best alignment is that of the best of all the alignments tried. This favours alignments covering more of the query, and scores then depend on which motif is the query too. Either way, score is the mean column correlation over the overlap.

Returns:

Type Description

ndarray | DataFrame

Without compare_to, an all-against-all numpy array with queries as rows and targets as columns, holding the mean column correlation, P-value or q-value as chosen by matrix_out. With compare_to, given as 0-based motif indices or as names, a DataFrame ordered by q-value with one row per match passing qvalue and the columns subject, subject_i, target, target_i, offset, strand, overlap, score, subject_consensus, target_consensus, pvalue and qvalue, where score is the mean column correlation.

Notes

Scoring alignments by summed Pearson correlations, with P-values from a null of database columns, is close to what TOMTOM does with its Pearson distance, so the two largely agree. The script scripts/compare_tomtom.py in the repository runs both on the same motifs. With 72 queries, 60 from CIS-BP and the 12 example motifs, against the 805 plant motifs of JASPAR 2024:

	`complete=True`, TOMTOM default	default, TOMTOM `-incomplete-scores`
log10 P-values, correlation over all 57,960 pairs	0.999	0.988
P-value rank correlation within queries, median	0.999	0.927
Each query's five best matches shared	95%	88%
Same best match	94%	78%
Strong matches on the same strand	99%	99%

Where the best match differs, the two choices are near ties, a tenth to a fifth of an order of magnitude apart.

Examples:

>>> import universalmotif as um
>>> motifs = um.example_motifs()[:4]
>>> um.compare_motifs(motifs).round(2).tolist()
[[1.0, 0.33, 0.11, 0.19], [0.33, 1.0, 0.46, 0.25], [0.11, 0.46, 1.0, 0.58], [0.19, 0.25, 0.58, 1.0]]

With compare_to, each match passing qvalue gets a row:

>>> wbox = um.example_motifs()[5]
>>> site = um.create_motif(["CTTTGACCAA", "GTTTGACCAA", "CTTTGACTAA"], name="site")
>>> hits = um.compare_motifs([wbox, site], compare_to="W-box")
>>> hits["target"].to_list(), hits["score"].round(2).to_list()
(['W-box', 'site'], [1.0, 0.65])

compare_columns ¶

compare_columns(
    x: Sequence[float],
    y: Sequence[float],
    method: str,
    bkg1: Sequence[float] | None = None,
    bkg2: Sequence[float] | None = None,
    nsites1: float = 100.0,
    nsites2: float = 100.0,
) -> float

Compare two columns with one of the comparison metrics: "PCC", "EUCL", "SW", "KL", "WEUCL", "ALLR", "BHAT", "HELL", "WPCC", "SEUCL", "MAN" or "ALLR_LL". The backgrounds are uniform by default and are used by the weighted metrics and ALLR, which also uses the site counts.

Examples:

>>> import universalmotif as um
>>> x = [0.7, 0.1, 0.1, 0.1]
>>> y = [0.6, 0.2, 0.1, 0.1]
>>> round(um.compare_columns(x, y, "PCC"), 3)
0.98
>>> round(um.compare_columns(x, y, "EUCL"), 3)
0.141

merge_motifs ¶

merge_motifs(
    motifs: Iterable[Motif],
    min_overlap: int = 5,
    rc: bool = True,
    new_name: str | None = None,
    weighted: bool = False,
    nthreads: int = 1,
) -> Motif

Merge motifs into one PPM.

Motifs of any alphabet can be merged. The motif with the highest total IC is the anchor. The others are aligned to it, with DNA and RNA motifs reverse complemented under rc where that fits better, and the aligned probabilities are averaged, weighted by nsites if weighted. Motifs that cannot overlap the anchor by min_overlap columns are left out, with a warning. The default name joins the merged motifs' names.

Examples:

>>> import universalmotif as um
>>> v1 = um.create_motif(["CTTTGACCAA", "GTTTGACCAA", "CTTTGACTAA"], name="v1")
>>> v2 = um.create_motif(["TTTGACCAAT", "TTTGACTAAT", "TTTGACCAAA"], name="v2")
>>> merged = um.merge_motifs([v1, v2], new_name="W-box merged")
>>> merged.name, merged.consensus, merged.type
('W-box merged', 'STTTGACYAAW', 'PPM')

merge_similar ¶

merge_similar(
    motifs: Iterable[Motif],
    qvalue: float = 0.01,
    null: str = "empirical",
    min_overlap: int = 5,
    rc: bool = True,
    bkg: Sequence[float]
    | Mapping[str, float]
    | None = None,
    weighted: bool = False,
    return_clusters: bool = False,
    nthreads: int = 1,
)

Cluster similar motifs and merge each cluster.

Motifs are linked when compare_motifs() gives them a q-value of at most qvalue in either direction, clusters are the connected groups, and each cluster of several motifs is merged with merge_motifs(). Since a single link is enough, motifs that are only alike in part, such as those of one family of transcription factors, can chain into one large cluster; cutting a motif_tree() at a fixed height avoids that.

Returns one motif per cluster, or with return_clusters a summary DataFrame of the input motifs with 0-based motif_i and cluster columns.

Examples:

>>> import universalmotif as um
>>> v1 = um.create_motif(["CTTTGACCAA", "GTTTGACCAA", "CTTTGACTAA"], name="v1")
>>> v2 = um.create_motif(["TTTGACCAAT", "TTTGACTAAT", "TTTGACCAAA"], name="v2")
>>> gbox = um.example_motifs()[3]
>>> [m.name for m in um.merge_similar([v1, v2, gbox])]
['merged_v1+v2', 'G-box']
>>> um.merge_similar([v1, v2, gbox], return_clusters=True)["cluster"].to_list()
[0, 0, 1]

motif_tree ¶

motif_tree(
    motifs,
    min_overlap: int = 6,
    rc: bool = True,
    method: str = "complete",
    labels: str | Sequence[str] = "name",
    nthreads: int = 1,
) -> MotifTree

Cluster motifs into a tree.

Motifs are compared with compare_motifs(), mean column correlations become distances as (1 - score) / 2, and the distances are clustered hierarchically. The agglomeration method is one of "complete", "average", "single", "ward.D", "ward.D2", "mcquitty", "median" or "centroid"; "ward.D" applies Ward's criterion to the distances as they are and "ward.D2" to the squared distances.

motifs can also be a square distance matrix, of which only the lower triangle is used. labels is either a list of labels or a motif attribute to label tips with, such as "name", "altname" or "family". Duplicate motif names get .1, .2 and so on, while the first of each is left as it is.

Examples:

>>> import universalmotif as um
>>> tree = um.motif_tree(um.example_motifs()[:6])
>>> tree
MotifTree(labels=['TATA-box', 'AP-1', 'CRE', 'G-box', 'E-box', 'W-box'], method='complete')
>>> [tree.labels[i] for i in tree.order]
['TATA-box', 'AP-1', 'W-box', 'E-box', 'CRE', 'G-box']
>>> tree.height.round(3).tolist()
[0.096, 0.211, 0.354, 0.384, 0.447]

MotifTree `dataclass` ¶

MotifTree(
    labels: list[str],
    merge: ndarray,
    height: ndarray,
    order: list[int],
    linkage: ndarray,
    newick: str,
    method: str = "complete",
    motifs: list[Motif] | None = None,
)

A clustering tree of motifs.

Attributes:

Name	Type	Description
`labels`	`list[str]`	Tip labels, in input order.
`merge`	`ndarray`	Merge steps, one per row. Row `i` joins two items, where `-k` is the `k`-th input and `k` is the cluster formed in row `k`, both counting from 1.
`height`	`ndarray`	Height of each merge.
`order`	`list[int]`	0-based input indices in dendrogram order.
`linkage`	`ndarray`	The same tree as a SciPy linkage matrix, for use with `scipy.cluster.hierarchy` functions such as `dendrogram` and `fcluster`.
`newick`	`str`	The tree in Newick format. Branch lengths are half the differences in merge height, written with up to 10 significant digits, and whitespace in tip labels becomes `_` while Newick punctuation becomes `-`.
`method`	`str`	Agglomeration method.
`motifs`	`list[Motif] \| None`	The motifs at the tips, for trees made from motifs.

Examples:

>>> import universalmotif as um
>>> dist = [[0.0, 0.1, 0.8], [0.1, 0.0, 0.7], [0.8, 0.7, 0.0]]
>>> tree = um.motif_tree(dist, labels=["a", "b", "c"])
>>> tree.merge.tolist(), tree.height.tolist()
([[-1, -2], [-3, 1]], [0.1, 0.8])
>>> tree.newick
'(c:0.4,(a:0.05,b:0.05):0.35);'

write_newick ¶

write_newick(file) -> None

Write the Newick string to a file.

Examples:

>>> import tempfile
>>> from pathlib import Path
>>> import universalmotif as um
>>> tree = um.motif_tree(um.example_motifs()[:3])
>>> with tempfile.TemporaryDirectory() as tmp:
...     path = Path(tmp) / "motifs.nwk"
...     tree.write_newick(path)
...     path.read_text().strip() == tree.newick
True

cutree ¶

cutree(
    k: int | None = None, h: float | None = None
) -> np.ndarray

Cut the tree into k groups, or at height h.

Returns one group per motif, numbered from 0 in the order the groups first appear among the motifs. Cutting by height needs increasing merge heights, which the centroid and median methods do not guarantee.

Examples:

>>> import universalmotif as um
>>> dist = [[0.0, 0.1, 0.8], [0.1, 0.0, 0.7], [0.8, 0.7, 0.0]]
>>> tree = um.motif_tree(dist, labels=["a", "b", "c"])
>>> tree.cutree(k=2).tolist()
[0, 0, 1]
>>> tree.cutree(h=0.05).tolist()
[0, 1, 2]

plot ¶

plot(
    layout: str = "rectangular",
    labels: bool = True,
    groups=None,
    tip_size=None,
    tip_logos: bool = False,
    branch_length: bool = False,
    legend: bool = True,
    colours=None,
    open_angle: float = 90,
    tip_logo_width: float | None = None,
    tip_logo_height: float = 0.9,
    text_size: float = 10,
    ax=None,
)

Draw the tree with matplotlib.

Parameters:

Name	Type	Description	Default
`layout`	`str`	`"rectangular"`, `"slanted"`, `"circular"` or `"fan"`, which is a circle with a gap of `open_angle` degrees.	`'rectangular'`
`labels`	`bool`	Label the tips with `self.labels`.	`True`
`groups`		Colour branches by group: a motif attribute such as `"family"`, one value per motif in the order of the tree's labels, or a mapping from labels to values. A branch takes a group's colour when every motif below it is in that group; `colours` maps groups to colours.	`None`
`tip_size`		Size the tips by one number per motif, given as for `groups`, or by a motif attribute, either `"icscore"` or one of `"pval"`, `"qval"` and `"eval"`, which are taken as their -log10.	`None`
`tip_logos`	`bool`	Draw each motif's logo at its tip, aligned as in `view_motifs()`, in the rectangular and slanted layouts only. `tip_logo_width` is the width of one position in tree units and `tip_logo_height` the logo height in tip spacings.	`False`
`branch_length`	`bool`	Draw branch lengths, which are half the merge heights as in the Newick string; by default tips are aligned, as a cladogram.	`False`
`ax`		Axes to draw on; a new figure by default.	`None`

Returns:

Type	Description
`The matplotlib axes.`

Examples:

>>> import matplotlib.pyplot as plt
>>> import universalmotif as um
>>> tree = um.motif_tree(um.example_motifs())
>>> ax = tree.plot(groups=tree.cutree(k=3), legend=False)
>>> type(ax).__name__
'Axes'
>>> plt.close("all")

Comparison and clustering¶

compare_motifs ¶

compare_columns ¶

merge_motifs ¶

merge_similar ¶

motif_tree ¶

MotifTree dataclass ¶

write_newick ¶

cutree ¶

plot ¶

MotifTree `dataclass` ¶