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
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
use std::collections::HashSet;
use anyhow::Context;
use ndarray::{Array2, Axis, stack};
use ordered_float::OrderedFloat;
use rust_htslib::{self, bam::HeaderView};
pub mod plp_from_records;
pub static BASE2IDX: [u8; 256] = {
let mut table = [0u8; 256];
table[b'G' as usize] = 1;
table[b'A' as usize] = 2;
table[b'T' as usize] = 3;
table[b'C' as usize] = 4;
table
};
#[derive(Debug)]
pub struct PlpInfo {
pub normed_count: Array2<f32>, // (4, step), 4: GATC . ratio
pub major: Vec<usize>,
pub minor: Vec<usize>,
}
impl PlpInfo {
// snp 5%, nh_indel: 10%, homo_indel: 45%
pub fn modify_ratio(
&mut self,
seq: &[u8],
snp_thr: f32,
nh_indel_thr: f32,
homo_indel_thr: f32,
) {
self.normed_count
.axis_iter_mut(Axis(1))
.enumerate()
.for_each(|(tt, mut gatc_ratio)| {
let maj = self.major[tt];
let mio = self.minor[tt];
/* 因为是 gap left alignment。所以如果是比对的 ins region，major向后移动一下
|
|
consensus: A C G - T
smc1 : A C G G T
*/
let maj = maj + mio.min(1);
if maj >= seq.len() {
return;
}
let base = seq[maj];
let base_idx = BASE2IDX[base as usize] as usize - 1;
let mut is_homo = false;
if maj > 0 {
is_homo |= seq[maj] == seq[maj - 1];
}
if (maj + 1) < seq.len() {
is_homo |= seq[maj] == seq[maj + 1];
}
let gap_ratio = 1.0 - gatc_ratio.sum();
let ins_region = self.minor[tt] > 0;
if !is_homo {
// non homo
if !ins_region {
// non-homo non ins region
/*
non homo insertion
|
consensus: A C G T
smc1 : A C G -
*/
if gap_ratio < nh_indel_thr {
gatc_ratio[base_idx] += gap_ratio;
}
for iter_idx in 0_usize..4 {
if iter_idx != base_idx {
if gatc_ratio[iter_idx] < snp_thr {
gatc_ratio[base_idx] += gatc_ratio[iter_idx];
gatc_ratio[iter_idx] = 0.0;
}
}
}
} else {
// non-homo ins retion
/*
non homo insertion ins retion
|
consensus: A C G - T
smc1 : A C G A T
*/
if gap_ratio > (1.0 - nh_indel_thr) {
gatc_ratio.mapv_inplace(|_| 0.0);
}
}
} else {
// homo
if !ins_region {
// homo non del region
/*
homo non insertion
|
consensus: A C G G
smc1 : A C - G
*/
if gap_ratio < homo_indel_thr {
gatc_ratio[base_idx] += gap_ratio;
}
for iter_idx in 0_usize..4 {
if iter_idx != base_idx {
if gatc_ratio[iter_idx] < snp_thr {
gatc_ratio[base_idx] += gatc_ratio[iter_idx];
gatc_ratio[iter_idx] = 0.0;
}
}
}
} else {
// homo ins retion
/*
homo insertion ins retion
|
|
consensus: A C G - T T
smc1 : A C G A T T
*/
if gap_ratio > (1.0 - homo_indel_thr) {
gatc_ratio.mapv_inplace(|_| 0.0);
}
}
}
});
}
// 如果当前位点的 ACGT 的比例加起来<ratio_thr, 那么扔掉该位点
pub fn drop_low_ratio_ins_locus(self, ratio_thr: f32) -> Self {
let low_ratio_locus = self
.normed_count
.axis_iter(Axis(1))
.enumerate()
.filter(|(_, ratios)| {
ratios
.iter()
.max_by_key(|v| OrderedFloat(**v))
.copied()
.unwrap()
< ratio_thr
})
.map(|(idx, _)| idx)
.collect::<HashSet<_>>();
if low_ratio_locus.is_empty() {
return self;
}
let normed_count = self
.normed_count
.axis_iter(Axis(1))
.enumerate()
.filter(|(idx, _)| !low_ratio_locus.contains(idx))
.map(|(_, arr)| arr)
.collect::<Vec<_>>();
let normed_count = stack(Axis(1), &normed_count).unwrap();
let major = self
.major
.into_iter()
.enumerate()
.filter(|(idx, _)| !low_ratio_locus.contains(idx))
.map(|(_, v)| v)
.collect::<Vec<_>>();
let minor = self
.minor
.into_iter()
.enumerate()
.filter(|(idx, _)| !low_ratio_locus.contains(idx))
.map(|(_, v)| v)
.collect::<Vec<_>>();
assert_eq!(major.len(), minor.len());
assert_eq!(normed_count.shape()[1], major.len());
Self {
normed_count,
major,
minor,
}
}
}
pub struct BamHeaderSeqInfo {
pub name: String,
pub length: usize,
}
pub fn extract_seq_info_from_header(
header_view: &HeaderView,
) -> anyhow::Result<Vec<BamHeaderSeqInfo>> {
let header = rust_htslib::bam::Header::from_template(header_view);
let header_hashmap = header.to_hashmap();
if !header_hashmap.contains_key("SQ") {
anyhow::bail!("invalid bam header. SQ not found");
}
let target_seq_infos = header_hashmap.get("SQ").unwrap();
// if target_seq_infos.len() > 1 {
// anyhow::bail!("more than one target seqs. not supported now");
// }
let mut results = vec![];
for seq_info in target_seq_infos {
if !seq_info.contains_key("LN") {
anyhow::bail!("invalid bam header. LN not found");
}
if !seq_info.contains_key("SN") {
anyhow::bail!("invalid bam header. SN not found");
}
let ln_str = seq_info.get("LN").unwrap();
let length = ln_str
.parse::<usize>()
.context(format!("parse {} to usize error", ln_str))?;
let sn = seq_info.get("SN").unwrap();
results.push(BamHeaderSeqInfo {
name: sn.to_string(),
length,
});
}
Ok(results)
}
#[cfg(test)]
mod test {
use crate::pileup_counter::PlpInfo;
#[test]
fn test_plp_info() {
// let mut plp_info = PlpInfo{normed_count: }
}
}