Repository: bacpop/ska.rust (rust)

Issue:
Reference bias in map with repeats
Proposed fix here: https://github.com/bacpop/ska.rust/compare/master...nickjcroucher:ska.rust:master

Reference copy can overwrite ambiguous bases in some cases:
> There is a reference sequence
> Query A contains the sequences: `TTTTATCTTGATTTTCT`
> Query B contains the sequences:
> ```
> TTTTATCTTGATTTTCT
> TTTATCTTAATTTTCTT
> ```
> Hence Query A contains the split k-mer: `AAGAAAAT AAGATAAA C`
> And Query B contains the split k-mer: `AAGAAAAT AAGATAAA Y`
> However, when mapped, a SNP is called distinguishing the sequences:
> ```
> >Query B
> TTTTAGTTTTATCTTAATTTTCTTA
> >Query A
> -------TTTATCTTGATTTTCTT
> ```
> The expected Y is converted to A because A is in the reference and matches the split k-mer arm.

Pull request description:
Fix ambiguous bases in pseudoalignment
Write them in (in some cases rewrite them in) before writing seq out

Closes #28 

Patch:
diff --git a/src/ska_dict.rs b/src/ska_dict.rs
index 2ecf737..fd8f6a4 100644
--- a/src/ska_dict.rs
+++ b/src/ska_dict.rs
@@ -257,7 +257,6 @@ where
     }
 }
 
-
 #[cfg(test)]
 mod tests {
     use super::*;
@@ -301,4 +300,3 @@ mod tests {
         test_obj_panic.add_palindrome_to_dict(555, 1);
     }
 }
-
diff --git a/src/ska_ref/aln_writer.rs b/src/ska_ref/aln_writer.rs
index 8e32fa0..853425a 100644
--- a/src/ska_ref/aln_writer.rs
+++ b/src/ska_ref/aln_writer.rs
@@ -11,6 +11,8 @@
 //! All sequence is stored in memory. Calling [`AlnWriter::get_seq()`] finalises
 //! and writes to the end of the last contig.
 
+use crate::ska_dict::bit_encoding::is_ambiguous;
+
 /// Stores indexes which keep track of the writing, the output sequence and a
 /// reference to the reference sequences which are written in to flanks of
 /// each match.
@@ -35,6 +37,8 @@ pub struct AlnWriter<'a> {
     /// Whether the finalise function has been run, filling to the end of the
     /// final contig.
     finalised: bool,
+    /// Buffer for amibguous bases, which are written in at the end.
+    _ambig_out: Vec<(u8, usize)>,
 }
 
 impl<'a> AlnWriter<'a> {
@@ -53,6 +57,7 @@ impl<'a> AlnWriter<'a> {
             seq_out: vec![b'-'; total_size],
             half_split_len,
             finalised: false,
+            _ambig_out: Vec::new(),
         }
     }
 
@@ -110,6 +115,11 @@ impl<'a> AlnWriter<'a> {
             self.fill_contig();
         }
         if mapped_pos < self.next_pos {
+            // Ambiguous bases may clash with the flanks, which are copied from
+            // reference. Deal with these in `finalise`.
+            if is_ambiguous(base) {
+                self._ambig_out.push((base, mapped_pos + self.chrom_offset));
+            }
             self.last_mapped = mapped_pos;
         } else {
             // Write bases between last match and this one
@@ -139,6 +149,10 @@ impl<'a> AlnWriter<'a> {
             while self.curr_chrom < self.ref_seq.len() {
                 self.fill_contig();
             }
+            // Make sure any ambiguous bases are correct
+            for (ambig_base, ambig_pos) in &self._ambig_out {
+                self.seq_out[*ambig_pos] = *ambig_base;
+            }
             self.finalised = true;
         }
     }


Test patch:
diff --git a/tests/map.rs b/tests/map.rs
index b87cd38..8cb21fc 100644
--- a/tests/map.rs
+++ b/tests/map.rs
@@ -58,6 +58,27 @@ fn map_aln() {
         .arg(sandbox.file_string("indel_test.fa", TestDir::Input))
         .assert()
         .stdout_eq_path(sandbox.file_string("map_aln_indels.stdout", TestDir::Correct));
+
+    Command::new(cargo_bin("ska"))
+        .current_dir(sandbox.get_wd())
+        .arg("build")
+        .arg("-k")
+        .arg("17")
+        .arg(sandbox.file_string("ambig_test_1.fa", TestDir::Input))
+        .arg(sandbox.file_string("ambig_test_2.fa", TestDir::Input))
+        .arg("-o")
+        .arg("ambig_map")
+        .arg("--single-strand")
+        .assert()
+        .success();
+
+    Command::new(cargo_bin("ska"))
+        .current_dir(sandbox.get_wd())
+        .arg("map")
+        .arg(sandbox.file_string("ambig_test_ref.fa", TestDir::Input))
+        .arg("ambig_map.skf")
+        .assert()
+        .stdout_eq_path(sandbox.file_string("map_aln_ambig.stdout", TestDir::Correct));
 }
 
 #[test]
