extended_htslib/bam/

pileup.rs

// Copyright 2014 Johannes Köster.
// Licensed under the MIT license (http://opensource.org/licenses/MIT)
// This file may not be copied, modified, or distributed
// except according to those terms.

use bio_types::sequence::SequenceRead;
#[cfg(feature = "pileuprayon")]
use rayon::prelude::*;
use std::cmp::Ordering;
use std::collections::HashMap;
use std::collections::hash_map::Entry::Occupied;
use std::collections::hash_map::Entry::Vacant;
use std::fmt;
use std::fmt::Display;
use std::iter;
use std::path::Path;
use std::slice;
use std::thread::available_parallelism;

use crate::bam::FetchDefinition;
use crate::bam::FetchDefinition::RegionString;
use crate::bam::Read;
use crate::bam::ext::BamRecordExtensions;
use crate::bam::ext::IterAlignedPairsFullCigar;
use crate::bam::record::Cigar;
use crate::faidx;
use crate::htslib;

use crate::bam;
use crate::bam::record;
use crate::errors::{Error, Result};

#[derive(Debug, Clone, Copy, PartialEq, Eq)]
#[non_exhaustive]
#[allow(non_snake_case)]
/// Config for RustPileUp. Number of arguments might change over time and therefore is not exhaustive, please use functions.
pub struct RustPileupConfig {
    pub onlyprimary: bool,
    recalculateBAM: bool,
    pub minqual: u8,
    pub minscore: u8,
    pub reversedel: bool,
    pub outputends: bool,
}
impl Default for RustPileupConfig {
    fn default() -> Self {
        Self::new(false, 0, 0, true, true)
    }
}
fn getreader<T>(
    reader: T,
    assemblyreader: Option<T>,
    index: Option<T>,
    assemblyreaderindex: Option<T>,
) -> bam::Result<(bam::IndexedReader, Option<faidx::Reader>)>
where
    T: AsRef<Path>,
{
    let (reader, assemblyreader) = match (reader, assemblyreader, index, assemblyreaderindex) {
        (a, Some(b), None, None) => (bam::IndexedReader::from_path_and_index(a, b)?, None),
        (a, None, None, None) => (bam::IndexedReader::from_path(a)?, None),
        (a, None, Some(c), None) => (
            bam::IndexedReader::from_path(a)?,
            Some(faidx::Reader::from_path(c)?),
        ),
        (a, None, Some(c), Some(d)) => (
            bam::IndexedReader::from_path(a)?,
            Some(faidx::Reader::from_path_and_index(c, d)?),
        ),
        (a, Some(b), Some(c), None) => (
            bam::IndexedReader::from_path_and_index(a, b)?,
            Some(faidx::Reader::from_path(c)?),
        ),
        (a, Some(b), Some(c), Some(d)) => (
            bam::IndexedReader::from_path_and_index(a, b)?,
            Some(faidx::Reader::from_path_and_index(c, d)?),
        ),
        (.., None, Some(a)) => {
            return Err(bam::Error::BamInvalidIndex {
                target: a.as_ref().display().to_string(),
            });
        }
    };
    Ok((reader, assemblyreader))
}
impl RustPileupConfig {
    /// Configuration for pileup. Recalculate BAM must be off (else panics)
    pub fn new(
        onlyprimary: bool,
        minqual: u8,
        minscore: u8,
        reversedel: bool,
        outputends: bool,
    ) -> Self {
        Self {
            onlyprimary,
            recalculateBAM: false,
            minqual,
            minscore,
            reversedel,
            outputends,
        }
    }
    /// Get BAM recalculation, is not implemented so cannot be changed and must be false
    #[allow(non_snake_case, unused)]
    pub fn getBAMrecalculate(&self) -> bool {
        self.recalculateBAM
    }
    pub fn isonlyprimary(&self) -> bool {
        self.onlyprimary
    }
    #[allow(unused)]
    pub fn setonlyprimary(&mut self, primary: bool) {
        self.onlyprimary = primary;
    }
    pub fn getminqual(&self) -> u8 {
        self.minqual
    }
    #[allow(unused)]
    pub fn setminqual(&mut self, minqual: u8) {
        self.minqual = minqual;
    }
    pub fn getminscore(&self) -> u8 {
        self.minscore
    }
    #[allow(unused)]
    pub fn setminscore(&mut self, minscore: u8) {
        self.minscore = minscore;
    }
    pub fn hasreversedel(&self) -> bool {
        self.reversedel
    }
    #[allow(unused)]
    pub fn setreversedel(&mut self, reversedel: bool) {
        self.reversedel = reversedel;
    }
    pub fn hasoutputends(&self) -> bool {
        self.outputends
    }
    #[allow(unused)]
    pub fn setoutputends(&mut self, outputends: bool) {
        self.outputends = outputends;
    }
}
#[derive(Debug)]
pub struct RustPileups<'a> {
    pileup: HashMap<(String, i64), RustPileup>,
    pub config: RustPileupConfig,
    region: FetchDefinition<'a>,
    assemblyreader: Option<faidx::Reader>,
}
#[derive(Debug, Eq, PartialEq, Clone)]
pub struct RustPileup {
    nbam: u64,
    index: Option<u64>,
    pub chrom_name: String,
    pub pos: u64,
    pub base: char,
    pub nreads: Vec<u64>,
    pub rbases: Vec<String>,
    pub qualities: Vec<Vec<u8>>,
}
impl Ord for RustPileup {
    fn cmp(&self, other: &Self) -> Ordering {
        match self.getchrom().cmp(&other.getchrom()) {
            Ordering::Equal => self.getpos().cmp(&other.pos),
            ord => ord,
        }
    }
}
impl PartialOrd for RustPileup {
    fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
        Some(self.cmp(&other))
    }
}
macro_rules! zip {
    ($x: expr) => ($x);
    ($x: expr, $($y: expr), +) => (
        $x.iter().zip(
            zip!($($y), +))
    )
}

