Crate mzalign

Crate mzalign 

Source
Expand description

§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-alignment1.

§Library features

§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

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);
// 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

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

  1. Schulte, D.; Snijder, J. A Handle on Mass Coincidence Errors in De Novo Sequencing of Antibodies by Bottom-up Proteomics. https://doi.org/10.1021/acs.jproteome.4c00188. ↩

  2. Ehrenmann, F.; Kaas, Q.; Lefranc, M.-P. IMGT/3Dstructure-DB and IMGT/DomainGapAlign: A Database and a Tool for Immunoglobulins or Antibodies, T Cell Receptors, MHC, IgSF and MhcSF. https://doi.org/10.1093/nar/gkp946. ↩

Modules§

matrix
Different scoring matrices that can be used. Matrices from: https://www.ncbi.nlm.nih.gov/IEB/ToolBox/CPP_DOC/lxr/source/src/util/tables/ and https://www.ncbi.nlm.nih.gov/IEB/ToolBox/C_DOC/lxr/source/data/
prelude
A subset of the types and traits that are envisioned to be used the most, importing this is a good starting point for working with the crate

Structs§

AlignIndex
An alignment index to store a collection of sequences to do alignment against with one or more other sequences. This improves performance because it can cache some intermediate results.
AlignScoring
The scoring parameters for the mass alignment.
AlignType
The type of alignment to perform
Alignment
An alignment of two reads. It has either a reference to the two sequences to prevent overzealous use of memory, or if needed use Self::to_owned to get a variant that clones the sequences and so can be used in more places.
ConsecutiveAlignment
A consecutive alignment, which align one sequence to multiple sequences.
Piece
A piece in an alignment, determining what step was taken in the alignment and how this impacted the score
Score
The score of an alignment
Stats
Statistics for an alignment with some helper functions to easily retrieve the number of interest.

Enums§

MatchType
The type of a single match step
PairMode
The pair mode of an alignment, which side is the peptidoform and which the database.
Side
The alignment specification for a single side

Functions§

align
Create an alignment of two peptides based on mass and homology. The substitution matrix is in the exact same order as the definition of AminoAcid. The AlignScoring sets the rules and exact scores while scoring. The AlignType controls the alignment behaviour, global/local or anything in between.
consecutive_align
Only available if features align and imgt are turned on. Align one sequence to multiple consecutive genes. Each gene can be controlled to be global to the left or free to allow unmatched residues between it and the previous gene. If the sequence is too short to cover all genes only the genes that could be matched are returned.
par_consecutive_align
Only available with if features align, rayon, and imgt are turned on. Align one sequence to multiple consecutive genes. Each gene can be controlled to be global to the left or free to allow unmatched residues between it and the previous gene. If the sequence is too short to cover all genes only the genes that could be matched are returned.