/// Finds aminoacid substitutions and deletions in query peptides relative to reference peptides, in one gene
///
/// ## Precondition
/// Nucleotide sequences and peptides are required to be stripped from insertions
///
///
/// ## Implementation details
/// We compare reference and query peptides (extracted by the preceding call to Nextalign),
/// one aminoacid at at time, and deduce changes. We then report the change and relevant nucleotide context surrounding
/// this change.
/// Previously we reported one-to-one mapping of aminoacid changes to corresponding nucleotide changes. However, it
/// was not always accurate, because if there are multiple nucleotide changes in a codon, the direct correspondence
/// might not always be established without knowing the order in which nucleotide changes have occurred. And in the
/// context of Nextclade we don't have this information.
Finds aminoacid substitutions and deletions in query peptides relative to reference peptides, in one gene
Finds aminoacid substitutions and deletions in query peptides relative to reference peptides in one gene
fn find_aa_changes_for_gene(
  qry_seq: &[Nuc],
  ref_seq: &[Nuc],
  ref_peptide: &[Aa],
  qry_peptide: &[Aa],
  gene: &Gene,
  alignment_range: &Range,
) -> FindAaChangesOutput {
  assert_eq!(ref_peptide.len(), qry_peptide.len());
  assert_eq!(qry_seq.len(), ref_seq.len());
  let mut aa_substitutions = Vec::<AaSub>::new();
  let mut aa_deletions = Vec::<AaDel>::new();
  let num_nucs = qry_seq.len();
  let num_codons = qry_peptide.len();
  for codon in 0..num_codons {
    let ref_aa = ref_peptide[codon];
    let qry_aa = qry_peptide[codon];
    // Find where the codon is in nucleotide sequences
    let codon_begin = if gene.strand != GeneStrand::Reverse {
      gene.start + codon * 3
    } else {
      gene.end - (codon + 1) * 3
    };
    let codon_end = codon_begin + 3;
    // If the codon is outside of nucleotide alignment, there is nothing to do
    if !alignment_range.contains(codon_begin) || !alignment_range.contains(codon_end) {
      continue;
    }
    // Provide surrounding context in nucleotide sequences: 1 codon to the left and 1 codon to the right
    let context_begin = clamp(codon_begin - 3, 0, num_nucs);
    let context_end = clamp(codon_end + 3, 0, num_nucs);
    let mut ref_context = (&ref_seq[context_begin..context_end]).to_vec();
    let mut query_context = (&qry_seq[context_begin..context_end]).to_vec();
    if gene.strand == GeneStrand::Reverse {
      reverse_complement_in_place(&mut ref_context);
      reverse_complement_in_place(&mut query_context);
    }
    let codon_nuc_range = Range {
      begin: codon_begin,
      end: codon_end,
    };
    let context_nuc_range = Range {
      begin: context_begin,
      end: context_end,
    };
    if qry_aa.is_gap() {
      // Gap in the ref sequence means that this is a deletion in the query sequence
      aa_deletions.push(AaDel {
        gene: gene.gene_name.clone(),
        reff: ref_aa,
        pos: codon,
        codon_nuc_range,
        ref_context: from_nuc_seq(&ref_context),
        query_context: from_nuc_seq(&query_context),
        context_nuc_range,
      });
    }
    // TODO: we might account for ambiguous amino acids in this condition
    else if qry_aa != ref_aa && qry_aa != Aa::X {
      // If not a gap and the state has changed, than it's a substitution
      aa_substitutions.push(AaSub {
        gene: gene.gene_name.clone(),
        reff: ref_aa,
        pos: codon,
        qry: qry_aa,
        codon_nuc_range,
        ref_context: from_nuc_seq(&ref_context),
        query_context: from_nuc_seq(&query_context),
        context_nuc_range,
      });
    }
  }
  FindAaChangesOutput {
    aa_substitutions,
    aa_deletions,
  }
}
fn find_aa_changes_for_gene ( qry_seq : & [ Nuc ] , ref_seq : & [ Nuc ] , ref_peptide : & [ Aa ] , qry_peptide : & [ Aa ] , gene : & Gene , alignment_range : & Range , ) -> FindAaChangesOutput { assert_eq ! ( ref_peptide len ( ) qry_peptide len ( ) ) ; assert_eq ! ( qry_seq len ( ) ref_seq len ( ) ) ; let mut aa_substitutions = Vec :: < AaSub > :: new ( ) ; let mut aa_deletions = Vec :: < AaDel > :: new ( ) ; let num_nucs = qry_seq . len ( ) ; let num_codons = qry_peptide . len ( ) ; for codon in 0 .. num_codons { let ref_aa = ref_peptide [ codon ] ; let qry_aa = qry_peptide [ codon ] ; // Find where the codon is in nucleotide sequences let codon_begin = if gene . strand != GeneStrand :: Reverse { gene . start + codon * 3 } else { gene . end - ( codon + 1 ) * 3 } ; let codon_end = codon_begin + 3 ; // If the codon is outside of nucleotide alignment, there is nothing to do if ! alignment_range . contains ( codon_begin ) || ! alignment_range . contains ( codon_end ) { continue ; } // Provide surrounding context in nucleotide sequences: 1 codon to the left and 1 codon to the right let context_begin = clamp ( codon_begin - 3 , 0 , num_nucs ) ; let context_end = clamp ( codon_end + 3 , 0 , num_nucs ) ; let mut ref_context = ( & ref_seq [ context_begin .. context_end ] ) . to_vec ( ) ; let mut query_context = ( & qry_seq [ context_begin .. context_end ] ) . to_vec ( ) ; if gene . strand == GeneStrand :: Reverse { reverse_complement_in_place ( & mut ref_context ) ; reverse_complement_in_place ( & mut query_context ) ; } let codon_nuc_range = Range { begin : codon_begin , end : codon_end , } ; let context_nuc_range = Range { begin : context_begin , end : context_end , } ; if qry_aa . is_gap ( ) { // Gap in the ref sequence means that this is a deletion in the query sequence aa_deletions . push ( AaDel { gene : gene . gene_name . clone ( ) , reff : ref_aa , pos : codon , codon_nuc_range , ref_context : from_nuc_seq ( & ref_context ) , query_context : from_nuc_seq ( & query_context ) , context_nuc_range , } ) ; } // TODO: we might account for ambiguous amino acids in this condition else if qry_aa != ref_aa && qry_aa != Aa :: X { // If not a gap and the state has changed, than it's a substitution aa_substitutions . push ( AaSub { gene : gene . gene_name . clone ( ) , reff : ref_aa , pos : codon , qry : qry_aa , codon_nuc_range , ref_context : from_nuc_seq ( & ref_context ) , query_context : from_nuc_seq ( & query_context ) , context_nuc_range , } ) ; } } FindAaChangesOutput { aa_substitutions , aa_deletions , } }
Finds aminoacid substitutions and deletions in query peptides relative to reference peptides, in one gene
Finds aminoacid substitutions and deletions in query peptides relative to reference peptides in one gene
qry_seq
&[Nuc]
ref_seq
&[Nuc]
ref_peptide
&[Aa]
qry_peptide
&[Aa]
gene
&Gene
alignment_range
&Range
precondition
Nucleotide sequences and peptides are required to be stripped from insertions
Nucleotide sequences and peptides are required to be stripped from insertions
implementation details
We compare reference and query peptides (extracted by the preceding call to Nextalign),
one aminoacid at at time, and deduce changes. We then report the change and relevant nucleotide context surrounding
this change.
Previously we reported one-to-one mapping of aminoacid changes to corresponding nucleotide changes. However, it
was not always accurate, because if there are multiple nucleotide changes in a codon, the direct correspondence
might not always be established without knowing the order in which nucleotide changes have occurred. And in the
context of Nextclade we don't have this information.
We compare reference and query peptides ( extracted by the preceding call to Nextalign ) one aminoacid at at time and deduce changes . We then report the change and relevant nucleotide context surrounding this change . Previously we reported one - to - one mapping of aminoacid changes to corresponding nucleotide changes . However it was not always accurate because if there are multiple nucleotide changes in a codon the direct correspondence might not always be established without knowing the order in which nucleotide changes have occurred . And in the context of Nextclade we don ' t have this information .