m

positions.tsv

@@ -1,2 +1,2 @@
-chr1	47239294	47239298
-chr1	47239294	47239298
+chr1	47239294	47239298	EH38E1702181
+chr1	47239294	47239298	EH38E1702181

src/main.rs

@@ -30,7 +30,8 @@ fn main() {
                     &mut bam, 
                     split.next().unwrap(), 
                     split.next().unwrap().parse().unwrap(),
-                    split.next().unwrap().parse().unwrap()
+                    split.next().unwrap().parse().unwrap(),
+                    split.next().unwrap().to_string()
                 )
             );
         }
@@ -49,10 +50,11 @@ struct Result {
     seqname: String,
     ref_pos: i64,
     inserted_seq: String,
-    depth: u32
+    depth: u32,
+    id: String
 }
 
-fn fetch_ins(bam: &mut IndexedReader, seqname: &str, start: u32, end: u32) -> Vec<Result> {
+fn fetch_ins(bam: &mut IndexedReader, seqname: &str, start: u32, end: u32, id: String) -> Vec<Result> {
     bam.fetch((seqname, start, end)).unwrap();
 
     // pileup over all covered sites
@@ -77,8 +79,14 @@ fn fetch_ins(bam: &mut IndexedReader, seqname: &str, start: u32, end: u32) -> Ve
                                 query_pos += l;
                             },
                             bam::record::Cigar::Ins(l) => {
+                                // let l = if query_pos + l <= seq.len() as u32 { *l } else { query_pos - seq.len() as u32 };
                                 let r = (query_pos as usize)..((query_pos + l) as usize);
-                                inserted_seq = String::from_utf8_lossy(&seq[r]).to_string();
+                                inserted_seq = if let Some(s) = seq.get(r.clone()) {
+                                     String::from_utf8_lossy(s).to_string()
+                                } else {
+                                    println!("{:?} {:?}", seq, r);
+                                    "NN".to_string()
+                                };
                                 break;
                             },
                             bam::record::Cigar::Del(l) => {
@@ -110,7 +118,8 @@ fn fetch_ins(bam: &mut IndexedReader, seqname: &str, start: u32, end: u32) -> Ve
                         seqname: seqname.to_string(),
                         ref_pos, 
                         inserted_seq,
-                        depth: pileup.depth()
+                        depth: pileup.depth(),
+                        id: id.clone()
                     });
                 },
                 _ => ()
@@ -127,6 +136,6 @@ impl Result {
     }
 
     fn tsv(&self) -> String {
-        format!("{}\t{}\t{}", self.read_name, self.hgvs(), self.depth)
+        format!("{}\t{}\t{}\t{}", self.id, self.read_name, self.hgvs(), self.depth)
     }
 }