Add tcrblosum support to TCRdist

[felixpetschko]

So far, the TCRdist metric used a distance matrix derived from the blosum62 substitution matrix. This PR extends TCRdistDistanceCalculator with a new base_matrix="tcrblosum" option alongside the existing default blosum62 behavior. This way, distance matrices based on the tcrblosum substitution matrices (different matrices for alpha and beta chain) are used for the TCRdist metric calculation.

I try to illustrate how I derived the tcrblosum based distance matrices in this google colab notebook.

The usage of the tcrblosum matrices was already discussed in #591.

[codecov]

Codecov Report

❌ Patch coverage is 20.00000% with 28 lines in your changes missing coverage. Please review.
✅ Project coverage is 19.60%. Comparing base (38f05cd) to head (86f4514).

Files with missing lines	Patch %	Lines
src/scirpy/ir_dist/metrics.py	18.75%	26 Missing ⚠️
src/scirpy/ir_dist/__init__.py	33.33%	2 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #685      +/-   ##
==========================================
- Coverage   19.75%   19.60%   -0.16%     
==========================================
  Files          51       51              
  Lines        4581     4607      +26     
==========================================
- Hits          905      903       -2     
- Misses       3676     3704      +28     

Files with missing lines	Coverage Δ
src/scirpy/ir_dist/__init__.py	21.21% <33.33%> (ø)
src/scirpy/ir_dist/metrics.py	13.58% <18.75%> (-0.99%)	⬇️

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[grst]

Implementation-wise this looks great!

What's still missing is

• changelog update
• Documentation-update of the user-facing (pp.ir_dist) method. Probably best to add a new metric tcrblosum or tcrdist_tcrblosum.
• Reference to the TCRblosum paper in the documentation
• Maybe tutorial update?

[grst]

Did you make any modifications to the original TCRblosum matrix? e.g. the cap at 4?

[felixpetschko]

First of all, you always have this conversion from the substitution matrix (blosum62, tcrblosum) to a distance matrix.
In the tcrblosum paper , they just menttion "Firstly, we transformed the tcrBLOSUM similarity matrix into a distance matrix according to the rules of TCRdist [12]".

In the original TCRdist paper they write:
"The mismatch distance is defined based on the BLOSUM62 (ref. 37) substitution matrix as follows: distance (a, a)=0; distance (a, b)=min (4, 4-BLOSUM62 (a, b)), where 4 is 1 unit greater than the most favourable BLOSUM62 score for a mismatch, and a and b are amino acids".

Now there would be two options:

1. also use the constant cap 4 like in the original tcrdist paper
2. use a cap that is one unit greater than the most favourable score

However, the most favourable score would be 1 for tcrblosum_alpha and 2 for tcrblosum_beta which would result in a quite low cap (see matrices in notebook). Therefore I just went with the fixed cap of 4 for the transformation formula from substitution matrix to distance matrix of the original tcrdist implementation:
distance(a, a) = 0 and distance(a, b) = min(4, 4 - score)

Unfortunately, I wasn't able to find out how they did it exactly in the tcrblosum paper. Another thing to consider is, that the distance values for alpha and beta chain might be compared later during the clonotype clustering. Therefore having two different caps might cause problems.

What do you think would be the best option?

[grst]

Dear @apostovskaya,

we are working on integrating your tcrBLOSUM matrix into scirpy, which contains a reimplementation of the TCRdist algorithm.

We are unsure how you intended the matrix to be used with TCRdist. Could you please clarify what you mean with

Firstly, we transformed the tcrBLOSUM similarity matrix into a distance matrix according to the rules of TCRdist [12]

in your paper and how you would suggest to set the cap?

[grst]

"The mismatch distance is defined based on the BLOSUM62 (ref. 37) substitution matrix as follows: distance (a, a)=0; distance (a, b)=min (4, 4-BLOSUM62 (a, b)), where 4 is 1 unit greater than the most favourable BLOSUM62 score for a mismatch, and a and b are amino acids".

tbh, I never understood this part about TCRdist. It skews the matrix quite a bit, basically assigning a distance of 4 to all mismatches that have a negative distance in BLOSUM62 (even if it's just -1). So the only way to get a score $1 \leq s \leq 3$ is if you have one of the rare pairs with a positive mismatch score according to BLOSUM62.

Having a cap of 2 with TCRblosum doesn't make sense to me, because then the strong negative effect of mismatches with C would be completely gone which defeats the purpose of TCRblosum. Using 4 is just as arbitrary.

[felixpetschko]

I updated the approach now according to what was discussed via email. The distance matrices are now computed from the original substitution matrices (blosum62, tcrblosum_alpha, tcrblosum_beta) rather than stored as precomputed distance tables.

In general, the calculation is now done as follows:
distance(a, a) = 0 and distance(a, b) = min(distance_cap, distance_offset - score)

The base_matrix(Literal["blosum62", "tcrblosum"] = "blosum62") and the distance_cap(int | None | Literal["default"] = "default") parameters can be set in the calling ir_dist method. For blosum62 distance_cap=4 by default and distance_offset=4 as always. For the tcrblosum matrices, distance_cap=None (no cap) by default and distance_offset=3=(global maximum off-diagonal substitution score of both alpha and beta matrices +1). This distance offset maps the most favorable off-diagonal score across both tcrblosum matrices (2) to distance 1. That way, the blosum62 defaults follow the original TCRdist paper while the tcrblosum defaults should preserve the signal present in the more sparse tcrblosum matrices.

My intention for the user call would be the following rather than creating an additional tcrdist_tcrblosum metric:

tcrdist with blosum62

ir.pp.ir_dist(
    mdata,
    metric="tcrdist",
    sequence="aa",
    cutoff=15,
    #base_matrix="blosum62", #set by default
    #distance_cap=4, #set by default
)

tcrdist with tcrblosum matrices

ir.pp.ir_dist(
    mdata,
    metric="tcrdist",
    sequence="aa",
    cutoff=15,
    base_matrix="tcrblosum",
    #distance_cap=None, #set by default
)

[grst]

Does this imply that letters in the alphabet, but not in the matrix have a distance of 0 (e.g. X vs. A)?

Is this what we want? Or maybe rather Inf?

FWIW, there's also variants of BLOSUM62 that include these extended amino acids:

[felixpetschko]

Does this imply that letters in the alphabet, but not in the matrix have a distance of 0 (e.g. X vs. A)?

Yes - this behavior comes from the original pwseqdist base implementation. Maybe this was the pragmatic approach to treat unknown residues as "uninformative" rather than as evidence of dissimilarity. One could see capturing true matches as more important than avoiding false matches.

The question is if we want to change this.

FWIW, there's also variants of BLOSUM62 that include these extended amino acids

However, probably we don't have something like that for the tcrblosum matrices.

[grst]

Implementation-wise this looks good, and I also prefer to make the tcrdist implementation generic such that it works with arbitrary matrices rather than introducing another metric.

We just need to make sure to advertise it appropriately in the tutorial, otherwise the functionality is a bit hidden.

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[felixpetschko]

We just need to make sure to advertise it appropriately in the tutorial, otherwise the functionality is a bit hidden.

I made a short tutorial update to mention tcrblosum usage.

[grst]

LGTM now!