twitcher 0.6.0

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
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157

use async_stream::stream;
use rust_htslib::bcf;
use tokio_stream::Stream;
use tracing::error;

use crate::{
    common::{
        coords::GenomeRegion,
        list_of_regions::{TargetQueryAnswer, Targets},
    },
    counter,
    vcf::pipeline::{
        clusterizer::{cluster::ClusteringSettings, phasing::GtAndPhase},
        reader::VCFReader,
    },
};

struct State {
    strategy: ClusteringSettings,
    record_buf: bcf::Record,
    buffer: Vec<bcf::Record>,
    last_rid: Option<u32>,
}

impl State {
    const fn new(strategy: ClusteringSettings, record_buf: bcf::Record) -> Self {
        Self {
            strategy,
            record_buf,
            buffer: vec![],
            last_rid: None,
        }
    }

    fn consume_record(&mut self, record: bcf::Record) -> Option<ClusterOrRecords> {
        let rid = record.rid();
        let record_belongs_to_current_buffer = self
            .strategy
            .belongs(self.buffer.last(), &record)
            .unwrap_or(false);
        if rid != self.last_rid {
            self.last_rid = rid;
            let ret = self.flush();
            if record_belongs_to_current_buffer {
                self.buffer.push(record);
            }
            ret
        } else if record_belongs_to_current_buffer {
            self.buffer.push(record);
            None
        } else {
            let ret = self.flush();
            self.buffer.push(record);
            ret
        }
    }

    fn flush(&mut self) -> Option<ClusterOrRecords> {
        if self.strategy.is_cluster(&self.buffer, |r| r, |_| None) {
            let buf = std::mem::take(&mut self.buffer);
            Some(ClusterOrRecords::Cluster(buf))
        } else if !self.buffer.is_empty() {
            let buf = std::mem::take(&mut self.buffer);
            Some(ClusterOrRecords::Records(buf))
        } else {
            None
        }
    }
}

pub enum ClusterOrRecords {
    Cluster(Vec<bcf::Record>),
    Records(Vec<bcf::Record>),
    SequenceDone,
}

impl std::fmt::Debug for ClusterOrRecords {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            Self::Records(records) => write!(f, "{} records", records.len()),
            Self::Cluster(records) => write!(f, "cluster of {}", records.len()),
            Self::SequenceDone => write!(f, "sequence done"),
        }
    }
}

pub(super) fn find_clusters(
    strategy: ClusteringSettings,
    mut reader: VCFReader,
    mut targets: Option<Targets>,
) -> impl Stream<Item = ClusterOrRecords> {
    let record_buf = reader.empty_record();

    let mut state = State::new(strategy, record_buf);

    stream! {
        loop {
            let read =
                tokio::task::block_in_place(|| reader.read(&mut state.record_buf));

            let event = match read {
                Some(Ok(())) => {
                    let reg = GenomeRegion::try_from(&state.record_buf);
                    let targets_include = if let (Some(targets), Ok(reg)) = (targets.as_mut(), reg) {
                        match targets.query(&reg) {
                            TargetQueryAnswer::Contained => true,
                            TargetQueryAnswer::NotContained => false,
                            TargetQueryAnswer::NoMoreTargets => break,
                        }
                    } else {
                        true // default: include
                    };
                    if targets_include {
                        if !state.record_buf.has_filter(&bcf::header::Id(0)) {
                            counter!("records.read.vcf_filter").inc(1);
                            None // Hole because filter
                        } else if !GtAndPhase::from(&state.record_buf).has_alt()
                        {
                            // Genotyped record with no alt allele on either haplotype: not a
                            // variant for this sample, so it must not participate in clustering.
                            counter!("records.read.no_alt").inc(1);
                            None // Hole because the genotype carries no alt allele
                        } else {
                            // filter PASS
                            counter!("records.read.pass").inc(1);
                            state.consume_record(state.record_buf.clone())
                        }
                    } else {
                        counter!("records.read.target_filter").inc(1);
                        None // Hole because excluded by targets
                    }
                }
                None if state.buffer.is_empty() => {
                    if reader.done() {
                        // the iterator is finished
                        break;
                    }
                    // there might be more to get, e.g. the next region.
                    Some(ClusterOrRecords::SequenceDone)
                }
                None => state.flush(),
                Some(Err(e)) => {
                    counter!("records.read.error").inc(1);
                    error!("Error reading record: {e}");
                    None // we have no entry here, but continue reading
                }
            };

            // Ignore holes
            let Some(event) = event else {
                continue;
            };

            yield event;
        }
    }
}