Examples#

Four runnable scripts in examples/. Each is self-contained: the only data is the 250k control repertoire bundled with the package, and the only import is seqtree plus the standard library.

python examples/01_gapped_search.py
python examples/02_sequence_dendrogram.py
python examples/03_indel_positions.py
python examples/04_island_profile.py

1. Gapped search, with a cutoff you can defend#

examples/01_gapped_search.py plants twelve relatives of CASSLGQAYEQYF in a slice of the control – some substituted, some a residue shorter, some a residue longer – and looks for them.

A Hamming ball can only ever reach references of the query’s own length:

Hamming ball (max_penalty 40): 10 hits, 0 of them length-different
gap-block ball  (max_penalty 60): 44 hits, 12 of them length-different

An indel costs gap_open = 28, so gapped hits never top a score-sorted list; the script surfaces them separately. Then it derives the cutoff instead of choosing one:

N = 2,012 references, M = 250,000 control. Rule of three: no E* below 3N/M = 0.024.
  E* = 1.0    theta per query:  CASSLG..=32  CSARDI..=50
  E* = 0.05   theta per query:  CASSLG..=16  CSARDI..=50
  E* = 0.001  theta per query:  CASSLG..=unreachable  CSARDI..=unreachable

     E*  theta   hits  planted  background  precision
    1.0     32      5        4           1       0.80
   0.05     16      2        2           0       1.00

Two things to take from this. The two queries get different cutoffs: CASSLGQAYEQYF sits in a dense, near-germline neighbourhood where the same score buys far more chance neighbours, so it is held to a stricter bar (32 against 50 at E* = 1.0). And E* = 0.001 comes back -1 – with 250k control sequences no cutoff can certify a bar below 3N/M. That is the rule of three, not a bug.

2. A dendrogram over an island#

examples/02_sequence_dendrogram.py builds a family with internal structure – a substitution branch and a branch defined by one deletion – and clusters it by average linkage on the gap-block distance. Thirty lines of stdlib; seqtree has no runtime dependencies and the example keeps it that way. Output is Newick plus an ASCII rendering.

merge                                           height
CASSLGAYEQYF        + CASSVGAYEQYF               6.0  -----
CASSWGQAYEQYF       + CASSWGQSYEQYF              6.0  -----
...
(...)               + (...)                     52.3  ------------------------------------

The deletion branch is the last to join: the root sits at 52.3 (average linkage, so the mean over cross-branch pairs), of which 28 is the irreducible gap_open every cross-length pair must pay — that floor is what pushes the one-deletion branch to the top. Then the tree is cut where the control says to cut it:

N = 20,012 (the annotation set), M = 250,000 (the control)
  E* = 1.0   mean theta 24.4   8/11 merges survive the cut   -> 4 clusters
  E* = 0.1   mean theta 15.1   6/11 merges survive the cut   -> 6 clusters

Warning

Do not build one of these over a whole repertoire: a dense condensed distance matrix over the ~90k CDR3 of a VDJdb shortlist is 32 GB. Cluster first, draw trees inside islands.

And do not read anything into merges above the calibrated cutoff. Distinct sequence islands for one epitope share no motif that similarity can find – only ~0.5% of cross-island co-specific pairs share even a rare central 4-mer, and none share a 6-mer. Above theta the tree is arithmetic, not biology.

3. Where the indel goes#

examples/03_indel_positions.py plants a deletion at a known position, adds k substitutions, and asks each rule to find it again.

=== deletion planted in the CORE (positions 4 .. L-5): exact recovery ===
  rule                         k=0 subs   k=1 subs   k=2 subs   k=3 subs
  fixed centre                    27.3%      25.0%      23.4%      24.4%
  central prior lam=21            28.3%      25.8%      24.0%      25.0%
  flat (score alone)             100.0%      92.4%      81.8%      76.5%

This is easy to misread. When two sequences really are one deletion plus a little noise, the score finds the gap. So why constrain it at all? Because that ground truth is a construction. Against TCR-pMHC crystal structures – where “correct” means “the residues that actually superpose” – the score agrees about a tenth of the time and a central prior 30–42%. In retrieval at a matched false-positive rate, letting the score choose cuts precision from 0.414 to 0.176.

The prior is an assumption, and the script prices it. Move the deletion to the anchors, where V(D)J recombination does not put it, and the prior is confidently wrong while the score is unbothered:

=== deletion planted at the EDGE (within 2 of an anchor), k = 1 sub ===
  rule                          exact   within 1
  fixed centre                   0.0%       0.0%
  central prior lam=21           0.0%       0.0%
  flat (score alone)            89.7%      98.6%

Finally, the marginal: where does each rule land on ordinary, unrelated pairs?

rule                        i/L mean     sd   histogram of i/L (0 = Cys end, 1 = Phe end)
fixed centre                   0.496  0.026   ....##....
central prior lam=21           0.496  0.027   ....##....
flat (score alone)             0.515  0.260   #+########

Left free, the score puts the block within two residues of an anchor on 25.9% of unrelated pairs. That is where a spurious low score hides: shift the whole sequence by one and the conserved CASS / EQYF frame still lines up. It is the same mechanism by which unrestricted affine alignment undercuts the block score by a median of 106 penalty units on unrelated pairs. It does not find a better alignment. It manufactures one.

4. An island profile, and the cutoff that decides whether to build one#

examples/04_island_profile.py fits an IslandProfile to a planted island and scores twenty held-out members against all 250,000 control junctions — twice, once per scorer. The point is not which wins on one island. It is that the E-value decides which question you are asking:

regime                         N     k        FPR   note
per-epitope island            90   138    0.0552%   E*=0.05
repertoire annotation     20,000     3    0.0012%   E*=0.05 unreachable (3N/M=0.240); using E*=0.240

k = floor(e_target · M / N) is how many control neighbours the calibrated cutoff may admit, so the false-positive rate is k / M — and it moves by a factor of fifty depending on whether you are building islands inside one epitope group or annotating a whole repertoire against them. At the first cutoff the two scorers are indistinguishable; at the second, over 108 real VDJdb islands, the profile recovers 48.5 % of held-out members against 37.6 %.

The script says so out loud rather than reporting its own twenty-member table as a result. Note also that E* = 0.05 is simply unreachable at N = 20,000: the rule of three certifies nothing below 3N/M, and thetas_from_scores() returns -1 instead of a cutoff it cannot defend.

It also demonstrates score_matrix(): min-over-members is a row minimum of the 250,000 × 40 gap-block matrix, which the C++ kernel produces in milliseconds.

See also#