|
|
@@ -1,7 +1,7 @@
|
|
|
use anyhow::{Context, Ok, Result};
|
|
|
use rust_htslib::{
|
|
|
bam,
|
|
|
- bam::{ext::BamRecordExtensions, record, Read},
|
|
|
+ bam::{ext::BamRecordExtensions, record, Read, Record},
|
|
|
};
|
|
|
|
|
|
use average::Mean;
|
|
|
@@ -178,12 +178,51 @@ pub fn get_start_end_qual(
|
|
|
*results.get_mut(index as usize).unwrap() += 1;
|
|
|
}
|
|
|
}
|
|
|
-
|
|
|
}
|
|
|
|
|
|
Ok(results)
|
|
|
}
|
|
|
|
|
|
+pub fn get_start_end_qual_rec(
|
|
|
+ bam: &mut rust_htslib::bam::IndexedReader,
|
|
|
+ chr: &str,
|
|
|
+ start: i32,
|
|
|
+ stop: i32,
|
|
|
+ mapq: u8,
|
|
|
+) -> Result<Vec<Vec<Record>>> {
|
|
|
+ bam.fetch((chr, start, stop))?;
|
|
|
+ let length = stop - start;
|
|
|
+ let mut results: Vec<Vec<Record>> = Vec::new();
|
|
|
+ results.resize(length as usize, vec![]);
|
|
|
+
|
|
|
+ for read in bam.records() {
|
|
|
+ let record = read.context(format!("Error while parsing record"))?;
|
|
|
+ let rstart = record.pos() as i32;
|
|
|
+ let rend = record.reference_end() as i32;
|
|
|
+
|
|
|
+ if rstart >= start && rstart < stop {
|
|
|
+ if record.mapq() >= mapq {
|
|
|
+ let index = rstart - start;
|
|
|
+ results
|
|
|
+ .get_mut(index as usize)
|
|
|
+ .unwrap()
|
|
|
+ .push(record.clone());
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ if rend >= start && rend < stop {
|
|
|
+ if record.mapq() >= mapq {
|
|
|
+ let index = rend - start;
|
|
|
+ results
|
|
|
+ .get_mut(index as usize)
|
|
|
+ .unwrap()
|
|
|
+ .push(record.clone());
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+ Ok(results)
|
|
|
+}
|
|
|
+
|
|
|
pub fn range_depths(
|
|
|
bam: &mut rust_htslib::bam::IndexedReader,
|
|
|
chr: &str,
|