Add Needleman-Wunsch distance as a faster alignment metric

[felixpetschko]

Summary

This PR adds a Numba-optimized Needleman-Wunsch sequence distance metric via metric="needleman_wunsch". Conceptually, this metric is very similar to the existing alignment metric, with two main differences: it is much faster for large datasets, and it uses a linear gap model, where every gap position receives the same penalty.

The existing alignment metric allows gap openings and gap extensions to be penalized differently. However, this distinction is not used with the current default parameters anyway, where gap_open == gap_extend. If the alignment metric is configured such that gap_open == gap_extend == gap_penalty, both metrics should return equal results.

The reason for using a linear gap model is mostly practical. Supporting separate gap-open and gap-extension penalties requires handling three dynamic programming matrices instead of one for each sequence pair, which adds substantial computational overhead. For CDR3 sequences, which are usually quite short, I would generally expect users to choose parameters that allow only a limited number of gap positions. For example, in my test runs I used cutoff=10 and gap_penalty=4, which allows up to two gap positions. In this setting, distinguishing between two separate single-position gaps and one adjacent two-position gap is probably not worth the additional runtime. Therefore, I do not think the restriction to a single gap penalty is a major limitation.

The problem with the current alignment metric is that it is too slow for large datasets. The fastalignment metric improves runtime, but it can also suffer from performance limitations and does not always return exact results. Besides that, the current default parameters do not allow gaps (cutoff=10, gap_open=gap_extend=11). In that default setting, only equal-length sequences can fall below the cutoff, so the metric effectively computes alignment distances between equal-length sequences without making use of meaningful gap placement.

Therefore, I implemented the Needleman-Wunsch distance in a similar style to the Numba-optimized TCRdist CPU implementation. It supports BLOSUM62 by default and can also use TCRBLOSUM alpha/beta matrices through base_matrix="tcrblosum". With the improved performance, I was able to run the 8 million-cell Omniscope COVID dataset with 64 CPU cores using cutoff=10 and gap_penalty=4 within around 12 hours.

Overall, I think this metric is a useful addition because it provides an exact alignment-based distance that can still handle larger datasets. It is more flexible than the Hamming distance, but substantially faster than the existing alignment metric. In contrast to TCRdist, it is also not specific to TCR CDR3 sequences and can therefore be used for BCR analyses as well.

The most useful default values for gap_penalty and cutoff are open for discussion, since the performance improvements make less restrictive parameter choices feasible. It could also be discussed whether the existing alignment metric is still needed if its additional gap parameterization and support for additional Parasail substitution matrices are not required.

Main changes

• Add NeedlemanWunschDistanceCalculator and expose it via metric="needleman_wunsch" in sequence_dist and ir_dist.
• Support BLOSUM62 and TCRBLOSUM alpha/beta substitution matrices in a similar way as the TCRdist metric.
• Implement a Numba-based Needleman-Wunsch dynamic programming algorithm with cutoff-based optimizations.
• Add test cases.

[codecov]

Codecov Report

❌ Patch coverage is 49.13295% with 88 lines in your changes missing coverage. Please review.
✅ Project coverage is 77.32%. Comparing base (79a8344) to head (01ceae6).

Files with missing lines	Patch %	Lines
src/scirpy/ir_dist/metrics.py	32.03%	87 Missing ⚠️
src/scirpy/ir_dist/_substitution_matrices.py	97.67%	1 Missing ⚠️

Additional details and impacted files

@@ Coverage Diff @@## main #725 +/- ##
==========================================
- Coverage 78.44% 77.32% -1.13% 
==========================================
Files 51 52 +1 Lines 4635 4785 +150 ==========================================
+ Hits 3636 3700 +64 - Misses 999 1085 +86 

Files with missing lines	Coverage Δ
src/scirpy/ir_dist/__init__.py	93.13% <100.00%> (+0.13%)	⬆️
src/scirpy/ir_dist/_substitution_matrices.py	97.67% <97.67%> (ø)
src/scirpy/ir_dist/metrics.py	50.86% <32.03%> (-5.83%)	⬇️

... and 1 file with indirect coverage changes

🚀 New features to boost your workflow:

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

[grst]

Hi @felixpetschko,

thanks for working on this! From a technical perspective this looks all great.

My main concern here is that we keep adding metrics without providing guidelines (and evidence for) which metric to use.

In contrast to TCRdist, it is also not specific to TCR CDR3 sequences and can therefore be used for BCR analyses as well.

In what way is TCRdist specific to TCR sequences that would prevent it from using it for BCR? Or rather in what way is the alignment distance superior for BCR? My understanding would be that the main difference is that TCRdist allows only for a single gap position, while alignment allows for multiple, but does this make a big difference in practice?

Also how does it compare to TCRdist in terms of speed? E.g. how long would TCRdist take on the omniscope dataset on the same hardware?

[felixpetschko]

Hi @grst

My main concern here is that we keep adding metrics without providing guidelines (and evidence for) which metric to use.

My intention was rather to provide a faster implementation of the existing alignment distance metric than introducing a new metric. It is probably hard to provide real evidence for which metric to use. However, Needleman-Wunsch is a standard alignment algorithm and is also used in related sequence-analysis tools, so I think it is reasonable to use it.

In what way is TCRdist specific to TCR sequences that would prevent it from using it for BCR? Or rather in what way is the alignment distance superior for BCR? My understanding would be that the main difference is that TCRdist allows only for a single gap position, while alignment allows for multiple, but does this make a big difference in practice?

Actually, I was mainly focusing on performance, and I do not have proof for which metric is better in which case. My reasoning was that TCRdist's approach, with trimming from the N and C terminus and a single gap region, seems more targeted towards the TCR model. Needleman-Wunsch might be easier to justify for BCR CDR3 comparisons because it does not impose TCRdist's trimming and single-gap-region assumptions, and can handle length differences with a general global alignment.

Also how does it compare to TCRdist in terms of speed? E.g. how long would TCRdist take on the omniscope dataset on the same hardware?

TCRdist can run the full Omniscope COVID dataset on the same hardware in around 2.3 hours with default parameters, which makes it around 5 times faster than Needleman-Wunsch in my test run. The main reason is that, with Scirpy's default parameters, the gap position is computed by a formula and it is not necessary to try different gap positions. In contrast, Needleman-Wunsch computes the optimal global alignment that minimizes the distance.

[grst]

Do you think we could implement this without any user-facing changes?
I.e. integrate this into the AlignmentDistanceCalculator and fall back to the parasail implementation when gap_open != gap_extend?

