potential bug: band boundary hit when trailing in-dels are present at the end of a simple stripe alignment

[mmolari]

Hello! Reporting a small potential bug that I ran into while using Nextalign on another project.

In my case the code was slightly different, but I believe that the issue could be reproduced by performing a banded alignment of these two sequences

    let ref_str = "------------------------------TTGGCCCCGGTGCTGTCCGTCAACACGTCGTCGTCCGGCGACCTACCTGGTCTCAAAGGAGGTTTTGTTAAATGAATTAGATGGGTAAGGTTACCACGTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA";
    let qry_str = "GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGTTGGCCCCGGTGCTGTCCGTCAACACGTCGTCGTCCGGCGACCTACCTGGTCTCAAAGGAGGTTTTGTTAAATGAATTAGATGGGTAAGGTTACCACGTCA------------------------------";

with band shift = -30 and band width = 1, with a simple stripe alignment (e.g. using align_aa, interpreting them as amino-acid alignments) results in a hit_boundary warning.

let qry_str = qry_str.replace('-', "");
let ref_str = ref_str.replace('-', "");

let qry_seq = to_aa_seq(&qry_str).unwrap();
let ref_seq = to_aa_seq(&ref_str).unwrap();

let gap_open_close = get_gap_open_close_scores_flat(&ref_seq, &params);
let band_width = 1;
let shift = -30;
let alignment = align_aa(&qry_seq, &ref_seq, &gap_open_close, band_width, shift, &params).unwrap();
assert!(!alignment.hit_boundary);

It might be desirable to avoid this, and not consider trailing indels as a boundary hit. If so, in my case this was solved by applying the following changes to the score_matrix function.

[mmolari]

It's not exactly how it was meant to be used, but I could for example reproduce the issue by adding the following test to align.rs:

  #[rstest]
  fn banded_alignment_trailing_indels_should_not_hit_boundary(ctx: Context) -> Result<(), Report> {
    use crate::align::gap_open::get_gap_open_close_scores_flat;
    use crate::alphabet::aa::to_aa_seq;

    // Test case from the issue: sequences with trailing indels that should not cause boundary hits
    #[rustfmt::skip]
    let ref_str = "------------------------------TTGGCCCCGGTGCTGTCCGTCAACACGTCGTCGTCCGGCGACCTACCTGGTCTCAAAGGAGGTTTTGTTAAATGAATTAGATGGGTAAGGTTACCACGTCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA";
    #[rustfmt::skip]
    let qry_str = "GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGTTGGCCCCGGTGCTGTCCGTCAACACGTCGTCGTCCGGCGACCTACCTGGTCTCAAAGGAGGTTTTGTTAAATGAATTAGATGGGTAAGGTTACCACGTCA------------------------------";

    // Remove gaps to get the actual sequences
    let ref_str = ref_str.replace('-', "");
    let qry_str = qry_str.replace('-', "");

    let qry_seq = to_aa_seq(&qry_str)?;
    let ref_seq = to_aa_seq(&ref_str)?;

    let nucl_ref_seq = to_nuc_seq(&ref_str)?;

    let gap_open_close = get_gap_open_close_scores_flat(&nucl_ref_seq, &ctx.params);
    let band_width = 1;
    let mean_shift = -30;

    let alignment = align_aa(&qry_seq, &ref_seq, &gap_open_close, &ctx.params, band_width, mean_shift);

    // The alignment should not hit the boundary despite having trailing indels
    assert!(
      !alignment.hit_boundary,
      "Alignment with trailing indels should not hit boundary"
    );
    Ok(())
  }

This fails with:

---- align::align::tests::banded_alignment_trailing_indels_should_not_hit_boundary stdout ----
The application panicked (crashed).
Message:  Alignment with trailing indels should not hit boundary
Location: packages/nextclade/src/align/align.rs:644

I think this could maybe also be reproduced directly on nucleotide alignments, by appropriate tuning of band and shift parameters.