impl RustPileup {
    /// Get number of bam
    pub fn getbam(&self) -> u64 {
        self.nbam
    }
    /// Get current index
    fn getindex(&self) -> Option<u64> {
        self.index
    }
    /// Get current index unchecked
    fn getindexunchecked(&self) -> u64 {
        self.index.unwrap()
    }
    fn setindex(&mut self, index: u64) -> bam::Result<()> {
        if index >= self.nbam {
            return Err(bam::Error::BamPileup);
        } else {
            self.index = Some(index);
            return Ok(());
        }
    }
    /// Set number of bam
    pub fn setbam(&mut self, num: u64, alignindex: bool) {
        self.nbam = num;
        if alignindex {
            self.setindex(num.saturating_sub(1))
                .unwrap_or_else(|_| unreachable!("Not reachable"));
        }
    }
    /// Get chromosome name
    pub fn getchrom(&self) -> &str {
        &self.chrom_name
    }
    /// Set chromosome name
    pub fn setchro(&mut self, chro: String) {
        self.chrom_name = chro;
    }
    /// Get position
    pub fn getpos(&self) -> u64 {
        self.pos
    }
    /// Set position
    pub fn setpos(&mut self, pos: u64) {
        self.pos = pos;
    }
    /// Get base (N if no assembly provided)
    pub fn getbase(&self) -> char {
        self.base
    }
    /// Set base
    pub fn setbase(&mut self, char: char) {
        self.base = char;
    }
    /// Get number of reads
    pub fn getnreads(&self) -> &[u64] {
        &self.nreads
    }
    /// Set number of reads, return error if index does not exist
    pub fn setnreads(&mut self, index: Option<usize>, reads: u64) -> Result<()> {
        match index {
            Some(a) if let Some(b) = self.nreads.get_mut(a) => {
                *b = reads;
                Ok(())
            }
            Some(b) if b == self.getnreads().len() => {
                self.nreads.push(reads);
                Ok(())
            }
            None => {
                self.nreads.push(reads);
                Ok(())
            }
            Some(_) => Err(crate::errors::Error::BamInvalidRecord),
        }
    }
    /// Get base reads information
    pub fn getbasereads(&self) -> &[String] {
        &self.rbases
    }
    /// Set base reads information
    pub fn setbasereads(&mut self, index: Option<usize>, reads: String) -> Result<()> {
        match index {
            Some(a) if let Some(b) = self.rbases.get_mut(a) => {
                *b = reads;
                Ok(())
            }
            Some(b) if b == self.rbases.len() => {
                self.rbases.push(reads);
                Ok(())
            }
            Some(_) => Err(crate::errors::Error::BamInvalidRecord),
            None => {
                self.rbases.push(reads);
                Ok(())
            }
        }
    }
    /// Get qualities
    pub fn getqualities(&self) -> &[Vec<u8>] {
        &self.qualities
    }
    /// Add qualities
    pub fn addqualities(&mut self, index: Option<usize>, qual: u8) -> Result<()> {
        match index {
            Some(a) if let Some(b) = self.qualities.get_mut(a) => {
                b.push(qual);
                Ok(())
            }
            Some(b) if b == self.getqualities().len() => {
                self.qualities.push(vec![qual]);
                Ok(())
            }
            Some(_) => Err(crate::errors::Error::BamInvalidRecord),
            None => {
                self.qualities.push(vec![qual]);
                Ok(())
            }
        }
    }
    /// Set qualities
    pub fn setqualities(&mut self, index: Option<usize>, qual: Vec<u8>) -> Result<()> {
        match index {
            Some(a) if let Some(b) = self.qualities.get_mut(a) => {
                *b = qual;
                Ok(())
            }
            Some(b) if b == self.getqualities().len() => {
                self.qualities.push(qual);
                Ok(())
            }
            Some(_) => Err(crate::errors::Error::BamInvalidRecord),
            None => {
                self.qualities.push(qual);
                Ok(())
            }
        }
    }
    /// Reset qualities
    pub fn resetqualities(&mut self, index: usize) -> Result<()> {
        if let Some(b) = self.qualities.get_mut(index) {
            b.clear();
            b.shrink_to_fit();
            Ok(())
        } else {
            Err(crate::errors::Error::BamInvalidRecord)
        }
    }
}
impl Default for RustPileup {
    fn default() -> Self {
        Self {
            nbam: 0,
            index: None,
            chrom_name: "N/A".to_string(),
            pos: 0,
            base: 'N',
            nreads: Vec::new(),
            rbases: Vec::new(),
            qualities: Vec::new(),
        }
    }
}
impl Display for RustPileup {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        let mut text = String::new();
        text.push_str(&format!("{}\t", self.getchrom()));
        text.push_str(&format!("{}\t", self.getpos()));
        text.push_str(&format!("{}\t", self.getbase()));
        for (bread, (nreads, qual)) in
            zip!(self.getbasereads(), self.getnreads(), self.getqualities())
        {
            let qual = if !qual.is_empty() {
                qual.iter().fold(String::new(), |mut acc, qual| {
                    let val = if *qual == u8::MAX { 0 } else { *qual };
                    let val = char::from_u32(u32::from(val).saturating_add(33)).unwrap_or('!');
                    acc.push_str(&format!("{}", val));
                    acc
                })
            } else {
                "*".to_string()
            };
            let bread = if bread.trim().len() == 0 {
                "*".to_string()
            } else {
                bread.to_string()
            };
            text.push_str(&format!("{}\t{}\t{}\t", nreads, bread, qual));
        }
        text = text.trim().to_string();
        write!(f, "{}", text)
    }
}
impl<'a> Display for RustPileups<'a> {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        let mut text = String::new();
        for pile in self.into_iter() {
            text.push_str(&format!("{}\n", pile.to_string()));
        }
        write!(f, "{}", text.trim())
    }
}
#[must_use]
fn casereverse<T>(record: &bam::Record, text: T) -> String
where
    T: AsRef<str>,
{
    if record.is_reverse() {
        text.as_ref().to_ascii_lowercase()
    } else {
        text.as_ref().to_ascii_uppercase()
    }
}
#[must_use]
fn casecharreverse(record: &bam::Record, text: char) -> char {
    if record.is_reverse() {
        text.to_ascii_lowercase()
    } else {
        text.to_ascii_uppercase()
    }
}
/// Splits [start, max) into up to n_chunks contiguous, non-overlapping sub-ranges.
#[cfg(feature = "pileuprayon")]
fn split_range(start: i64, max: i64, n_chunks: usize) -> Vec<(i64, i64)> {
    let n_chunks = n_chunks.max(1);
    let total = max.saturating_sub(start).max(0) as usize;
    if total == 0 {
        return vec![(start, max)];
    }
    let chunk_size = ((total + n_chunks - 1) / n_chunks) as i64; // div_ceil
    let mut ranges = Vec::new();
    let mut cur = start;
    while cur < max {
        let end = (cur + chunk_size).min(max);
        ranges.push((cur, end));
        cur = end;
    }
    ranges
}

