Skip to content

Commit fbabc7b

Browse files
author
Your Name
committed
made the binary version optional with the -usedump this will close #42 and #43 
1 parent 80b8e65 commit fbabc7b

4 files changed

Lines changed: 93 additions & 53 deletions

File tree

‎ngsLCA.cpp‎

Lines changed: 1 addition & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -752,9 +752,7 @@ int main_lca(int argc, char **argv) {
752752
int2int *i2i = NULL;
753753
// fprintf(stderr,"p->header: %p\n",p->header);
754754
if (p->header)
755-
i2i = bamRefId2tax(p->header, p->acc2taxfile, p->htsfile, errmap, p->tempfolder, p->reallyDump, p->filteredAcc2taxfile,NULL);
756-
757-
755+
i2i = bamRefId2tax(p->header, p->acc2taxfile, p->htsfile, errmap, p->tempfolder, p->useDump, p->filteredAcc2taxfile,NULL);
758756
// map of taxid -> taxid
759757
int2int parent;
760758
// map of taxid -> rank

‎ngsLCA_cli.cpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -44,7 +44,7 @@ pars *pars_init() {
4444
p->weighttype = 0;
4545
p->tempfolder = strdup("");
4646
p->ignore_errors = 0;
47-
p->reallyDump = 0;
47+
p->useDump = 0;
4848
p->maxreads = -1;
4949
return p;
5050
}
@@ -248,8 +248,8 @@ pars *get_pars(int argc, char **argv) {
248248
p->nthreads = atoi(val);
249249
else if (!strcasecmp("--weight_type", key))
250250
p->weighttype = atoi(val);
251-
else if (!strcasecmp("--reallyDump", key))
252-
p->reallyDump = atoi(val);
251+
else if (!strcasecmp("--useDump", key))
252+
p->useDump = atoi(val);
253253
else if (!strcasecmp("--ignore_errors", key)||!strcasecmp("-i", key))
254254
p->ignore_errors++;
255255
else if (!strcasecmp("--temp", key)) {

‎ngsLCA_cli.h‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -41,7 +41,7 @@ typedef struct {
4141
int weighttype;
4242
char *tempfolder;
4343
int ignore_errors;
44-
int reallyDump;
44+
int useDump;
4545
long maxreads;
4646
} pars;
4747

‎shared.cpp‎

Lines changed: 88 additions & 46 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,14 @@
1515
#include <map> // for operator!=, __map_iterator, operator==
1616
#include <vector> // for vector
1717

18+
19+
int file_is_older(const char *file1, const char *file2) {
20+
struct stat st1, st2;
21+
if (stat(file1, &st1) != 0) return -1;
22+
if (stat(file2, &st2) != 0) return -1;
23+
return st1.st_mtime < st2.st_mtime; // return 1 hvis file1 er ældre
24+
}
25+
1826
BGZF *getbgzf(const char *str1, const char *mode, int nthreads) {
1927
BGZF *fp = NULL;
2028
fp = bgzf_open(str1, mode);
@@ -44,6 +52,15 @@ BGZF *getbgzf3(const char *str1, const char *str2, const char *str3, const char
4452
return getbgzf(tmp, mode, nthreads);
4553
}
4654

55+
char *getfilename4(const char *str1, const char *str2, const char *str3, const char *str4) {
56+
unsigned tmp_l = strlen(str1) + strlen(str2) + strlen(str3) + strlen(str4) + 5;
57+
char *tmp=(char*)calloc(tmp_l,1);
58+
snprintf(tmp, tmp_l, "%s%s%s%s", str1, str2, str3, str4);
59+
return tmp;//getbgzf(tmp, mode, nthreads);
60+
}
61+
62+
63+
4764
BGZF *getbgzf4(const char *str1, const char *str2, const char *str3, const char *str4, const char *mode, int nthreads) {
4865
unsigned tmp_l = strlen(str1) + strlen(str2) + strlen(str3) + strlen(str4) + 5;
4966
char tmp[tmp_l];
@@ -229,102 +246,127 @@ void parse_nodes2(int2int &parent, int2intvec &child) {
229246
}
230247

231248
int SIG_COND = 1;
232-
int2int *bamRefId2tax(bam_hdr_t *hdr, char *acc2taxfile, char *bamfile, int2int &errmap, char *tempfolder, int forceDump, char *filteredAcc2taxfile,char2int *acc2taxidmap) {
233-
fprintf(stderr, "\t-> Starting to extract (acc->taxid) from binary file: \'%s\'\n", acc2taxfile);
249+
int2int *bamRefId2tax(bam_hdr_t *hdr, char *acc2taxfile, char *bamfile, int2int &errmap,
250+
char *tempfolder, int usedump, char *filteredAcc2taxfile,
251+
char2int *acc2taxidmap) {
252+
fprintf(stderr, "\t-> Starting to extract (acc->taxid) from binary file: '%s'\n", acc2taxfile);
234253
fflush(stderr);
235-
int dodump = !fexists4(tempfolder, basename(acc2taxfile), basename(bamfile), ".bin");
236-
if(acc2taxidmap!=NULL)
237-
forceDump = 1;
238-
dodump += forceDump;
239-
240-
fprintf(stderr, "\t-> Checking if need to reload acc2tax and dump. dodump=%d forcedump=%d acc2taxidmap: %p\n", dodump,forceDump,acc2taxidmap);
254+
255+
int binExists = fexists4(tempfolder, basename(acc2taxfile), basename(bamfile), ".bin");
256+
int useBinary = 0;
257+
258+
if (acc2taxidmap != NULL) {
259+
useBinary = 0;
260+
} else if (binExists && usedump) {
261+
useBinary = 1;
262+
} else if (!binExists && usedump) {
263+
useBinary = 0;
264+
} else if (binExists && !usedump) {
265+
useBinary = 1;
266+
} else {
267+
useBinary = 0;
268+
}
269+
270+
fprintf(stderr, "\t-> Binary exists: %d, usedump: %d, acc2taxidmap: %p => useBinary: %d\n",
271+
binExists, usedump, acc2taxidmap, useBinary);
241272

242273
time_t t = time(NULL);
243274
BGZF *fp = NULL;
244-
if (dodump)
245-
fp = getbgzf4(tempfolder, basename(acc2taxfile), basename(bamfile), ".bin", "wb", 4);
246-
else
275+
if (!useBinary && usedump)
276+
fp = getbgzf4(tempfolder, basename(acc2taxfile), basename(bamfile), ".bin", "wb", 4);
277+
else if (useBinary){
278+
char *binfile = getfilename4(tempfolder, basename(acc2taxfile), basename(bamfile), ".bin"); //not cleaned up
279+
280+
int binIsOlder = 0;
281+
if (binExists && usedump) {
282+
binIsOlder = file_is_older(binfile, bamfile);
283+
if (binIsOlder == 1) {
284+
fprintf(stderr, "\n\t-> WARNING: Binary cache file '%s' is older than BAM file '%s'.\n", binfile, bamfile);
285+
fprintf(stderr, "\t-> This may indicate out-of-date accession-taxid mapping. Refusing to use it.\n");
286+
exit(1);
287+
}
288+
}
247289
fp = getbgzf4(tempfolder, basename(acc2taxfile), basename(bamfile), ".bin", "rb", 4);
248-
249-
// This contains refname(as int) -> taxid
290+
}
250291
int2int *am = new int2int;
251292

252-
// Open the filtered output file if specified
253293
gzFile filteredFile = Z_NULL;
254-
if (filteredAcc2taxfile != NULL && dodump) {
294+
if (filteredAcc2taxfile != NULL && !useBinary) {
255295
filteredFile = gzopen(filteredAcc2taxfile, "wb");
256296
if (!filteredFile) {
257297
fprintf(stderr, "Error opening file '%s' for writing\n", filteredAcc2taxfile);
258298
exit(1);
259299
}
260-
gzprintf(filteredFile, "BAM_Reference\tTaxID\tBam_hdr_index\n"); // Write header
300+
gzprintf(filteredFile, "BAM_Reference\tTaxID\tBam_hdr_index\n");
261301
}
262302

263-
if (dodump) {
264-
char buf[4096];
303+
if (!useBinary) {
265304
int at = 0;
266-
267305
kstring_t *kstr = (kstring_t *)malloc(sizeof(kstring_t));
268306
kstr->l = kstr->m = 0;
269307
kstr->s = NULL;
270308
BGZF *fp2 = getbgzf(acc2taxfile, "rb", 2);
271-
bgzf_getline(fp2, '\n', kstr); // Skip header
309+
bgzf_getline(fp2, '\n', kstr); // skip header
272310
kstr->l = 0;
273-
int nprocs =0;
311+
int nprocs = 0;
274312
while (SIG_COND && bgzf_getline(fp2, '\n', kstr)) {
275-
if (kstr->l == 0)
276-
break;
313+
if (kstr->l == 0) break;
277314
if (!((at++ % 100000)) && isatty(fileno(stderr)))
278-
fprintf(stderr, "\r\t-> At linenr: %d in \'%s\' ", at, acc2taxfile);
315+
fprintf(stderr, "\r\t-> At linenr: %d in '%s' ", at, acc2taxfile);
316+
279317
char *tok = strtok(kstr->s, "\t\n ");
280318
char *key = strtok(NULL, "\t\n ");
281319
tok = strtok(NULL, "\t\n ");
282320
int val = atoi(tok);
283-
if(acc2taxidmap!=NULL){
284-
if(!acc2taxidmap->insert(std::pair<char*,int>(strdup(key),val)).second){
285-
fprintf(stderr,"\t-> Problem inserting key: %s with value: %d\n",key,val);
286-
exit(0);
287-
}
288-
}
289-
int valinbam = bam_name2id(hdr, key);
321+
if (!key) continue;
290322

291-
if (valinbam == -1)
292-
continue;
293-
nprocs++;
294-
if (fp != NULL) {
295-
assert(bgzf_write(fp, &valinbam, sizeof(int)) == sizeof(int));
296-
assert(bgzf_write(fp, &val, sizeof(int)) == sizeof(int));
323+
if (acc2taxidmap != NULL) {
324+
if (!acc2taxidmap->insert(std::pair<char*, int>(strdup(key), val)).second) {
325+
fprintf(stderr, "\t-> Problem inserting key: %s with value: %d\n", key, val);
326+
exit(1);
327+
}
297328
}
298329

330+
int valinbam = bam_name2id(hdr, key);
331+
if (valinbam == -1) continue;
332+
nprocs++;
333+
334+
if (fp != NULL)
335+
assert(bgzf_write(fp, &valinbam, sizeof(int)) == sizeof(int) &&
336+
bgzf_write(fp, &val, sizeof(int)) == sizeof(int));
337+
299338
if (am->find(valinbam) != am->end())
300-
fprintf(stderr, "\t-> Duplicate entries found \'%s\'\n", key);
339+
fprintf(stderr, "\t-> Duplicate entries found '%s'\n", key);
301340
(*am)[valinbam] = val;
302341

303-
// Write to the filtered file if specified
304342
if (filteredFile != Z_NULL)
305-
gzprintf(filteredFile, "%s\t%d\t%d\n", key, val,valinbam);
343+
gzprintf(filteredFile, "%s\t%d\t%d\n", key, val, valinbam);
306344
kstr->l = 0;
307345
}
308346
bgzf_close(fp2);
309-
fprintf(stderr,"\t-> Number of items from acc2tax file that is relevant from the bamheader: %d\n",nprocs);
310-
fflush(stderr);
347+
free(kstr->s);
348+
free(kstr);
349+
350+
fprintf(stderr, "\t-> Number of items from acc2tax file that is relevant from the bamheader: %d\n", nprocs);
351+
fflush(stderr);
311352
} else {
312353
int valinbam, val;
313354
while (bgzf_read(fp, &valinbam, sizeof(int))) {
314355
assert(bgzf_read(fp, &val, sizeof(int)) == sizeof(int));
315356
(*am)[valinbam] = val;
316-
317357
}
318358
}
319359

320-
// Close the filtered file if it was opened
321360
if (filteredFile != Z_NULL) {
322361
gzclose(filteredFile);
323362
fprintf(stderr, "\n\t-> Filtered acc2taxfile saved to '%s'\n", filteredAcc2taxfile);
324363
}
325364

326-
bgzf_close(fp);
327-
fprintf(stderr, "\t-> Number of entries to use from accession to taxid: %lu, time taken: %.2f sec acc2taxidmap.size(): %lu\n", am->size(), (float)(time(NULL) - t),acc2taxidmap!=NULL?acc2taxidmap->size():0);
365+
if (fp != NULL)
366+
bgzf_close(fp);
367+
368+
fprintf(stderr, "\t-> Number of entries to use from accession to taxid: %lu, time taken: %.2f sec acc2taxidmap.size(): %lu\n",
369+
am->size(), (float)(time(NULL) - t), acc2taxidmap != NULL ? acc2taxidmap->size() : 0);
328370
return am;
329371
}
330372

0 commit comments

Comments
 (0)