| 1 |
/* |
| 2 |
* prep_reads.cpp |
| 3 |
* TopHat |
| 4 |
* |
| 5 |
* Created by Cole Trapnell on 9/2/08. |
| 6 |
* Copyright 2008 Cole Trapnell. All rights reserved. |
| 7 |
* Derived from maq "catfilter", by Heng Li at Sanger |
| 8 |
*/ |
| 9 |
|
| 10 |
#ifdef HAVE_CONFIG_H |
| 11 |
#include <config.h> |
| 12 |
#endif |
| 13 |
|
| 14 |
#include <stdio.h> |
| 15 |
#include <cassert> |
| 16 |
#include <vector> |
| 17 |
#include <cstring> |
| 18 |
#include <cstdlib> |
| 19 |
|
| 20 |
#include "common.h" |
| 21 |
#include "reads.h" |
| 22 |
#include "tokenize.h" |
| 23 |
#include "qual.h" |
| 24 |
|
| 25 |
using namespace std; |
| 26 |
|
| 27 |
void format_qual_string(string& qual_str) |
| 28 |
{ |
| 29 |
for (size_t i = 0; i < qual_str.size(); ++i) |
| 30 |
{ |
| 31 |
qual_str[i] = charToPhred33(qual_str[i], |
| 32 |
solexa_quals, |
| 33 |
phred64_quals); |
| 34 |
} |
| 35 |
} |
| 36 |
|
| 37 |
vector<string> flt_reads_fnames; |
| 38 |
bool readmap_loaded=false; |
| 39 |
vector<bool> readmap; //for filter_multihits |
| 40 |
vector<bool> mate_readmap; //for filter_multihits |
| 41 |
|
| 42 |
void load_readmap(string& flt_fname, vector<bool>& rmap) { |
| 43 |
//readmap_loaded=false; |
| 44 |
if (flt_fname.empty()) return; |
| 45 |
ReadStream rdstream(flt_fname, NULL, true); |
| 46 |
Read read; |
| 47 |
while (rdstream.get_direct(read, reads_format)) { |
| 48 |
uint32_t rnum=(uint32_t)atol(read.name.c_str()); |
| 49 |
if (rnum>=(uint32_t) rmap.size()) |
| 50 |
rmap.resize(rnum+1, false); |
| 51 |
rmap[rnum] = true; |
| 52 |
} |
| 53 |
readmap_loaded=true; |
| 54 |
} |
| 55 |
|
| 56 |
bool check_readmap(vector<bool>& rmap, uint32_t read_id) { |
| 57 |
if (read_id>=rmap.size()) return false; |
| 58 |
else return rmap[read_id]; |
| 59 |
} |
| 60 |
|
| 61 |
void flt_reads_and_hits(vector<string>& reads_files) { |
| 62 |
if (!readmap_loaded) |
| 63 |
err_die("Error: filtering reads not enabled, aborting."); |
| 64 |
if (aux_outfile.empty()) |
| 65 |
err_die("Error: auxiliary output file not provided."); |
| 66 |
// -- filter mappings: |
| 67 |
string fext=getFext(flt_mappings); |
| 68 |
if (fext=="bam") { |
| 69 |
samfile_t* fbam=samopen(flt_mappings.c_str(), "rb", 0); |
| 70 |
if (fbam==NULL) |
| 71 |
err_die("Error opening BAM file %s!\n", flt_mappings.c_str()); |
| 72 |
bam1_t *b = bam_init1(); |
| 73 |
string aux_index(aux_outfile); |
| 74 |
aux_index+=".index"; |
| 75 |
GBamWriter wbam(aux_outfile.c_str(), fbam->header, aux_index); |
| 76 |
while (samread(fbam, b) > 0) { |
| 77 |
char* rname = bam1_qname(b); |
| 78 |
uint32_t rid=(uint32_t)atol(rname); |
| 79 |
if (check_readmap(readmap, rid)) { |
| 80 |
//write this mapping into the output file |
| 81 |
wbam.write(b, rid); |
| 82 |
} |
| 83 |
} |
| 84 |
bam_destroy1(b); |
| 85 |
samclose(fbam); |
| 86 |
} |
| 87 |
else { |
| 88 |
bool is_sam=false; |
| 89 |
string s(flt_mappings); |
| 90 |
for(size_t i=0; i!=s.length(); ++i) |
| 91 |
s[i] = std::tolower(s[i]); |
| 92 |
if (fext=="sam" || s.rfind(".sam.")+fext.length()+5==s.length()) is_sam=true; |
| 93 |
string unzcmd=getUnpackCmd(flt_mappings); |
| 94 |
FZPipe hitsfile(flt_mappings, unzcmd); |
| 95 |
FLineReader fh(hitsfile); |
| 96 |
FZPipe outfile; |
| 97 |
outfile.openWrite(aux_outfile.c_str(), zpacker); |
| 98 |
if (outfile.file==NULL) err_die("Error: cannot create file %s", aux_outfile.c_str()); |
| 99 |
const char* line; |
| 100 |
while ((line=fh.nextLine())) { |
| 101 |
if (is_sam && line[0]=='@') { |
| 102 |
//copy the header |
| 103 |
fprintf(outfile.file, "%s\n", line); |
| 104 |
continue; |
| 105 |
} |
| 106 |
char* tab=strchr((char*)line, '\t'); |
| 107 |
if (tab==NULL) |
| 108 |
err_die("Error: cannot find tab character in %s mappings line:\n%s", |
| 109 |
flt_mappings.c_str(),line); |
| 110 |
*tab=0; |
| 111 |
uint32_t rid = (uint32_t) atol(line); |
| 112 |
if (rid==0) |
| 113 |
err_die("Error: invalid read ID (%s) parsed from mapping file %s", |
| 114 |
line, flt_mappings.c_str()); |
| 115 |
*tab='\t'; |
| 116 |
if (check_readmap(readmap, rid)) { |
| 117 |
fprintf(outfile.file, "%s\n", line); |
| 118 |
} |
| 119 |
} |
| 120 |
outfile.close(); |
| 121 |
hitsfile.close(); |
| 122 |
} |
| 123 |
// -- now filter reads |
| 124 |
//FILE* findex = NULL; |
| 125 |
GBamWriter* wbam=NULL; |
| 126 |
FILE* fout=NULL; |
| 127 |
if (std_outfile.empty()) { |
| 128 |
fout=stdout; |
| 129 |
} |
| 130 |
else { //output file name explicitely given |
| 131 |
if (getFext(std_outfile)=="bam") { |
| 132 |
if (sam_header.empty()) err_die("Error: sam header file not provided.\n"); |
| 133 |
wbam = new GBamWriter(std_outfile.c_str(), sam_header.c_str(), index_outfile); |
| 134 |
//wbam = new GBamWriter(std_outfile, index_outfile); |
| 135 |
} |
| 136 |
else { |
| 137 |
fout = fopen(std_outfile.c_str(), "w"); |
| 138 |
if (fout==NULL) |
| 139 |
err_die("Error: cannot create file %s\n", std_outfile.c_str()); |
| 140 |
} |
| 141 |
} |
| 142 |
/* |
| 143 |
if (wbam==NULL && !index_outfile.empty()) { |
| 144 |
findex = fopen(index_outfile.c_str(), "w"); |
| 145 |
if (findex == NULL) |
| 146 |
err_die("Error: cannot create file %s\n", index_outfile.c_str()); |
| 147 |
} |
| 148 |
*/ |
| 149 |
for (size_t fi = 0; fi < reads_files.size(); ++fi) { |
| 150 |
//only one file expected here, this is not the initial prep_reads |
| 151 |
Read read; |
| 152 |
ReadStream readstream(reads_files[fi], NULL, true); |
| 153 |
//skip_lines(fr); |
| 154 |
while (readstream.get_direct(read)) { |
| 155 |
uint32_t rnum=(uint32_t)atol(read.name.c_str()); |
| 156 |
if (check_readmap(readmap, rnum)) { |
| 157 |
if (wbam) { |
| 158 |
GBamRecord bamrec(read.name.c_str(), -1, 0, false, |
| 159 |
read.seq.c_str(), NULL, read.qual.c_str()); |
| 160 |
wbam->write(bamrec.get_b(), rnum); |
| 161 |
} |
| 162 |
else { |
| 163 |
fprintf(fout, "@%s\n%s\n+%s\n%s\n", |
| 164 |
read.name.c_str(), |
| 165 |
read.seq.c_str(), |
| 166 |
read.alt_name.c_str(), |
| 167 |
read.qual.c_str()); |
| 168 |
} |
| 169 |
} |
| 170 |
} |
| 171 |
} |
| 172 |
if (wbam) delete wbam; |
| 173 |
if (fout && fout!=stdout) fclose(fout); |
| 174 |
} |
| 175 |
|
| 176 |
|
| 177 |
void writePrepBam(GBamWriter* wbam, Read& read, uint32_t rid, char trashcode=0, int matenum=0) { |
| 178 |
if (wbam==NULL) return; |
| 179 |
string rnum; |
| 180 |
str_appendUInt(rnum,rid); |
| 181 |
string rname(read.name); |
| 182 |
size_t pl=rname.length()-1; |
| 183 |
GBamRecord bamrec(rnum.c_str(), -1, 0, false, read.seq.c_str(), NULL, read.qual.c_str()); |
| 184 |
if (matenum) { |
| 185 |
if (matenum==1) bamrec.set_flag(BAM_FREAD1); |
| 186 |
else if (matenum==2) bamrec.set_flag(BAM_FREAD2); |
| 187 |
} |
| 188 |
/* |
| 189 |
int matenum=0; |
| 190 |
if (rname[pl-1]=='/') { |
| 191 |
if (rname[pl]=='1') { |
| 192 |
bamrec.set_flag(BAM_FREAD1); |
| 193 |
matenum=1; |
| 194 |
} |
| 195 |
else if (rname[pl]=='2') { |
| 196 |
bamrec.set_flag(BAM_FREAD2); |
| 197 |
matenum=2; |
| 198 |
} |
| 199 |
*/ |
| 200 |
if (matenum && rname[pl-1]=='/') |
| 201 |
rname.resize(pl-1); |
| 202 |
// } |
| 203 |
|
| 204 |
bamrec.add_aux("ZN", 'Z', rname.length(), (uint8_t*)rname.c_str()); |
| 205 |
|
| 206 |
if (trashcode) { |
| 207 |
bamrec.set_flag(BAM_FQCFAIL); |
| 208 |
bamrec.add_aux("ZT", 'A', 1, (uint8_t*)&trashcode); |
| 209 |
} |
| 210 |
wbam->write(bamrec.get_b(), rid); |
| 211 |
} |
| 212 |
|
| 213 |
bool processRead(int matenum, Read& read, uint32_t next_id, int& num_reads_chucked, |
| 214 |
int& multimap_chucked, GBamWriter* wbam, FILE* fout, FILE* fqindex, |
| 215 |
int& min_read_len, int& max_read_len, uint64_t& fout_offset, vector<bool>& rmap) { |
| 216 |
if (read.seq.length()<12) { |
| 217 |
++num_reads_chucked; |
| 218 |
writePrepBam(wbam, read, next_id, 'S', matenum); |
| 219 |
return false; |
| 220 |
} |
| 221 |
if ((int)read.seq.length()<min_read_len) |
| 222 |
min_read_len=read.seq.length(); |
| 223 |
if ((int)read.seq.length()>max_read_len) |
| 224 |
max_read_len=read.seq.length(); |
| 225 |
|
| 226 |
// daehwan - check this later, it's due to bowtie |
| 227 |
if (color && read.seq[1] == '4') { |
| 228 |
++num_reads_chucked; |
| 229 |
writePrepBam(wbam, read, next_id, 'c', matenum); |
| 230 |
return false; |
| 231 |
} |
| 232 |
|
| 233 |
if (readmap_loaded && check_readmap(rmap, next_id)) { |
| 234 |
++num_reads_chucked; |
| 235 |
++multimap_chucked; |
| 236 |
writePrepBam(wbam, read, next_id, 'M', matenum); |
| 237 |
return false; |
| 238 |
} |
| 239 |
format_qual_string(read.qual); |
| 240 |
std::transform(read.seq.begin(), read.seq.end(), read.seq.begin(), ::toupper); |
| 241 |
char counts[256]; |
| 242 |
memset(counts, 0, sizeof(counts)); |
| 243 |
// Count up the bad characters |
| 244 |
for (unsigned int i = 0; i != read.seq.length(); ++i) |
| 245 |
{ |
| 246 |
char c = (char)toupper(read.seq[i]); |
| 247 |
counts[(size_t)c]++; |
| 248 |
} |
| 249 |
|
| 250 |
double percent_A = (double)(counts[(size_t)'A']) / read.seq.length(); |
| 251 |
double percent_C = (double)(counts[(size_t)'C']) / read.seq.length(); |
| 252 |
double percent_G = (double)(counts[(size_t)'G']) / read.seq.length(); |
| 253 |
double percent_T = (double)(counts[(size_t)'T']) / read.seq.length(); |
| 254 |
double percent_N = (double)(counts[(size_t)'N']) / read.seq.length(); |
| 255 |
double percent_4 = (double)(counts[(size_t)'4']) / read.seq.length(); |
| 256 |
|
| 257 |
// Chuck the read if there are at least 5 'N's or if it's mostly |
| 258 |
// (>90%) 'N's and 'A's |
| 259 |
char trash_code=0; |
| 260 |
if (percent_A > 0.9 || |
| 261 |
percent_C > 0.9 || |
| 262 |
percent_G > 0.9 || |
| 263 |
percent_T > 0.9) |
| 264 |
trash_code='L'; |
| 265 |
else if (percent_N >= 0.1 || percent_4 >=0.1) |
| 266 |
trash_code='N'; |
| 267 |
if (trash_code) { |
| 268 |
++num_reads_chucked; |
| 269 |
writePrepBam(wbam, read, next_id, trash_code, matenum); |
| 270 |
return false; |
| 271 |
} |
| 272 |
|
| 273 |
if (wbam) { |
| 274 |
if (reads_format == FASTA) |
| 275 |
{ |
| 276 |
if (color) |
| 277 |
read.qual = string(read.seq.length()-1, 'I').c_str(); |
| 278 |
else |
| 279 |
read.qual = string(read.seq.length(), 'I').c_str(); |
| 280 |
} |
| 281 |
writePrepBam(wbam, read, next_id, 0, matenum); |
| 282 |
} |
| 283 |
else { |
| 284 |
// daehwan - we should not use buf in printf function |
| 285 |
// because it may contain some control characters such as "\" from quality values. |
| 286 |
// Here, buf is only used for calculating the file offset |
| 287 |
char buf[2048] = {0}; |
| 288 |
if (reads_format == FASTQ or (reads_format == FASTA && quals)) |
| 289 |
{ |
| 290 |
sprintf(buf, "@%u\n%s\n+%s\n%s\n", |
| 291 |
next_id, |
| 292 |
read.seq.c_str(), |
| 293 |
read.name.c_str(), |
| 294 |
read.qual.c_str()); |
| 295 |
|
| 296 |
fprintf(fout, "@%u\n%s\n+%s\n%s\n", |
| 297 |
next_id, |
| 298 |
read.seq.c_str(), |
| 299 |
read.name.c_str(), |
| 300 |
read.qual.c_str()); |
| 301 |
} |
| 302 |
else if (reads_format == FASTA) |
| 303 |
{ |
| 304 |
string qual; |
| 305 |
if (color) |
| 306 |
qual = string(read.seq.length()-1, 'I').c_str(); |
| 307 |
else |
| 308 |
qual = string(read.seq.length(), 'I').c_str(); |
| 309 |
|
| 310 |
sprintf(buf, "@%u\n%s\n+%s\n%s\n", |
| 311 |
next_id, |
| 312 |
read.seq.c_str(), |
| 313 |
read.name.c_str(), |
| 314 |
qual.c_str()); |
| 315 |
|
| 316 |
fprintf(fout, "@%u\n%s\n+%s\n%s\n", |
| 317 |
next_id, |
| 318 |
read.seq.c_str(), |
| 319 |
read.name.c_str(), |
| 320 |
qual.c_str()); |
| 321 |
} |
| 322 |
else |
| 323 |
{ |
| 324 |
assert(0); |
| 325 |
} |
| 326 |
|
| 327 |
if (fqindex != NULL) |
| 328 |
{ |
| 329 |
if ((next_id - num_reads_chucked) % INDEX_REC_COUNT == 0) |
| 330 |
fprintf(fqindex, "%d\t%lu\n", next_id, fout_offset); |
| 331 |
} |
| 332 |
|
| 333 |
fout_offset += strlen(buf); |
| 334 |
} |
| 335 |
return true; |
| 336 |
} //validate read |
| 337 |
|
| 338 |
|
| 339 |
const char* ERR_FILE_CREATE="Error: cannot create file %s\n"; |
| 340 |
|
| 341 |
void process_reads(vector<string>& reads_fnames, vector<FZPipe>& quals_files, |
| 342 |
vector<string>& mate_fnames, vector<FZPipe>& mate_quals_files) |
| 343 |
{ |
| 344 |
//TODO: add the option to write the garbage reads into separate file(s) |
| 345 |
int num_reads_chucked = 0; |
| 346 |
int multimap_chucked = 0; |
| 347 |
int mate_num_reads_chucked = 0; |
| 348 |
int mate_multimap_chucked = 0; |
| 349 |
int min_read_len = 20000000; |
| 350 |
int max_read_len = 0; |
| 351 |
int mate_min_read_len = 20000000; |
| 352 |
int mate_max_read_len = 0; |
| 353 |
uint32_t next_id = 0; |
| 354 |
FILE* fw=NULL; //aux output file |
| 355 |
string outfname; //std_outfile after instancing template |
| 356 |
string mate_outfname; |
| 357 |
string idxfname; //index_outfile after instancing template |
| 358 |
string mate_idxfname; |
| 359 |
bool have_mates = (mate_fnames.size() > 0); |
| 360 |
|
| 361 |
if (!aux_outfile.empty()) { |
| 362 |
fw=fopen(aux_outfile.c_str(), "w"); |
| 363 |
if (fw==NULL) |
| 364 |
err_die(ERR_FILE_CREATE,aux_outfile.c_str()); |
| 365 |
} |
| 366 |
|
| 367 |
FILE* fqindex = NULL; //fastq index |
| 368 |
FILE* mate_fqindex = NULL; |
| 369 |
GBamWriter* wbam=NULL; |
| 370 |
GBamWriter* mate_wbam=NULL; |
| 371 |
FILE* fout=NULL; |
| 372 |
FILE* mate_fout=NULL; |
| 373 |
uint64_t fout_offset = 0; |
| 374 |
uint64_t mate_fout_offset = 0; |
| 375 |
if (std_outfile.empty()) { |
| 376 |
fout=stdout; |
| 377 |
//for PE reads, flt_side will decide which side is printed (can't be both) |
| 378 |
if (have_mates && flt_side==2) |
| 379 |
err_die("Error: --flt-side option required for PE reads directed to stdout!\n"); |
| 380 |
} |
| 381 |
else { //output file name explicitely given |
| 382 |
//could be a template |
| 383 |
if (std_outfile.find("%side%") != string::npos) { |
| 384 |
outfname=str_replace(std_outfile, "%side%", "left"); |
| 385 |
if (have_mates) |
| 386 |
mate_outfname=str_replace(std_outfile, "%side%", "right"); |
| 387 |
} |
| 388 |
else { |
| 389 |
outfname=std_outfile; |
| 390 |
} |
| 391 |
if (index_outfile.find("%side%") != string::npos) { |
| 392 |
idxfname=str_replace(index_outfile, "%side%", "left"); |
| 393 |
if (have_mates) |
| 394 |
mate_idxfname=str_replace(index_outfile, "%side%", "right"); |
| 395 |
} |
| 396 |
else { |
| 397 |
idxfname=index_outfile; |
| 398 |
} |
| 399 |
if (getFext(outfname)=="bam") { |
| 400 |
if (sam_header.empty()) err_die("Error: sam header file not provided.\n"); |
| 401 |
wbam = new GBamWriter(outfname.c_str(), sam_header.c_str(), idxfname); |
| 402 |
if (!mate_outfname.empty()) { |
| 403 |
mate_wbam = new GBamWriter(mate_outfname.c_str(), sam_header.c_str(), mate_idxfname); |
| 404 |
} |
| 405 |
} |
| 406 |
else { //fastq output |
| 407 |
fout = fopen(outfname.c_str(), "w"); |
| 408 |
if (fout==NULL) |
| 409 |
err_die(ERR_FILE_CREATE, outfname.c_str()); |
| 410 |
mate_fout = fopen(mate_outfname.c_str(), "w"); |
| 411 |
if (mate_fout==NULL) |
| 412 |
err_die(ERR_FILE_CREATE, mate_outfname.c_str()); |
| 413 |
} |
| 414 |
} |
| 415 |
if (wbam==NULL && !idxfname.empty()) { |
| 416 |
//fastq file output, indexed |
| 417 |
fqindex = fopen(idxfname.c_str(), "w"); |
| 418 |
if (fqindex == NULL) |
| 419 |
err_die(ERR_FILE_CREATE, idxfname.c_str()); |
| 420 |
if (!mate_idxfname.empty()) { |
| 421 |
mate_fqindex = fopen(mate_idxfname.c_str(), "w"); |
| 422 |
if (mate_fqindex == NULL) |
| 423 |
err_die(ERR_FILE_CREATE, mate_idxfname.c_str()); |
| 424 |
} |
| 425 |
|
| 426 |
} |
| 427 |
for (size_t fi = 0; fi < reads_fnames.size(); ++fi) |
| 428 |
{ |
| 429 |
Read read; |
| 430 |
Read mate_read; |
| 431 |
FZPipe* fq=NULL; |
| 432 |
FZPipe* mate_fq=NULL; |
| 433 |
if (quals) |
| 434 |
fq = & quals_files[fi]; |
| 435 |
ReadStream reads(reads_fnames[fi], fq, true); |
| 436 |
ReadStream* mate_reads=NULL; |
| 437 |
if (have_mates) { |
| 438 |
if (quals) |
| 439 |
mate_fq = & mate_quals_files[fi]; |
| 440 |
mate_reads=new ReadStream(mate_fnames[fi], mate_fq, true); |
| 441 |
} |
| 442 |
while (reads.get_direct(read, reads_format)) { |
| 443 |
// Get the next read from the file |
| 444 |
//if (!next_fastx_read(fr, read, reads_format, ((quals) ? &frq : NULL))) |
| 445 |
// break; |
| 446 |
++next_id; |
| 447 |
//IMPORTANT: to keep paired reads in sync, this must be |
| 448 |
//incremented BEFORE any reads are chucked ! |
| 449 |
int matenum=0; |
| 450 |
if (have_mates) { |
| 451 |
if (!mate_reads->get_direct(mate_read, reads_format)) |
| 452 |
err_die("Error: could not retrieve mate for read '%s' !\n",read.name.c_str()); |
| 453 |
//check if reads are paired correctly |
| 454 |
matenum=1; |
| 455 |
size_t nl=read.name.length()-1; |
| 456 |
size_t mate_nl=mate_read.name.length()-1; |
| 457 |
bool mate_match=(nl==mate_nl); |
| 458 |
if (mate_match && read.name[nl-1]=='/' ) --nl; |
| 459 |
if (mate_match && strncmp(read.name.data(), mate_read.name.data(), nl+1)!=0) |
| 460 |
mate_match=false; |
| 461 |
if (!mate_match) { |
| 462 |
err_die("Error: read #%d (%s) does not have a matching mate in the same order (%s found instead)\n", |
| 463 |
next_id, read.name.c_str(), mate_read.name.c_str()); |
| 464 |
} |
| 465 |
} |
| 466 |
processRead(matenum, read, next_id, num_reads_chucked, multimap_chucked, |
| 467 |
wbam, fout, fqindex, min_read_len, max_read_len, fout_offset, readmap); |
| 468 |
if (mate_reads) { |
| 469 |
matenum=2; |
| 470 |
processRead(matenum, mate_read, next_id, mate_num_reads_chucked, mate_multimap_chucked, |
| 471 |
mate_wbam, mate_fout, mate_fqindex, mate_min_read_len, mate_max_read_len, |
| 472 |
mate_fout_offset, mate_readmap); |
| 473 |
} |
| 474 |
|
| 475 |
} //while !fr.isEof() |
| 476 |
} //for each input file |
| 477 |
fprintf(stderr, "%u out of %u reads have been filtered out\n", |
| 478 |
num_reads_chucked, next_id); |
| 479 |
if (readmap_loaded) |
| 480 |
fprintf(stderr, "\t(%u filtered out due to %s)\n", |
| 481 |
multimap_chucked, flt_reads_fnames[0].c_str()); |
| 482 |
if (have_mates) |
| 483 |
fprintf(stderr, "%u out of %u read mates have been filtered out\n", |
| 484 |
mate_num_reads_chucked, next_id); |
| 485 |
if (readmap_loaded && have_mates && mate_multimap_chucked) |
| 486 |
fprintf(stderr, "\t(%u mates filtered out due to %s)\n", |
| 487 |
mate_multimap_chucked, flt_reads_fnames[1].c_str()); |
| 488 |
|
| 489 |
if (wbam) { delete wbam; } |
| 490 |
if (mate_wbam) { delete mate_wbam; } |
| 491 |
if (fout && fout!=stdout) fclose(fout); |
| 492 |
if (mate_fout) fclose(mate_fout); |
| 493 |
if (fw!=NULL) { |
| 494 |
string side(""); |
| 495 |
if (have_mates) |
| 496 |
side="left_"; |
| 497 |
fprintf(fw, "%smin_read_len=%d\n", side.c_str(), min_read_len - (color ? 1 : 0)); |
| 498 |
fprintf(fw, "%smax_read_len=%d\n", side.c_str(), max_read_len - (color ? 1 : 0)); |
| 499 |
fprintf(fw, "%sreads_in =%d\n", side.c_str(), next_id); |
| 500 |
fprintf(fw, "%sreads_out=%d\n", side.c_str(), next_id-num_reads_chucked); |
| 501 |
if (have_mates) { |
| 502 |
side="right_"; |
| 503 |
fprintf(fw, "%smin_read_len=%d\n", side.c_str(), mate_min_read_len - (color ? 1 : 0)); |
| 504 |
fprintf(fw, "%smax_read_len=%d\n", side.c_str(), mate_max_read_len - (color ? 1 : 0)); |
| 505 |
fprintf(fw, "%sreads_in =%d\n", side.c_str(), next_id); |
| 506 |
fprintf(fw, "%sreads_out=%d\n", side.c_str(), next_id-mate_num_reads_chucked); |
| 507 |
} |
| 508 |
fclose(fw); |
| 509 |
} |
| 510 |
|
| 511 |
if (fqindex) fclose(fqindex); |
| 512 |
if (mate_fqindex) fclose(mate_fqindex); |
| 513 |
} |
| 514 |
|
| 515 |
void print_usage() |
| 516 |
{ |
| 517 |
fprintf(stderr, "Usage:\n prep_reads [--filter-multi <multi.fq>] <reads1.fa/fq,...,readsN.fa/fq> [<read_qual_files>,..] \\" |
| 518 |
"[<mate_reads1.fa/fq,...,mate_readsN.fa/fq> [<mate_qual_files>,..]\n"); |
| 519 |
} |
| 520 |
|
| 521 |
void open_qual_files(vector<FZPipe>& quals_files, string& quals_file_list) { |
| 522 |
vector<string> quals_file_names; |
| 523 |
tokenize(quals_file_list, ",", quals_file_names); |
| 524 |
for (size_t i = 0; i < quals_file_names.size(); ++i) |
| 525 |
{ |
| 526 |
FZPipe seg_file(quals_file_names[i], true); |
| 527 |
if (seg_file.file == NULL) |
| 528 |
{ |
| 529 |
fprintf(stderr, "Error: cannot open qual. file %s\n", |
| 530 |
quals_file_names[i].c_str()); |
| 531 |
exit(1); |
| 532 |
} |
| 533 |
quals_files.push_back(seg_file); |
| 534 |
} |
| 535 |
} |
| 536 |
|
| 537 |
int main(int argc, char *argv[]) |
| 538 |
{ |
| 539 |
fprintf(stderr, "prep_reads v%s (%s)\n", PACKAGE_VERSION, SVN_REVISION); |
| 540 |
fprintf(stderr, "---------------------------\n"); |
| 541 |
|
| 542 |
int parse_ret = parse_options(argc, argv, print_usage); |
| 543 |
if (parse_ret) |
| 544 |
return parse_ret; |
| 545 |
if(optind >= argc) |
| 546 |
{ |
| 547 |
print_usage(); |
| 548 |
return 1; |
| 549 |
} |
| 550 |
|
| 551 |
string reads_file_list(argv[optind++]); |
| 552 |
vector<string> reads_filenames; |
| 553 |
tokenize(reads_file_list, ",",reads_filenames); |
| 554 |
vector<FZPipe> quals_files; |
| 555 |
if (quals) |
| 556 |
{ |
| 557 |
if (optind>=argc) { |
| 558 |
err_die("Error: quality value file(s) not provided !\n"); |
| 559 |
} |
| 560 |
string quals_file_list = argv[optind++]; |
| 561 |
open_qual_files(quals_files, quals_file_list); |
| 562 |
if (quals_files.size()!=reads_filenames.size()) |
| 563 |
err_die("Error: number of quality value files must much the number of read files!\n"); |
| 564 |
} |
| 565 |
string mate_file_list; |
| 566 |
vector<string> mate_filenames; |
| 567 |
vector<FZPipe> mate_quals_files; |
| 568 |
if (optind<argc) |
| 569 |
mate_file_list = argv[optind++]; |
| 570 |
|
| 571 |
if (!mate_file_list.empty()) { |
| 572 |
tokenize(mate_file_list, ",",mate_filenames); |
| 573 |
if (reads_filenames.size()!=mate_filenames.size()) { |
| 574 |
err_die("Error: same number of input files required for paired reads!\n"); |
| 575 |
} |
| 576 |
if (quals) |
| 577 |
{ |
| 578 |
if (optind>=argc) { |
| 579 |
err_die("Error: mate quality value file(s) not provided !\n"); |
| 580 |
} |
| 581 |
string mate_quals_file_list = argv[optind++]; |
| 582 |
open_qual_files(mate_quals_files, mate_quals_file_list); |
| 583 |
if (mate_quals_files.size()!=mate_filenames.size()) |
| 584 |
err_die("Error: number of quality value files must much the number of read files!\n"); |
| 585 |
} |
| 586 |
} |
| 587 |
if (!flt_reads.empty()) { |
| 588 |
//for multi-mapped prefiltering usage |
| 589 |
readmap_loaded = false; |
| 590 |
tokenize(flt_reads, ",", flt_reads_fnames); |
| 591 |
load_readmap(flt_reads_fnames[0], readmap); |
| 592 |
if (flt_reads_fnames.size()==2) |
| 593 |
load_readmap(flt_reads_fnames[1], mate_readmap); |
| 594 |
} |
| 595 |
if (flt_mappings.empty()) |
| 596 |
process_reads(reads_filenames, quals_files, mate_filenames, mate_quals_files); |
| 597 |
else //special use case: filter previous mappings (when prefiltering) |
| 598 |
flt_reads_and_hits(reads_filenames); |
| 599 |
|
| 600 |
return 0; |
| 601 |
} |