/// Processes one sub-region: opens its own reader, fetches overlapping reads,
/// and builds a local pileup map for positions strictly inside [chunk_start, chunk_max).
#[cfg(feature = "pileuprayon")]
fn process_chunk<'a, T>(
    bam_path: T,
    bam_index: Option<T>,
    assembly_path: Option<T>,
    assembly_path_index: Option<T>,
    config: Option<RustPileupConfig>,
    actualpileup: Option<&mut RustPileups>,
    region: FetchDefinition<'a>,
) -> bam::Result<Option<RustPileups<'a>>>
where
    T: AsRef<Path> + Send + Sync + Clone,
{
    let (contigname, chunk_start, chunk_max) = match region.clone() {
        RegionString(a, b, c) => (String::from_utf8_lossy(a), b, c),
        _ => return Err(bam::Error::Fetch),
    };
    let (reader, assemblyreader) =
        getreader(bam_path, bam_index, assembly_path, assembly_path_index)?;
    match (config, actualpileup) {
        (Some(config), None) => {
            let mut pileup = RustPileups {
                pileup: HashMap::new(),
                region,
                config,
                assemblyreader: None,
            };
            for i in chunk_start..chunk_max {
                pileup.pileup.insert(
                    (contigname.to_string(), i.saturating_sub(1)),
                    RustPileup::default(),
                );
            }
            pileup.goonrecord(reader, assemblyreader)?;
            Ok(Some(pileup))
        }
        (None, Some(b)) => {
            b.goonrecord(reader, assemblyreader)?;
            Ok(None)
        }
        _ => return Err(bam::Error::BamInvalidRecord),
    }
}
impl<'a> RustPileups<'a> {
    /// Add another BAM to existing Pileup config
    pub fn addextrabam<T>(&mut self, reader: T, index: Option<T>) -> bam::Result<()>
    where
        T: AsRef<Path>,
    {
        let (name, start, max) = match &self.region {
            RegionString(name, start, end) => (name, *start, *end),
            _ => return Err(bam::Error::Fetch),
        };
        //#[cfg(not(feature = "pileuprayon"))]
        {
            let mut reader = getreader(reader, index, None, None).map(|(a, _)| a)?;
            for i in start..=max {
                if let Some(a) = self.pileup.get_mut(&(
                    String::from_utf8_lossy(name).to_string(),
                    i.saturating_sub(1),
                )) {
                    a.setbam(a.getbam() + 1, true);
                    a.setnreads(a.getindex().and_then(|f| f.try_into().ok()), 0)?;
                } else {
                    return Err(Error::BamInvalidRecord);
                }
            }
            let _ = reader.set_threads(
                available_parallelism()
                    .unwrap_or(std::num::NonZero::new(1).unwrap())
                    .get(),
            );
            self.goonrecord(reader, None)?;
        }
        /*
        Cannot share faidx
        #[cfg(feature = "pileuprayon")]
        {
            let pileup = {
                let n_chunks = rayon::current_num_threads();
                let chunk_ranges = split_range(start, max, n_chunks);
                let chunk_pileups = chunk_ranges
                    .par_iter()
                    .map(|(chunk_start, chunk_max)| {
                        let region = FetchDefinition::RegionString(name, *chunk_start, *chunk_max);
                        let val = process_chunk(
                            reader.as_ref().clone(),
                            index.map(|b| b.as_ref()),
                            None,
                            None,
                            None,
                            Some(self),
                            region,
                        );
                        val
                    })
                    .collect::<Result<Vec<_>, _>>()?;
                let mut row = chunk_pileups.into_iter();
                let mut pileup: RustPileups = row.next().unwrap();
                pileup.region = region;
                for chunk_map in row {
                    pileup.extend(chunk_map);
                }
                pileup
            };
        }
        */
        Ok(())
    }
    fn goonrecord(
        &mut self,
        mut reader: bam::IndexedReader,
        assemblyreader: Option<faidx::Reader>,
    ) -> bam::Result<()> {
        let pileup = &mut self.pileup;
        let config = &self.config;
        let mut record = bam::Record::new();
        reader.fetch(self.region.clone())?;
        match (self.assemblyreader.is_some(), assemblyreader) {
            (_, Some(a)) => {
                self.assemblyreader = Some(a);
            }
            (true, _) => (),
            (false, None) => {
                self.assemblyreader = None;
            }
        };
        let (name, start, max) = match self.region.clone() {
            RegionString(a, b, c) => (String::from_utf8_lossy(a), b, c),
            _ => return Err(bam::Error::Fetch),
        };
        while let Some(v) = reader.read(&mut record) {
            if v.is_err() {
                continue;
            }
            if (config.isonlyprimary() && !record.is_primary())
                || (record.mapq() < config.getminscore())
            {
                continue;
            }
            let records: Vec<([Option<i64>; 2], Cigar, u8)> =
                IterAlignedPairsFullCigar::new(record.aligned_pairs_full(), &record).collect();
            let mut currentgread = 0;
            for (rangeindex, ([rread, gread], cigar, qual)) in records.iter().enumerate() {
                let gread = match gread {
                    Some(a) => *a,
                    None if currentgread != 0 => currentgread,
                    _ => continue,
                };
                if gread < start.saturating_sub(1) {
                    continue;
                }
                if gread > max.saturating_sub(1) {
                    break;
                }
                let rread = match rread {
                    Some(a) => Some(*a),
                    None => records
                        .iter()
                        .skip(rangeindex)
                        .find(|p| &p.1 != cigar)
                        .and_then(|([rread, _], _, _)| match rread {
                            Some(a) => Some(*a),
                            None => None,
                        }),
                };
                currentgread = gread;
                let base = match &self.assemblyreader {
                    Some(reader) => reader
                        .fetch_seq_string(
                            &name,
                            gread.try_into().unwrap_or_default(),
                            gread.try_into().unwrap_or_default(),
                        )
                        .map_or('n', |d| d.chars().next().unwrap_or('n')),
                    _ => 'n',
                };
                let entry = pileup
                    .entry((name.to_string(), gread))
                    .or_insert(RustPileup::default());
                if entry.getindex().is_none() && entry.getbam() == 0 {
                    entry.setbam(1, true);
                    entry.setchro(name.to_string());
                    entry.setbase(base);
                    entry.setpos(gread.saturating_add(1).try_into().unwrap_or(u64::MIN));
                }
                // hit is now per-entry instead of function-scoped:
                let hit = entry.getindex().map(|a| a as usize);
                let index = entry.getindexunchecked() as usize;
                if entry.getnreads().get(index).is_none_or(|p| p == &0) {
                    entry.setbasereads(Some(index), String::new())?;
                    entry.setqualities(Some(index), vec![])?;
                    if entry.getnreads().get(index).is_none() {
                        entry.setnreads(Some(index), 0)?;
                    }
                }
                let initial: u64 = match hit {
                    Some(a) => entry.getnreads().get(a).copied().unwrap_or(0),
                    None => 0,
                };
                let mut readbase = if let Some(rread) = rread {
                    record
                        .seq()
                        .rangeextract(rread as usize..rread.saturating_add(1) as usize)
                        .unwrap_or(std::borrow::Cow::Borrowed("n"))
                        .chars()
                        .next()
                        .unwrap_or('n')
                } else {
                    'n'
                };
                readbase = casecharreverse(&record, readbase);
                if *qual <= config.getminqual() {
                    continue;
                }
                if matches!(cigar, Cigar::Ins(_))
                    && let Some(_rread) = rread
                {
                    () //Would be counted afterwards
                } else {
                    entry.setnreads(hit, initial.saturating_add(1))?;
                    entry.addqualities(hit, *qual)?;
                }
                if record.reference_start() == gread && config.hasoutputends() {
                    let mut info = entry
                        .getbasereads()
                        .get(index)
                        .map_or(String::new(), |f| f.to_string());
                    info.push_str(&format!(
                        "^{}",
                        char::from_u32(u32::from(record.mapq().saturating_add(33))).unwrap_or('!')
                    ));
                    entry.setbasereads(hit, info)?;
                } else if record.reference_end() == gread && config.hasoutputends() {
                    let mut info = entry
                        .getbasereads()
                        .get(index)
                        .map_or(String::new(), |f| f.to_string());
                    info.push_str(&format!(
                        "{}$",
                        char::from_u32(u32::from(record.mapq().saturating_add(33))).unwrap_or('!')
                    ));
                    entry.setbasereads(hit, info)?;
                }
                match cigar {
                    Cigar::Equal(_) | Cigar::Match(_) if readbase.eq_ignore_ascii_case(&base) => {
                        let mut info = entry
                            .getbasereads()
                            .get(index)
                            .map_or(String::new(), |f| f.to_string());
                        info.push_str(if record.is_reverse() { "," } else { "." });
                        entry.setbasereads(Some(index), info)?;
                    }
                    /* Cigar::Equal(_) => {
                        return Err(bam::Error::BamParseCigar {
                            msg: format!(
                                "Equal in cigar does not match reference at position {} for read {}",
                                gread,
                                String::from_utf8_lossy(record.name())
                            ),
                        });
                    } */
                    Cigar::Diff(_) | Cigar::Match(_) | Cigar::Equal(_) => {
                        let mut info = entry
                            .getbasereads()
                            .get(index)
                            .map_or(String::new(), |f| f.to_string());
                        info.push_str(&readbase.to_string());
                        entry.setbasereads(Some(index), info)?;
                    }
                    Cigar::Ins(n) if let Some(_rread) = rread => {
                        //Would duplicate
                        ()
                        /* let mut info = entry
                            .getbasereads()
                            .get(index)
                            .map_or(String::new(), |f| f.to_string());
                        let recordseq = record.seq();
                        let val = recordseq
                            .rangeextract(rread as usize..rread.saturating_add(*n as i64) as usize)
                            .unwrap_or(std::borrow::Cow::Owned(
                                "N".repeat((*n).try_into().unwrap_or_default()),
                            ));
                        let val: std::borrow::Cow<str> =
                            std::borrow::Cow::Owned(casereverse(&record, val));
                        info.push_str(&format!("+{}{}", n, val));
                        entry.setbasereads(Some(index), info)?; */
                    }
                    Cigar::Ins(_) => {
                        return Err(bam::Error::BamParseCigar {
                            msg: format!(
                                "Insertion in cigar does not match reference at position {} for read {}",
                                gread,
                                String::from_utf8_lossy(record.name())
                            ),
                        });
                    }
                    Cigar::Del(_) => {
                        let mut info = entry
                            .getbasereads()
                            .get(index)
                            .map_or(String::new(), |f| f.to_string());
                        if record.is_reverse() && config.hasreversedel() {
                            info.push_str("#");
                        } else {
                            info.push_str("*");
                        }
                        entry.setbasereads(Some(index), info)?;
                    }
                    _ => {
                        let mut info = entry
                            .getbasereads()
                            .get(index)
                            .map_or(String::new(), |f| f.to_string());
                        info.push_str("N");
                        entry.setbasereads(Some(index), info)?;
                    }
                }
                if let Some(([rstart, gstart], cigar, _)) =
                    records.iter().skip(rangeindex).find(|p| &p.1 != cigar)
                    && (gstart.is_some_and(|f| f.abs_diff(gread) <= 1)
                        || rstart.is_some_and(|f| f.abs_diff(rread.unwrap_or_default()) <= 1))
                //Last range before next record
                {
                    match cigar {
                        Cigar::Del(n) => {
                            let val = match (&self.assemblyreader, gstart) {
                                (Some(reader), Some(gstart)) => reader
                                    .fetch_seq_string(
                                        &name,
                                        (*gstart).try_into().unwrap_or_default(),
                                        gstart
                                            .saturating_add(
                                                (*n).saturating_sub(1)
                                                    .try_into()
                                                    .unwrap_or_default(),
                                            )
                                            .try_into()
                                            .unwrap(),
                                    )
                                    .unwrap_or(
                                        "N".repeat((*n).try_into().unwrap_or_default()).to_string(),
                                    ),
                                _ => "N".repeat((*n).try_into().unwrap_or_default()).to_string(),
                            };
                            let val: std::borrow::Cow<str> =
                                std::borrow::Cow::Owned(casereverse(&record, val));
                            let mut info = entry
                                .getbasereads()
                                .get(index)
                                .map_or(String::new(), |f| f.to_string());
                            info.push_str(&format!("-{}{}", n, val));
                            entry.setbasereads(Some(index), info)?;
                        }
                        Cigar::Ins(n) => {
                            let mut info = entry
                                .getbasereads()
                                .get(index)
                                .map_or(String::new(), |f| f.to_string());
                            let recordseq = record.seq();
                            let val = match rstart {
                                Some(rstart) => recordseq
                                    .rangeextract(
                                        *rstart as usize
                                            ..rstart.saturating_add((*n) as i64) as usize,
                                    )
                                    .unwrap_or(std::borrow::Cow::Owned(
                                        "N".repeat((*n).try_into().unwrap_or_default()),
                                    )),
                                _ => std::borrow::Cow::Owned(
                                    "N".repeat((*n).try_into().unwrap_or_default()),
                                ),
                            };
                            let val: std::borrow::Cow<str> =
                                std::borrow::Cow::Owned(casereverse(&record, val));
                            info.push_str(&format!("+{}{}", n, val));
                            entry.setbasereads(Some(index), info)?;
                        }
                        _ => (),
                    }
                }
            }
        }
        Ok(())
    }
    /// Generate a new pileup from BAM and assembly. Recalculate BAM is not supported and would panic.
    /// **In benchmark this function does 400x worse/slower than samtools and is, at this stage, not recommended.**
    /// # Examples
    /// ```
    /// use extended_htslib::bam::pileup::{Pileup, RustPileup, RustPileups,RustPileupConfig};
    /// use extended_htslib::bam::{FetchDefinition, Read};
    /// use extended_htslib::{bam, faidx};
    /// use std::fs;
    ///
    ///
    /// let mut config = RustPileupConfig::default();
    /// let pos = FetchDefinition::RegionString("chr14".as_bytes(), 99812480, 99815480);
    /// config.setminqual(0);
    /// assert_eq!(config.getminqual(),0);
    /// let po = RustPileups::new("test/locus.bam", None, Some("test/locus.fasta"), None, pos, config).unwrap();
    /// # println!("{}",po.to_string());
    /// assert_eq!(po.into_iter().next().unwrap().to_string(),"chr14	99812480	G	1	^],	E");
    /// ```
    /// # Errors
    /// Bad region.
    pub fn new<T>(
        #[allow(unused_mut)] mut reader: T,
        #[allow(unused_mut)] mut index: Option<T>,
        assemblyreader: Option<T>,
        #[allow(unused_mut)] mut assemblyreaderindex: Option<T>,
        region: bam::FetchDefinition<'a>,
        config: RustPileupConfig,
    ) -> bam::Result<Self>
    where
        T: AsRef<Path> + Send + Sync + Clone,
    {
        if config.recalculateBAM {
            unimplemented!("Recalculate BAM is not implemented");
        }
        let (name, start, max) = match &region {
            RegionString(name, start, end) => (name, *start, *end),
            _ => return Err(bam::Error::Fetch),
        };
        #[cfg(not(feature = "pileuprayon"))]
        let pileup = {
            let (reader, assemblyindex) =
                getreader(reader, index, assemblyreader, assemblyreaderindex)?;
            let mut pileup = RustPileups {
                pileup: HashMap::new(),
                region: region.clone(),
                config,
                assemblyreader: None,
            };
            for i in start..=max {
                pileup.pileup.insert(
                    (
                        String::from_utf8_lossy(name).to_string(),
                        i.saturating_sub(1),
                    ),
                    RustPileup::default(),
                );
            }
            pileup.goonrecord(reader, assemblyindex)?;
            pileup
        };
        #[cfg(feature = "pileuprayon")]
        let pileup = {
            let n_chunks = rayon::current_num_threads();
            let chunk_ranges = split_range(start, max, n_chunks);
            let chunk_pileups = chunk_ranges
                .par_iter()
                .map(|(chunk_start, chunk_max)| {
                    let region = FetchDefinition::RegionString(name, *chunk_start, *chunk_max);
                    let val = process_chunk(
                        reader.clone(),
                        index.clone(),
                        assemblyreader.clone(),
                        assemblyreaderindex.clone(),
                        Some(config),
                        None,
                        region,
                    );
                    val
                })
                .collect::<Result<Vec<_>, _>>()?
                .into_iter()
                .collect::<Option<Vec<_>>>()
                .ok_or_else(|| bam::Error::BamInvalidRecord)?;
            let mut row = chunk_pileups.into_iter();
            let mut pileup: RustPileups = row.next().unwrap();
            pileup.region = region;
            for chunk_map in row {
                pileup.extend(chunk_map);
            }
            pileup
        };
        Ok(pileup)
    }
}
impl<'a> Extend<RustPileup> for RustPileups<'a> {
    fn extend<T: IntoIterator<Item = RustPileup>>(&mut self, iter: T) {
        for elem in iter {
            match self.pileup.entry((
                elem.chrom_name.clone(),
                elem.pos.saturating_sub(1).try_into().unwrap_or_default(),
            )) {
                Vacant(a) => {
                    a.insert(elem);
                }
                Occupied(mut b) => {
                    let entry = b.get_mut();
                    if *entry == elem {
                        continue;
                    }
                    let index = entry.getindex().and_then(|a| usize::try_from(a).ok());
                    if let Some(newindex) = index
                        && let (Some(nread), Some(bread), Some(mut qual)) = (
                            elem.getnreads().first(),
                            elem.getbasereads().first(),
                            elem.qualities.clone().get_mut(0),
                        )
                    {
                        let _ = entry.setnreads(
                            index,
                            entry.getnreads().get(newindex).map_or(0, |f| *f) + (*nread),
                        );
                        let _ = entry.setbasereads(
                            index,
                            format!(
                                "{}{}",
                                entry
                                    .getbasereads()
                                    .get(newindex)
                                    .map_or(String::new(), |f: &String| f.to_string()),
                                bread
                            ),
                        );
                        let mut new = Vec::new();
                        let veco = entry.qualities.get_mut(newindex).unwrap_or(&mut new);
                        veco.append(&mut qual);
                    }
                }
            }
        }
    }
}
impl<'a> IntoIterator for RustPileups<'a> {
    type IntoIter = std::vec::IntoIter<Self::Item>;
    type Item = RustPileup;
    fn into_iter(self) -> Self::IntoIter {
        let mut vec: Vec<RustPileup> = self.pileup.into_values().collect();
        if !vec.is_sorted() {
            vec.sort_unstable();
        }
        vec.into_iter()
    }
}
impl<'a> IntoIterator for &RustPileups<'a> {
    type IntoIter = std::vec::IntoIter<Self::Item>;
    type Item = &'a RustPileup;
    fn into_iter(self) -> Self::IntoIter {
        let mut vec: Vec<&'a RustPileup> = self
            .pileup
            .values()
            .map(|v| unsafe { &*(v as *const RustPileup) })
            .collect();
        if !vec.is_sorted() {
            vec.sort_unstable();
        }
        vec.into_iter()
    }
}
/// Iterator over alignments of a pileup.
pub type Alignments<'a> = iter::Map<
    slice::Iter<'a, htslib::bam_pileup1_t>,
    fn(&'a htslib::bam_pileup1_t) -> Alignment<'a>,
>;

/// A pileup over one genomic position.
#[derive(Debug)]
pub struct Pileup {
    inner: *const htslib::bam_pileup1_t,
    depth: u32,
    tid: u32,
    pos: u32,
}

impl Pileup {
    pub fn tid(&self) -> u32 {
        self.tid
    }

    pub fn pos(&self) -> u32 {
        self.pos
    }

    pub fn depth(&self) -> u32 {
        self.depth
    }

    pub fn alignments(&self) -> Alignments<'_> {
        self.inner().iter().map(Alignment::new)
    }

    fn inner(&self) -> &[htslib::bam_pileup1_t] {
        unsafe {
            slice::from_raw_parts(
                self.inner as *mut htslib::bam_pileup1_t,
                self.depth as usize,
            )
        }
    }
}

