changes

src/lib.rs

@@ -1,6 +1,6 @@
-use anyhow::{Result, Ok};
+use anyhow::{Ok, Result};
 use log::info;
-use minimap2::Mapping;
+use minimap2::{Aligner, Mapping};
 use rust_htslib::bam::{self, Record};
 use std::{
     collections::{HashMap, VecDeque},
@@ -20,7 +20,9 @@ pub struct Chromosome {
 #[derive(Debug, Clone)]
 pub struct Contig {
     pub id: String,
+    // contig seq on ref
     pub mappings: Vec<Mapping>,
+    // reads on ref
     pub supporting_records: Vec<Record>,
     pub sequence: String,
     pub contig_ref: ContigRef,
@@ -314,10 +316,28 @@ impl Chromosome {
 
 impl Contig {
     pub fn sort(&mut self) {
-        self.mappings.sort_by(|a, b| a.target_start.cmp(&b.target_start));
+        self.mappings
+            .sort_by(|a, b| a.target_start.cmp(&b.target_start));
     }
 
     pub fn write_bam(&self, path: &str) -> Result<()> {
+        let aligner = Aligner::builder()
+            .asm5()
+            .with_threads(8)
+            .with_cigar()
+            .with_seq(self.sequence.as_bytes())
+            .expect("Unable to build index");
+
+        let mut mappings = Vec::new();
+        for record in self.supporting_records.iter() {
+            let seq = record.seq().as_bytes();
+            let alignment = aligner
+                .map(&seq, false, false, None, None)
+                .expect("Unable to align");
+            mappings.push(alignment);
+        }
+        let mappings: Vec<_> = mappings.into_iter().flatten().collect();
+
         let mut header = bam::Header::new();
         let mut sq_record = rust_htslib::bam::header::HeaderRecord::new(b"SQ");
         sq_record.push_tag(b"SN", self.id.clone());
@@ -327,8 +347,14 @@ impl Contig {
         let mut out = bam::Writer::from_path(path, &header, bam::Format::Bam).unwrap();
 
         // copy reverse reads to new BAM file
-        for mapping in self.mappings.iter() {
-            let record = minimap2::htslib::mapping_to_record(Some(mapping), self.sequence.as_bytes(), header.clone(), None, Some(self.id.as_bytes()));
+        for mapping in mappings.iter() {
+            let record = minimap2::htslib::mapping_to_record(
+                Some(mapping),
+                self.sequence.as_bytes(),
+                header.clone(),
+                None,
+                Some(self.id.as_bytes()),
+            );
             out.write(&record)?;
         }
         rust_htslib::bam::index::build(path, None, rust_htslib::bam::index::Type::Bai, 5)?;