120 |
|
if (fline==NULL) return false; |
121 |
|
} |
122 |
|
//must be on '+' line here |
123 |
< |
if (fline==NULL || (reads_format == FASTQ && fline[0] != '+') || (reads_format == FASTA && quals && fline[0] != '>')) { |
123 |
> |
if (fline==NULL || (reads_format == FASTQ && fline[0] != '+') || |
124 |
> |
(reads_format == FASTA && quals && fline[0] != '>')) { |
125 |
|
err_exit("Error: '+' not found for fastq record %s\n",fline); |
126 |
|
return false; |
127 |
|
} |
132 |
|
while ((fline=fr.nextLine())!=NULL) { |
133 |
|
if (integer_quals) |
134 |
|
{ |
135 |
< |
vector<string> integer_qual_values; |
136 |
< |
tokenize(string(fline), " ", integer_qual_values); |
135 |
> |
vector<string> integer_qual_values; |
136 |
> |
tokenize(string(fline), " ", integer_qual_values); |
137 |
|
|
138 |
< |
string temp_qual; |
139 |
< |
for (size_t i = 0; i < integer_qual_values.size(); ++i) |
140 |
< |
{ |
141 |
< |
int qual_value = atoi(integer_qual_values[i].c_str()); |
142 |
< |
if (qual_value < 0) qual_value = 0; |
143 |
< |
temp_qual.push_back((char)(qual_value + 33)); |
144 |
< |
} |
138 |
> |
string temp_qual; |
139 |
> |
for (size_t i = 0; i < integer_qual_values.size(); ++i) |
140 |
> |
{ |
141 |
> |
int qual_value = atoi(integer_qual_values[i].c_str()); |
142 |
> |
if (qual_value < 0) qual_value = 0; |
143 |
> |
temp_qual.push_back((char)(qual_value + 33)); |
144 |
> |
} |
145 |
|
|
146 |
< |
qual.append(temp_qual); |
146 |
> |
qual.append(temp_qual); |
147 |
|
} |
148 |
|
else |
149 |
|
qual.append(fline); |
150 |
< |
|
150 |
< |
if ((!color && qual.length()>=seq.length()) || (color && qual.length()+1>=seq.length())) break; |
150 |
> |
if (qual.length()>=seq.length()-1) break; |
151 |
|
} |
152 |
|
// final check |
153 |
< |
if ((!color && seq.length()!=qual.length()) || (color && seq.length()!=qual.length()+1)) { |
154 |
< |
err_exit("Error: qual length (%d) differs from seq length (%d) for fastq record %s!\n", |
153 |
> |
if (color) { |
154 |
> |
if (seq.length()==qual.length()) { |
155 |
> |
//discard first qv |
156 |
> |
qual=qual.substr(1); |
157 |
> |
} |
158 |
> |
if (seq.length()!=qual.length()+1) { |
159 |
> |
err_exit("Error: length of quality string does not match seq length (%d) for color read %s!\n", |
160 |
> |
seq.length(), alt_name.c_str()); |
161 |
> |
} |
162 |
> |
} |
163 |
> |
else { |
164 |
> |
if (seq.length()!=qual.length()) { |
165 |
> |
err_exit("Error: qual string length (%d) differs from seq length (%d) for read %s!\n", |
166 |
|
qual.length(), seq.length(), alt_name.c_str()); |
167 |
< |
return false; |
167 |
> |
//return false; |
168 |
|
} |
169 |
+ |
} |
170 |
|
// |
159 |
– |
|
171 |
|
return !(qual.empty()); |
172 |
|
} |
173 |
|
|
174 |
|
bool next_fastx_read(FLineReader& fr, Read& read, ReadFormat reads_format, |
175 |
|
FLineReader* frq) { |
176 |
+ |
if (fr.pushed_read) |
177 |
+ |
{ |
178 |
+ |
read = fr.last_read; |
179 |
+ |
fr.pushed_read = false; |
180 |
+ |
return true; |
181 |
+ |
} |
182 |
+ |
|
183 |
|
read.clear(); |
184 |
|
char* buf=NULL; |
185 |
|
while ((buf=fr.nextLine())!=NULL) { |
237 |
|
read.qual.append(temp_qual); |
238 |
|
} |
239 |
|
else { |
240 |
< |
// if (color && read.qual.length()==0 && buf[0]=='!') |
223 |
< |
// read.qual.append(&(buf[1])); //some color qual strings start with '!' for the adaptor |
224 |
< |
// else |
225 |
< |
read.qual.append(buf); |
240 |
> |
read.qual.append(buf); |
241 |
|
} |
242 |
< |
if ((!color && read.qual.length()>=read.seq.length()) |
243 |
< |
|| (color && read.qual.length()+1>=read.seq.length())) break; |
242 |
> |
if (read.qual.length()>=read.seq.length()-1) |
243 |
> |
break; |
244 |
|
} //while qv lines |
245 |
|
|
246 |
|
// final check |
247 |
< |
if ((!color && read.seq.length()!=read.qual.length()) || (color && read.seq.length()!=read.qual.length()+1)) { |
247 |
> |
if (color) { |
248 |
> |
if (read.seq.length()==read.qual.length()) { |
249 |
> |
//discard first qv |
250 |
> |
read.qual=read.qual.substr(1); |
251 |
> |
} |
252 |
> |
if (read.seq.length()!=read.qual.length()+1) { |
253 |
> |
err_exit("Error: length of quality string does not match sequence length (%d) for color read %s!\n", |
254 |
> |
read.seq.length(), read.alt_name.c_str()); |
255 |
> |
} |
256 |
> |
} |
257 |
> |
else { |
258 |
> |
if (read.seq.length()!=read.qual.length()) { |
259 |
|
err_exit("Error: qual length (%d) differs from seq length (%d) for fastq record %s!\n", |
260 |
|
read.qual.length(), read.seq.length(), read.alt_name.c_str()); |
261 |
|
return false; |
262 |
|
} |
263 |
< |
// |
264 |
< |
return !(read.qual.empty()); |
263 |
> |
} |
264 |
> |
|
265 |
> |
fr.last_read = read; |
266 |
> |
return !(read.seq.empty()); |
267 |
|
} |
268 |
|
|
269 |
|
// This could be faster. |
593 |
|
} |
594 |
|
} |
595 |
|
|
596 |
+ |
|
597 |
+ |
bool ReadStream::next_read(Read& r, ReadFormat read_format) { |
598 |
+ |
FLineReader fr(fstream.file); |
599 |
+ |
while (read_pq.size()<500000 && !r_eof) { |
600 |
+ |
//keep the queue topped off |
601 |
+ |
Read rf; |
602 |
+ |
if (!next_fastx_read(fr, rf, read_format)) { |
603 |
+ |
r_eof=true; |
604 |
+ |
break; |
605 |
+ |
} |
606 |
+ |
//Read read=Read(rf); |
607 |
+ |
uint64_t id = (uint64_t)atol(rf.name.c_str()); |
608 |
+ |
read_pq.push(make_pair(id, rf)); |
609 |
+ |
} |
610 |
+ |
if (read_pq.size()==0) |
611 |
+ |
return false; |
612 |
+ |
const pair<uint64_t, Read>& t = read_pq.top(); |
613 |
+ |
r=t.second; //copy strings |
614 |
+ |
//free(t.second); |
615 |
+ |
read_pq.pop(); |
616 |
+ |
return true; |
617 |
+ |
} |
618 |
+ |
|
619 |
+ |
// reads must ALWAYS requested in increasing order of their ID |
620 |
+ |
bool ReadStream::getRead(uint64_t r_id, |
621 |
+ |
Read& read, |
622 |
+ |
ReadFormat read_format, |
623 |
+ |
bool strip_slash, |
624 |
+ |
uint64_t begin_id, |
625 |
+ |
uint64_t end_id, |
626 |
+ |
FILE* um_out, //unmapped reads output |
627 |
+ |
bool um_write_found //write the found ones |
628 |
+ |
) { |
629 |
+ |
if (!fstream.file) |
630 |
+ |
err_die("Error: calling ReadStream::getRead() with no file handle!"); |
631 |
+ |
if (r_id<last_id) |
632 |
+ |
err_die("Error: ReadStream::getRead() called with out-of-order id#!"); |
633 |
+ |
last_id=r_id; |
634 |
+ |
bool found=false; |
635 |
+ |
while (!found) { |
636 |
+ |
read.clear(); |
637 |
+ |
// Get the next read from the file |
638 |
+ |
if (!next_read(read, read_format)) |
639 |
+ |
break; |
640 |
+ |
|
641 |
+ |
if (strip_slash) { |
642 |
+ |
string::size_type slash = read.name.rfind("/"); |
643 |
+ |
if (slash != string::npos) |
644 |
+ |
read.name.resize(slash); |
645 |
+ |
} |
646 |
+ |
uint64_t id = (uint64_t)atol(read.name.c_str()); |
647 |
+ |
if (id >= end_id) |
648 |
+ |
return false; |
649 |
+ |
|
650 |
+ |
if (id < begin_id) |
651 |
+ |
continue; |
652 |
+ |
|
653 |
+ |
if (id == r_id) |
654 |
+ |
{ |
655 |
+ |
found=true; |
656 |
+ |
} |
657 |
+ |
else if (id > r_id) |
658 |
+ |
{ |
659 |
+ |
// this can't be right, read_pq should be maintained by next_read() |
660 |
+ |
//read_pq.push(make_pair(id, read)); |
661 |
+ |
fprintf(stderr, "WARNING: unexpected ID order encountered!(%lul > %lu)\n", |
662 |
+ |
id, r_id); |
663 |
+ |
break; |
664 |
+ |
} |
665 |
+ |
|
666 |
+ |
if (um_out && (um_write_found || !found)) { |
667 |
+ |
//write unmapped reads |
668 |
+ |
fprintf(um_out, "@%s\n%s\n+\n%s\n", read.alt_name.c_str(), |
669 |
+ |
read.seq.c_str(), read.qual.c_str()); |
670 |
+ |
} |
671 |
+ |
} //while reads |
672 |
+ |
return found; |
673 |
+ |
} |
674 |
+ |
|
675 |
+ |
|
676 |
|
bool get_read_from_stream(uint64_t insert_id, |
677 |
< |
FILE* reads_file, |
677 |
> |
FLineReader& fr, |
678 |
|
ReadFormat reads_format, |
679 |
|
bool strip_slash, |
680 |
< |
char read_name [], |
681 |
< |
char read_seq [], |
682 |
< |
char read_alt_name [], |
575 |
< |
char read_qual []) |
680 |
> |
Read& read, |
681 |
> |
FILE* um_out, |
682 |
> |
bool um_write_found) |
683 |
|
{ |
684 |
< |
Read read; |
685 |
< |
FLineReader fr(reads_file); |
686 |
< |
while(!fr.isEof()) |
687 |
< |
{ |
581 |
< |
read.clear(); |
582 |
< |
|
583 |
< |
// Get the next read from the file |
584 |
< |
if (!next_fastx_read(fr, read, reads_format)) |
585 |
< |
break; |
684 |
> |
bool found=false; |
685 |
> |
while(!found && !fr.isEof()) |
686 |
> |
{ |
687 |
> |
read.clear(); |
688 |
|
|
689 |
< |
if (strip_slash) |
690 |
< |
{ |
691 |
< |
string::size_type slash = read.name.rfind("/"); |
692 |
< |
if (slash != string::npos) |
693 |
< |
read.name.resize(slash); |
694 |
< |
} |
695 |
< |
|
696 |
< |
if ((uint64_t)atoi(read.name.c_str()) == insert_id) |
697 |
< |
{ |
698 |
< |
if (read_name) strcpy(read_name, read.name.c_str()); |
699 |
< |
if (read_seq) strcpy(read_seq, read.seq.c_str()); |
700 |
< |
if (read_alt_name) strcpy(read_alt_name, read.alt_name.c_str()); |
701 |
< |
if (read_qual) strcpy(read_qual, read.qual.c_str()); |
702 |
< |
return true; |
703 |
< |
} |
704 |
< |
|
705 |
< |
//rt.get_id(read.name, ref_str); |
706 |
< |
} |
707 |
< |
|
708 |
< |
return false; |
689 |
> |
// Get the next read from the file |
690 |
> |
if (!next_fastx_read(fr, read, reads_format)) |
691 |
> |
break; |
692 |
> |
if (strip_slash) |
693 |
> |
{ |
694 |
> |
string::size_type slash = read.name.rfind("/"); |
695 |
> |
if (slash != string::npos) |
696 |
> |
read.name.resize(slash); |
697 |
> |
} |
698 |
> |
uint64_t read_id = (uint64_t)atoi(read.name.c_str()); |
699 |
> |
if (read_id == insert_id) |
700 |
> |
{ |
701 |
> |
found=true; |
702 |
> |
} |
703 |
> |
else if (read_id > insert_id) |
704 |
> |
{ |
705 |
> |
fr.pushBack_read(); |
706 |
> |
break; |
707 |
> |
} |
708 |
> |
|
709 |
> |
if (um_out && (um_write_found || !found)) { |
710 |
> |
//write unmapped reads |
711 |
> |
fprintf(um_out, "@%s\n%s\n+\n%s\n", read.alt_name.c_str(), |
712 |
> |
read.seq.c_str(), read.qual.c_str()); |
713 |
> |
} |
714 |
> |
//rt.get_id(read.name, ref_str); |
715 |
> |
} //while reads |
716 |
> |
return found; |
717 |
|
} |