/// An aligned read in a pileup.
pub struct Alignment<'a> {
    inner: &'a htslib::bam_pileup1_t,
}

impl fmt::Debug for Alignment<'_> {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "Alignment")
    }
}

impl<'a> Alignment<'a> {
    pub fn new(inner: &'a htslib::bam_pileup1_t) -> Self {
        Alignment { inner }
    }

    /// Position within the read. None if either `is_del` or `is_refskip`.
    pub fn qpos(&self) -> Option<usize> {
        if self.is_del() || self.is_refskip() {
            // there is no alignment position in such a case
            None
        } else {
            Some(self.inner.qpos as usize)
        }
    }

    /// Insertion, deletion (with length) if indel starts at next base or None otherwise.
    pub fn indel(&self) -> Indel {
        match self.inner.indel {
            len if len < 0 => Indel::Del(-len as u32),
            len if len > 0 => Indel::Ins(len as u32),
            _ => Indel::None,
        }
    }

    /// Whether there is a deletion in the alignment at this position.
    pub fn is_del(&self) -> bool {
        self.inner.is_del() != 0
    }

    /// Whether the alignment starts at this position.
    pub fn is_head(&self) -> bool {
        self.inner.is_head() != 0
    }

    /// Whether the alignment ends at this position.
    pub fn is_tail(&self) -> bool {
        self.inner.is_tail() != 0
    }

