twitcher 0.6.9

Find template switch mutations in genomic data
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117

use std::ops::Deref;

use anyhow::Context as _;
use rust_htslib::bcf;

use crate::{common::coords::GenomeRegion, vcf::pipeline::clusterizer::phasing::GtAndPhase};

/// A record as it was read from the input file, together with the genotype that this program
/// considers effective for it.
///
/// The wrapped record is bound to the *input* file's header, which this program does not own
/// and must not extend: a FORMAT field the input does not declare cannot be written onto such
/// a record. Therefore [`InputRecord`] hands out shared access only — `Deref` to
/// [`bcf::Record`] exposes every reading accessor, while the `&mut` needed by the `push_*`
/// family is unreachable. Records are made writable exclusively by the writer, which rebinds
/// them to the output header (see `OutputWriter::adopt`).
///
/// Consequently, genotype corrections found by local read-based phasing are not written back
/// onto the record; they are carried in `local_gt` and take effect through [`Self::gt`].
/// [`Self::input_gt`] keeps reporting what the input file said, so the writer can tell whether
/// the emitted phasing differs from the input.
#[derive(Debug, Clone)]
pub struct InputRecord {
    record: bcf::Record,
    input_gt: GtAndPhase,
    local_gt: Option<GtAndPhase>,
}

impl InputRecord {
    /// Wrap a record read from the input file, parsing its genotype once.
    pub fn new(record: bcf::Record) -> Self {
        let input_gt = GtAndPhase::from_record(&record);
        Self {
            record,
            input_gt,
            local_gt: None,
        }
    }

    /// Wrap a record whose genotype has already been parsed, avoiding a second parse.
    pub const fn with_gt(record: bcf::Record, input_gt: GtAndPhase) -> Self {
        Self {
            record,
            input_gt,
            local_gt: None,
        }
    }

    /// The genotype to act on: the one resolved by local phasing if there is one, else the
    /// input file's.
    pub const fn gt(&self) -> GtAndPhase {
        match self.local_gt {
            Some(gt) => gt,
            None => self.input_gt,
        }
    }

    /// The genotype as declared by the input file, regardless of local phasing.
    pub const fn input_gt(&self) -> GtAndPhase {
        self.input_gt
    }

    /// A copy of this record with the genotype resolved by local read-based phasing attached.
    pub fn with_local_gt(&self, local_gt: GtAndPhase) -> Self {
        Self {
            record: self.record.clone(),
            input_gt: self.input_gt,
            local_gt: Some(local_gt),
        }
    }

    /// The underlying record. Shared access only, see the type-level documentation.
    pub const fn record(&self) -> &bcf::Record {
        &self.record
    }

    /// Hand the record over to the writer, which rebinds it to the output header.
    pub(super) fn into_record(self) -> bcf::Record {
        self.record
    }
}

impl Deref for InputRecord {
    type Target = bcf::Record;

    fn deref(&self) -> &Self::Target {
        &self.record
    }
}

impl TryFrom<&InputRecord> for GenomeRegion {
    type Error = anyhow::Error;

    fn try_from(rec: &InputRecord) -> Result<Self, Self::Error> {
        Self::try_from(rec.record())
    }
}

impl TryFrom<&[InputRecord]> for GenomeRegion {
    type Error = anyhow::Error;

    fn try_from(recs: &[InputRecord]) -> Result<Self, Self::Error> {
        if recs.is_empty() {
            anyhow::bail!("empty list of records");
        }
        let min_pos = recs.iter().map(|r| r.pos()).min().context("not empty")?;
        let max_end = recs.iter().map(|r| r.end()).max().context("not empty")?;
        let rec = recs.first().context("not empty")?;
        let chr = rec
            .header()
            .rid2name(rec.rid().context("no rid for record.")?)
            .context("read error")?;
        let gpos =
            crate::common::coords::GenomePosition::new_0(chr.into(), usize::try_from(min_pos)?);
        Ok(Self::new_bounded(gpos, usize::try_from(max_end - min_pos)?))
    }
}