|
|
@@ -128,7 +128,7 @@ pub fn get_n_start_end_qual(
|
|
|
chr: &str,
|
|
|
start: i32,
|
|
|
stop: i32,
|
|
|
- mapq: u8
|
|
|
+ mapq: u8,
|
|
|
) -> Result<usize> {
|
|
|
bam.fetch((chr, start, stop))?;
|
|
|
|
|
|
@@ -138,7 +138,7 @@ pub fn get_n_start_end_qual(
|
|
|
let record = read.context(format!("eRR"))?;
|
|
|
let rstart = record.pos() as i32;
|
|
|
let rend = record.reference_end() as i32;
|
|
|
- if rstart >= start && rstart < stop || rend >= start && rend < stop {
|
|
|
+ if rstart >= start && rstart < stop || rend >= start && rend < stop {
|
|
|
if record.mapq() >= mapq {
|
|
|
n_start += 1;
|
|
|
}
|
|
|
@@ -146,9 +146,43 @@ pub fn get_n_start_end_qual(
|
|
|
}
|
|
|
|
|
|
Ok(n_start)
|
|
|
- // Ok(start_positions.len())
|
|
|
}
|
|
|
|
|
|
+pub fn get_start_end_qual(
|
|
|
+ bam: &mut rust_htslib::bam::IndexedReader,
|
|
|
+ chr: &str,
|
|
|
+ start: i32,
|
|
|
+ stop: i32,
|
|
|
+ mapq: u8,
|
|
|
+) -> Result<Vec<i32>> {
|
|
|
+ bam.fetch((chr, start, stop))?;
|
|
|
+ let length = stop - start;
|
|
|
+ let mut results = Vec::new();
|
|
|
+ results.resize(length as usize, 0);
|
|
|
+
|
|
|
+ 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() += 1;
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ if rend >= start && rend < stop {
|
|
|
+ if record.mapq() >= mapq {
|
|
|
+ let index = rend - start;
|
|
|
+ *results.get_mut(index as usize).unwrap() += 1;
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ }
|
|
|
+
|
|
|
+ Ok(results)
|
|
|
+}
|
|
|
|
|
|
pub fn range_depths(
|
|
|
bam: &mut rust_htslib::bam::IndexedReader,
|