    /// Whether this position is marked as refskip in the CIGAR string.
    pub fn is_refskip(&self) -> bool {
        self.inner.is_refskip() != 0
    }

    /// The corresponding record.
    pub fn record(&self) -> record::Record {
        record::Record::from_inner(self.inner.b)
    }
}

#[derive(PartialEq, Eq, Debug, Copy, Clone, Hash)]
pub enum Indel {
    Ins(u32),
    Del(u32),
    None,
}

/// Iterator over pileups.
#[derive(Debug)]
pub struct Pileups<'a, R: bam::Read> {
    #[allow(dead_code)]
    reader: &'a mut R,
    itr: htslib::bam_plp_t,
}

impl<'a, R: bam::Read> Pileups<'a, R> {
    pub fn new(reader: &'a mut R, itr: htslib::bam_plp_t) -> Self {
        Pileups { reader, itr }
    }

    /// Warning: because htslib internally uses signed integer for depth this method
    /// will panic if `depth` exceeds `i32::MAX`.
    pub fn set_max_depth(&mut self, depth: u32) {
        if depth > i32::MAX as u32 {
            panic!(
                "Maximum value for pileup depth is {} but {} was provided",
                i32::MAX,
                depth
            )
        }
        let intdepth = depth as i32;
        unsafe {
            htslib::bam_plp_set_maxcnt(self.itr, intdepth);
        }
    }
}