diff --git a/tests/test_files_in/ambig_test_1.fa b/tests/test_files_in/ambig_test_1.fa
new file mode 100644
index 0000000..4deb034
--- /dev/null
+++ b/tests/test_files_in/ambig_test_1.fa
@@ -0,0 +1,2 @@
+>query_A
+TTTTATCTTGATTTTCTT
diff --git a/tests/test_files_in/ambig_test_2.fa b/tests/test_files_in/ambig_test_2.fa
new file mode 100644
index 0000000..953e404
--- /dev/null
+++ b/tests/test_files_in/ambig_test_2.fa
@@ -0,0 +1,4 @@
+>query_B_1
+TTTTATCTTGATTTTCTT
+>query_B_2
+TTTTATCTTAATTTTCTT
diff --git a/tests/test_files_in/ambig_test_ref.fa b/tests/test_files_in/ambig_test_ref.fa
new file mode 100644
index 0000000..ac65dbb
--- /dev/null
+++ b/tests/test_files_in/ambig_test_ref.fa
@@ -0,0 +1,2 @@
+>ref
+TTTTAGTTTTATCTTAATTTTCTTA
diff --git a/tests/test_results_correct/map_aln_ambig.stdout b/tests/test_results_correct/map_aln_ambig.stdout
new file mode 100644
index 0000000..eba6c9c
--- /dev/null
+++ b/tests/test_results_correct/map_aln_ambig.stdout
@@ -0,0 +1,4 @@
+>ambig_test_1
+-------TTTATCTTGATTTTCTT-
+>ambig_test_2
+------TTTTATCTTRATTTTCTT-
diff --git a/tests/test_results_correct/map_aln_k9.stdout b/tests/test_results_correct/map_aln_k9.stdout
index f90ac74..8888733 100644
--- a/tests/test_results_correct/map_aln_k9.stdout
+++ b/tests/test_results_correct/map_aln_k9.stdout
@@ -1,4 +1,4 @@
 >test_1
-------------GATTTGAAGGCGTTTATGATATTTAGCTGTACGAGAGTCTTTTAAGAGTSTTTTGATGGTTT------------
+------------GATTTGAAGGCGTTTATGATATTTAGCTGTACGAGAGTSTTTTAAGAGTSTTTTGATGGTTT------------
 >test_2
 ------------GATTTGAAGGCGTTTATGATATATAGCTGTACGAGAGTCTTTTAAAAGTGTTTTGATGGTTT------------
diff --git a/tests/test_results_correct/map_ss.stdout b/tests/test_results_correct/map_ss.stdout
index f512e77..83a9fd0 100644
--- a/tests/test_results_correct/map_ss.stdout
+++ b/tests/test_results_correct/map_ss.stdout
@@ -1,4 +1,4 @@
 >test_1
-------------GATTTGAAGGCGTTTATGATATTTAGCTGTACGAGAGTCTTTTAAGAGTSTTTTGATGGTTT------------
+------------GATTTGAAGGCGTTTATGATATTTAGCTGTACGAGAGTSTTTTAAGAGTSTTTTGATGGTTT------------
 >test_2_rc
 --------------------------------------------------------------------------------------
diff --git a/tests/test_results_correct/map_vcf_ss.stdout b/tests/test_results_correct/map_vcf_ss.stdout
index 777ceb8..a1d39ef 100644
--- a/tests/test_results_correct/map_vcf_ss.stdout
+++ b/tests/test_results_correct/map_vcf_ss.stdout
@@ -51,7 +51,7 @@ fake_ref	47	.	G	.	.	.	.	GT	0	.
 fake_ref	48	.	A	.	.	.	.	GT	0	.
 fake_ref	49	.	G	.	.	.	.	GT	0	.
 fake_ref	50	.	T	.	.	.	.	GT	0	.
-fake_ref	51	.	C	.	.	.	.	GT	0	.
+fake_ref	51	.	C	N	.	.	.	GT	1	.
 fake_ref	52	.	T	.	.	.	.	GT	0	.
 fake_ref	53	.	T	.	.	.	.	GT	0	.
 fake_ref	54	.	T	.	.	.	.	GT	0	.