Skip to content

Commit

Permalink
Added maximal score to alignment, fixed local alignment bug, preferre…
Browse files Browse the repository at this point in the history
…d smaller Isobaric/Rotated sets
  • Loading branch information
douweschulte committed Dec 19, 2023
1 parent b4ed90c commit fbaefae
Show file tree
Hide file tree
Showing 3 changed files with 22 additions and 11 deletions.
6 changes: 4 additions & 2 deletions src/align/alignment.rs
Original file line number Diff line number Diff line change
Expand Up @@ -15,9 +15,11 @@ use crate::{LinearPeptide, MolecularFormula};
pub struct Alignment {
/// The absolute score of this alignment
pub absolute_score: isize,
/// The maximal score of this alignment: the average score of the sequence slices on sequence a and b if they were aligned to themself, rounded down.
/// Think of it like this: `align(sequence_a.sequence[start_a..len_a], sequence_a.sequence[start_a..len_a])`.
pub maximal_score: isize,
/// The normalised score, normalised for the alignment length and for the used alphabet.
/// The normalisation is as follows `score / max_score` where max score is the average score of the sequence slices on sequence a and b if they were aligned to themself.
/// Think of it like this: `align(sequence_a.sequence[start_a..len_a], sequence_a.sequence[start_a..len_a])`
/// The normalisation is as follows `absolute_score / max_score`.
pub normalised_score: f64,
/// The path or steps taken for the alignment
pub path: Vec<Piece>,
Expand Down
22 changes: 15 additions & 7 deletions src/align/mass_alignment.rs
Original file line number Diff line number Diff line change
Expand Up @@ -140,7 +140,7 @@ pub fn align<const STEPS: usize>(
let high_score = high.0;
while ty == Type::Global || !(high.1 == 0 && high.2 == 0) {
let value = matrix[high.1][high.2].clone();
if value.step_a == 0 && value.step_b == 0 {
if value.step_a == 0 && value.step_b == 0 || ty == Type::Local && value.score < 0 {
break;
}
high = (
Expand All @@ -152,18 +152,20 @@ pub fn align<const STEPS: usize>(
}

let path: Vec<Piece> = path.into_iter().rev().collect();
let max_score = seq_a.sequence
let max_score = (seq_a.sequence
[high.1..high.1 + path.iter().map(|p| p.step_a as usize).sum::<usize>()]
.iter()
.map(|a| alphabet[a.aminoacid as usize][a.aminoacid as usize] as isize)
.sum::<isize>()
+ seq_b.sequence[high.2..high.2 + path.iter().map(|p| p.step_b as usize).sum::<usize>()]
.iter()
.map(|a| alphabet[a.aminoacid as usize][a.aminoacid as usize] as isize)
.sum::<isize>();
.sum::<isize>())
/ 2;
Alignment {
absolute_score: high_score,
normalised_score: high_score as f64 / max_score as f64 * 2.0,
normalised_score: high_score as f64 / max_score as f64,
maximal_score: max_score,
path,
start_a: high.1,
start_b: high.2,
Expand Down Expand Up @@ -199,7 +201,13 @@ fn score_pair(
1,
)
}
(false, true) => Piece::new(score + ISOMASS as isize, ISOMASS, MatchType::Isobaric, 1, 1), // TODO: I/L/J is now also scored as isobaric, which is correct but the score is considerably lower then in previous iterations
(false, true) => Piece::new(
score + ISOBARIC as isize,
ISOBARIC,
MatchType::Isobaric,
1,
1,
), // TODO: I/L/J is now also scored as isobaric, which is correct but the score is considerably lower then in previous iterations
(false, false) => Piece::new(
score + MISMATCH as isize,
MISMATCH,
Expand Down Expand Up @@ -231,9 +239,9 @@ fn score(
});
#[allow(clippy::cast_possible_wrap)]
let local = if rotated {
SWITCHED * a.len() as i8
BASE_SPECIAL + ROTATED * a.len() as i8
} else {
ISOMASS * (a.len() + b.len()) as i8 / 2
BASE_SPECIAL + ISOBARIC * (a.len() + b.len()) as i8 / 2
};
Some(Piece::new(
score + local as isize,
Expand Down
5 changes: 3 additions & 2 deletions src/align/scoring.rs
Original file line number Diff line number Diff line change
Expand Up @@ -22,8 +22,9 @@ pub enum MatchType {

pub const MISMATCH: i8 = -1;
pub const MASS_MISMATCH_PENALTY: i8 = -1;
pub const SWITCHED: i8 = 3;
pub const ISOMASS: i8 = 2;
pub const BASE_SPECIAL: i8 = 1;
pub const ROTATED: i8 = 3;
pub const ISOBARIC: i8 = 2;
pub const GAP_START_PENALTY: i8 = -5;
pub const GAP_EXTEND_PENALTY: i8 = -1;

Expand Down

0 comments on commit fbaefae

Please sign in to comment.