impl<R: bam::Read> Iterator for Pileups<'_, R> {
    type Item = Result<Pileup>;

    #[allow(clippy::match_bool)]
    fn next(&mut self) -> Option<Result<Pileup>> {
        let (mut tid, mut pos, mut depth) = (0i32, 0i32, 0i32);
        let inner = unsafe { htslib::bam_plp_auto(self.itr, &mut tid, &mut pos, &mut depth) };

        match inner.is_null() {
            true if depth == -1 => Some(Err(Error::BamPileup)),
            true => None,
            false => Some(Ok(Pileup {
                inner,
                depth: depth as u32,
                tid: tid as u32,
                pos: pos as u32,
            })),
        }
    }
}

impl<R: bam::Read> Drop for Pileups<'_, R> {
    fn drop(&mut self) {
        unsafe {
            htslib::bam_plp_reset(self.itr);
            htslib::bam_plp_destroy(self.itr);
        }
    }
}

#[cfg(test)]
mod tests {

    use std::fs::{self};
    use std::path::PathBuf;

    use crate::bam;
    use crate::bam::pileup::{RustPileupConfig, RustPileups};
    use crate::bam::{FetchDefinition, Read};
    #[test]
    fn testpileup() {
        let bam = PathBuf::from("test/locus.bam");
        let fasta = PathBuf::from("test/locus.fasta");
        // samtools mpileup --reverse-del -a -B -Q 0 -q 0 -r "chr14:99812480-99815480" -f test/locus.fasta test/locus.bam > test/locus.pileup
        let pos = FetchDefinition::RegionString("chr14".as_bytes(), 99812480, 99815480);
        let config = RustPileupConfig::default();
        let po = RustPileups::new(bam, None, Some(fasta), None, pos, config).unwrap();
        let val = fs::read_to_string("test/locus.pileup").unwrap();
        assert_eq!(val.trim(), po.to_string(), "Pileup does not match");
    }
    #[test]
    fn test_multi_pileup() {
        let bam = PathBuf::from("test/locus.bam");
        let bam2 = "test/locusmd.bam";
        let fasta = PathBuf::from("test/locus.fasta");
        // samtools mpileup --reverse-del -a -B -Q 0 -q 0 -r "chr14:99812480-99815480" -f test/locus.fasta test/locus.bam test/locusmd.bam > test/locusmulti.pileup
        let pos = FetchDefinition::RegionString("chr14".as_bytes(), 99812480, 99815480);
        let config = RustPileupConfig::default();
        let mut po = RustPileups::new(bam, None, Some(fasta), None, pos, config).unwrap();
        po.addextrabam(bam2, None).unwrap();
        assert_eq!(
            fs::read_to_string("test/locusmulti.pileup").unwrap().trim(),
            po.to_string(),
            "Pileup does not match"
        );
    }
    #[test]
    fn test_max_pileup() {
        let mut bam = bam::Reader::from_path("test/test.bam").unwrap();
        let mut p = bam.pileup();
        p.set_max_depth(0u32);
        p.set_max_depth(800u32);
    }

    #[test]
    #[should_panic]
    fn test_max_pileup_to_high() {
        let mut bam = bam::Reader::from_path("test/test.bam").unwrap();
        let mut p = bam.pileup();
        p.set_max_depth((i32::MAX as u32) + 1);
    }
}