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
use crate::alignment::{PileupAlignment, CIGAR_STATE_UNINIT};
use log::error;
use rust_htslib::bam::record::Cigar;
#[inline(always)]
/// This is a port of htslib's cigar_resolver2 from sam.c. I didn't try to reinvent the wheel by
/// implementing it in highly-idiomatic Rust; the algorithm is delicate enough as is. Shouldn't be
/// changed unless you have a very good reason.
pub fn resolve_cigar(plp: &mut PileupAlignment, pos: i64) {
let cs = &mut plp.cstate;
let cig = &cs.cig;
let ncig = cig.len();
let mut k: usize = 0;
let mut op: Cigar;
//////////////////////////////////////////////////////////////
// PHASE 1: SEEK TO CIGAR OPERATOR CONTAINING QUERY COORDINATE
//////////////////////////////////////////////////////////////
// never processed
if cs.icig == CIGAR_STATE_UNINIT {
plp.qpos = 0;
if ncig == 1 {
match cig[0] {
Cigar::Match(_) | Cigar::Equal(_) | Cigar::Diff(_) => {
cs.icig = 0;
cs.bam_pos = plp.rec.pos();
cs.iseq = 0;
}
_ => (),
}
} else {
cs.icig = 0;
cs.bam_pos = plp.rec.pos();
cs.iseq = 0;
for idx in 0..ncig {
k = idx;
match cig[k] {
Cigar::Match(_) | Cigar::Del(_) | Cigar::RefSkip(_) | Cigar::Equal(_) | Cigar::Diff(_) => break,
Cigar::Ins(l) | Cigar::SoftClip(l) => {
cs.iseq += l as usize;
}
// pad and hardclip, they don't consume anything, so we move on.
_ => (),
}
}
assert!(k < ncig);
cs.icig = k;
}
} else {
op = cig[cs.icig];
// the position we want is not in this operator, but one downstream.
if pos - cs.bam_pos >= i64::from(op.len()) {
assert!(cs.icig < ncig);
// Now we peek the next operator after this one. If it is reference consuming, we jump
// to it.
let op2 = cig[cs.icig + 1];
match op2 {
Cigar::Match(_) | Cigar::Del(_) | Cigar::RefSkip(_) | Cigar::Equal(_) | Cigar::Diff(_) => {
// if the old cigar we choose to move past is a read or reference consuming,
// update the indexes accordingly.
let op = cig[cs.icig];
match op {
Cigar::Match(l) | Cigar::Equal(l) | Cigar::Diff(l) => {
cs.iseq += l as usize;
}
_ => (),
};
cs.bam_pos += op.len() as i64;
cs.icig += 1;
}
// the next operator is not reference consuming, so we move on until we find one.
_ => {
// 1. update state by moving past current cigar operator.
let op = cig[cs.icig];
match op {
Cigar::Match(l) | Cigar::Equal(l) | Cigar::Diff(l) => {
cs.iseq += l as usize;
}
_ => (),
}
cs.bam_pos += op.len() as i64;
// search all following cigar strings for a ref consumer
for idx in (cs.icig + 1)..ncig {
k = idx;
let next_op = cig[k];
match next_op {
Cigar::Match(_) | Cigar::Del(_) | Cigar::RefSkip(_) | Cigar::Equal(_) | Cigar::Diff(_) => {
break
} // found it!
// didn't find it, but need to up query pos...
Cigar::Ins(l) | Cigar::SoftClip(l) => cs.iseq += l as usize,
_ => (),
}
}
cs.icig = k;
}
}
assert!(cs.icig < ncig);
}
}
////////////////////////////////////////////////////
// PHASE 2: EXTRACTING PILEUP INFORMATION FROM CIGAR
////////////////////////////////////////////////////
// At this stage, we have hit the cigar operator that contains our queried position. Now look
// at the alignment and extract needed info.
op = cig[cs.icig];
// TODO: RESET OPERATORS HERE
plp.indel = 0;
plp.refskip = false;
plp.del = false;
// our position is right at the edge of an operation, so peek the next one.
if (cs.bam_pos + op.len() as i64) as usize - 1 == pos as usize && cs.icig + 1 < ncig {
let op2 = cig[cs.icig + 1];
match op2 {
Cigar::Del(l) if !matches!(op, Cigar::Del(_)) => {
plp.indel = -(l as i32);
for idx in (cs.icig + 2)..ncig {
k = idx;
// check that this is correct
if matches!(cs.cig[k], Cigar::Del(_)) {
plp.indel -= cs.cig[k].len() as i32;
} else {
break;
}
}
}
Cigar::Ins(l) => {
plp.indel = l as i32;
for idx in (cs.icig + 2)..ncig {
k = idx;
match cig[k] {
Cigar::Ins(l) => {
plp.indel += l as i32;
}
Cigar::Pad(_) => (),
_ => break,
};
}
}
Cigar::Pad(_) if cs.icig + 2 < ncig => {
let mut next_op: Cigar;
let mut total_insertion_length = 0;
for idx in (cs.icig + 2)..ncig {
k = idx;
next_op = cig[k];
match next_op {
Cigar::Ins(l) => {
total_insertion_length += l;
}
Cigar::Pad(_) => (),
_ => break,
}
}
if total_insertion_length > 0 {
plp.indel = total_insertion_length as i32;
}
}
// CHECK THAT THIS IS CORRECT
_ => (),
}
}
match op {
Cigar::Match(_) | Cigar::Equal(_) | Cigar::Diff(_) => {
plp.qpos = cs.iseq + (pos - cs.bam_pos) as usize;
}
Cigar::Del(_) | Cigar::RefSkip(_) => {
plp.del = true;
plp.refskip = matches!(op, Cigar::RefSkip(_));
plp.qpos = cs.iseq;
}
_ => {
error!("Bad cigar resolution");
}
}
plp.cigar_index = plp.cstate.icig;
plp.head = pos == plp.rec.pos();
plp.tail = pos == plp.rec.pos() + plp.cstate.read_len_from_cigar - 1;
}