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¹.

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² 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:

Identity: 0.500 (2/4), Similarity: 0.750 (3/4), Gaps: 0.000 (0/4), Score: 0.706 (12/17), Equal mass, Tolerance: 10 ppm, Alignment: global
Start: A 0 B 0, Path: 1=1:2i1=
AN·A A
AGGA B
 ╶╴

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:

1. 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.
2. Rotation: the same amino acids but in a different order, example: AHK on KAH.
3. 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.

Identity: 0.176 (3/17), Mass similarity: 0.353 (6/17), Similarity: 0.176 (3/17), Gaps: 0.118 (2/17), Score: 0.100 (9/90), Mass difference: 240.101 Da 120.775 ‰
Start: A 0 B 0, Path: 1i3X2D2=1:2i1=6X
Tolerance: 10 ppm, Alignment: global (==), Maximal isobaric step: 4

IPIWKTLKN·HGHDDWF A
LFSA--LKGGHEIYFER B
─xxx++  ╶╴ xxxxxx

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

Identity: 0.600 (3/5), Mass similarity: 1.000 (5/5), Similarity: 0.600 (3/5), Gaps: 0.000 (0/5), Score: 0.808 (21/26), Equal mass
Start: A 6 B 4, Path: 2=1:2i1=
Tolerance: 10 ppm, Alignment: local (⁐==⁐), Maximal isobaric step: 4

IPIWKTLKN·HGHDDWF A
  LFSALKGGHEIYFER B
        ╶╴

Local alignment, the best parts of the peptidoforms are aligned

Identity: 0.357 (5/14), Mass similarity: 0.429 (6/14), Similarity: 0.357 (5/14), Gaps: 0.071 (1/14), Score: 0.250 (20/80), Mass difference: 315.108 Da 177.239 ‰
Start: A 2 B 0, Path: 1i3X2=1X1D2=3X1=
Tolerance: 10 ppm, Alignment: either global (-==-), Maximal isobaric step: 4

IPIWKTLKNHGHDDWF   A
  LFSALKG-GHEIYFER B
  ─xxx  x+  xxx

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.

Selected: Mus musculus House mouse IGHV1-18-26*01 / IgV1-18-26*01
Identity: 0.847 (83/98), Mass similarity: 0.888 (87/98), Similarity: 0.847 (83/98), Gaps: 0.000 (0/98), Score: 0.848 (440/519), Mass difference: 127.979 Da 11.715 ‰
Start: IGHV1-18-26*01 0 Query 0, Path: 26=1X4=3X6=1X4=1X1=1X6=4r3=1X2=1X8=1X2=1X21=
Tolerance: 10 ppm, Alignment: global left, either global right (==-), Maximal isobaric step: 4

FR1      10       20     CDR1    FR2  40        50
EVQLQQSGPELVKPGASVKISCKASGYTFTDYNMHWVKQSHGKSLEWIGY IGHV1-18-26*01
EVQLQQSGPELVKPGASVKISCKASGFTFTDFSIHWVKQSQGKSLDWVGY Query
                          x    xxx      x    x x
CDR2    FR3       70        80        90      C3
IYPYNGGTGYNQKFKSKATLTVDNSSSTAYMELRSLTSEDSAVYYCAR   IGHV1-18-26*01
IYPYTGGNGYNLKFQSKATLTVDTSSTTAYMELRSLTSEDSAVYYCARRE Query
    ╶──╴   x  x        x  x
GNFVGAMDYWGQGTSVTVSSAKTTAPSVYPLAPVCGDTTGSSVTLGCLVK Query
GYFPEPVTLTWNSGSLSSGVHTFPAVLQSDLYTLSSSVTVTSSTWPSQSI Query
TCNVAHPASSTKVDKKIEPRGPTIKPCPPCKCPAPNLLGGPSVFIFPPKI Query
KDVLMISLSPIVTCVVVDVSEDDPDVQISWFVNNVEVHTAQTQTHREDYN Query
STLRVVSALPIQHQDWMSGKEFKCKVNNKDLPAPIERTISKPKGSVRAPQ Query
VYVLPPPEEEMTKKQVTLTCMVTDFMPEDIYVEWTNNGKTELNYKNTEPV Query
LDSDGSYFMYSKLRVEKKNWVERNSYSCSVVHEGLHNHHTTKSFSRTPGK Query

Aligning a V gene with a full antibody sequence

Selected: Mus musculus House mouse IGHV1-18-26*01 / IgV1-18-26*01
Identity: 0.759 (63/83), Mass similarity: 0.807 (67/83), Similarity: 0.759 (63/83), Gaps: 0.060 (5/83), Score: 0.770 (331/430), Mass difference: 725.292 Da 78.644 ‰
Start: IGHV1-18-26*01 0 Query 0, Path: 5D21=1X4=3X6=1X4=1X1=1X6=4r3=1X2=1X8=1X2=1X6=
Tolerance: 10 ppm, Alignment: global left, either global right (==-), Maximal isobaric step: 4

FR1      10       20     CDR1    FR2  40        50
EVQLQQSGPELVKPGASVKISCKASGYTFTDYNMHWVKQSHGKSLEWIGY IGHV1-18-26*01
-----QSGPELVKPGASVKISCKASGFTFTDFSIHWVKQSQGKSLDWVGY Query
                          x    xxx      x    x x
CDR2    FR3       70        80        90      C3
IYPYNGGTGYNQKFKSKATLTVDNSSSTAYMELRSLTSEDSAVYYCAR   IGHV1-18-26*01
IYPYTGGNGYNLKFQSKATLTVDTSSTTAYMEL                  Query
    ╶──╴   x  x        x  x

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.