(would be interesting to know if in the history of scirpy anyone has ever changed these default parameters. I'd guess not).

[felixpetschko]

Do you think we could implement this without any user-facing changes?

Yes, I will do that!

(would be interesting to know if in the history of scirpy anyone has ever changed these default parameters. I'd guess not).

However, I think we should definitely change the default params and set them in a way such that gaps are allowed. Otherwise there is no alignment done at all. I would allow at least 2 gaps such that it's worth to even run the dynamic programming alignment algorithm.

[grst]

Just changing defaults is also not very good practice... so we'd at least have to warn about it.
But maybe it is better to keep needleman_wunsch as a separate alignment metric and deprecate "alignment" and "fastalignment" distances.

We could complement it with a "metrics" guide in the documentation that explains the pros/cons and usecases of the different metrics.

[felixpetschko]

Alright, then let's do it like that 👍

[felixpetschko]

Now I marked alignment and fastalignment as deprecated in favor of needleman_wunsch and added a guide for choosing a sequence distance metric via #728.

[grst]

I have a few more comments, mostly on improving the organization of the metrics submodule.

[grst]

This would fail on joblib backends that do not support return_as="generator", e.g. dask.

See scirpy.utils._parallelize_with_joblib for a helper that addresses this.

[grst]

I'd move all these definitions to the top of the file, or maybe even better, a separate submodule within the ir_dist package.

[grst]

Also, maybe worth renaming this simply to AA_ALPHABET? Or is this still specific to parasail in any way?

[grst]

... or maybe make them dataclasses, that match together alphabet and substitution matrix

@dataclassclassSubstitutionMatrix:
alphabet: strmatrix: np.ndarrayBLOSUM62=SubstitutionMatrix(alphabet="ARN...", matrix=np.array([...]))
TCRBLOSUM_ALPHA=SubstitutionMatrix(...)

[grst]

Is there any reason for defining these as class variables instead of directly referencing the constants?

[grst]

Again, wouldn't it be easier to read if the constants were directly used everywhere?

[grst]

Could this become a metric-agnostic helper function? I think we already have similar code in other metrics...

[grst]

Suggested change

	If`gap_open == gap_extend`, useNeedlemanWunschDistanceCalculatorinstead.
	If`gap_open == gap_extend` (thedefault), useNeedlemanWunschDistanceCalculatorinstead, whichprovidesidenticalresultswhilebeingmuchfaster. Ifyouactuallyhaveause-caseforaffinegappenalties, pleaseletusknowbyopeninganissueonGitHub.

[grst]

Suggested change

	If`gap_open == gap_extend`, useNeedlemanWunschDistanceCalculatorinstead.
	If`gap_open == gap_extend` (thedefault), useNeedlemanWunschDistanceCalculatorinstead, whichprovidesidenticalresultswhilebeingmuchfaster. Ifyouactuallyhaveause-caseforaffinegappenalties, pleaseletusknowbyopeninganissueonGitHub.

[grst]

How has this been derived? Parasail?

[felixpetschko]

Yes, it was calculated with the deprecated AlignmentDistanceCalculator and the WU3k dataset (that only contains the 20 canonical amino acids).

[grst]

regarding deprecations, take a look at #735 please that switches to the decorators provided by scverse-misc.

[felixpetschko]

Hi @grst! I have implemented the requested changes - in particular, substitution matrices and the corresponding alphabets are now grouped in a dedicated _substitution_matrices submodule using a small dataclass, and the matrix mapping and substitution-to-distance conversion are provided as metric-independent helper functions. I also switched the block-wise parallelization to Scirpy’s joblib helper and updated the deprecations to use scverse-misc.

Besides that, I encountered a bug that was already present in the old AlignmentDistanceCalculator: With the ambiguous symbols in the blosum62 matrix (especially X) it is possible to get negative distance values, e.g. AAAXX vs AAAAA would result in a distance of -2. This situation occurs if there are values that score higher than self alignment. In the old implementation this causes datatype overflow or in the case of a distance of -1 ( +1 for storing = 0) the value is not stored at all due to the sparse format. In theory, one could just set negative result distance values to 0 (1 for storing). But logic-wise it seems a little bit odd to me that self-alignment cannot be interpreted as the maximum achievable alignment score. It also leads to problems with my performance optimization and setting the cutoff, because you cannot infer the number of allowed gaps by the cutoff and the gap penalty. This is due to the fact that parts of the sequence that score higher than self alignment can compensate for gap penalties. The symbols affected by this are X of the blosum62 matrix with ambiguous symbols and T from the TCRblosum beta matrix.

Therefore my suggestion would be to restrict the Needleman-Wunsch distance to the blosum62 matrix with only the 20 canonical amino acids (this is also the blosum62 matrix that can be chosen for TCRdist). In this case the matrix entries of the (unknown) ambiguous symbols would just be automatically set to 0 and the problems mentioned would be resolved. Currently, it is implemented in this PR like that. What do you think about that?

[felixpetschko]

Summary

This PR adds a Numba-optimized Needleman-Wunsch sequence distance metric via metric="needleman_wunsch". Conceptually, this metric is very similar to the existing alignment metric, with two main differences: it is much faster for large datasets, and it uses a linear gap model, where every gap position receives the same penalty.

The existing alignment metric allows gap openings and gap extensions to be penalized differently. However, this distinction is not used with the current default parameters anyway, where gap_open == gap_extend. If the alignment metric is configured such that gap_open == gap_extend == gap_penalty, both metrics should return equal results.

The reason for using a linear gap model is mostly practical. Supporting separate gap-open and gap-extension penalties requires handling three dynamic programming matrices instead of one for each sequence pair, which adds substantial computational overhead. For CDR3 sequences, which are usually quite short, I would generally expect users to choose parameters that allow only a limited number of gap positions. For example, in my test runs I used cutoff=10 and gap_penalty=4, which allows up to two gap positions. In this setting, distinguishing between two separate single-position gaps and one adjacent two-position gap is probably not worth the additional runtime. Therefore, I do not think the restriction to a single gap penalty is a major limitation.

The problem with the current alignment metric is that it is too slow for large datasets. The fastalignment metric improves runtime, but it can also suffer from performance limitations and does not always return exact results. Besides that, the current default parameters do not allow gaps (cutoff=10, gap_open=gap_extend=11). In that default setting, only equal-length sequences can fall below the cutoff, so the metric effectively computes alignment distances between equal-length sequences without making use of meaningful gap placement.

Therefore, I implemented the Needleman-Wunsch distance in a similar style to the Numba-optimized TCRdist CPU implementation. It supports BLOSUM62 by default and can also use TCRBLOSUM alpha/beta matrices through base_matrix="tcrblosum". With the improved performance, I was able to run the 8 million-cell Omniscope COVID dataset with 64 CPU cores using cutoff=10 and gap_penalty=4 within around 12 hours.

Overall, I think this metric is a useful addition because it provides an exact alignment-based distance that can still handle larger datasets. It is more flexible than the Hamming distance, but substantially faster than the existing alignment metric. In contrast to TCRdist, it is also not specific to TCR CDR3 sequences and can therefore be used for BCR analyses as well.

The most useful default values for gap_penalty and cutoff are open for discussion, since the performance improvements make less restrictive parameter choices feasible. It could also be discussed whether the existing alignment metric is still needed if its additional gap parameterization and support for additional Parasail substitution matrices are not required.

Main changes

• Add NeedlemanWunschDistanceCalculator and expose it via metric="needleman_wunsch" in sequence_dist and ir_dist.
• Support BLOSUM62 and TCRBLOSUM alpha/beta substitution matrices in a similar way as the TCRdist metric.
• Implement a Numba-based Needleman-Wunsch dynamic programming algorithm with cutoff-based optimizations.
• Add test cases.

[codecov]

Codecov Report

❌ Patch coverage is 49.13295% with 88 lines in your changes missing coverage. Please review.
✅ Project coverage is 77.32%. Comparing base (79a8344) to head (01ceae6).

Files with missing lines	Patch %	Lines
src/scirpy/ir_dist/metrics.py	32.03%	87 Missing ⚠️
src/scirpy/ir_dist/_substitution_matrices.py	97.67%	1 Missing ⚠️

Additional details and impacted files

@@ Coverage Diff @@## main #725 +/- ##
==========================================
- Coverage 78.44% 77.32% -1.13% 
==========================================
Files 51 52 +1 Lines 4635 4785 +150 ==========================================
+ Hits 3636 3700 +64 - Misses 999 1085 +86 

Files with missing lines	Coverage Δ
src/scirpy/ir_dist/__init__.py	93.13% <100.00%> (+0.13%)	⬆️
src/scirpy/ir_dist/_substitution_matrices.py	97.67% <97.67%> (ø)
src/scirpy/ir_dist/metrics.py	50.86% <32.03%> (-5.83%)	⬇️

... and 1 file with indirect coverage changes

🚀 New features to boost your workflow:

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

[grst]

Hi @felixpetschko,

thanks for working on this! From a technical perspective this looks all great.

My main concern here is that we keep adding metrics without providing guidelines (and evidence for) which metric to use.

In contrast to TCRdist, it is also not specific to TCR CDR3 sequences and can therefore be used for BCR analyses as well.

In what way is TCRdist specific to TCR sequences that would prevent it from using it for BCR? Or rather in what way is the alignment distance superior for BCR? My understanding would be that the main difference is that TCRdist allows only for a single gap position, while alignment allows for multiple, but does this make a big difference in practice?

Also how does it compare to TCRdist in terms of speed? E.g. how long would TCRdist take on the omniscope dataset on the same hardware?

[felixpetschko]

Hi @grst

My main concern here is that we keep adding metrics without providing guidelines (and evidence for) which metric to use.

My intention was rather to provide a faster implementation of the existing alignment distance metric than introducing a new metric. It is probably hard to provide real evidence for which metric to use. However, Needleman-Wunsch is a standard alignment algorithm and is also used in related sequence-analysis tools, so I think it is reasonable to use it.

In what way is TCRdist specific to TCR sequences that would prevent it from using it for BCR? Or rather in what way is the alignment distance superior for BCR? My understanding would be that the main difference is that TCRdist allows only for a single gap position, while alignment allows for multiple, but does this make a big difference in practice?

Actually, I was mainly focusing on performance, and I do not have proof for which metric is better in which case. My reasoning was that TCRdist's approach, with trimming from the N and C terminus and a single gap region, seems more targeted towards the TCR model. Needleman-Wunsch might be easier to justify for BCR CDR3 comparisons because it does not impose TCRdist's trimming and single-gap-region assumptions, and can handle length differences with a general global alignment.

Also how does it compare to TCRdist in terms of speed? E.g. how long would TCRdist take on the omniscope dataset on the same hardware?

TCRdist can run the full Omniscope COVID dataset on the same hardware in around 2.3 hours with default parameters, which makes it around 5 times faster than Needleman-Wunsch in my test run. The main reason is that, with Scirpy's default parameters, the gap position is computed by a formula and it is not necessary to try different gap positions. In contrast, Needleman-Wunsch computes the optimal global alignment that minimizes the distance.

[grst]

Do you think we could implement this without any user-facing changes?
I.e. integrate this into the AlignmentDistanceCalculator and fall back to the parasail implementation when gap_open != gap_extend?

(would be interesting to know if in the history of scirpy anyone has ever changed these default parameters. I'd guess not).

[felixpetschko]

Do you think we could implement this without any user-facing changes?

Yes, I will do that!

(would be interesting to know if in the history of scirpy anyone has ever changed these default parameters. I'd guess not).

However, I think we should definitely change the default params and set them in a way such that gaps are allowed. Otherwise there is no alignment done at all. I would allow at least 2 gaps such that it's worth to even run the dynamic programming alignment algorithm.

[grst]

Just changing defaults is also not very good practice... so we'd at least have to warn about it.
But maybe it is better to keep needleman_wunsch as a separate alignment metric and deprecate "alignment" and "fastalignment" distances.

We could complement it with a "metrics" guide in the documentation that explains the pros/cons and usecases of the different metrics.

[felixpetschko]

Alright, then let's do it like that 👍

[felixpetschko]

Now I marked alignment and fastalignment as deprecated in favor of needleman_wunsch and added a guide for choosing a sequence distance metric via #728.

[grst]

I have a few more comments, mostly on improving the organization of the metrics submodule.

[grst]

This would fail on joblib backends that do not support return_as="generator", e.g. dask.

See scirpy.utils._parallelize_with_joblib for a helper that addresses this.

[grst]

I'd move all these definitions to the top of the file, or maybe even better, a separate submodule within the ir_dist package.

[grst]

Also, maybe worth renaming this simply to AA_ALPHABET? Or is this still specific to parasail in any way?

[grst]

... or maybe make them dataclasses, that match together alphabet and substitution matrix

@dataclassclassSubstitutionMatrix:
alphabet: strmatrix: np.ndarrayBLOSUM62=SubstitutionMatrix(alphabet="ARN...", matrix=np.array([...]))
TCRBLOSUM_ALPHA=SubstitutionMatrix(...)

[grst]

Is there any reason for defining these as class variables instead of directly referencing the constants?

[grst]

Again, wouldn't it be easier to read if the constants were directly used everywhere?

[grst]

Could this become a metric-agnostic helper function? I think we already have similar code in other metrics...

[grst]

Suggested change

	If`gap_open == gap_extend`, useNeedlemanWunschDistanceCalculatorinstead.
	If`gap_open == gap_extend` (thedefault), useNeedlemanWunschDistanceCalculatorinstead, whichprovidesidenticalresultswhilebeingmuchfaster. Ifyouactuallyhaveause-caseforaffinegappenalties, pleaseletusknowbyopeninganissueonGitHub.

[grst]

Suggested change

	If`gap_open == gap_extend`, useNeedlemanWunschDistanceCalculatorinstead.
	If`gap_open == gap_extend` (thedefault), useNeedlemanWunschDistanceCalculatorinstead, whichprovidesidenticalresultswhilebeingmuchfaster. Ifyouactuallyhaveause-caseforaffinegappenalties, pleaseletusknowbyopeninganissueonGitHub.

[grst]

How has this been derived? Parasail?

[felixpetschko]

Yes, it was calculated with the deprecated AlignmentDistanceCalculator and the WU3k dataset (that only contains the 20 canonical amino acids).

[grst]

regarding deprecations, take a look at #735 please that switches to the decorators provided by scverse-misc.

[felixpetschko]

Hi @grst! I have implemented the requested changes - in particular, substitution matrices and the corresponding alphabets are now grouped in a dedicated _substitution_matrices submodule using a small dataclass, and the matrix mapping and substitution-to-distance conversion are provided as metric-independent helper functions. I also switched the block-wise parallelization to Scirpy’s joblib helper and updated the deprecations to use scverse-misc.

Besides that, I encountered a bug that was already present in the old AlignmentDistanceCalculator: With the ambiguous symbols in the blosum62 matrix (especially X) it is possible to get negative distance values, e.g. AAAXX vs AAAAA would result in a distance of -2. This situation occurs if there are values that score higher than self alignment. In the old implementation this causes datatype overflow or in the case of a distance of -1 ( +1 for storing = 0) the value is not stored at all due to the sparse format. In theory, one could just set negative result distance values to 0 (1 for storing). But logic-wise it seems a little bit odd to me that self-alignment cannot be interpreted as the maximum achievable alignment score. It also leads to problems with my performance optimization and setting the cutoff, because you cannot infer the number of allowed gaps by the cutoff and the gap penalty. This is due to the fact that parts of the sequence that score higher than self alignment can compensate for gap penalties. The symbols affected by this are X of the blosum62 matrix with ambiguous symbols and T from the TCRblosum beta matrix.

Therefore my suggestion would be to restrict the Needleman-Wunsch distance to the blosum62 matrix with only the 20 canonical amino acids (this is also the blosum62 matrix that can be chosen for TCRdist). In this case the matrix entries of the (unknown) ambiguous symbols would just be automatically set to 0 and the problems mentioned would be resolved. Currently, it is implemented in this PR like that. What do you think about that?

[felixpetschko]

Summary

This PR adds a Numba-optimized Needleman-Wunsch sequence distance metric via metric="needleman_wunsch". Conceptually, this metric is very similar to the existing alignment metric, with two main differences: it is much faster for large datasets, and it uses a linear gap model, where every gap position receives the same penalty.

The existing alignment metric allows gap openings and gap extensions to be penalized differently. However, this distinction is not used with the current default parameters anyway, where gap_open == gap_extend. If the alignment metric is configured such that gap_open == gap_extend == gap_penalty, both metrics should return equal results.

The reason for using a linear gap model is mostly practical. Supporting separate gap-open and gap-extension penalties requires handling three dynamic programming matrices instead of one for each sequence pair, which adds substantial computational overhead. For CDR3 sequences, which are usually quite short, I would generally expect users to choose parameters that allow only a limited number of gap positions. For example, in my test runs I used cutoff=10 and gap_penalty=4, which allows up to two gap positions. In this setting, distinguishing between two separate single-position gaps and one adjacent two-position gap is probably not worth the additional runtime. Therefore, I do not think the restriction to a single gap penalty is a major limitation.

The problem with the current alignment metric is that it is too slow for large datasets. The fastalignment metric improves runtime, but it can also suffer from performance limitations and does not always return exact results. Besides that, the current default parameters do not allow gaps (cutoff=10, gap_open=gap_extend=11). In that default setting, only equal-length sequences can fall below the cutoff, so the metric effectively computes alignment distances between equal-length sequences without making use of meaningful gap placement.

Therefore, I implemented the Needleman-Wunsch distance in a similar style to the Numba-optimized TCRdist CPU implementation. It supports BLOSUM62 by default and can also use TCRBLOSUM alpha/beta matrices through base_matrix="tcrblosum". With the improved performance, I was able to run the 8 million-cell Omniscope COVID dataset with 64 CPU cores using cutoff=10 and gap_penalty=4 within around 12 hours.

Overall, I think this metric is a useful addition because it provides an exact alignment-based distance that can still handle larger datasets. It is more flexible than the Hamming distance, but substantially faster than the existing alignment metric. In contrast to TCRdist, it is also not specific to TCR CDR3 sequences and can therefore be used for BCR analyses as well.

The most useful default values for gap_penalty and cutoff are open for discussion, since the performance improvements make less restrictive parameter choices feasible. It could also be discussed whether the existing alignment metric is still needed if its additional gap parameterization and support for additional Parasail substitution matrices are not required.

Main changes

• Add NeedlemanWunschDistanceCalculator and expose it via metric="needleman_wunsch" in sequence_dist and ir_dist.
• Support BLOSUM62 and TCRBLOSUM alpha/beta substitution matrices in a similar way as the TCRdist metric.
• Implement a Numba-based Needleman-Wunsch dynamic programming algorithm with cutoff-based optimizations.
• Add test cases.

[codecov]

Codecov Report

❌ Patch coverage is 49.13295% with 88 lines in your changes missing coverage. Please review.
✅ Project coverage is 77.32%. Comparing base (79a8344) to head (01ceae6).

Files with missing lines	Patch %	Lines
src/scirpy/ir_dist/metrics.py	32.03%	87 Missing ⚠️
src/scirpy/ir_dist/_substitution_matrices.py	97.67%	1 Missing ⚠️

Additional details and impacted files

@@ Coverage Diff @@## main #725 +/- ##
==========================================
- Coverage 78.44% 77.32% -1.13% 
==========================================
Files 51 52 +1 Lines 4635 4785 +150 ==========================================
+ Hits 3636 3700 +64 - Misses 999 1085 +86 

Files with missing lines	Coverage Δ
src/scirpy/ir_dist/__init__.py	93.13% <100.00%> (+0.13%)	⬆️
src/scirpy/ir_dist/_substitution_matrices.py	97.67% <97.67%> (ø)
src/scirpy/ir_dist/metrics.py	50.86% <32.03%> (-5.83%)	⬇️

... and 1 file with indirect coverage changes

🚀 New features to boost your workflow:

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

[grst]

Hi @felixpetschko,

thanks for working on this! From a technical perspective this looks all great.

My main concern here is that we keep adding metrics without providing guidelines (and evidence for) which metric to use.

In contrast to TCRdist, it is also not specific to TCR CDR3 sequences and can therefore be used for BCR analyses as well.

In what way is TCRdist specific to TCR sequences that would prevent it from using it for BCR? Or rather in what way is the alignment distance superior for BCR? My understanding would be that the main difference is that TCRdist allows only for a single gap position, while alignment allows for multiple, but does this make a big difference in practice?

Also how does it compare to TCRdist in terms of speed? E.g. how long would TCRdist take on the omniscope dataset on the same hardware?

[felixpetschko]

Hi @grst

My main concern here is that we keep adding metrics without providing guidelines (and evidence for) which metric to use.

My intention was rather to provide a faster implementation of the existing alignment distance metric than introducing a new metric. It is probably hard to provide real evidence for which metric to use. However, Needleman-Wunsch is a standard alignment algorithm and is also used in related sequence-analysis tools, so I think it is reasonable to use it.

In what way is TCRdist specific to TCR sequences that would prevent it from using it for BCR? Or rather in what way is the alignment distance superior for BCR? My understanding would be that the main difference is that TCRdist allows only for a single gap position, while alignment allows for multiple, but does this make a big difference in practice?

Actually, I was mainly focusing on performance, and I do not have proof for which metric is better in which case. My reasoning was that TCRdist's approach, with trimming from the N and C terminus and a single gap region, seems more targeted towards the TCR model. Needleman-Wunsch might be easier to justify for BCR CDR3 comparisons because it does not impose TCRdist's trimming and single-gap-region assumptions, and can handle length differences with a general global alignment.

Also how does it compare to TCRdist in terms of speed? E.g. how long would TCRdist take on the omniscope dataset on the same hardware?

TCRdist can run the full Omniscope COVID dataset on the same hardware in around 2.3 hours with default parameters, which makes it around 5 times faster than Needleman-Wunsch in my test run. The main reason is that, with Scirpy's default parameters, the gap position is computed by a formula and it is not necessary to try different gap positions. In contrast, Needleman-Wunsch computes the optimal global alignment that minimizes the distance.

[grst]

Do you think we could implement this without any user-facing changes?
I.e. integrate this into the AlignmentDistanceCalculator and fall back to the parasail implementation when gap_open != gap_extend?

(would be interesting to know if in the history of scirpy anyone has ever changed these default parameters. I'd guess not).

[felixpetschko]

Do you think we could implement this without any user-facing changes?

Yes, I will do that!

(would be interesting to know if in the history of scirpy anyone has ever changed these default parameters. I'd guess not).

However, I think we should definitely change the default params and set them in a way such that gaps are allowed. Otherwise there is no alignment done at all. I would allow at least 2 gaps such that it's worth to even run the dynamic programming alignment algorithm.

[grst]

Just changing defaults is also not very good practice... so we'd at least have to warn about it.
But maybe it is better to keep needleman_wunsch as a separate alignment metric and deprecate "alignment" and "fastalignment" distances.

We could complement it with a "metrics" guide in the documentation that explains the pros/cons and usecases of the different metrics.

[felixpetschko]

Alright, then let's do it like that 👍

[felixpetschko]

Now I marked alignment and fastalignment as deprecated in favor of needleman_wunsch and added a guide for choosing a sequence distance metric via #728.

[grst]

I have a few more comments, mostly on improving the organization of the metrics submodule.

[grst]

This would fail on joblib backends that do not support return_as="generator", e.g. dask.

See scirpy.utils._parallelize_with_joblib for a helper that addresses this.

[grst]

I'd move all these definitions to the top of the file, or maybe even better, a separate submodule within the ir_dist package.

[grst]

Also, maybe worth renaming this simply to AA_ALPHABET? Or is this still specific to parasail in any way?

[grst]

... or maybe make them dataclasses, that match together alphabet and substitution matrix

@dataclassclassSubstitutionMatrix:
alphabet: strmatrix: np.ndarrayBLOSUM62=SubstitutionMatrix(alphabet="ARN...", matrix=np.array([...]))
TCRBLOSUM_ALPHA=SubstitutionMatrix(...)

[grst]

Is there any reason for defining these as class variables instead of directly referencing the constants?

[grst]

Again, wouldn't it be easier to read if the constants were directly used everywhere?

[grst]

Could this become a metric-agnostic helper function? I think we already have similar code in other metrics...

[grst]

Suggested change

	If`gap_open == gap_extend`, useNeedlemanWunschDistanceCalculatorinstead.
	If`gap_open == gap_extend` (thedefault), useNeedlemanWunschDistanceCalculatorinstead, whichprovidesidenticalresultswhilebeingmuchfaster. Ifyouactuallyhaveause-caseforaffinegappenalties, pleaseletusknowbyopeninganissueonGitHub.

[grst]

Suggested change

	If`gap_open == gap_extend`, useNeedlemanWunschDistanceCalculatorinstead.
	If`gap_open == gap_extend` (thedefault), useNeedlemanWunschDistanceCalculatorinstead, whichprovidesidenticalresultswhilebeingmuchfaster. Ifyouactuallyhaveause-caseforaffinegappenalties, pleaseletusknowbyopeninganissueonGitHub.

[grst]

How has this been derived? Parasail?

[felixpetschko]

Yes, it was calculated with the deprecated AlignmentDistanceCalculator and the WU3k dataset (that only contains the 20 canonical amino acids).

[grst]

regarding deprecations, take a look at #735 please that switches to the decorators provided by scverse-misc.

[felixpetschko]

Hi @grst! I have implemented the requested changes - in particular, substitution matrices and the corresponding alphabets are now grouped in a dedicated _substitution_matrices submodule using a small dataclass, and the matrix mapping and substitution-to-distance conversion are provided as metric-independent helper functions. I also switched the block-wise parallelization to Scirpy’s joblib helper and updated the deprecations to use scverse-misc.

Besides that, I encountered a bug that was already present in the old AlignmentDistanceCalculator: With the ambiguous symbols in the blosum62 matrix (especially X) it is possible to get negative distance values, e.g. AAAXX vs AAAAA would result in a distance of -2. This situation occurs if there are values that score higher than self alignment. In the old implementation this causes datatype overflow or in the case of a distance of -1 ( +1 for storing = 0) the value is not stored at all due to the sparse format. In theory, one could just set negative result distance values to 0 (1 for storing). But logic-wise it seems a little bit odd to me that self-alignment cannot be interpreted as the maximum achievable alignment score. It also leads to problems with my performance optimization and setting the cutoff, because you cannot infer the number of allowed gaps by the cutoff and the gap penalty. This is due to the fact that parts of the sequence that score higher than self alignment can compensate for gap penalties. The symbols affected by this are X of the blosum62 matrix with ambiguous symbols and T from the TCRblosum beta matrix.

Therefore my suggestion would be to restrict the Needleman-Wunsch distance to the blosum62 matrix with only the 20 canonical amino acids (this is also the blosum62 matrix that can be chosen for TCRdist). In this case the matrix entries of the (unknown) ambiguous symbols would just be automatically set to 0 and the problems mentioned would be resolved. Currently, it is implemented in this PR like that. What do you think about that?

[felixpetschko]

Summary

This PR adds a Numba-optimized Needleman-Wunsch sequence distance metric via metric="needleman_wunsch". Conceptually, this metric is very similar to the existing alignment metric, with two main differences: it is much faster for large datasets, and it uses a linear gap model, where every gap position receives the same penalty.

The existing alignment metric allows gap openings and gap extensions to be penalized differently. However, this distinction is not used with the current default parameters anyway, where gap_open == gap_extend. If the alignment metric is configured such that gap_open == gap_extend == gap_penalty, both metrics should return equal results.

The reason for using a linear gap model is mostly practical. Supporting separate gap-open and gap-extension penalties requires handling three dynamic programming matrices instead of one for each sequence pair, which adds substantial computational overhead. For CDR3 sequences, which are usually quite short, I would generally expect users to choose parameters that allow only a limited number of gap positions. For example, in my test runs I used cutoff=10 and gap_penalty=4, which allows up to two gap positions. In this setting, distinguishing between two separate single-position gaps and one adjacent two-position gap is probably not worth the additional runtime. Therefore, I do not think the restriction to a single gap penalty is a major limitation.

The problem with the current alignment metric is that it is too slow for large datasets. The fastalignment metric improves runtime, but it can also suffer from performance limitations and does not always return exact results. Besides that, the current default parameters do not allow gaps (cutoff=10, gap_open=gap_extend=11). In that default setting, only equal-length sequences can fall below the cutoff, so the metric effectively computes alignment distances between equal-length sequences without making use of meaningful gap placement.

Therefore, I implemented the Needleman-Wunsch distance in a similar style to the Numba-optimized TCRdist CPU implementation. It supports BLOSUM62 by default and can also use TCRBLOSUM alpha/beta matrices through base_matrix="tcrblosum". With the improved performance, I was able to run the 8 million-cell Omniscope COVID dataset with 64 CPU cores using cutoff=10 and gap_penalty=4 within around 12 hours.

Overall, I think this metric is a useful addition because it provides an exact alignment-based distance that can still handle larger datasets. It is more flexible than the Hamming distance, but substantially faster than the existing alignment metric. In contrast to TCRdist, it is also not specific to TCR CDR3 sequences and can therefore be used for BCR analyses as well.

The most useful default values for gap_penalty and cutoff are open for discussion, since the performance improvements make less restrictive parameter choices feasible. It could also be discussed whether the existing alignment metric is still needed if its additional gap parameterization and support for additional Parasail substitution matrices are not required.

Main changes

• Add NeedlemanWunschDistanceCalculator and expose it via metric="needleman_wunsch" in sequence_dist and ir_dist.
• Support BLOSUM62 and TCRBLOSUM alpha/beta substitution matrices in a similar way as the TCRdist metric.
• Implement a Numba-based Needleman-Wunsch dynamic programming algorithm with cutoff-based optimizations.
• Add test cases.

[codecov]

Codecov Report

❌ Patch coverage is 49.13295% with 88 lines in your changes missing coverage. Please review.
✅ Project coverage is 77.32%. Comparing base (79a8344) to head (01ceae6).

Files with missing lines	Patch %	Lines
src/scirpy/ir_dist/metrics.py	32.03%	87 Missing ⚠️
src/scirpy/ir_dist/_substitution_matrices.py	97.67%	1 Missing ⚠️

Additional details and impacted files

@@ Coverage Diff @@## main #725 +/- ##
==========================================
- Coverage 78.44% 77.32% -1.13% 
==========================================
Files 51 52 +1 Lines 4635 4785 +150 ==========================================
+ Hits 3636 3700 +64 - Misses 999 1085 +86 

Files with missing lines	Coverage Δ
src/scirpy/ir_dist/__init__.py	93.13% <100.00%> (+0.13%)	⬆️
src/scirpy/ir_dist/_substitution_matrices.py	97.67% <97.67%> (ø)
src/scirpy/ir_dist/metrics.py	50.86% <32.03%> (-5.83%)	⬇️

... and 1 file with indirect coverage changes

🚀 New features to boost your workflow:

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

[grst]

Hi @felixpetschko,

thanks for working on this! From a technical perspective this looks all great.

My main concern here is that we keep adding metrics without providing guidelines (and evidence for) which metric to use.

In contrast to TCRdist, it is also not specific to TCR CDR3 sequences and can therefore be used for BCR analyses as well.

In what way is TCRdist specific to TCR sequences that would prevent it from using it for BCR? Or rather in what way is the alignment distance superior for BCR? My understanding would be that the main difference is that TCRdist allows only for a single gap position, while alignment allows for multiple, but does this make a big difference in practice?

Also how does it compare to TCRdist in terms of speed? E.g. how long would TCRdist take on the omniscope dataset on the same hardware?

[felixpetschko]

Hi @grst

My main concern here is that we keep adding metrics without providing guidelines (and evidence for) which metric to use.

My intention was rather to provide a faster implementation of the existing alignment distance metric than introducing a new metric. It is probably hard to provide real evidence for which metric to use. However, Needleman-Wunsch is a standard alignment algorithm and is also used in related sequence-analysis tools, so I think it is reasonable to use it.

In what way is TCRdist specific to TCR sequences that would prevent it from using it for BCR? Or rather in what way is the alignment distance superior for BCR? My understanding would be that the main difference is that TCRdist allows only for a single gap position, while alignment allows for multiple, but does this make a big difference in practice?

Actually, I was mainly focusing on performance, and I do not have proof for which metric is better in which case. My reasoning was that TCRdist's approach, with trimming from the N and C terminus and a single gap region, seems more targeted towards the TCR model. Needleman-Wunsch might be easier to justify for BCR CDR3 comparisons because it does not impose TCRdist's trimming and single-gap-region assumptions, and can handle length differences with a general global alignment.

Also how does it compare to TCRdist in terms of speed? E.g. how long would TCRdist take on the omniscope dataset on the same hardware?

TCRdist can run the full Omniscope COVID dataset on the same hardware in around 2.3 hours with default parameters, which makes it around 5 times faster than Needleman-Wunsch in my test run. The main reason is that, with Scirpy's default parameters, the gap position is computed by a formula and it is not necessary to try different gap positions. In contrast, Needleman-Wunsch computes the optimal global alignment that minimizes the distance.

[grst]

Do you think we could implement this without any user-facing changes?
I.e. integrate this into the AlignmentDistanceCalculator and fall back to the parasail implementation when gap_open != gap_extend?

(would be interesting to know if in the history of scirpy anyone has ever changed these default parameters. I'd guess not).

[felixpetschko]

Do you think we could implement this without any user-facing changes?

Yes, I will do that!

(would be interesting to know if in the history of scirpy anyone has ever changed these default parameters. I'd guess not).

However, I think we should definitely change the default params and set them in a way such that gaps are allowed. Otherwise there is no alignment done at all. I would allow at least 2 gaps such that it's worth to even run the dynamic programming alignment algorithm.

[grst]

Just changing defaults is also not very good practice... so we'd at least have to warn about it.
But maybe it is better to keep needleman_wunsch as a separate alignment metric and deprecate "alignment" and "fastalignment" distances.

We could complement it with a "metrics" guide in the documentation that explains the pros/cons and usecases of the different metrics.

[felixpetschko]

Alright, then let's do it like that 👍

[felixpetschko]

Now I marked alignment and fastalignment as deprecated in favor of needleman_wunsch and added a guide for choosing a sequence distance metric via #728.

[grst]

I have a few more comments, mostly on improving the organization of the metrics submodule.

[grst]

This would fail on joblib backends that do not support return_as="generator", e.g. dask.

See scirpy.utils._parallelize_with_joblib for a helper that addresses this.

[grst]

I'd move all these definitions to the top of the file, or maybe even better, a separate submodule within the ir_dist package.

[grst]

Also, maybe worth renaming this simply to AA_ALPHABET? Or is this still specific to parasail in any way?

[grst]

... or maybe make them dataclasses, that match together alphabet and substitution matrix

@dataclassclassSubstitutionMatrix:
alphabet: strmatrix: np.ndarrayBLOSUM62=SubstitutionMatrix(alphabet="ARN...", matrix=np.array([...]))
TCRBLOSUM_ALPHA=SubstitutionMatrix(...)

[grst]

Is there any reason for defining these as class variables instead of directly referencing the constants?

[grst]

Again, wouldn't it be easier to read if the constants were directly used everywhere?

[grst]

Could this become a metric-agnostic helper function? I think we already have similar code in other metrics...

[grst]

Suggested change

	If`gap_open == gap_extend`, useNeedlemanWunschDistanceCalculatorinstead.
	If`gap_open == gap_extend` (thedefault), useNeedlemanWunschDistanceCalculatorinstead, whichprovidesidenticalresultswhilebeingmuchfaster. Ifyouactuallyhaveause-caseforaffinegappenalties, pleaseletusknowbyopeninganissueonGitHub.

[grst]

Suggested change

	If`gap_open == gap_extend`, useNeedlemanWunschDistanceCalculatorinstead.
	If`gap_open == gap_extend` (thedefault), useNeedlemanWunschDistanceCalculatorinstead, whichprovidesidenticalresultswhilebeingmuchfaster. Ifyouactuallyhaveause-caseforaffinegappenalties, pleaseletusknowbyopeninganissueonGitHub.

[grst]

How has this been derived? Parasail?

[felixpetschko]

Yes, it was calculated with the deprecated AlignmentDistanceCalculator and the WU3k dataset (that only contains the 20 canonical amino acids).

[grst]

regarding deprecations, take a look at #735 please that switches to the decorators provided by scverse-misc.

[felixpetschko]

Hi @grst! I have implemented the requested changes - in particular, substitution matrices and the corresponding alphabets are now grouped in a dedicated _substitution_matrices submodule using a small dataclass, and the matrix mapping and substitution-to-distance conversion are provided as metric-independent helper functions. I also switched the block-wise parallelization to Scirpy’s joblib helper and updated the deprecations to use scverse-misc.

Besides that, I encountered a bug that was already present in the old AlignmentDistanceCalculator: With the ambiguous symbols in the blosum62 matrix (especially X) it is possible to get negative distance values, e.g. AAAXX vs AAAAA would result in a distance of -2. This situation occurs if there are values that score higher than self alignment. In the old implementation this causes datatype overflow or in the case of a distance of -1 ( +1 for storing = 0) the value is not stored at all due to the sparse format. In theory, one could just set negative result distance values to 0 (1 for storing). But logic-wise it seems a little bit odd to me that self-alignment cannot be interpreted as the maximum achievable alignment score. It also leads to problems with my performance optimization and setting the cutoff, because you cannot infer the number of allowed gaps by the cutoff and the gap penalty. This is due to the fact that parts of the sequence that score higher than self alignment can compensate for gap penalties. The symbols affected by this are X of the blosum62 matrix with ambiguous symbols and T from the TCRblosum beta matrix.

Therefore my suggestion would be to restrict the Needleman-Wunsch distance to the blosum62 matrix with only the 20 canonical amino acids (this is also the blosum62 matrix that can be chosen for TCRdist). In this case the matrix entries of the (unknown) ambiguous symbols would just be automatically set to 0 and the problems mentioned would be resolved. Currently, it is implemented in this PR like that. What do you think about that?

[felixpetschko]

Summary

This PR adds a Numba-optimized Needleman-Wunsch sequence distance metric via metric="needleman_wunsch". Conceptually, this metric is very similar to the existing alignment metric, with two main differences: it is much faster for large datasets, and it uses a linear gap model, where every gap position receives the same penalty.

The existing alignment metric allows gap openings and gap extensions to be penalized differently. However, this distinction is not used with the current default parameters anyway, where gap_open == gap_extend. If the alignment metric is configured such that gap_open == gap_extend == gap_penalty, both metrics should return equal results.

The reason for using a linear gap model is mostly practical. Supporting separate gap-open and gap-extension penalties requires handling three dynamic programming matrices instead of one for each sequence pair, which adds substantial computational overhead. For CDR3 sequences, which are usually quite short, I would generally expect users to choose parameters that allow only a limited number of gap positions. For example, in my test runs I used cutoff=10 and gap_penalty=4, which allows up to two gap positions. In this setting, distinguishing between two separate single-position gaps and one adjacent two-position gap is probably not worth the additional runtime. Therefore, I do not think the restriction to a single gap penalty is a major limitation.

The problem with the current alignment metric is that it is too slow for large datasets. The fastalignment metric improves runtime, but it can also suffer from performance limitations and does not always return exact results. Besides that, the current default parameters do not allow gaps (cutoff=10, gap_open=gap_extend=11). In that default setting, only equal-length sequences can fall below the cutoff, so the metric effectively computes alignment distances between equal-length sequences without making use of meaningful gap placement.

Therefore, I implemented the Needleman-Wunsch distance in a similar style to the Numba-optimized TCRdist CPU implementation. It supports BLOSUM62 by default and can also use TCRBLOSUM alpha/beta matrices through base_matrix="tcrblosum". With the improved performance, I was able to run the 8 million-cell Omniscope COVID dataset with 64 CPU cores using cutoff=10 and gap_penalty=4 within around 12 hours.

Overall, I think this metric is a useful addition because it provides an exact alignment-based distance that can still handle larger datasets. It is more flexible than the Hamming distance, but substantially faster than the existing alignment metric. In contrast to TCRdist, it is also not specific to TCR CDR3 sequences and can therefore be used for BCR analyses as well.

The most useful default values for gap_penalty and cutoff are open for discussion, since the performance improvements make less restrictive parameter choices feasible. It could also be discussed whether the existing alignment metric is still needed if its additional gap parameterization and support for additional Parasail substitution matrices are not required.

Main changes

• Add NeedlemanWunschDistanceCalculator and expose it via metric="needleman_wunsch" in sequence_dist and ir_dist.
• Support BLOSUM62 and TCRBLOSUM alpha/beta substitution matrices in a similar way as the TCRdist metric.
• Implement a Numba-based Needleman-Wunsch dynamic programming algorithm with cutoff-based optimizations.
• Add test cases.

[codecov]

Codecov Report

❌ Patch coverage is 49.13295% with 88 lines in your changes missing coverage. Please review.
✅ Project coverage is 77.32%. Comparing base (79a8344) to head (01ceae6).

Files with missing lines	Patch %	Lines
src/scirpy/ir_dist/metrics.py	32.03%	87 Missing ⚠️
src/scirpy/ir_dist/_substitution_matrices.py	97.67%	1 Missing ⚠️

Additional details and impacted files

@@ Coverage Diff @@## main #725 +/- ##
==========================================
- Coverage 78.44% 77.32% -1.13% 
==========================================
Files 51 52 +1 Lines 4635 4785 +150 ==========================================
+ Hits 3636 3700 +64 - Misses 999 1085 +86 

Files with missing lines	Coverage Δ
src/scirpy/ir_dist/__init__.py	93.13% <100.00%> (+0.13%)	⬆️
src/scirpy/ir_dist/_substitution_matrices.py	97.67% <97.67%> (ø)
src/scirpy/ir_dist/metrics.py	50.86% <32.03%> (-5.83%)	⬇️

... and 1 file with indirect coverage changes

🚀 New features to boost your workflow:

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

[grst]

Hi @felixpetschko,

thanks for working on this! From a technical perspective this looks all great.

My main concern here is that we keep adding metrics without providing guidelines (and evidence for) which metric to use.

In contrast to TCRdist, it is also not specific to TCR CDR3 sequences and can therefore be used for BCR analyses as well.

In what way is TCRdist specific to TCR sequences that would prevent it from using it for BCR? Or rather in what way is the alignment distance superior for BCR? My understanding would be that the main difference is that TCRdist allows only for a single gap position, while alignment allows for multiple, but does this make a big difference in practice?

Also how does it compare to TCRdist in terms of speed? E.g. how long would TCRdist take on the omniscope dataset on the same hardware?

[felixpetschko]

Hi @grst

My main concern here is that we keep adding metrics without providing guidelines (and evidence for) which metric to use.

My intention was rather to provide a faster implementation of the existing alignment distance metric than introducing a new metric. It is probably hard to provide real evidence for which metric to use. However, Needleman-Wunsch is a standard alignment algorithm and is also used in related sequence-analysis tools, so I think it is reasonable to use it.

In what way is TCRdist specific to TCR sequences that would prevent it from using it for BCR? Or rather in what way is the alignment distance superior for BCR? My understanding would be that the main difference is that TCRdist allows only for a single gap position, while alignment allows for multiple, but does this make a big difference in practice?

Actually, I was mainly focusing on performance, and I do not have proof for which metric is better in which case. My reasoning was that TCRdist's approach, with trimming from the N and C terminus and a single gap region, seems more targeted towards the TCR model. Needleman-Wunsch might be easier to justify for BCR CDR3 comparisons because it does not impose TCRdist's trimming and single-gap-region assumptions, and can handle length differences with a general global alignment.

Also how does it compare to TCRdist in terms of speed? E.g. how long would TCRdist take on the omniscope dataset on the same hardware?

TCRdist can run the full Omniscope COVID dataset on the same hardware in around 2.3 hours with default parameters, which makes it around 5 times faster than Needleman-Wunsch in my test run. The main reason is that, with Scirpy's default parameters, the gap position is computed by a formula and it is not necessary to try different gap positions. In contrast, Needleman-Wunsch computes the optimal global alignment that minimizes the distance.

[grst]

Do you think we could implement this without any user-facing changes?
I.e. integrate this into the AlignmentDistanceCalculator and fall back to the parasail implementation when gap_open != gap_extend?

(would be interesting to know if in the history of scirpy anyone has ever changed these default parameters. I'd guess not).

[felixpetschko]

Do you think we could implement this without any user-facing changes?

Yes, I will do that!

(would be interesting to know if in the history of scirpy anyone has ever changed these default parameters. I'd guess not).

However, I think we should definitely change the default params and set them in a way such that gaps are allowed. Otherwise there is no alignment done at all. I would allow at least 2 gaps such that it's worth to even run the dynamic programming alignment algorithm.

[grst]

Just changing defaults is also not very good practice... so we'd at least have to warn about it.
But maybe it is better to keep needleman_wunsch as a separate alignment metric and deprecate "alignment" and "fastalignment" distances.

We could complement it with a "metrics" guide in the documentation that explains the pros/cons and usecases of the different metrics.

[felixpetschko]

Alright, then let's do it like that 👍

[felixpetschko]

Now I marked alignment and fastalignment as deprecated in favor of needleman_wunsch and added a guide for choosing a sequence distance metric via #728.

[grst]

I have a few more comments, mostly on improving the organization of the metrics submodule.

[grst]

This would fail on joblib backends that do not support return_as="generator", e.g. dask.

See scirpy.utils._parallelize_with_joblib for a helper that addresses this.

[grst]

I'd move all these definitions to the top of the file, or maybe even better, a separate submodule within the ir_dist package.

[grst]

Also, maybe worth renaming this simply to AA_ALPHABET? Or is this still specific to parasail in any way?

[grst]

... or maybe make them dataclasses, that match together alphabet and substitution matrix

@dataclassclassSubstitutionMatrix:
alphabet: strmatrix: np.ndarrayBLOSUM62=SubstitutionMatrix(alphabet="ARN...", matrix=np.array([...]))
TCRBLOSUM_ALPHA=SubstitutionMatrix(...)

[grst]

Is there any reason for defining these as class variables instead of directly referencing the constants?

[grst]

Again, wouldn't it be easier to read if the constants were directly used everywhere?

[grst]

Could this become a metric-agnostic helper function? I think we already have similar code in other metrics...

[grst]

Suggested change

	If`gap_open == gap_extend`, useNeedlemanWunschDistanceCalculatorinstead.
	If`gap_open == gap_extend` (thedefault), useNeedlemanWunschDistanceCalculatorinstead, whichprovidesidenticalresultswhilebeingmuchfaster. Ifyouactuallyhaveause-caseforaffinegappenalties, pleaseletusknowbyopeninganissueonGitHub.

[grst]

Suggested change

	If`gap_open == gap_extend`, useNeedlemanWunschDistanceCalculatorinstead.
	If`gap_open == gap_extend` (thedefault), useNeedlemanWunschDistanceCalculatorinstead, whichprovidesidenticalresultswhilebeingmuchfaster. Ifyouactuallyhaveause-caseforaffinegappenalties, pleaseletusknowbyopeninganissueonGitHub.

[grst]

How has this been derived? Parasail?

[felixpetschko]

Yes, it was calculated with the deprecated AlignmentDistanceCalculator and the WU3k dataset (that only contains the 20 canonical amino acids).

[grst]

regarding deprecations, take a look at #735 please that switches to the decorators provided by scverse-misc.

[felixpetschko]

Hi @grst! I have implemented the requested changes - in particular, substitution matrices and the corresponding alphabets are now grouped in a dedicated _substitution_matrices submodule using a small dataclass, and the matrix mapping and substitution-to-distance conversion are provided as metric-independent helper functions. I also switched the block-wise parallelization to Scirpy’s joblib helper and updated the deprecations to use scverse-misc.

Besides that, I encountered a bug that was already present in the old AlignmentDistanceCalculator: With the ambiguous symbols in the blosum62 matrix (especially X) it is possible to get negative distance values, e.g. AAAXX vs AAAAA would result in a distance of -2. This situation occurs if there are values that score higher than self alignment. In the old implementation this causes datatype overflow or in the case of a distance of -1 ( +1 for storing = 0) the value is not stored at all due to the sparse format. In theory, one could just set negative result distance values to 0 (1 for storing). But logic-wise it seems a little bit odd to me that self-alignment cannot be interpreted as the maximum achievable alignment score. It also leads to problems with my performance optimization and setting the cutoff, because you cannot infer the number of allowed gaps by the cutoff and the gap penalty. This is due to the fact that parts of the sequence that score higher than self alignment can compensate for gap penalties. The symbols affected by this are X of the blosum62 matrix with ambiguous symbols and T from the TCRblosum beta matrix.

Therefore my suggestion would be to restrict the Needleman-Wunsch distance to the blosum62 matrix with only the 20 canonical amino acids (this is also the blosum62 matrix that can be chosen for TCRdist). In this case the matrix entries of the (unknown) ambiguous symbols would just be automatically set to 0 and the problems mentioned would be resolved. Currently, it is implemented in this PR like that. What do you think about that?

[felixpetschko]

Summary

This PR adds a Numba-optimized Needleman-Wunsch sequence distance metric via metric="needleman_wunsch". Conceptually, this metric is very similar to the existing alignment metric, with two main differences: it is much faster for large datasets, and it uses a linear gap model, where every gap position receives the same penalty.

The existing alignment metric allows gap openings and gap extensions to be penalized differently. However, this distinction is not used with the current default parameters anyway, where gap_open == gap_extend. If the alignment metric is configured such that gap_open == gap_extend == gap_penalty, both metrics should return equal results.

The reason for using a linear gap model is mostly practical. Supporting separate gap-open and gap-extension penalties requires handling three dynamic programming matrices instead of one for each sequence pair, which adds substantial computational overhead. For CDR3 sequences, which are usually quite short, I would generally expect users to choose parameters that allow only a limited number of gap positions. For example, in my test runs I used cutoff=10 and gap_penalty=4, which allows up to two gap positions. In this setting, distinguishing between two separate single-position gaps and one adjacent two-position gap is probably not worth the additional runtime. Therefore, I do not think the restriction to a single gap penalty is a major limitation.

The problem with the current alignment metric is that it is too slow for large datasets. The fastalignment metric improves runtime, but it can also suffer from performance limitations and does not always return exact results. Besides that, the current default parameters do not allow gaps (cutoff=10, gap_open=gap_extend=11). In that default setting, only equal-length sequences can fall below the cutoff, so the metric effectively computes alignment distances between equal-length sequences without making use of meaningful gap placement.

Therefore, I implemented the Needleman-Wunsch distance in a similar style to the Numba-optimized TCRdist CPU implementation. It supports BLOSUM62 by default and can also use TCRBLOSUM alpha/beta matrices through base_matrix="tcrblosum". With the improved performance, I was able to run the 8 million-cell Omniscope COVID dataset with 64 CPU cores using cutoff=10 and gap_penalty=4 within around 12 hours.

Overall, I think this metric is a useful addition because it provides an exact alignment-based distance that can still handle larger datasets. It is more flexible than the Hamming distance, but substantially faster than the existing alignment metric. In contrast to TCRdist, it is also not specific to TCR CDR3 sequences and can therefore be used for BCR analyses as well.

The most useful default values for gap_penalty and cutoff are open for discussion, since the performance improvements make less restrictive parameter choices feasible. It could also be discussed whether the existing alignment metric is still needed if its additional gap parameterization and support for additional Parasail substitution matrices are not required.

Main changes

• Add NeedlemanWunschDistanceCalculator and expose it via metric="needleman_wunsch" in sequence_dist and ir_dist.
• Support BLOSUM62 and TCRBLOSUM alpha/beta substitution matrices in a similar way as the TCRdist metric.
• Implement a Numba-based Needleman-Wunsch dynamic programming algorithm with cutoff-based optimizations.
• Add test cases.

[codecov]

Codecov Report

❌ Patch coverage is 49.13295% with 88 lines in your changes missing coverage. Please review.
✅ Project coverage is 77.32%. Comparing base (79a8344) to head (01ceae6).

Files with missing lines	Patch %	Lines
src/scirpy/ir_dist/metrics.py	32.03%	87 Missing ⚠️
src/scirpy/ir_dist/_substitution_matrices.py	97.67%	1 Missing ⚠️

Additional details and impacted files

@@ Coverage Diff @@## main #725 +/- ##
==========================================
- Coverage 78.44% 77.32% -1.13% 
==========================================
Files 51 52 +1 Lines 4635 4785 +150 ==========================================
+ Hits 3636 3700 +64 - Misses 999 1085 +86 

Files with missing lines	Coverage Δ
src/scirpy/ir_dist/__init__.py	93.13% <100.00%> (+0.13%)	⬆️
src/scirpy/ir_dist/_substitution_matrices.py	97.67% <97.67%> (ø)
src/scirpy/ir_dist/metrics.py	50.86% <32.03%> (-5.83%)	⬇️

... and 1 file with indirect coverage changes

🚀 New features to boost your workflow:

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

[grst]

Hi @felixpetschko,

thanks for working on this! From a technical perspective this looks all great.

My main concern here is that we keep adding metrics without providing guidelines (and evidence for) which metric to use.

In contrast to TCRdist, it is also not specific to TCR CDR3 sequences and can therefore be used for BCR analyses as well.

In what way is TCRdist specific to TCR sequences that would prevent it from using it for BCR? Or rather in what way is the alignment distance superior for BCR? My understanding would be that the main difference is that TCRdist allows only for a single gap position, while alignment allows for multiple, but does this make a big difference in practice?

Also how does it compare to TCRdist in terms of speed? E.g. how long would TCRdist take on the omniscope dataset on the same hardware?

[felixpetschko]

Hi @grst

My main concern here is that we keep adding metrics without providing guidelines (and evidence for) which metric to use.

My intention was rather to provide a faster implementation of the existing alignment distance metric than introducing a new metric. It is probably hard to provide real evidence for which metric to use. However, Needleman-Wunsch is a standard alignment algorithm and is also used in related sequence-analysis tools, so I think it is reasonable to use it.

In what way is TCRdist specific to TCR sequences that would prevent it from using it for BCR? Or rather in what way is the alignment distance superior for BCR? My understanding would be that the main difference is that TCRdist allows only for a single gap position, while alignment allows for multiple, but does this make a big difference in practice?

Actually, I was mainly focusing on performance, and I do not have proof for which metric is better in which case. My reasoning was that TCRdist's approach, with trimming from the N and C terminus and a single gap region, seems more targeted towards the TCR model. Needleman-Wunsch might be easier to justify for BCR CDR3 comparisons because it does not impose TCRdist's trimming and single-gap-region assumptions, and can handle length differences with a general global alignment.

Also how does it compare to TCRdist in terms of speed? E.g. how long would TCRdist take on the omniscope dataset on the same hardware?

TCRdist can run the full Omniscope COVID dataset on the same hardware in around 2.3 hours with default parameters, which makes it around 5 times faster than Needleman-Wunsch in my test run. The main reason is that, with Scirpy's default parameters, the gap position is computed by a formula and it is not necessary to try different gap positions. In contrast, Needleman-Wunsch computes the optimal global alignment that minimizes the distance.

[grst]

Do you think we could implement this without any user-facing changes?
I.e. integrate this into the AlignmentDistanceCalculator and fall back to the parasail implementation when gap_open != gap_extend?

(would be interesting to know if in the history of scirpy anyone has ever changed these default parameters. I'd guess not).

[felixpetschko]

Do you think we could implement this without any user-facing changes?

Yes, I will do that!

(would be interesting to know if in the history of scirpy anyone has ever changed these default parameters. I'd guess not).

However, I think we should definitely change the default params and set them in a way such that gaps are allowed. Otherwise there is no alignment done at all. I would allow at least 2 gaps such that it's worth to even run the dynamic programming alignment algorithm.

[grst]

Just changing defaults is also not very good practice... so we'd at least have to warn about it.
But maybe it is better to keep needleman_wunsch as a separate alignment metric and deprecate "alignment" and "fastalignment" distances.

We could complement it with a "metrics" guide in the documentation that explains the pros/cons and usecases of the different metrics.

[felixpetschko]

Alright, then let's do it like that 👍

[felixpetschko]

Now I marked alignment and fastalignment as deprecated in favor of needleman_wunsch and added a guide for choosing a sequence distance metric via #728.

[grst]

I have a few more comments, mostly on improving the organization of the metrics submodule.

[grst]

This would fail on joblib backends that do not support return_as="generator", e.g. dask.

See scirpy.utils._parallelize_with_joblib for a helper that addresses this.

[grst]

I'd move all these definitions to the top of the file, or maybe even better, a separate submodule within the ir_dist package.

[grst]

Also, maybe worth renaming this simply to AA_ALPHABET? Or is this still specific to parasail in any way?

[grst]

... or maybe make them dataclasses, that match together alphabet and substitution matrix

@dataclassclassSubstitutionMatrix:
alphabet: strmatrix: np.ndarrayBLOSUM62=SubstitutionMatrix(alphabet="ARN...", matrix=np.array([...]))
TCRBLOSUM_ALPHA=SubstitutionMatrix(...)

[grst]

Is there any reason for defining these as class variables instead of directly referencing the constants?

[grst]

Again, wouldn't it be easier to read if the constants were directly used everywhere?

[grst]

Could this become a metric-agnostic helper function? I think we already have similar code in other metrics...

[grst]

Suggested change

	If`gap_open == gap_extend`, useNeedlemanWunschDistanceCalculatorinstead.
	If`gap_open == gap_extend` (thedefault), useNeedlemanWunschDistanceCalculatorinstead, whichprovidesidenticalresultswhilebeingmuchfaster. Ifyouactuallyhaveause-caseforaffinegappenalties, pleaseletusknowbyopeninganissueonGitHub.

[grst]

Suggested change

	If`gap_open == gap_extend`, useNeedlemanWunschDistanceCalculatorinstead.
	If`gap_open == gap_extend` (thedefault), useNeedlemanWunschDistanceCalculatorinstead, whichprovidesidenticalresultswhilebeingmuchfaster. Ifyouactuallyhaveause-caseforaffinegappenalties, pleaseletusknowbyopeninganissueonGitHub.

[grst]

How has this been derived? Parasail?

[felixpetschko]

Yes, it was calculated with the deprecated AlignmentDistanceCalculator and the WU3k dataset (that only contains the 20 canonical amino acids).

[grst]

regarding deprecations, take a look at #735 please that switches to the decorators provided by scverse-misc.

[felixpetschko]

Hi @grst! I have implemented the requested changes - in particular, substitution matrices and the corresponding alphabets are now grouped in a dedicated _substitution_matrices submodule using a small dataclass, and the matrix mapping and substitution-to-distance conversion are provided as metric-independent helper functions. I also switched the block-wise parallelization to Scirpy’s joblib helper and updated the deprecations to use scverse-misc.

Besides that, I encountered a bug that was already present in the old AlignmentDistanceCalculator: With the ambiguous symbols in the blosum62 matrix (especially X) it is possible to get negative distance values, e.g. AAAXX vs AAAAA would result in a distance of -2. This situation occurs if there are values that score higher than self alignment. In the old implementation this causes datatype overflow or in the case of a distance of -1 ( +1 for storing = 0) the value is not stored at all due to the sparse format. In theory, one could just set negative result distance values to 0 (1 for storing). But logic-wise it seems a little bit odd to me that self-alignment cannot be interpreted as the maximum achievable alignment score. It also leads to problems with my performance optimization and setting the cutoff, because you cannot infer the number of allowed gaps by the cutoff and the gap penalty. This is due to the fact that parts of the sequence that score higher than self alignment can compensate for gap penalties. The symbols affected by this are X of the blosum62 matrix with ambiguous symbols and T from the TCRblosum beta matrix.

Therefore my suggestion would be to restrict the Needleman-Wunsch distance to the blosum62 matrix with only the 20 canonical amino acids (this is also the blosum62 matrix that can be chosen for TCRdist). In this case the matrix entries of the (unknown) ambiguous symbols would just be automatically set to 0 and the problems mentioned would be resolved. Currently, it is implemented in this PR like that. What do you think about that?

[felixpetschko]

Summary

This PR adds a Numba-optimized Needleman-Wunsch sequence distance metric via metric="needleman_wunsch". Conceptually, this metric is very similar to the existing alignment metric, with two main differences: it is much faster for large datasets, and it uses a linear gap model, where every gap position receives the same penalty.

The existing alignment metric allows gap openings and gap extensions to be penalized differently. However, this distinction is not used with the current default parameters anyway, where gap_open == gap_extend. If the alignment metric is configured such that gap_open == gap_extend == gap_penalty, both metrics should return equal results.

The reason for using a linear gap model is mostly practical. Supporting separate gap-open and gap-extension penalties requires handling three dynamic programming matrices instead of one for each sequence pair, which adds substantial computational overhead. For CDR3 sequences, which are usually quite short, I would generally expect users to choose parameters that allow only a limited number of gap positions. For example, in my test runs I used cutoff=10 and gap_penalty=4, which allows up to two gap positions. In this setting, distinguishing between two separate single-position gaps and one adjacent two-position gap is probably not worth the additional runtime. Therefore, I do not think the restriction to a single gap penalty is a major limitation.

The problem with the current alignment metric is that it is too slow for large datasets. The fastalignment metric improves runtime, but it can also suffer from performance limitations and does not always return exact results. Besides that, the current default parameters do not allow gaps (cutoff=10, gap_open=gap_extend=11). In that default setting, only equal-length sequences can fall below the cutoff, so the metric effectively computes alignment distances between equal-length sequences without making use of meaningful gap placement.

Therefore, I implemented the Needleman-Wunsch distance in a similar style to the Numba-optimized TCRdist CPU implementation. It supports BLOSUM62 by default and can also use TCRBLOSUM alpha/beta matrices through base_matrix="tcrblosum". With the improved performance, I was able to run the 8 million-cell Omniscope COVID dataset with 64 CPU cores using cutoff=10 and gap_penalty=4 within around 12 hours.

Overall, I think this metric is a useful addition because it provides an exact alignment-based distance that can still handle larger datasets. It is more flexible than the Hamming distance, but substantially faster than the existing alignment metric. In contrast to TCRdist, it is also not specific to TCR CDR3 sequences and can therefore be used for BCR analyses as well.

The most useful default values for gap_penalty and cutoff are open for discussion, since the performance improvements make less restrictive parameter choices feasible. It could also be discussed whether the existing alignment metric is still needed if its additional gap parameterization and support for additional Parasail substitution matrices are not required.

Main changes

• Add NeedlemanWunschDistanceCalculator and expose it via metric="needleman_wunsch" in sequence_dist and ir_dist.
• Support BLOSUM62 and TCRBLOSUM alpha/beta substitution matrices in a similar way as the TCRdist metric.
• Implement a Numba-based Needleman-Wunsch dynamic programming algorithm with cutoff-based optimizations.
• Add test cases.

[codecov]

Codecov Report

❌ Patch coverage is 49.13295% with 88 lines in your changes missing coverage. Please review.
✅ Project coverage is 77.32%. Comparing base (79a8344) to head (01ceae6).

Files with missing lines	Patch %	Lines
src/scirpy/ir_dist/metrics.py	32.03%	87 Missing ⚠️
src/scirpy/ir_dist/_substitution_matrices.py	97.67%	1 Missing ⚠️

Additional details and impacted files

@@ Coverage Diff @@## main #725 +/- ##
==========================================
- Coverage 78.44% 77.32% -1.13% 
==========================================
Files 51 52 +1 Lines 4635 4785 +150 ==========================================
+ Hits 3636 3700 +64 - Misses 999 1085 +86 

Files with missing lines	Coverage Δ
src/scirpy/ir_dist/__init__.py	93.13% <100.00%> (+0.13%)	⬆️
src/scirpy/ir_dist/_substitution_matrices.py	97.67% <97.67%> (ø)
src/scirpy/ir_dist/metrics.py	50.86% <32.03%> (-5.83%)	⬇️

... and 1 file with indirect coverage changes

🚀 New features to boost your workflow:

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

[grst]

Hi @felixpetschko,

thanks for working on this! From a technical perspective this looks all great.

My main concern here is that we keep adding metrics without providing guidelines (and evidence for) which metric to use.

In contrast to TCRdist, it is also not specific to TCR CDR3 sequences and can therefore be used for BCR analyses as well.

In what way is TCRdist specific to TCR sequences that would prevent it from using it for BCR? Or rather in what way is the alignment distance superior for BCR? My understanding would be that the main difference is that TCRdist allows only for a single gap position, while alignment allows for multiple, but does this make a big difference in practice?

Also how does it compare to TCRdist in terms of speed? E.g. how long would TCRdist take on the omniscope dataset on the same hardware?

[felixpetschko]

Hi @grst

My main concern here is that we keep adding metrics without providing guidelines (and evidence for) which metric to use.

My intention was rather to provide a faster implementation of the existing alignment distance metric than introducing a new metric. It is probably hard to provide real evidence for which metric to use. However, Needleman-Wunsch is a standard alignment algorithm and is also used in related sequence-analysis tools, so I think it is reasonable to use it.

In what way is TCRdist specific to TCR sequences that would prevent it from using it for BCR? Or rather in what way is the alignment distance superior for BCR? My understanding would be that the main difference is that TCRdist allows only for a single gap position, while alignment allows for multiple, but does this make a big difference in practice?

Actually, I was mainly focusing on performance, and I do not have proof for which metric is better in which case. My reasoning was that TCRdist's approach, with trimming from the N and C terminus and a single gap region, seems more targeted towards the TCR model. Needleman-Wunsch might be easier to justify for BCR CDR3 comparisons because it does not impose TCRdist's trimming and single-gap-region assumptions, and can handle length differences with a general global alignment.

Also how does it compare to TCRdist in terms of speed? E.g. how long would TCRdist take on the omniscope dataset on the same hardware?

TCRdist can run the full Omniscope COVID dataset on the same hardware in around 2.3 hours with default parameters, which makes it around 5 times faster than Needleman-Wunsch in my test run. The main reason is that, with Scirpy's default parameters, the gap position is computed by a formula and it is not necessary to try different gap positions. In contrast, Needleman-Wunsch computes the optimal global alignment that minimizes the distance.

[grst]

Do you think we could implement this without any user-facing changes?
I.e. integrate this into the AlignmentDistanceCalculator and fall back to the parasail implementation when gap_open != gap_extend?

(would be interesting to know if in the history of scirpy anyone has ever changed these default parameters. I'd guess not).

[felixpetschko]

Do you think we could implement this without any user-facing changes?

Yes, I will do that!

(would be interesting to know if in the history of scirpy anyone has ever changed these default parameters. I'd guess not).

However, I think we should definitely change the default params and set them in a way such that gaps are allowed. Otherwise there is no alignment done at all. I would allow at least 2 gaps such that it's worth to even run the dynamic programming alignment algorithm.

[grst]

Just changing defaults is also not very good practice... so we'd at least have to warn about it.
But maybe it is better to keep needleman_wunsch as a separate alignment metric and deprecate "alignment" and "fastalignment" distances.

We could complement it with a "metrics" guide in the documentation that explains the pros/cons and usecases of the different metrics.

[felixpetschko]

Alright, then let's do it like that 👍

[felixpetschko]

Now I marked alignment and fastalignment as deprecated in favor of needleman_wunsch and added a guide for choosing a sequence distance metric via #728.

[grst]

I have a few more comments, mostly on improving the organization of the metrics submodule.

[grst]

This would fail on joblib backends that do not support return_as="generator", e.g. dask.

See scirpy.utils._parallelize_with_joblib for a helper that addresses this.

[grst]

I'd move all these definitions to the top of the file, or maybe even better, a separate submodule within the ir_dist package.

[grst]

Also, maybe worth renaming this simply to AA_ALPHABET? Or is this still specific to parasail in any way?

[grst]

... or maybe make them dataclasses, that match together alphabet and substitution matrix

@dataclassclassSubstitutionMatrix:
alphabet: strmatrix: np.ndarrayBLOSUM62=SubstitutionMatrix(alphabet="ARN...", matrix=np.array([...]))
TCRBLOSUM_ALPHA=SubstitutionMatrix(...)

[grst]

Is there any reason for defining these as class variables instead of directly referencing the constants?

[grst]

Again, wouldn't it be easier to read if the constants were directly used everywhere?

[grst]

Could this become a metric-agnostic helper function? I think we already have similar code in other metrics...

[grst]

Suggested change

	If`gap_open == gap_extend`, useNeedlemanWunschDistanceCalculatorinstead.
	If`gap_open == gap_extend` (thedefault), useNeedlemanWunschDistanceCalculatorinstead, whichprovidesidenticalresultswhilebeingmuchfaster. Ifyouactuallyhaveause-caseforaffinegappenalties, pleaseletusknowbyopeninganissueonGitHub.

[grst]

Suggested change

	If`gap_open == gap_extend`, useNeedlemanWunschDistanceCalculatorinstead.
	If`gap_open == gap_extend` (thedefault), useNeedlemanWunschDistanceCalculatorinstead, whichprovidesidenticalresultswhilebeingmuchfaster. Ifyouactuallyhaveause-caseforaffinegappenalties, pleaseletusknowbyopeninganissueonGitHub.

[grst]

How has this been derived? Parasail?

[felixpetschko]

Yes, it was calculated with the deprecated AlignmentDistanceCalculator and the WU3k dataset (that only contains the 20 canonical amino acids).

[grst]

regarding deprecations, take a look at #735 please that switches to the decorators provided by scverse-misc.

[felixpetschko]

Hi @grst! I have implemented the requested changes - in particular, substitution matrices and the corresponding alphabets are now grouped in a dedicated _substitution_matrices submodule using a small dataclass, and the matrix mapping and substitution-to-distance conversion are provided as metric-independent helper functions. I also switched the block-wise parallelization to Scirpy’s joblib helper and updated the deprecations to use scverse-misc.

Besides that, I encountered a bug that was already present in the old AlignmentDistanceCalculator: With the ambiguous symbols in the blosum62 matrix (especially X) it is possible to get negative distance values, e.g. AAAXX vs AAAAA would result in a distance of -2. This situation occurs if there are values that score higher than self alignment. In the old implementation this causes datatype overflow or in the case of a distance of -1 ( +1 for storing = 0) the value is not stored at all due to the sparse format. In theory, one could just set negative result distance values to 0 (1 for storing). But logic-wise it seems a little bit odd to me that self-alignment cannot be interpreted as the maximum achievable alignment score. It also leads to problems with my performance optimization and setting the cutoff, because you cannot infer the number of allowed gaps by the cutoff and the gap penalty. This is due to the fact that parts of the sequence that score higher than self alignment can compensate for gap penalties. The symbols affected by this are X of the blosum62 matrix with ambiguous symbols and T from the TCRblosum beta matrix.

Therefore my suggestion would be to restrict the Needleman-Wunsch distance to the blosum62 matrix with only the 20 canonical amino acids (this is also the blosum62 matrix that can be chosen for TCRdist). In this case the matrix entries of the (unknown) ambiguous symbols would just be automatically set to 0 and the problems mentioned would be resolved. Currently, it is implemented in this PR like that. What do you think about that?

[felixpetschko]

Summary

This PR adds a Numba-optimized Needleman-Wunsch sequence distance metric via metric="needleman_wunsch". Conceptually, this metric is very similar to the existing alignment metric, with two main differences: it is much faster for large datasets, and it uses a linear gap model, where every gap position receives the same penalty.

The existing alignment metric allows gap openings and gap extensions to be penalized differently. However, this distinction is not used with the current default parameters anyway, where gap_open == gap_extend. If the alignment metric is configured such that gap_open == gap_extend == gap_penalty, both metrics should return equal results.

The reason for using a linear gap model is mostly practical. Supporting separate gap-open and gap-extension penalties requires handling three dynamic programming matrices instead of one for each sequence pair, which adds substantial computational overhead. For CDR3 sequences, which are usually quite short, I would generally expect users to choose parameters that allow only a limited number of gap positions. For example, in my test runs I used cutoff=10 and gap_penalty=4, which allows up to two gap positions. In this setting, distinguishing between two separate single-position gaps and one adjacent two-position gap is probably not worth the additional runtime. Therefore, I do not think the restriction to a single gap penalty is a major limitation.

The problem with the current alignment metric is that it is too slow for large datasets. The fastalignment metric improves runtime, but it can also suffer from performance limitations and does not always return exact results. Besides that, the current default parameters do not allow gaps (cutoff=10, gap_open=gap_extend=11). In that default setting, only equal-length sequences can fall below the cutoff, so the metric effectively computes alignment distances between equal-length sequences without making use of meaningful gap placement.

Therefore, I implemented the Needleman-Wunsch distance in a similar style to the Numba-optimized TCRdist CPU implementation. It supports BLOSUM62 by default and can also use TCRBLOSUM alpha/beta matrices through base_matrix="tcrblosum". With the improved performance, I was able to run the 8 million-cell Omniscope COVID dataset with 64 CPU cores using cutoff=10 and gap_penalty=4 within around 12 hours.

Overall, I think this metric is a useful addition because it provides an exact alignment-based distance that can still handle larger datasets. It is more flexible than the Hamming distance, but substantially faster than the existing alignment metric. In contrast to TCRdist, it is also not specific to TCR CDR3 sequences and can therefore be used for BCR analyses as well.

The most useful default values for gap_penalty and cutoff are open for discussion, since the performance improvements make less restrictive parameter choices feasible. It could also be discussed whether the existing alignment metric is still needed if its additional gap parameterization and support for additional Parasail substitution matrices are not required.

Main changes

• Add NeedlemanWunschDistanceCalculator and expose it via metric="needleman_wunsch" in sequence_dist and ir_dist.
• Support BLOSUM62 and TCRBLOSUM alpha/beta substitution matrices in a similar way as the TCRdist metric.
• Implement a Numba-based Needleman-Wunsch dynamic programming algorithm with cutoff-based optimizations.
• Add test cases.

[codecov]

Codecov Report

❌ Patch coverage is 49.13295% with 88 lines in your changes missing coverage. Please review.
✅ Project coverage is 77.32%. Comparing base (79a8344) to head (01ceae6).

Files with missing lines	Patch %	Lines
src/scirpy/ir_dist/metrics.py	32.03%	87 Missing ⚠️
src/scirpy/ir_dist/_substitution_matrices.py	97.67%	1 Missing ⚠️

Additional details and impacted files

@@ Coverage Diff @@## main #725 +/- ##
==========================================
- Coverage 78.44% 77.32% -1.13% 
==========================================
Files 51 52 +1 Lines 4635 4785 +150 ==========================================
+ Hits 3636 3700 +64 - Misses 999 1085 +86 

Files with missing lines	Coverage Δ
src/scirpy/ir_dist/__init__.py	93.13% <100.00%> (+0.13%)	⬆️
src/scirpy/ir_dist/_substitution_matrices.py	97.67% <97.67%> (ø)
src/scirpy/ir_dist/metrics.py	50.86% <32.03%> (-5.83%)	⬇️

... and 1 file with indirect coverage changes

🚀 New features to boost your workflow:

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

[grst]

Hi @felixpetschko,

thanks for working on this! From a technical perspective this looks all great.

My main concern here is that we keep adding metrics without providing guidelines (and evidence for) which metric to use.

In contrast to TCRdist, it is also not specific to TCR CDR3 sequences and can therefore be used for BCR analyses as well.

In what way is TCRdist specific to TCR sequences that would prevent it from using it for BCR? Or rather in what way is the alignment distance superior for BCR? My understanding would be that the main difference is that TCRdist allows only for a single gap position, while alignment allows for multiple, but does this make a big difference in practice?

Also how does it compare to TCRdist in terms of speed? E.g. how long would TCRdist take on the omniscope dataset on the same hardware?

[felixpetschko]

Hi @grst

My main concern here is that we keep adding metrics without providing guidelines (and evidence for) which metric to use.

My intention was rather to provide a faster implementation of the existing alignment distance metric than introducing a new metric. It is probably hard to provide real evidence for which metric to use. However, Needleman-Wunsch is a standard alignment algorithm and is also used in related sequence-analysis tools, so I think it is reasonable to use it.

In what way is TCRdist specific to TCR sequences that would prevent it from using it for BCR? Or rather in what way is the alignment distance superior for BCR? My understanding would be that the main difference is that TCRdist allows only for a single gap position, while alignment allows for multiple, but does this make a big difference in practice?

Actually, I was mainly focusing on performance, and I do not have proof for which metric is better in which case. My reasoning was that TCRdist's approach, with trimming from the N and C terminus and a single gap region, seems more targeted towards the TCR model. Needleman-Wunsch might be easier to justify for BCR CDR3 comparisons because it does not impose TCRdist's trimming and single-gap-region assumptions, and can handle length differences with a general global alignment.

Also how does it compare to TCRdist in terms of speed? E.g. how long would TCRdist take on the omniscope dataset on the same hardware?

TCRdist can run the full Omniscope COVID dataset on the same hardware in around 2.3 hours with default parameters, which makes it around 5 times faster than Needleman-Wunsch in my test run. The main reason is that, with Scirpy's default parameters, the gap position is computed by a formula and it is not necessary to try different gap positions. In contrast, Needleman-Wunsch computes the optimal global alignment that minimizes the distance.

[grst]

Do you think we could implement this without any user-facing changes?
I.e. integrate this into the AlignmentDistanceCalculator and fall back to the parasail implementation when gap_open != gap_extend?

(would be interesting to know if in the history of scirpy anyone has ever changed these default parameters. I'd guess not).

[felixpetschko]

Do you think we could implement this without any user-facing changes?

Yes, I will do that!

(would be interesting to know if in the history of scirpy anyone has ever changed these default parameters. I'd guess not).

However, I think we should definitely change the default params and set them in a way such that gaps are allowed. Otherwise there is no alignment done at all. I would allow at least 2 gaps such that it's worth to even run the dynamic programming alignment algorithm.

[grst]

Just changing defaults is also not very good practice... so we'd at least have to warn about it.
But maybe it is better to keep needleman_wunsch as a separate alignment metric and deprecate "alignment" and "fastalignment" distances.

We could complement it with a "metrics" guide in the documentation that explains the pros/cons and usecases of the different metrics.

[felixpetschko]

Alright, then let's do it like that 👍

[felixpetschko]

Now I marked alignment and fastalignment as deprecated in favor of needleman_wunsch and added a guide for choosing a sequence distance metric via #728.

[grst]

I have a few more comments, mostly on improving the organization of the metrics submodule.

[grst]

This would fail on joblib backends that do not support return_as="generator", e.g. dask.

See scirpy.utils._parallelize_with_joblib for a helper that addresses this.

[grst]

I'd move all these definitions to the top of the file, or maybe even better, a separate submodule within the ir_dist package.

[grst]

Also, maybe worth renaming this simply to AA_ALPHABET? Or is this still specific to parasail in any way?

[grst]

... or maybe make them dataclasses, that match together alphabet and substitution matrix

@dataclassclassSubstitutionMatrix:
alphabet: strmatrix: np.ndarrayBLOSUM62=SubstitutionMatrix(alphabet="ARN...", matrix=np.array([...]))
TCRBLOSUM_ALPHA=SubstitutionMatrix(...)

[grst]

Is there any reason for defining these as class variables instead of directly referencing the constants?

[grst]

Again, wouldn't it be easier to read if the constants were directly used everywhere?

[grst]

Could this become a metric-agnostic helper function? I think we already have similar code in other metrics...

[grst]

Suggested change

	If`gap_open == gap_extend`, useNeedlemanWunschDistanceCalculatorinstead.
	If`gap_open == gap_extend` (thedefault), useNeedlemanWunschDistanceCalculatorinstead, whichprovidesidenticalresultswhilebeingmuchfaster. Ifyouactuallyhaveause-caseforaffinegappenalties, pleaseletusknowbyopeninganissueonGitHub.

[grst]

Suggested change

	If`gap_open == gap_extend`, useNeedlemanWunschDistanceCalculatorinstead.
	If`gap_open == gap_extend` (thedefault), useNeedlemanWunschDistanceCalculatorinstead, whichprovidesidenticalresultswhilebeingmuchfaster. Ifyouactuallyhaveause-caseforaffinegappenalties, pleaseletusknowbyopeninganissueonGitHub.

[grst]

How has this been derived? Parasail?

[felixpetschko]

Yes, it was calculated with the deprecated AlignmentDistanceCalculator and the WU3k dataset (that only contains the 20 canonical amino acids).

[grst]

regarding deprecations, take a look at #735 please that switches to the decorators provided by scverse-misc.

[felixpetschko]

Hi @grst! I have implemented the requested changes - in particular, substitution matrices and the corresponding alphabets are now grouped in a dedicated _substitution_matrices submodule using a small dataclass, and the matrix mapping and substitution-to-distance conversion are provided as metric-independent helper functions. I also switched the block-wise parallelization to Scirpy’s joblib helper and updated the deprecations to use scverse-misc.

Besides that, I encountered a bug that was already present in the old AlignmentDistanceCalculator: With the ambiguous symbols in the blosum62 matrix (especially X) it is possible to get negative distance values, e.g. AAAXX vs AAAAA would result in a distance of -2. This situation occurs if there are values that score higher than self alignment. In the old implementation this causes datatype overflow or in the case of a distance of -1 ( +1 for storing = 0) the value is not stored at all due to the sparse format. In theory, one could just set negative result distance values to 0 (1 for storing). But logic-wise it seems a little bit odd to me that self-alignment cannot be interpreted as the maximum achievable alignment score. It also leads to problems with my performance optimization and setting the cutoff, because you cannot infer the number of allowed gaps by the cutoff and the gap penalty. This is due to the fact that parts of the sequence that score higher than self alignment can compensate for gap penalties. The symbols affected by this are X of the blosum62 matrix with ambiguous symbols and T from the TCRblosum beta matrix.

Therefore my suggestion would be to restrict the Needleman-Wunsch distance to the blosum62 matrix with only the 20 canonical amino acids (this is also the blosum62 matrix that can be chosen for TCRdist). In this case the matrix entries of the (unknown) ambiguous symbols would just be automatically set to 0 and the problems mentioned would be resolved. Currently, it is implemented in this PR like that. What do you think about that?