Browse Source

qual for qnames at base

Thomas 1 year ago
parent
commit
896ce97f87
1 changed files with 8 additions and 3 deletions
  1. 8 3
      src/lib.rs

+ 8 - 3
src/lib.rs

@@ -445,6 +445,7 @@ pub fn qnames_at_base(
     chr: &str,
     chr: &str,
     start: i32,
     start: i32,
     with_next_ins: bool,
     with_next_ins: bool,
+    min_qual: u8,
 ) -> Result<Vec<(String, u8)>> {
 ) -> Result<Vec<(String, u8)>> {
     let stop = start + 1;
     let stop = start + 1;
     let mut bases = Vec::new();
     let mut bases = Vec::new();
@@ -455,14 +456,18 @@ pub fn qnames_at_base(
             "Can't pilup bam at position {}:{}-{}",
             "Can't pilup bam at position {}:{}-{}",
             chr, start, stop
             chr, start, stop
         ))?;
         ))?;
+        
         let position = pileup.pos() as i32;
         let position = pileup.pos() as i32;
         if position == start {
         if position == start {
             for alignment in pileup.alignments() {
             for alignment in pileup.alignments() {
+                let record = alignment.record();
+                if record.mapq() < min_qual {
+                    continue;
+                }
                 match alignment.indel() {
                 match alignment.indel() {
                     bam::pileup::Indel::Ins(_len) => bam_pileup.push(b'I'),
                     bam::pileup::Indel::Ins(_len) => bam_pileup.push(b'I'),
                     bam::pileup::Indel::Del(_len) => bam_pileup.push(b'D'),
                     bam::pileup::Indel::Del(_len) => bam_pileup.push(b'D'),
                     _ => {
                     _ => {
-                        let record = alignment.record();
                         let qname = String::from_utf8(record.qname().to_vec())?;
                         let qname = String::from_utf8(record.qname().to_vec())?;
                         if record.seq_len() > 0 {
                         if record.seq_len() > 0 {
                             if let Some(b) = hts_base_at(&record, start as u32, with_next_ins)? {
                             if let Some(b) = hts_base_at(&record, start as u32, with_next_ins)? {
@@ -870,7 +875,7 @@ mod tests {
         let ins_pos = ("chr7", 511_746);
         let ins_pos = ("chr7", 511_746);
         let mut bam = IndexedReader::from_path(diag_bam_path).unwrap();
         let mut bam = IndexedReader::from_path(diag_bam_path).unwrap();
 
 
-        let reads = qnames_at_base(&mut bam, ins_pos.0, ins_pos.1, true).unwrap();
+        let reads = qnames_at_base(&mut bam, ins_pos.0, ins_pos.1, true, 50).unwrap();
 
 
         let n_ins = reads.iter().filter(|(_, b)| *b == b'I').count();
         let n_ins = reads.iter().filter(|(_, b)| *b == b'I').count();
         assert_eq!(n_ins, 5);
         assert_eq!(n_ins, 5);
@@ -885,7 +890,7 @@ mod tests {
         let ins_pos = ("chr7", 511_886);
         let ins_pos = ("chr7", 511_886);
         let mut bam = IndexedReader::from_path(diag_bam_path).unwrap();
         let mut bam = IndexedReader::from_path(diag_bam_path).unwrap();
 
 
-        let reads = qnames_at_base(&mut bam, ins_pos.0, ins_pos.1, false).unwrap();
+        let reads = qnames_at_base(&mut bam, ins_pos.0, ins_pos.1, false, 50).unwrap();
 
 
         let n_del = reads.iter().filter(|(_, b)| *b == b'D').count();
         let n_del = reads.iter().filter(|(_, b)| *b == b'D').count();
         assert_eq!(n_del, 1);
         assert_eq!(n_del, 1);