Mass alignment

If you want to know how two peptidoforms/proteoforms are related and want to keep all mass-based mistakes in mind you need to align the peptidoforms with mass-alignment[^1].

Library features

Align two peptidoforms
With custom scoring
Global, local, either-global and mixtures of these all
Indexed alignment to speed up repeated alignments against the same database
Consecutive alignment (only when the imgt feature is turned on) to create a DomainGapAlign[^2] like alignment with the mass-based alignment

The alignment itself

A mass-based alignment handles the case in which multiple amino acids are wrong, but the total mass of this set of amino acids is equal to the mass of a set of different amino acids on the other peptide. This is quite common in mass spectrometry where mistakes based on mass coincidences are very common. For example N has the same mass as GG, so if we want to make a mass spectrometry faithful alignment of ANA with AGGA the result should reflect this fact:

Generated using this algorithm bound to a cli tool: https://github.com/snijderlab/align-cli

use mzcore::{prelude::*, sequence::SimpleLinear};
use mzalign::*;
// For other use cases Ontologies::init() makes more sense.
let ontologies = &mzcore::ontology::STATIC_ONTOLOGIES;
let a = Peptidoform::pro_forma("ANA", ontologies).unwrap().0.into_simple_linear().unwrap();
let b = Peptidoform::pro_forma("AGGA", ontologies).unwrap().0.into_simple_linear().unwrap();
let alignment = align::<4, &Peptidoform<SimpleLinear>, &Peptidoform<SimpleLinear>>(&a, &b, AlignScoring::default(), AlignType::GLOBAL);
assert_eq!(alignment.short(), "1=1:2i1=");
assert_eq!(alignment.ppm().value, 0.0);

This extends the Smith-Waterman/Needleman-Wunsch alignment states (Match, Mismatch, Insertion, and Deletion) with the following additional states:

IndentityMassMismatch: the same amino acid but the mass of the sequence element is different, this happens when either of the options has a modification that is not present on the other, this could be an error or could be fine depending on the use of the alignment, use PairMode to control this behaviour.
Rotation: the same amino acids but in a different order, example: AHK on KAH.
Isobaric: different amino acids (or different modifications) that have the same mass, example N on GG, or M[Oxidation] on F.

As is visible in the examples the Rotation errors can be longer than just one amino acids, and Isobaric errors can even be of different lengths for the two peptidoforms. This means that these errors need special care to be visualised properly. Also because these need to handle these different lengths the main loop of the alignment needs to loop over all combinations of lengths from 1 to the maximum chosen length. This means that algorithmically speaking this alignment is slower than SW/NW.

AlignTypes

Alignments can be done global (as is Needleman-Wunsch) or local (as is Smith-Waterman). Global means that both peptidoforms have to be fully matched, local means that both peptidoforms can have unmatched sequences on both sides of the match. Either global is a variation that enforces that at least one of the two peptidoforms has to be fully matched or said inversely only one peptidoform can have unmatched sequence. This can be used when it is known that sequences are related but neither is known a priori to be a full superset of the other sequence. For example if you are matching a peptide to a database and the database does not contain the full protein sequence but only a part of the sequence. As is common in antibodies, those are built from three/four separate genes spliced together.

Global alignment, the full peptidoforms are matched even if that match is quite poor

Local alignment, the best parts of the peptidoforms are aligned

Either global, for both sides (left and right) only one peptidoform can have unmatched sequence

In this crate AlignType can be used to set these types of alignments. Note that this allows setting one of these types per side and per peptidoform. This allows for encoding as much knowledge into the alignment as possible. For example if one would build a tool to align an antibody sequence to germline V genes you could make it fixed to a global left alignment (meaning both sequences have to start together) and end with an either global right alignment if it is not known if the user will supply full antibody sequences or might even supply truncated V gene sequences as shown below.

Aligning a V gene with a full antibody sequence

Aligning a V gene with a beginning and end truncated antibody sequence, showing that the missing start is scored negatively as a leading deletion but the missing end is scored neutral as an allowed truncation

Example usage

# fn main() -> Result<(), Vec<context_error::BoxedError<'static, context_error::BasicKind>>> {
use mzcore::{prelude::*, sequence::SimpleLinear};
use mzalign::{prelude::*, align};
// For other use cases Ontologies::init() makes more sense.
let ontologies = &mzcore::ontology::STATIC_ONTOLOGIES;
// Check how this peptide compares to a similar peptide (using the feature `align`)
let first_peptide = Peptidoform::pro_forma("IVQEVT", ontologies)?.0.into_simple_linear().unwrap();
let second_peptide = Peptidoform::pro_forma("LVQVET", ontologies)?.0.into_simple_linear().unwrap();
// Align the two peptides using mass based alignment
// IVQEVT A
// LVQVET B
// ─  ╶╴
let alignment = align::<4, &Peptidoform<SimpleLinear>, &Peptidoform<SimpleLinear>>(
                  &first_peptide,
                  &second_peptide,
                  AlignScoring::default(),
                  AlignType::GLOBAL);
# dbg!(&alignment);
// Calculate some more statistics on this alignment
let stats = alignment.stats();
assert_eq!(stats.mass_similar, 6); // 6 out of the 6 positions are mass similar
# Ok(()) }

Compilation features

imgt - enables access to the IMGT database of antibodies germline sequences, with annotations. This also turns on the use of consecutive alignments.
rayon - enables parallel iterators using rayon, this also turns on rayon in imgt.