1 /* str_finder.c -- Short Tandem Repeat finder.
2 Originally from Crumble (https://github.com/jkbonfield/crumble)
3
4 Copyright (C) 2015-2016, 2021 Genome Research Ltd.
5
6 Author: James Bonfield <jkb@sanger.ac.uk>
7
8 Permission is hereby granted, free of charge, to any person obtaining a copy
9 of this software and associated documentation files (the "Software"), to deal
10 in the Software without restriction, including without limitation the rights
11 to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12 copies of the Software, and to permit persons to whom the Software is
13 furnished to do so, subject to the following conditions:
14
15 The above copyright notice and this permission notice shall be included in
16 all copies or substantial portions of the Software.
17
18 THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20 FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
21 THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
23 FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
24 DEALINGS IN THE SOFTWARE. */
25
26 #include <stdio.h>
27 #include <stdlib.h>
28 #include <string.h>
29 #include <stdint.h>
30 #include <ctype.h>
31
32 #include "str_finder.h"
33 #include "utlist.h"
34
35 #define MAX(a,b) ((a)>(b)?(a):(b))
36 #define MIN(a,b) ((a)<(b)?(a):(b))
37
38 typedef unsigned char uc;
39
add_rep(rep_ele ** list,char * cons,int clen,int pos,int rlen,int lower_only,unsigned int w)40 static void add_rep(rep_ele **list, char *cons, int clen, int pos, int rlen,
41 int lower_only, unsigned int w) {
42 rep_ele *el, *tmp, *prev;
43 char *cp1, *cp2, *cp_end;
44 int i;
45
46 // Already handled this in previous overlap?
47 if (*list) {
48 tmp = DL_TAIL(*list);
49 if (tmp->start <= pos-rlen*2+1 && tmp->end >= pos)
50 return;
51 }
52
53 // Find current and last occurence of repeated word.
54
55 cp2 = &cons[pos+1];
56 // If unpadded, this is quicker: cp1 = &cons[pos+1-rlen];
57
58 for (cp1 = &cons[pos], i = 1; i < rlen; cp1--) // compensate for pads
59 if (*cp1 == '*')
60 continue;
61 else
62 i++;
63 while (*cp1 == '*')
64 cp1--;
65
66
67 // Scan ahead to see how much further it goes.
68 cp_end = &cons[clen];
69 while (cp2 < cp_end) {
70 if (*cp1 != *cp2)
71 break;
72
73 w<<=2;
74 w|=*cp2;
75 cp1++;
76 cp2++;
77 }
78
79 if (!(el = malloc(sizeof(*el))))
80 return;
81
82 el->end = pos + cp2-&cons[pos+1];
83 el->rep_len = rlen;
84 pos++;
85 while (rlen--) {
86 while (cons[--pos] == '*');
87 while (cons[--pos] == '*');
88 }
89 //pos++;
90 while (pos > 1 && cons[pos-1] == '*') pos--;
91 el->start = pos;
92
93 // Check it meets the lower-case only criteria
94 if (lower_only) {
95 int lc = 0;
96 for (i = el->start; i <= el->end; i++) {
97 if (islower(cons[i])) {
98 lc = 1;
99 break;
100 }
101 }
102
103 if (!lc) {
104 free(el);
105 return;
106 }
107 }
108
109 // Remove any older items on the list that are entirely contained within el
110 if (*list) {
111 tmp = DL_TAIL(*list);
112 do {
113 prev = tmp->prev;
114 if (tmp->end < el->start)
115 break;
116
117 if (tmp->start >= el->start) {
118 DL_DELETE(*list, tmp);
119 free(tmp);
120 }
121
122 if (tmp == DL_HEAD(*list))
123 break;
124 tmp = prev;
125 } while (*list);
126 }
127
128 DL_APPEND(*list, el);
129
130 return;
131 }
132
133 /*
134 * Finds repeated homopolymers up to 8-mers.
135 * Note this assumes cons is 0-3, so N of 4 may rarely give false hits.
136 *
137 * Returns a list of rep_ele structs holding the start,end tuples of repeats;
138 * NULL on failure.
139 */
find_STR(char * cons,int len,int lower_only)140 rep_ele *find_STR(char *cons, int len, int lower_only) {
141 int i, j;
142 uint32_t w = 0;
143 rep_ele *reps = NULL;
144
145 for (i = j = 0; i < len && j < 15; i++) {
146 if (cons[i] == '*') continue;
147
148 w <<= 2;
149 w |= cons[i];
150 //printf("%3d %c w=%08x\n", i, cons[i], w);
151 if (j>= 1 && (w&0x0003) == ((w>> 2)&0x0003))
152 add_rep(&reps, cons, len, i, 1, lower_only, w);
153 if (j>= 3 && (w&0x000f) == ((w>> 4)&0x000f))
154 add_rep(&reps, cons, len, i, 2, lower_only, w);
155 if (j>= 5 && (w&0x003f) == ((w>> 6)&0x003f))
156 add_rep(&reps, cons, len, i, 3, lower_only, w);
157 if (j>= 7 && (w&0x00ff) == ((w>> 8)&0x00ff))
158 add_rep(&reps, cons, len, i, 4, lower_only, w);
159 if (j>= 9 && (w&0x03ff) == ((w>>10)&0x03ff))
160 add_rep(&reps, cons, len, i, 5, lower_only, w);
161 if (j>=11 && (w&0x0fff) == ((w>>12)&0x0fff))
162 add_rep(&reps, cons, len, i, 6, lower_only, w);
163 if (j>=13 && (w&0x3fff) == ((w>>14)&0x3fff))
164 add_rep(&reps, cons, len, i, 7, lower_only, w);
165
166 j++;
167 }
168
169 for (; i < len; i++) {
170 if (cons[i] == '*') continue;
171
172 w <<= 2;
173 w |= cons[i];
174 //printf("%3d %c w=%08x\n", i, cons[i], w);
175 if ((w&0xffff) == ((w>>16)&0xffff))
176 add_rep(&reps, cons, len, i, 8, lower_only, w);
177 else if ((w&0x3fff) == ((w>>14)&0x3fff))
178 add_rep(&reps, cons, len, i, 7, lower_only, w);
179 else if ((w&0x0fff) == ((w>>12)&0x0fff))
180 add_rep(&reps, cons, len, i, 6, lower_only, w);
181 else if ((w&0x03ff) == ((w>>10)&0x03ff))
182 add_rep(&reps, cons, len, i, 5, lower_only, w);
183 else if ((w&0x00ff) == ((w>> 8)&0x00ff))
184 add_rep(&reps, cons, len, i, 4, lower_only, w);
185 else if ((w&0x003f) == ((w>> 6)&0x003f))
186 add_rep(&reps, cons, len, i, 3, lower_only, w);
187 else if ((w&0x000f) == ((w>> 4)&0x000f))
188 add_rep(&reps, cons, len, i, 2, lower_only, w);
189 else if ((w&0x0003) == ((w>> 2)&0x0003))
190 add_rep(&reps, cons, len, i, 1, lower_only, w);
191 }
192
193 return reps;
194 }
195
196 /* -----------------------------------------------------------------------------
197 * Computes repeat regions in the consensus and then provides a bit mask
198 * indicating the extend of the STRs.
199 *
200 * The purpose of this is to identify where a read needs to span the entire
201 * region in order to validate how many copies of a repeat word are present.
202 * This only really has a major impact when indels are involved.
203 *
204 * For example, given this multiple alignment:
205 *
206 * S1 GATCGGACGAGAG
207 * S2 GATCGGACGAGAGAGAGAGAGT
208 * S3 GATCGGACGAGAGAGAGAG**TCGGAC
209 * S4 GGACGAGAGAGAGAGAGTCGGAC
210 * S5 CGAGAGAGAGAG**TCGGAC
211 * S6 AGAGAGAGTCGGAC
212 *
213 * We have subseq of GAGAGAGAGAG** vs GAGAGAGAGAGAG. The first and last
214 * (S1 and S6) sequences do not span and so we do not know which allele they
215 * match. Specifically as the pad is at the right hand end, the alignment of
216 * S6 gives incorrect weight to the consensus as it is stating AG when it
217 * may actually be ** at that point.
218 *
219 * By identifying the repeats we can soft clip as follows:
220 *
221 * S1 GATCGGACgagag
222 * S2 GATCGGACGAGAGAGAGAGAGT
223 * S3 GATCGGACGAGAGAGAGAG**TCGGAC
224 * S4 GGACGAGAGAGAGAGAGTCGGAC
225 * S5 CGAGAGAGAGAG**TCGGAC
226 * S6 agagagagTCGGAC
227 *
228 * Returns an array of STR vs no-STR values.
229 * 0 => non repetitive.
230 * 1+ => repeat with consecutive bit-number for repeat size.
231 *
232 * Eg: AGGGGAGGAGAAGAC
233 * 1111 1111
234 * 2222222
235 * 444444
236 * => 011331137754440
237 */
cons_mark_STR(char * cons,int len,int lower_only)238 char *cons_mark_STR(char *cons, int len, int lower_only) {
239 rep_ele *reps, *elt, *tmp;
240 char *str;
241
242 str = calloc(1, len);
243 reps = find_STR(cons, len, lower_only);
244
245 DL_FOREACH_SAFE(reps, elt, tmp) {
246 int i, v = 0;
247
248 //printf("%2d .. %2d %.*s\n", elt->start, elt->end,
249 // elt->end - elt->start+1, &cons[elt->start]);
250
251 // What is there?
252 for (i = MAX(elt->start-1,0); i <= MIN(elt->end+1,len-1); i++)
253 v |= str[i];
254
255 for (i = 0; i < 8; i++) {
256 if (!(v&(1<<i)))
257 break;
258 }
259 v = (i == 8) ? 1 : (1<<i);
260
261 // Add new if available, or just overload 1 if not
262 for (i = elt->start; i <= elt->end; i++)
263 str[i] |= v;
264
265 DL_DELETE(reps, elt);
266 free(elt);
267 }
268
269 return str;
270 }
271