1 |
#include "GArgs.h" |
2 |
#include "GStr.h" |
3 |
#include "GHash.hh" |
4 |
#include "GList.hh" |
5 |
|
6 |
#define USAGE "Usage:\n\ |
7 |
fqtrim [-5 <5'adapter>] [-3 <3'adapter>] [-l <minlen>] [-C] [-D] [-Q] \\\n\ |
8 |
[-n <rename_prefix>] [-o <trimmed.fq>] [-r <discarded.lst>] <input.fq>\n\ |
9 |
\n\ |
10 |
Trim low quality bases at 3' end, optional adapter, filter for low complexity and\n\ |
11 |
optionally collapse duplicates in the input fastq data.\n\ |
12 |
\n\ |
13 |
Options:\n\ |
14 |
-n rename all the reads using the <prefix> followed by a read counter;\n\ |
15 |
if -C option was given, the suffix \"_x<N>\" is appended, with <N> being\n\ |
16 |
the read duplication count\n\ |
17 |
-o write the trimmed/filtered fastq into <trimmed.fq>(instead of stdout)\n\ |
18 |
-5 trim the given adapter or primer sequence at the 5' end of each read\n\ |
19 |
(e.g. -5 CGACAGGTTCAGAGTTCTACAGTCCGACGATC)\n\ |
20 |
-3 trim the given adapter sequence at the 3' end of each read\n\ |
21 |
(e.g. -3 TCGTATGCCGTCTTCTGCTTG)\n\ |
22 |
-l minimum \"clean\" length at the high quality 5'end that a read must\n\ |
23 |
have in order to pass the filter (default: 16)\n\ |
24 |
-r write a simple \"trash report\" file listing the discarded reads with a\n\ |
25 |
one letter code for the reason why they were trashed\n\ |
26 |
-C collapse duplicate reads and add a prefix with their count to the read\n\ |
27 |
name\n\ |
28 |
-D apply a low-complexity (dust) filter and discard any read that has at\n\ |
29 |
least half of it masked by this\n\ |
30 |
-Q quality values in the input data are interpreted as phred64, convert\n\ |
31 |
them to phred33\n\ |
32 |
" |
33 |
// example 3' adapter for miRNAs: TCGTATGCCGTCTTCTGCTTG |
34 |
FILE* f_out=NULL; //stdout if not provided |
35 |
FILE* f_in=NULL; //input fastq (stdin if not provided) |
36 |
FILE* freport=NULL; |
37 |
bool debug=false; |
38 |
bool doCollapse=false; |
39 |
bool doDust=false; |
40 |
bool trashReport=false; |
41 |
bool rawFormat=false; |
42 |
int min_read_len=16; |
43 |
int dust_cutoff=16; |
44 |
bool isfasta=false; |
45 |
bool phred64=false; |
46 |
GStr prefix; |
47 |
GStr adapter3; |
48 |
GStr adapter5; |
49 |
int a3len=0; |
50 |
int a5len=0; |
51 |
// adaptor matching metrics -- for extendMatch() function |
52 |
static const int a_m_score=2; //match score |
53 |
static const int a_mis_score=-3; //mismatch |
54 |
static const int a_dropoff_score=7; |
55 |
static const int a_min_score=7; |
56 |
|
57 |
// element in dhash: |
58 |
class FqDupRec { |
59 |
public: |
60 |
int count; //how many of these reads are the same |
61 |
int len; //length of qv |
62 |
char* firstname; //optional, only if we want to keep the original read names |
63 |
char* qv; |
64 |
FqDupRec(GStr* qstr=NULL, const char* rname=NULL) { |
65 |
len=0; |
66 |
qv=NULL; |
67 |
firstname=NULL; |
68 |
count=0; |
69 |
if (qstr!=NULL) { |
70 |
qv=Gstrdup(qstr->chars()); |
71 |
len=qstr->length(); |
72 |
count++; |
73 |
} |
74 |
if (rname!=NULL) firstname=Gstrdup(rname); |
75 |
} |
76 |
~FqDupRec() { |
77 |
GFREE(qv); |
78 |
GFREE(firstname); |
79 |
} |
80 |
void add(GStr& d) { //collapse another record into this one |
81 |
if (d.length()!=len) |
82 |
GError("Error at FqDupRec::add(): cannot collapse reads with different length!\n"); |
83 |
count++; |
84 |
for (int i=0;i<len;i++) |
85 |
qv[i]=(qv[i]+d[i])/2; |
86 |
} |
87 |
}; |
88 |
|
89 |
void openfw(FILE* &f, GArgs& args, char opt) { |
90 |
GStr s=args.getOpt(opt); |
91 |
if (!s.is_empty()) { |
92 |
if (s=='-') f=stdout; |
93 |
else { |
94 |
f=fopen(s,"w"); |
95 |
if (f==NULL) GError("Error creating file: %s\n", s.chars()); |
96 |
} |
97 |
} |
98 |
} |
99 |
|
100 |
#define FWCLOSE(fh) if (fh!=NULL && fh!=stdout) fclose(fh) |
101 |
#define FRCLOSE(fh) if (fh!=NULL && fh!=stdin) fclose(fh) |
102 |
|
103 |
GHash<FqDupRec> dhash; //hash to keep track of duplicates |
104 |
|
105 |
bool ntrim(GStr& rseq, GStr &rqv, int &l5, int &l3); //returns true if any trimming occured |
106 |
|
107 |
int dust(GStr& seq); |
108 |
bool trim_adapter3(GStr& seq, int &l5, int &l3); //returns true if any trimming occured |
109 |
bool trim_adapter5(GStr& seq, int &l5, int &l3); //returns true if any trimming occured |
110 |
|
111 |
void convertQ64(char* q, int len); |
112 |
void convertQ64(GStr& q); |
113 |
|
114 |
int main(int argc, char * const argv[]) { |
115 |
GArgs args(argc, argv, "QRDCl:d:3:5:n:r:o:"); |
116 |
int e; |
117 |
int icounter=0; //counter for input reads |
118 |
int outcounter=0; //counter for output reads |
119 |
if ((e=args.isError())>0) { |
120 |
GMessage("%s\nInvalid argument: %s\n", USAGE, argv[e]); |
121 |
exit(224); |
122 |
} |
123 |
debug=(args.getOpt('D')!=NULL); |
124 |
phred64=(args.getOpt('Q')!=NULL); |
125 |
doCollapse=(args.getOpt('C')!=NULL); |
126 |
doDust=(args.getOpt('D')!=NULL); |
127 |
rawFormat=(args.getOpt('R')!=NULL); |
128 |
if (rawFormat) { |
129 |
GError("Sorry, raw qseq format parsing is not implemented yet!\n"); |
130 |
} |
131 |
prefix=args.getOpt('n'); |
132 |
GStr s=args.getOpt('l'); |
133 |
if (!s.is_empty()) |
134 |
min_read_len=s.asInt(); |
135 |
s=args.getOpt('d'); |
136 |
if (!s.is_empty()) { |
137 |
dust_cutoff=s.asInt(); |
138 |
doDust=true; |
139 |
} |
140 |
|
141 |
if (args.getOpt('3')!=NULL) { |
142 |
adapter3=args.getOpt('3'); |
143 |
adapter3.upper(); |
144 |
a3len=adapter3.length(); |
145 |
} |
146 |
if (args.getOpt('5')!=NULL) { |
147 |
adapter5=args.getOpt('5'); |
148 |
adapter5.upper(); |
149 |
a5len=adapter5.length(); |
150 |
} |
151 |
trashReport= (args.getOpt('r')!=NULL); |
152 |
if (args.startNonOpt()==0) { |
153 |
GMessage(USAGE); |
154 |
exit(224); |
155 |
} |
156 |
|
157 |
openfw(f_out, args, 'o'); |
158 |
if (trashReport) |
159 |
openfw(freport, args, 'r'); |
160 |
char* infile=NULL; |
161 |
while ((infile=args.nextNonOpt())!=NULL) { |
162 |
GStr infname(infile); |
163 |
if (strcmp(infile,"-")==0) { |
164 |
f_in=stdin; infname="stdin"; } |
165 |
else { |
166 |
f_in=fopen(infile,"r"); |
167 |
if (f_in==NULL) |
168 |
GError("Cannot open input file %s!\n",infile); |
169 |
} |
170 |
GLineReader fq(f_in); |
171 |
char* l=NULL; |
172 |
while ((l=fq.getLine())!=NULL) { |
173 |
GStr rname; //current read name |
174 |
GStr rseq; //current read sequence |
175 |
GStr rqv; //current read quality values |
176 |
GStr s; |
177 |
if (rawFormat) { |
178 |
//TODO: implement qseq parsing here |
179 |
//if (raw type=N) then continue; //skip invalid/bad records |
180 |
|
181 |
} //raw qseq format |
182 |
else { // FASTQ or FASTA |
183 |
isfasta=(l[0]=='>'); |
184 |
if (!isfasta && l[0]!='@') GError("Error: fastq record marker not detected!\n"); |
185 |
s=l; |
186 |
rname=&(l[1]); |
187 |
icounter++; |
188 |
//GMessage("readname=%s\n",rname.chars()); |
189 |
for (int i=0;i<rname.length();i++) |
190 |
if (rname[i]<=' ') { rname.cut(i); break; } |
191 |
//now get the sequence |
192 |
if ((l=fq.getLine())==NULL) |
193 |
GError("Error: unexpected EOF after header for %s\n",rname.chars()); |
194 |
rseq=l; //this must be the DNA line |
195 |
if (isfasta) { |
196 |
while ((l=fq.getLine())!=NULL) { |
197 |
//fasta can have multiple sequence lines |
198 |
if (l[0]=='>') { |
199 |
fq.pushBack(); |
200 |
break; // |
201 |
} |
202 |
rseq+=l; |
203 |
} |
204 |
} //multi-line fasta file |
205 |
if (!isfasta) { |
206 |
if ((l=fq.getLine())==NULL) |
207 |
GError("Error: unexpected EOF after sequence for %s\n", rname.chars()); |
208 |
if (l[0]!='+') GError("Error: fastq qv header marker not detected!\n"); |
209 |
if ((l=fq.getLine())==NULL) |
210 |
GError("Error: unexpected EOF after qv header for %s\n", rname.chars()); |
211 |
rqv=l; |
212 |
if (rqv.length()!=rseq.length()) |
213 |
GError("Error: qv len != seq len for %s\n", rname.chars()); |
214 |
} |
215 |
} //<-- FASTA or FASTQ |
216 |
rseq.upper(); |
217 |
int l5=0; |
218 |
int l3=rseq.length()-1; |
219 |
if (l3-l5+1<min_read_len) { |
220 |
if (trashReport) { |
221 |
fprintf(freport, "%s\ts\t%s\n",rname.chars(), rseq.chars()); |
222 |
} |
223 |
continue; |
224 |
} |
225 |
if (ntrim(rseq, rqv, l5, l3)) { |
226 |
if (l3-l5+1<min_read_len) { |
227 |
if (trashReport) { |
228 |
fprintf(freport, "%s\tN\t%s\n", rname.chars(), rseq.chars()); |
229 |
} |
230 |
continue; //invalid read |
231 |
} |
232 |
//-- keep only the l5..l3 range |
233 |
rseq=rseq.substr(l5, l3-l5+1); |
234 |
if (!rqv.is_empty()) |
235 |
rqv=rqv.substr(l5, l3-l5+1); |
236 |
} |
237 |
|
238 |
if (a3len>0) { |
239 |
if (trim_adapter3(rseq, l5, l3)) { |
240 |
if (l3-l5+1<min_read_len) { |
241 |
if (trashReport) { |
242 |
fprintf(freport, "%s\tA3\t%s\n",rname.chars(), rseq.chars()); |
243 |
} |
244 |
continue; |
245 |
} |
246 |
//-- keep only the l5..l3 range |
247 |
rseq=rseq.substr(l5, l3-l5+1); |
248 |
if (!rqv.is_empty()) |
249 |
rqv=rqv.substr(l5, l3-l5+1); |
250 |
}//some adapter was trimmed |
251 |
} //adapter trimming |
252 |
if (a5len>0) { |
253 |
if (trim_adapter5(rseq, l5, l3)) { |
254 |
if (l3-l5+1<min_read_len) { |
255 |
if (trashReport) { |
256 |
fprintf(freport, "%s\tA5\t%s\n",rname.chars(), rseq.chars()); |
257 |
} |
258 |
continue; |
259 |
} |
260 |
//-- keep only the l5..l3 range |
261 |
rseq=rseq.substr(l5, l3-l5+1); |
262 |
if (!rqv.is_empty()) |
263 |
rqv=rqv.substr(l5, l3-l5+1); |
264 |
}//some adapter was trimmed |
265 |
} //adapter trimming |
266 |
if (doCollapse) { |
267 |
//keep read for later |
268 |
FqDupRec* dr=dhash.Find(rseq.chars()); |
269 |
if (dr==NULL) { //new entry |
270 |
//if (prefix.is_empty()) |
271 |
dhash.Add(rseq.chars(), |
272 |
new FqDupRec(&rqv, rname.chars())); |
273 |
//else dhash.Add(rseq.chars(), new FqDupRec(rqv.chars(),rqv.length())); |
274 |
} |
275 |
else |
276 |
dr->add(rqv); |
277 |
} //collapsing duplicates |
278 |
else { //not collapsing duplicates |
279 |
//do the dust filter now |
280 |
if (doDust) { |
281 |
int dustbases=dust(rseq); |
282 |
if (dustbases>(rseq.length()>>1)) { |
283 |
if (trashReport) { |
284 |
fprintf(freport, "%s\tD\t%s\n",rname.chars(),rseq.chars()); |
285 |
} |
286 |
continue; |
287 |
} |
288 |
} |
289 |
//print this record here |
290 |
outcounter++; |
291 |
if (isfasta) { |
292 |
if (prefix.is_empty()) |
293 |
fprintf(f_out, ">%s\n%s\n", rname.chars(), rseq.chars()); |
294 |
else |
295 |
fprintf(f_out, ">%s%08d\n%s\n", prefix.chars(), outcounter, |
296 |
rseq.chars()); |
297 |
} |
298 |
else { //fastq |
299 |
if (phred64) convertQ64(rqv); |
300 |
if (prefix.is_empty()) |
301 |
fprintf(f_out, "@%s\n%s\n+\n%s\n", rname.chars(), rseq.chars(),rqv.chars()); |
302 |
else |
303 |
fprintf(f_out, "@%s_%08d\n%s\n+\n%s\n", prefix.chars(), outcounter, |
304 |
rseq.chars(),rqv.chars() ); |
305 |
} |
306 |
} //not collapsing duplicates |
307 |
} //for each fastq record |
308 |
} //while each input file |
309 |
FRCLOSE(f_in); |
310 |
if (doCollapse) { |
311 |
outcounter=0; |
312 |
int maxdup_count=1; |
313 |
char* maxdup_seq=NULL; |
314 |
dhash.startIterate(); |
315 |
FqDupRec* qd=NULL; |
316 |
char* seq=NULL; |
317 |
while ((qd=dhash.NextData(seq))!=NULL) { |
318 |
GStr rseq(seq); |
319 |
//do the dusting here |
320 |
if (doDust) { |
321 |
int dustbases=dust(rseq); |
322 |
if (dustbases>(rseq.length()>>1)) { |
323 |
if (trashReport && qd->firstname!=NULL) { |
324 |
fprintf(freport, "%s_x%d\tD\t%s\n",qd->firstname, qd->count,seq); |
325 |
} |
326 |
continue; |
327 |
} |
328 |
} |
329 |
outcounter++; |
330 |
if (qd->count>maxdup_count) { |
331 |
maxdup_count=qd->count; |
332 |
maxdup_seq=seq; |
333 |
} |
334 |
if (isfasta) { |
335 |
if (prefix.is_empty()) { |
336 |
fprintf(f_out, "@%s_x%d\n%s\n", qd->firstname, qd->count, |
337 |
rseq.chars()); |
338 |
} |
339 |
else { //use custom read name |
340 |
fprintf(f_out, "@%s%08d_x%d\n%s\n", prefix.chars(), outcounter, |
341 |
qd->count, rseq.chars()); |
342 |
} |
343 |
} |
344 |
else { //fastq format |
345 |
if (phred64) convertQ64(qd->qv, qd->len); |
346 |
if (prefix.is_empty()) { |
347 |
fprintf(f_out, "@%s_x%d\n%s\n+\n%s\n", qd->firstname, qd->count, |
348 |
rseq.chars(), qd->qv); |
349 |
} |
350 |
else { //use custom read name |
351 |
fprintf(f_out, "@%s%08d_x%d\n%s\n+\n%s\n", prefix.chars(), outcounter, |
352 |
qd->count, rseq.chars(), qd->qv); |
353 |
} |
354 |
} |
355 |
}//for each element of dhash |
356 |
if (maxdup_count>1) { |
357 |
GMessage("Maximum read multiplicity: x %d (read: %s)\n",maxdup_count, maxdup_seq); |
358 |
} |
359 |
} //report collapsed dhash entries |
360 |
GMessage("Number of input reads: %9d\n", icounter); |
361 |
GMessage(" Output records: %9d\n", outcounter); |
362 |
if (trashReport) { |
363 |
FWCLOSE(freport); |
364 |
} |
365 |
|
366 |
FWCLOSE(f_out); |
367 |
} |
368 |
|
369 |
class NData { |
370 |
public: |
371 |
int NPos[1024]; //there should be no reads longer than 1K ? |
372 |
int NCount; |
373 |
int end5; |
374 |
int end3; |
375 |
int seqlen; |
376 |
double perc_N; //percentage of Ns in end5..end3 range only! |
377 |
const char* seq; |
378 |
bool valid; |
379 |
NData() { |
380 |
NCount=0; |
381 |
end5=0; |
382 |
end3=0; |
383 |
seq=NULL; |
384 |
perc_N=0; |
385 |
valid=true; |
386 |
} |
387 |
void init(GStr& rseq) { |
388 |
NCount=0; |
389 |
perc_N=0; |
390 |
seqlen=rseq.length(); |
391 |
seq=rseq.chars(); |
392 |
end5=1; |
393 |
end3=seqlen; |
394 |
for (int i=0;i<seqlen;i++) |
395 |
if (rseq[i]=='N') {// if (!ichrInStr(rseq[i], "ACGT") |
396 |
NPos[NCount]=i; |
397 |
NCount++; |
398 |
} |
399 |
perc_N=(NCount*100.0)/seqlen; |
400 |
valid=true; |
401 |
} |
402 |
void N_calc() { //only in the region after trimming |
403 |
int n=0; |
404 |
for (int i=end3-1;i<end5;i++) { |
405 |
if (seq[i]=='N') n++; |
406 |
} |
407 |
perc_N=(n*100.0)/(end5-end3+1); |
408 |
} |
409 |
}; |
410 |
|
411 |
static NData feat; |
412 |
int perc_lenN=12; // incremental distance from ends, in percentage of |
413 |
// sequence length, where N-trimming is done (default:12 %) (autolimited to 20) |
414 |
|
415 |
void N_analyze(int l5, int l3, int p5, int p3) { |
416 |
/* assumes feat was filled properly */ |
417 |
int old_dif, t5,t3,v; |
418 |
if (l3-l5<min_read_len || p5>p3 ) { |
419 |
feat.end5=l5+1; |
420 |
feat.end3=l3+1; |
421 |
return; |
422 |
} |
423 |
|
424 |
t5=feat.NPos[p5]-l5; |
425 |
t3=l3-feat.NPos[p3]; |
426 |
old_dif=p3-p5; |
427 |
v=(int)((((double)(l3-l5))*perc_lenN)/100); |
428 |
if (v>20) v=20; /* enforce N-search limit for very long reads */ |
429 |
if (t5 < v ) { |
430 |
l5=feat.NPos[p5]+1; |
431 |
p5++; |
432 |
} |
433 |
if (t3 < v) { |
434 |
l3=feat.NPos[p3]-1; |
435 |
p3--; |
436 |
} |
437 |
/* restNs=p3-p5; number of Ns in the new CLR */ |
438 |
if (p3-p5==old_dif) { /* no change, return */ |
439 |
/* don't trim if this may shorten the read too much */ |
440 |
//if (l5-l3<min_read_len) return; |
441 |
feat.end5=l5+1; |
442 |
feat.end3=l3+1; |
443 |
return; |
444 |
} |
445 |
else |
446 |
N_analyze(l5,l3, p5,p3); |
447 |
} |
448 |
|
449 |
|
450 |
bool ntrim(GStr& rseq, GStr &rqv, int &l5, int &l3) { |
451 |
//count Ns in the sequence |
452 |
feat.init(rseq); |
453 |
l5=feat.end5-1; |
454 |
l3=feat.end3-1; |
455 |
N_analyze(feat.end5-1, feat.end3-1, 0, feat.NCount-1); |
456 |
if (l5==feat.end5-1 && l3==feat.end3-1) |
457 |
return false; //nothing changed |
458 |
l5=feat.end5-1; |
459 |
l3=feat.end3-1; |
460 |
if (l3-l5+1<min_read_len) { |
461 |
feat.valid=false; |
462 |
return true; |
463 |
} |
464 |
feat.N_calc(); |
465 |
if (feat.perc_N>6.2) { //not valid if more than 1 N per 16 bases |
466 |
feat.valid=false; |
467 |
l3=l5+1; |
468 |
return true; |
469 |
} |
470 |
return true; |
471 |
} |
472 |
|
473 |
//--------------- dust functions ---------------- |
474 |
class DNADuster { |
475 |
public: |
476 |
int dustword; |
477 |
int dustwindow; |
478 |
int dustwindow2; |
479 |
int dustcutoff; |
480 |
int mv, iv, jv; |
481 |
int counts[32*32*32]; |
482 |
int iis[32*32*32]; |
483 |
DNADuster(int cutoff=16, int winsize=32, int wordsize=3) { |
484 |
dustword=wordsize; |
485 |
dustwindow=winsize; |
486 |
dustwindow2 = (winsize>>1); |
487 |
dustcutoff=cutoff; |
488 |
mv=0; |
489 |
iv=0; |
490 |
jv=0; |
491 |
} |
492 |
void setWindowSize(int value) { |
493 |
dustwindow = value; |
494 |
dustwindow2 = (dustwindow >> 1); |
495 |
} |
496 |
void setWordSize(int value) { |
497 |
dustword=value; |
498 |
} |
499 |
void wo1(int len, const char* s, int ivv) { |
500 |
int i, ii, j, v, t, n, n1, sum; |
501 |
int js, nis; |
502 |
n = 32 * 32 * 32; |
503 |
n1 = n - 1; |
504 |
nis = 0; |
505 |
i = 0; |
506 |
ii = 0; |
507 |
sum = 0; |
508 |
v = 0; |
509 |
for (j=0; j < len; j++, s++) { |
510 |
ii <<= 5; |
511 |
if (*s<=32) { |
512 |
i=0; |
513 |
continue; |
514 |
} |
515 |
ii |= *s - 'A'; //assume uppercase! |
516 |
ii &= n1; |
517 |
i++; |
518 |
if (i >= dustword) { |
519 |
for (js=0; js < nis && iis[js] != ii; js++) ; |
520 |
if (js == nis) { |
521 |
iis[nis] = ii; |
522 |
counts[ii] = 0; |
523 |
nis++; |
524 |
} |
525 |
if ((t = counts[ii]) > 0) { |
526 |
sum += t; |
527 |
v = 10 * sum / j; |
528 |
if (mv < v) { |
529 |
mv = v; |
530 |
iv = ivv; |
531 |
jv = j; |
532 |
} |
533 |
} |
534 |
counts[ii]++; |
535 |
} |
536 |
} |
537 |
} |
538 |
|
539 |
int wo(int len, const char* s, int* beg, int* end) { |
540 |
int i, l1; |
541 |
l1 = len - dustword + 1; |
542 |
if (l1 < 0) { |
543 |
*beg = 0; |
544 |
*end = len - 1; |
545 |
return 0; |
546 |
} |
547 |
mv = 0; |
548 |
iv = 0; |
549 |
jv = 0; |
550 |
for (i=0; i < l1; i++) { |
551 |
wo1(len-i, s+i, i); |
552 |
} |
553 |
*beg = iv; |
554 |
*end = iv + jv; |
555 |
return mv; |
556 |
} |
557 |
|
558 |
void dust(const char* seq, char* seqmsk, int seqlen, int cutoff=0) { //, maskFunc maskfn) { |
559 |
int i, j, l, a, b, v; |
560 |
if (cutoff==0) cutoff=dustcutoff; |
561 |
a=0;b=0; |
562 |
//GMessage("Dust cutoff=%d\n", cutoff); |
563 |
for (i=0; i < seqlen; i += dustwindow2) { |
564 |
l = (seqlen > i+dustwindow) ? dustwindow : seqlen-i; |
565 |
v = wo(l, seq+i, &a, &b); |
566 |
if (v > cutoff) { |
567 |
//for (j = a; j <= b && j < dustwindow2; j++) { |
568 |
for (j = a; j <= b; j++) { |
569 |
seqmsk[i+j]='N';//could be made lowercase instead |
570 |
} |
571 |
} |
572 |
} |
573 |
//return first; |
574 |
} |
575 |
}; |
576 |
|
577 |
static DNADuster duster; |
578 |
|
579 |
int dust(GStr& rseq) { |
580 |
char* seq=Gstrdup(rseq.chars()); |
581 |
duster.dust(rseq.chars(), seq, rseq.length(), dust_cutoff); |
582 |
//check the number of Ns: |
583 |
int ncount=0; |
584 |
for (int i=0;i<rseq.length();i++) { |
585 |
if (seq[i]=='N') ncount++; |
586 |
} |
587 |
GFREE(seq); |
588 |
return ncount; |
589 |
} |
590 |
|
591 |
// ------------------ adapter matching - triplet matching |
592 |
//look for matching triplets amongst the first 9 nucleotides of the 3' adaptor |
593 |
// or the last 9 nucleotides for the 5' adapter |
594 |
//when a triplet match is found, simply try to extend the alignment using a drop-off scheme |
595 |
// check minimum score and |
596 |
// for 3' adapter trimming: |
597 |
// require that the right end of the alignment for either the adaptor OR the read must be |
598 |
// < 3 distance from its right end |
599 |
// for 5' adapter trimming: |
600 |
// require that the left end of the alignment for either the adaptor OR the read must |
601 |
// be at coordinate 0 |
602 |
|
603 |
bool extendMatch(const char* a, int alen, int ai, |
604 |
const char* b, int blen, int bi, int mlen, int& l5, int& l3) { |
605 |
//so the alignment starts at ai in a, bi in b, with a perfect match of length mlen |
606 |
//GStr dbg(a); |
607 |
//GMessage(">> in %s\n\textending hit: %s at position %d\n", a, (dbg.substr(ai, mlen)).chars(), ai); |
608 |
int a_l=ai; //alignment coordinates on a |
609 |
int a_r=ai+mlen-1; |
610 |
int b_l=bi; //alignment coordinates on b |
611 |
int b_r=bi+mlen-1; |
612 |
int ai_maxscore=ai; |
613 |
int bi_maxscore=bi; |
614 |
int score=mlen*a_m_score; |
615 |
int maxscore=score; |
616 |
int mism5score=a_mis_score-2; //increase penalty for mismatches at 5' end |
617 |
//try to extend to the left first, if possible |
618 |
while (ai>0 && bi>0) { |
619 |
ai--; |
620 |
bi--; |
621 |
score+= (a[ai]==b[bi])? a_m_score : mism5score; |
622 |
if (score>maxscore) { |
623 |
ai_maxscore=ai; |
624 |
bi_maxscore=bi; |
625 |
maxscore=score; |
626 |
} |
627 |
else if (maxscore-score>a_dropoff_score) break; |
628 |
} |
629 |
a_l=ai_maxscore; |
630 |
b_l=bi_maxscore; |
631 |
//now extend to the right |
632 |
ai_maxscore=a_r; |
633 |
bi_maxscore=b_r; |
634 |
ai=a_r; |
635 |
bi=b_r; |
636 |
score=maxscore; |
637 |
//sometimes there are extra AAAAs at the end of the read, ignore those |
638 |
if (strcmp(&a[alen-4],"AAAA")==0) { |
639 |
alen-=3; |
640 |
while (a[alen-1]=='A' && alen>ai) alen--; |
641 |
} |
642 |
while (ai<alen-1 && bi<blen-1) { |
643 |
ai++; |
644 |
bi++; |
645 |
//score+= (a[ai]==b[bi])? a_m_score : a_mis_score; |
646 |
if (a[ai]==b[bi]) { //match |
647 |
score+=a_m_score; |
648 |
if (ai>=alen-2) { |
649 |
score+=a_m_score-(alen-ai-1); |
650 |
} |
651 |
} |
652 |
else { //mismatch |
653 |
score+=a_mis_score; |
654 |
} |
655 |
if (score>maxscore) { |
656 |
ai_maxscore=ai; |
657 |
bi_maxscore=bi; |
658 |
maxscore=score; |
659 |
} |
660 |
else if (maxscore-score>a_dropoff_score) break; |
661 |
} |
662 |
a_r=ai_maxscore; |
663 |
b_r=bi_maxscore; |
664 |
if (maxscore>=a_min_score && (alen-a_r<3 || blen-b_r<3 || a_l<2 || b_l<2)) { |
665 |
if (a_l<alen-a_r-1) { |
666 |
//adapter closer to the left end (typical for 5' adapter) |
667 |
l5=a_r+1; |
668 |
l3=alen-1; |
669 |
} |
670 |
else { |
671 |
//adapter matching at the right end (typical for 3' adapter) |
672 |
l5=0; |
673 |
l3=a_l-1; |
674 |
} |
675 |
return true; |
676 |
} |
677 |
return false; |
678 |
} |
679 |
|
680 |
bool trim_adapter3(GStr& seq, int&l5, int &l3) { |
681 |
int rlen=seq.length(); |
682 |
l5=0; |
683 |
l3=rlen-1; |
684 |
//first try a full match, we might get lucky |
685 |
int fi=-1; |
686 |
if ((fi=seq.index(adapter3))>=0) { |
687 |
if (fi<rlen-fi-a3len) {//match is closer to the right end |
688 |
l5=fi+a3len; |
689 |
l3=rlen-1; |
690 |
} |
691 |
else { |
692 |
l5=0; |
693 |
l3=fi-1; |
694 |
} |
695 |
return true; |
696 |
} |
697 |
//also, for fast detection of other adapter-only reads that start past |
698 |
// the beginning of the adapter sequence, try to see if the first a3len-4 |
699 |
// characters of the read are a substring of the adapter |
700 |
GStr rstart=seq.substr(0,a3len-4); |
701 |
if ((fi=adapter3.index(rstart))>=0) { |
702 |
l3=rlen-1; |
703 |
l5=a3len-4; |
704 |
while (fi+l5<a3len && l5<l3 && adapter3[fi+l5]==seq[l5]) l5++; |
705 |
return true; |
706 |
} |
707 |
|
708 |
//another easy case: first half of the adaptor matches |
709 |
int a3hlen=a3len>>1; |
710 |
GStr ahalf=adapter3.substr(0, a3hlen); |
711 |
if ((fi=seq.index(ahalf))>=0) { |
712 |
extendMatch(seq.chars(), rlen, fi, |
713 |
adapter3.chars(), a3len, 0, a3hlen, l5,l3); |
714 |
return true; |
715 |
} |
716 |
//no easy cases, so let's do the word hashing |
717 |
for (int iw=0;iw<6;iw++) { |
718 |
GStr aword=adapter3.substr(iw,3); |
719 |
if ((fi=seq.index(aword))>=0 && rlen-fi>3) { |
720 |
if (extendMatch(seq.chars(), rlen, fi, adapter3.chars(), |
721 |
a3len, iw, 3, l5,l3)) return true; |
722 |
} |
723 |
} |
724 |
return false; //no adapter parts found |
725 |
} |
726 |
|
727 |
bool trim_adapter5(GStr& seq, int&l5, int &l3) { |
728 |
int rlen=seq.length(); |
729 |
l5=0; |
730 |
l3=rlen-1; |
731 |
//try to see if adapter is fully included in the read |
732 |
int fi=-1; |
733 |
if ((fi=seq.index(adapter5))>=0) { |
734 |
if (fi<rlen-fi-a5len) {//match is closer to the right end |
735 |
l5=fi+a5len; |
736 |
l3=rlen-1; |
737 |
} |
738 |
else { |
739 |
l5=0; |
740 |
l3=fi-1; |
741 |
} |
742 |
return true; |
743 |
} |
744 |
//for fast detection of adapter-rich reads, check if the first 12 |
745 |
//characters of the read are a substring of the adapter |
746 |
GStr rstart=seq.substr(1,12); |
747 |
if ((fi=adapter5.index(rstart))>=0) { |
748 |
//l3=rlen-1; |
749 |
//l5=a5len-4; |
750 |
//while (fi+l5<a5len && l5<l3 && adapter5[fi+l5]==seq[l5]) l5++; |
751 |
//return true; |
752 |
extendMatch(seq.chars(), rlen, 1, |
753 |
adapter5.chars(), a5len, fi, 12, l5,l3); |
754 |
return true; |
755 |
} |
756 |
|
757 |
//another easy case: last 12 characters of the adaptor found as a substring of the read |
758 |
int aplen=12; |
759 |
int apstart=a5len-aplen-2; |
760 |
if (apstart<0) { apstart=0; aplen=a5len-1; } |
761 |
GStr bstr=adapter5.substr(apstart, aplen); |
762 |
if ((fi=seq.index(bstr))>=0) { |
763 |
extendMatch(seq.chars(), rlen, fi, |
764 |
adapter5.chars(), a5len, apstart, aplen, l5,l3); |
765 |
return true; |
766 |
} |
767 |
//no easy cases, find a triplet match as a seed for alignment extension |
768 |
//find triplets at the right end of the adapter |
769 |
for (int iw=0;iw<6;iw++) { |
770 |
apstart=a5len-iw-4; |
771 |
GStr aword=adapter5.substr(apstart,3); |
772 |
if ((fi=seq.index(aword))>=0) { |
773 |
if (extendMatch(seq.chars(), rlen, fi, adapter5.chars(), |
774 |
a5len, apstart, 3, l5,l3)) { |
775 |
return true; |
776 |
} |
777 |
} |
778 |
} |
779 |
return false; //no adapter parts found |
780 |
} |
781 |
|
782 |
//conversion of phred64 to phread33 |
783 |
void convertQ64(GStr& q) { |
784 |
for (int i=0;i<q.length();i++) q[i]-=31; |
785 |
} |
786 |
|
787 |
void convertQ64(char* q, int len) { |
788 |
for (int i=0;i<len;i++) q[i]-=31; |
789 |
} |