diff --git a/Makefile b/Makefile index 3d26d9a..0b1854a 100644 --- a/Makefile +++ b/Makefile @@ -128,6 +128,9 @@ misc: $(LIBHTS) $(OBJ) # --- Automatisk afhængighedshåndtering --- -include $(OBJ:.o=.d) +bfgs.o: bfgs.cpp $(LIBHTS) + $(CXX) -c $(CPPFLAGS) $(filter-out -Wall -Wextra,$(CXXFLAGS)) -w $< -o $@ + %.o: %.c $(LIBHTS) $(CC) -c $(CPPFLAGS) $(CFLAGS) $< -o $@ diff --git a/dfit_helppage.cpp b/dfit_helppage.cpp index 2331013..764e529 100644 --- a/dfit_helppage.cpp +++ b/dfit_helppage.cpp @@ -2,7 +2,6 @@ #include #include #include -#include #include #include #include diff --git a/dfit_optim.cpp b/dfit_optim.cpp index 71a3a71..2dbddc2 100644 --- a/dfit_optim.cpp +++ b/dfit_optim.cpp @@ -1,5 +1,4 @@ #include -#include #include #include #include @@ -12,7 +11,10 @@ double **read1_ugly_matrix(const char *fname){ fprintf(stderr,"\t-> Reading file: \'%s\'\n",fname); gzFile gz = Z_NULL; - assert((gz=gzopen(fname,"rb"))!=Z_NULL); + if((gz=gzopen(fname,"rb"))==Z_NULL){ + fprintf(stderr,"\t-> Problem opening file: %s will exit\n",fname); + exit(1); + } char buf[4096]; char *taxid = NULL; gzgets(gz,buf,4096); @@ -24,9 +26,12 @@ double **read1_ugly_matrix(const char *fname){ char *tok = strtok(buf,"\t\n "); if(taxid==NULL) taxid = strdup(tok); - else - assert(strcmp(taxid,tok)==0); - + else{ + if (strcmp(taxid, tok) != 0) { + fprintf(stderr, "\t-> Error: taxid mismatch (%s vs %s) will exit\n", taxid, tok); + exit(1); + } + } tok = strtok(NULL,"\t\n "); int is5 = 1; if(strcmp(tok,"3'")==0) diff --git a/main_aggregate_stat.cpp b/main_aggregate_stat.cpp index 27e61e7..78c8acb 100644 --- a/main_aggregate_stat.cpp +++ b/main_aggregate_stat.cpp @@ -13,7 +13,6 @@ #include #include #include // for std::sort -#include #include #include "profile.h" @@ -173,7 +172,11 @@ std::map read_dfit(char *fname){ BGZF *fpfpfp = NULL; fpfpfp = bgzf_open(fname, "rb"); - assert(fpfpfp!=NULL); + if(fpfpfp==NULL){ + fprintf(stderr,"\t-> Problem opening file: %s will exit\n",fname); + exit(1); + } + kstring_t *kstr = new kstring_t; kstr->s=NULL;kstr->l=kstr->m = 0; while(bgzf_getline(fpfpfp,'\n',kstr)){ @@ -200,25 +203,6 @@ std::map read_dfit(char *fname){ return ret; } -/* -void aggr_stat3000(std::map &stats,int2int &parent){ - std::map dasmap; - for(std::map::iterator it = stats.begin();it!=stats.end();it++) - dasmap[it->first] = it->second.nreads; - //fprintf(stderr,"dasmap.size(): %lu stats.size():%lu\n",dasmap.size(), stats.size()); - - for(std::map::iterator itt=dasmap.begin();itt!=dasmap.end();itt++){ - // for(int i=0;ifirst; - int2int::iterator it = parent.find(focal_taxid); - assert(it!=parent.end()); - int target = it->second; - to_root(focal_taxid,target,stats,parent,itt->second); - } - -} - -*/ void aggr_stat3000(std::map &stats,int2int &parent){ std::map dasmap; std::map datamap; @@ -234,7 +218,10 @@ void aggr_stat3000(std::map &stats,int2int &parent){ for(std::map::iterator itt=dasmap.begin();itt!=dasmap.end();itt++){ int focal_taxid = itt->first; int2int::iterator it = parent.find(focal_taxid); - assert(it!=parent.end()); + if (it == parent.end()) { + fprintf(stderr, "\t-> Error: iterator reached end (key:%d not found),will exit\n",focal_taxid); + exit(1); + } int target = it->second; if(target!=focal_taxid) @@ -276,10 +263,26 @@ int main_aggregate(int argc, char **argv) { fprintf(stderr,"\t-> --names file.txt.gz must be defined with --nodes is defined\n"); exit(1); } - fprintf(stderr,"aggregate infile_bdamage: %s infile_names: %s infile_nodes: %s infile_lcastat: %s infile_dfit: %s outfile_name: %s\n",infile_bdamage,infile_names,infile_nodes,infile_lcastat,infile_dfit,outfile_name); + fprintf(stderr, + "aggregate infile_bdamage: %s infile_names: %s infile_nodes: %s infile_lcastat: %s infile_dfit: %s outfile_name: %s\n", + infile_bdamage ? infile_bdamage : "NULL", + infile_names ? infile_names : "NULL", + infile_nodes ? infile_nodes : "NULL", + infile_lcastat ? infile_lcastat : "NULL", + infile_dfit ? infile_dfit : "NULL", + outfile_name ? outfile_name : "NULL" + ); if(outfile_name==NULL) outfile_name = strdup(infile_bdamage); - fprintf(stderr,"aggregate infile_bdamage: %s infile_names: %s infile_nodes: %s infile_lcastat: %s infile_dfit: %s outfile_name: %s\n",infile_bdamage,infile_names,infile_nodes,infile_lcastat,infile_dfit,outfile_name); + fprintf(stderr, + "aggregate infile_bdamage: %s infile_names: %s infile_nodes: %s infile_lcastat: %s infile_dfit: %s outfile_name: %s\n", + infile_bdamage ? infile_bdamage : "NULL", + infile_names ? infile_names : "NULL", + infile_nodes ? infile_nodes : "NULL", + infile_lcastat ? infile_lcastat : "NULL", + infile_dfit ? infile_dfit : "NULL", + outfile_name ? outfile_name : "NULL" + ); char buf[1024]; snprintf(buf, 1024, "%s.stat.gz", outfile_name); fprintf(stderr, "\t-> Dumping file: \'%s\'\n", buf); @@ -328,7 +331,10 @@ int main_aggregate(int argc, char **argv) { } if(dfit_int_char.size()>0){ std::map::iterator it = dfit_int_char.find(-1); - assert(it!=dfit_int_char.end()); + if (it == dfit_int_char.end()) { + fprintf(stderr, "\t-> Error: iterator reached end in dfit_int_char, will exit\n"); + exit(1); + } ksprintf(kstr,"\t%s",it->second); } ksprintf(kstr,"\n"); diff --git a/main_dfit.cpp b/main_dfit.cpp index c1a8103..3b605ac 100644 --- a/main_dfit.cpp +++ b/main_dfit.cpp @@ -79,6 +79,7 @@ void make_bootstrap_data(double **in,double **out,int howmany,mrand_t *rand_allo } void print_dat(double **dat, int howmany, int libprep) { + (void) libprep; int rows = 4; // [metadata, kvec, nvec, freq] int cols = 2 * howmany; printf("printing d:\n"); @@ -503,7 +504,13 @@ void slave_block(std::map &retmap,int howmany,sam_hdr_t *hdr,int2c } // Sigma and Z - double stats[2+2*(int)dat[0][0]]; + int n = 2 + 2 * (int)dat[0][0]; + + double *stats = (double *)malloc(n * sizeof(double)); + if (stats == NULL) { + fprintf(stderr, "\t-> Error: failed to allocate memory for stats, will exit\n"); + exit(1); + } getstat(dat,pars,stats); // stats contains standard deviation, then significance, then the calculated Dx for each position then the normalized @@ -626,7 +633,7 @@ void slave_block(std::map &retmap,int howmany,sam_hdr_t *hdr,int2c } ksprintf(kstr,"\n"); delete[] CI_val; // Free the array of pointers - + free(stats); } // Free the memory @@ -772,9 +779,9 @@ int main_dfit(int argc, char **argv) { snprintf(bootbuf, 1024, "%s.boot.stat.gz", outfile_name); bootfp = bgzf_open(bootbuf, "wb"); ksprintf(bootkstr,"taxid\tA_b\tq_b\tc_b\tphi_b\n"); - if(bgzf_write(bootfp,bootkstr->s,bootkstr->l)!=bootkstr->l){ + if(bgzf_write(bootfp,bootkstr->s,bootkstr->l)!=(ssize_t)bootkstr->l){ fprintf(stderr,"\t-> Problem writing boot file\n"); - exit(0); + exit(1); } bootkstr->l = 0; } @@ -884,7 +891,12 @@ int main_dfit(int argc, char **argv) { dings[i].bootkstr = bootkstr; dings[i].showfits = showfits; } - pthread_t mythreads[nthreads]; + + pthread_t *mythreads = (pthread_t *)malloc(nthreads * sizeof(pthread_t)); + if (mythreads == NULL) { + fprintf(stderr, "\t-> Error: failed to allocate memory for threads, will exit\n"); + exit(1); + } for(size_t i=0;is,kstr->l) != kstr->l){ + if(bgzf_write(fpfpfp,kstr->s,kstr->l) != (ssize_t) kstr->l){ fprintf(stderr, "\t-> Problems write to output BGZ file\n"); exit(1); } @@ -907,7 +919,7 @@ int main_dfit(int argc, char **argv) { if(doboot>0){ - if(bgzf_write(bootfp,bootkstr->s,bootkstr->l) != bootkstr->l){ + if(bgzf_write(bootfp,bootkstr->s,bootkstr->l) != (ssize_t) bootkstr->l){ fprintf(stderr, "\t-> Problemst write to output BGZ file\n"); exit(1); } diff --git a/main_pmd.cpp b/main_pmd.cpp index 7caef00..34ff71f 100644 --- a/main_pmd.cpp +++ b/main_pmd.cpp @@ -15,10 +15,6 @@ #include "profile.h" // for reconstructRefWithPosHTS -//extern htsFormat *dingding2; -//htsFormat *dingding3 = dingding2; - - int nproc = 0; // number of reads processed extern char refToChar[256]; @@ -34,7 +30,7 @@ double pmd_stat(char *seq, char *ref, int len, char *qs) { // make revcom static char compref1[512]; for (int i = 0; i < len; i++) - compref1[i] = revTable[ref[len - i - 1]]; + compref1[i] = revTable[(unsigned char)ref[len - i - 1]]; // do flipcount int count[2] = {0, 0}; @@ -45,17 +41,11 @@ double pmd_stat(char *seq, char *ref, int len, char *qs) { count[1] = count[1] + 1; } - char *myref = NULL; - if (count[0] > count[1]) - myref = compref1; - else - myref = ref; - // do stat double llh1 = 0; double llh2 = 0; for (int i = 0; i < len; i++) { - double eps = qsToProp[qs[i]]; + double eps = qsToProp[(unsigned char)qs[i]]; if (ref[i] == 'C') { double Dz = C0 + p * pow(1 - p, i + 1); double pm1 = (1 - ppi) * (1 - eps) * (1 - Dz) + ppi * eps / 3 * (1 - Dz) + ppi / 3 * (1 - eps) * Dz + (1 - ppi) * eps / 3 * Dz; @@ -83,7 +73,7 @@ double pmd_stat(char *seq, char *ref, int len, char *qs) { return llh1 - llh2; } -void wrapper(const bam1_t *b, const char *reconstructedReference, const std::vector &reconstructedReferencePos, const int &minQualBase, int MAXLENGTH, float **mm5p, float **mm3p, float incval, char myread[512], char myref[512]) { +void wrapper(const bam1_t *b, const char *reconstructedReference, const std::vector /*&reconstructedReferencePos*/, const int &minQualBase, int /*MAXLENGTH*/, float /***mm5p*/, float /***mm3p*/, float /*incval*/, char myread[512], char myref[512]) { const char *alphabetHTSLIB = "NACNGNNNTNNNNNNN"; char refeBase; char readBase; @@ -225,7 +215,7 @@ void parse_sequencingdata(char *refName, char *fname, int mapped_only, int se_on memset(myread, 'N', 512); memset(myrefe, 'N', 512); reconstructRefWithPosHTS(b, mypair, reconstructedRef); - wrapper(b, mypair.first->s, mypair.second, 0, 0, NULL, NULL, 1, myread, myrefe); + wrapper(b, mypair.first->s, mypair.second, 0, 0, 0, 0, 1, myread, myrefe); #if 0 fprintf(stderr,"\nmyread:\n%.*s\nmyReference:\n%.*s\n",b->core.l_qseq,myread,b->core.l_qseq,myrefe); fprintf(stderr,"---read[%d]----\n",nproc); @@ -237,7 +227,6 @@ void parse_sequencingdata(char *refName, char *fname, int mapped_only, int se_on } bam_destroy1(b); - free(dingding3); sam_hdr_destroy(hdr); sam_close(in); free(fname); @@ -359,6 +348,7 @@ int main_pmd(int argc, char **argv) { "\t[ALL done] cpu-time used = %.2f sec\n" "\t[ALL done] walltime used = %.2f sec\n", (float)(clock() - t) / CLOCKS_PER_SEC, (float)(time(NULL) - t2)); + (void)nthreads;//<- avoid compiler complain. return 0; } diff --git a/metaDMG.cpp b/metaDMG.cpp index 0c944d1..e6b960c 100644 --- a/metaDMG.cpp +++ b/metaDMG.cpp @@ -8,7 +8,6 @@ #include // for gzprintf, gzclose, gzgets, gzopen, Z_NULL #include -#include // for assert #include // for fprintf, NULL, stderr, stdout, fopen #include // for atoi, exit, free, atof, calloc, malloc #include // for strdup, strcmp, strtok, strlen, strncpy @@ -321,10 +320,13 @@ int main_getdamage(int argc, char **argv) { tmp2.push_back(mylen); seqlens[whichref] = tmp2; } else { - if(it->second.size()<1000000) + if(it->second.size()<1000000)//<- we dont need more than a million values do we? it->second.push_back(mygc); it = seqlens.find(whichref); - assert(it != seqlens.end()); + if (it == seqlens.end()) { + fprintf(stderr, "\t-> Error: iterator reached end in seqlens, will exit\n"); + exit(1); + } if(it->second.size()<1000000) it->second.push_back(mylen); } @@ -340,13 +342,22 @@ int main_getdamage(int argc, char **argv) { fprintf(stderr, "\t-> Outputting overall statistic in file: \"%s\"\n", buf); gzFile fpstat = NULL; - assert((fpstat = gzopen(buf, "wb")) != NULL); + if((fpstat = gzopen(buf, "wb")) == NULL){ + fprintf(stderr,"\t-> Error problem opening file %s will exit\n",buf); + exit(1); + } gzprintf(fpstat,"taxid\tnreads\tmean_len\tvar_len\tmean_gc\tvar_gc\tlca\trank\n"); for (std::map >::iterator it = gcconts.begin(); it != gcconts.end(); it++) { std::map::iterator it2 = dmg->assoc.find(it->first); - assert(it2 != dmg->assoc.end()); + if (it2 == dmg->assoc.end()) { + fprintf(stderr, "\t-> Error: iterator reached end in dmg->assoc, will exit\n"); + exit(1); + } std::map >::iterator it3 = seqlens.find(it->first); - assert(it3 != seqlens.end()); + if (it3 == seqlens.end()) { + fprintf(stderr, "\t-> Error: iterator reached end in seqlens (it3), will exit\n"); + exit(1); + } if (0) gzprintf(fpstat, "%d\t%lu\t%f\t%f\t%f\t%f\tNA\tNA\n", it->first, it2->second.nreads, mean(it3->second), var(it3->second), mean(it->second), var(it->second)); else{ @@ -407,9 +418,22 @@ int main_print(int argc, char **argv) { infile = strdup(*argv); } - fprintf(stderr, "infile: %s inbam: %s search: %d ctga: %d countout: %d nodes: %s names: %s howmany: %d\n", infile, inbam, search, ctga, countout, infile_nodes, infile_names, howmany); + fprintf(stderr, + "infile: %s inbam: %s search: %d ctga: %d countout: %d nodes: %s names: %s howmany: %d\n", + infile ? infile : "NULL", + inbam ? inbam : "NULL", + search, + ctga, + countout, + infile_nodes ? infile_nodes : "NULL", + infile_names ? infile_names : "NULL", + howmany + ); + if (!infile) { + fprintf(stderr, "\t-> Error: infile is NULL or could not be opened, will exit\n"); + exit(1); + } - assert(infile); int2char name_map; if (infile_names != NULL) name_map = parse_names(infile_names); @@ -460,7 +484,10 @@ int main_print(int argc, char **argv) { } int printlength; - assert(sizeof(int) == bgzf_read(bgfp, &printlength, sizeof(int))); + if (bgzf_read(bgfp, &printlength, sizeof(int)) != sizeof(int)) { + fprintf(stderr, "\t-> Error: failed to read expected number of bytes from bgzf file, will exit\n"); + exit(1); + } fprintf(stderr, "\t-> printlength(howmany) from inputfile: %d\n", printlength); int ref_nreads[2]; @@ -478,12 +505,23 @@ int main_print(int argc, char **argv) { float data[16]; while (1) { int nread = bgzf_read(bgfp, ref_nreads, 2 * sizeof(int)); - double ctgas[2 * printlength]; + double *ctgas = (double *)malloc(2 * printlength * sizeof(double)); + if (ctgas == NULL) { + fprintf(stderr, "\t-> Error: failed to allocate memory for ctgas, will exit\n"); + exit(1); + } if (nread == 0) break; - assert(nread == 2 * sizeof(int)); + if (nread != 2 * sizeof(int)) { + fprintf(stderr, "\t-> Error: unexpected number of bytes read (nread != 2*sizeof(int)), will exit\n"); + exit(1); + } for (int at = 0; at < printlength; at++) { - assert(16 * sizeof(float) == bgzf_read(bgfp, data, sizeof(float) * 16)); + if (bgzf_read(bgfp, data, sizeof(float) * 16) != 16 * sizeof(float)) { + fprintf(stderr, "\t-> Error: failed to read expected number of bytes (16 floats) from bgzf file, will exit\n"); + exit(1); + } + if (at >= howmany) continue; if (search == -1 || search == ref_nreads[0]) { @@ -538,7 +576,10 @@ int main_print(int argc, char **argv) { } for (int at = 0; at < printlength; at++) { - assert(16 * sizeof(float) == bgzf_read(bgfp, data, sizeof(float) * 16)); + if (bgzf_read(bgfp, data, sizeof(float) * 16) != 16 * sizeof(float)) { + fprintf(stderr, "\t-> Error: failed to read expected number of bytes (16 floats) from bgzf file, will exit\n"); + exit(1); + } if (at >= howmany) continue; if (search == -1 || search == ref_nreads[0]) { @@ -588,6 +629,7 @@ int main_print(int argc, char **argv) { fprintf(stdout, "\n"); } } + free(ctgas); } //cleanup for(int2char::iterator it=name_map.begin();it!=name_map.end();it++) @@ -647,8 +689,20 @@ int main_print2(int argc, char **argv) { infile = strdup(*argv); } - fprintf(stderr, "infile: %s inbam: %s names: %s search: %d ctga: %d countout: %d nodes: %s\n", infile, inbam, acc2tax, search, ctga, countout, infile_nodes); - assert(infile); + fprintf(stderr, + "infile: %s inbam: %s names: %s search: %d ctga: %d countout: %d nodes: %s\n", + infile ? infile : "NULL", + inbam ? inbam : "NULL", + acc2tax ? acc2tax : "NULL", + search, + ctga, + countout, + infile_nodes ? infile_nodes : "NULL" + ); + if (!infile) { + fprintf(stderr, "\t-> Error: infile is NULL or could not be opened, will exit\n"); + exit(1); + } int2char name_map; if (acc2tax != NULL) name_map = parse_names(acc2tax); @@ -699,7 +753,10 @@ int main_print2(int argc, char **argv) { } int printlength; - assert(sizeof(int) == bgzf_read(bgfp, &printlength, sizeof(int))); + if (bgzf_read(bgfp, &printlength, sizeof(int)) != sizeof(int)) { + fprintf(stderr, "\t-> Error: failed to read expected number of bytes for printlength, will exit\n"); + exit(1); + } fprintf(stderr, "\t-> printlength(howmany) from inputfile: %d\n", printlength); int ref_nreads[2]; @@ -726,13 +783,24 @@ int main_print2(int argc, char **argv) { while (1) { int nread = bgzf_read(bgfp, ref_nreads, 2 * sizeof(int)); - double ctgas[2 * printlength]; + double *ctgas = (double *)malloc(2 * printlength * sizeof(double)); + if (ctgas == NULL) { + fprintf(stderr, "\t-> Error: failed to allocate memory for ctgas, will exit\n"); + exit(1); + } + if (nread == 0) break; fprintf(stderr, "ref: %d nreads: %d\n", ref_nreads[0], ref_nreads[1]); - assert(nread == 2 * sizeof(int)); - for (int at = 0; at < printlength; at++) { - assert(16 * sizeof(float) == bgzf_read(bgfp, data, sizeof(float) * 16)); + if (nread != 2 * sizeof(int)) { + fprintf(stderr, "\t-> Error: unexpected number of bytes read (nread != 2*sizeof(int)), will exit\n"); + exit(1); + } + for (int at = 0; at < printlength; at++) { + if (bgzf_read(bgfp, data, sizeof(float) * 16) != 16 * sizeof(float)) { + fprintf(stderr, "\t-> Error: failed to read expected number of bytes (16 floats) from bgzf file, will exit\n"); + exit(1); + } if ((at + 1) > howmany) continue; if (search == -1 || search == ref_nreads[0]) { @@ -787,7 +855,10 @@ int main_print2(int argc, char **argv) { } for (int at = 0; at < printlength; at++) { - assert(16 * sizeof(float) == bgzf_read(bgfp, data, sizeof(float) * 16)); + if (bgzf_read(bgfp, data, sizeof(float) * 16) != 16 * sizeof(float)) { + fprintf(stderr, "\t-> Error: failed to read expected number of bytes (16 floats) from bgzf file, will exit\n"); + exit(1); + } if (at + 1 > howmany) continue; if (search == -1 || search == ref_nreads[0]) { @@ -838,6 +909,7 @@ int main_print2(int argc, char **argv) { fprintf(stdout, "\n"); } } + free(ctgas); } //clean up for(int2char::iterator it=name_map.begin();it!=name_map.end();it++) @@ -864,123 +936,29 @@ int main_print2(int argc, char **argv) { return 0; } -int main_merge(int argc, char **argv) { - - fprintf(stderr, "./metaDMG-cpp merge file.lca file.bdamage.gz [-names file.gz -bam file.bam -howmany 5 -nodes trestructure.gz]\n"); - if (argc <= 2) - return 0; - char *infile_lca = argv[1]; - char *infile_bdamage = argv[2]; - char *infile_nodes = NULL; - - char *inbam = NULL; - char *acc2tax = NULL; - int howmany = 5; - ++argv; - while (*(++argv)) { - if (strcasecmp("-names", *argv) == 0) - acc2tax = strdup(*(++argv)); - if (strcasecmp("-bam", *argv) == 0) - inbam = strdup(*(++argv)); - if (strcasecmp("-howmany", *argv) == 0) - howmany = atoi(*(++argv)); - if (strcasecmp("-nodes", *argv) == 0) - infile_nodes = strdup(*(++argv)); - } - - fprintf(stderr, "infile_lca: %s infile_bdamage: %s nodes: %s\n", infile_lca, infile_bdamage, infile_nodes); - - // map of taxid -> taxid - int2int parent; - // map of taxid -> rank - int2char rank; - // map of parent -> child taxids - int2intvec child; - - if (infile_nodes != NULL) - parse_nodes(infile_nodes, rank, parent, child, 1); - - std::map retmap = load_bdamage3(infile_bdamage, howmany); - // fprintf(stderr,"retmap.size():%lu\n",retmap.size()); - int2char name_map; - if (acc2tax != NULL) - name_map = parse_names(acc2tax); - - BGZF *bgfp = NULL; - samFile *samfp = NULL; - bam_hdr_t *hdr = NULL; - - if (inbam != NULL) { - if (((bgfp = bgzf_open(inbam, "r"))) == NULL) { - fprintf(stderr, "Could not open input BAM file: %s\n", inbam); - return 1; - } - } - - if (inbam != NULL) { - if (((samfp = sam_open_format(inbam, "r", NULL))) == NULL) { - fprintf(stderr, "Could not open input BAM file: %s\n", inbam); - return 1; - } - if (((hdr = sam_hdr_read(samfp))) == NULL) { - fprintf(stderr, "Could not read header for: %s\n", inbam); - return 1; - } - } - - gzFile fp = Z_NULL; - fp = gzopen(infile_lca, "rb"); - assert(fp != Z_NULL); - char buf[4096]; - char orig[4096]; - gzgets(fp, buf, 4096); // skipheader - buf[strlen(buf) - 1] = '\0'; - fprintf(stderr, "%s: VERSION:%s\n", buf, METADAMAGE_VERSION); - float presize = retmap.size(); - while (gzgets(fp, buf, 4096)) { - strncpy(orig, buf, 4096); - // fprintf(stderr,"buf: %s\n",buf); - char *tok = strtok(buf, "\t\n "); - int taxid = atoi(strtok(NULL, ":")); - // fprintf(stderr,"taxid: %d\n",taxid); - double *dbl = getval(retmap, child, taxid, howmany); - double dbldbl[3 * howmany + 1]; // 3 because ct,ga,other - dbldbl[0] = dbl[0]; - for (int i = 0; i < 3 * howmany; i++) - dbldbl[i + 1] = dbl[1 + i] / dbl[0]; - - orig[strlen(orig) - 1] = '\0'; - fprintf(stdout, "%s\t%d:%.0f", orig, taxid, dbldbl[0]); - for (int i = 0; i < 3 * howmany; i++) - fprintf(stdout, "\t%f", dbldbl[1 + i]); - fprintf(stdout, "\n"); - } - float postsize = retmap.size(); - fprintf(stderr, "\t-> pre: %f post:%f grownbyfactor: %f\n", presize, postsize, postsize / presize); - if (bgfp) - bgzf_close(bgfp); - if (hdr) - bam_hdr_destroy(hdr); - if (samfp) - sam_close(samfp); - if (fp != Z_NULL) - gzclose(fp); - return 0; -} - int2int getlcadist(char *fname) { // fprintf(stderr,"fname: %s\n",fname); int2int lcadist; - if(strlen(fname)>1000){ - fprintf(stderr,"\t-> Ridiculus long filename: \'%s\'\n",fname); + if (fname == NULL) { + fprintf(stderr, "\t-> Error: fname is NULL, will exit\n"); exit(1); } + char tmp[1024]; - snprintf(tmp, 1000, "%s.stat", fname); - fprintf(stderr, "tmp: %s\n", tmp); + + size_t len = strlen(fname); + if (len + 5 >= sizeof(tmp)) { // ".stat" + '\0' = 5 + fprintf(stderr,"\t-> Ridiculus long filename: \'%s\'\n",fname);//<- harry potter + exit(1); + } + snprintf(tmp, sizeof(tmp), "%s.stat", fname); + FILE *fp = NULL; fp = fopen(tmp, "rb"); - assert(fp != NULL); + if (fp == NULL) { + fprintf(stderr, "\t-> Error: failed to open file %s, will exit\n", tmp); + exit(1); + } char buf[4096]; while (fgets(buf, 4096, fp)) { char *tok1 = strtok(buf, "\t\n "); @@ -996,15 +974,19 @@ int2int getlcadist(char *fname) { return lcadist; } - - std::map getcsv(char *fname) { std::map ret; FILE *fp = NULL; fp = fopen(fname, "rb"); - assert(fp != NULL); + if (fp == NULL) { + fprintf(stderr, "\t-> Error: failed to open file %s, will exit\n", fname); + exit(1); + } char buf[4096]; - assert(fgets(buf, 4096, fp) != NULL); // skip header + if (fgets(buf, 4096, fp) == NULL) { + fprintf(stderr, "\t-> Error: failed to read header line from file, will exit\n"); + exit(1); + } while (fgets(buf, 4096, fp)) { int key = atoi(strtok(buf, "\t\n, ")); double *valval = new double[2]; @@ -1016,182 +998,116 @@ std::map getcsv(char *fname) { return ret; } -int main_merge2(int argc, char **argv) { - fprintf(stderr, "./metaDMG-cpp merge2 file.lca file.bdamage.gz christian.csv [-names file.gz -bam file.bam -howmany 5 -nodes trestructure.gz]\n"); - if (argc <= 3) - return 0; - char *infile_lca = argv[1]; - char *infile_bdamage = argv[2]; - char *infile_christian = argv[3]; - char *infile_nodes = NULL; - - char *acc2tax = NULL; - int howmany = 5; - fprintf(stdout, "#./metaDMG-cpp "); - for (int i = 0; i < argc; i++) - fprintf(stdout, "%s ", argv[i]); - fprintf(stdout, "\n"); - ++argv; - while (*(++argv)) { - if (strcasecmp("-names", *argv) == 0) - acc2tax = strdup(*(++argv)); - if (strcasecmp("-howmany", *argv) == 0) - howmany = atoi(*(++argv)); - if (strcasecmp("-nodes", *argv) == 0) - infile_nodes = strdup(*(++argv)); - } - - fprintf(stderr, "infile_lca: %s infile_bdamage: %s nodes: %s\n", infile_lca, infile_bdamage, infile_nodes); - - int2int lcadist = getlcadist(infile_lca); - std::map chris = getcsv(infile_christian); - // map of taxid -> taxid - int2int parent; - // map of taxid -> rank - int2char rank; - // map of parent -> child taxids - int2intvec child; - - if (infile_nodes != NULL) - parse_nodes(infile_nodes, rank, parent, child, 1); - - std::map retmap = load_bdamage3(infile_bdamage, howmany); - fprintf(stderr, "\t-> Number of entries in damage pattern file: %lu\n", retmap.size()); - - int2char name_map; - if (acc2tax != NULL) - name_map = parse_names(acc2tax); - - BGZF *bgfp = NULL; - - gzFile fp = Z_NULL; - fp = gzopen(infile_lca, "rb"); - assert(fp != Z_NULL); - char buf[4096]; - char orig[4096]; - gzgets(fp, buf, 4096); // skipheader - buf[strlen(buf) - 1] = '\0'; - fprintf(stdout, "%s: VERSION:%s\n", buf, METADAMAGE_VERSION); - float presize = retmap.size(); - double *rawval = new double[2]; - rawval[0] = rawval[1] = -1.0; - while (gzgets(fp, buf, 4096)) { - strncpy(orig, buf, 4096); - // fprintf(stderr,"buf: %s\n",buf); - char *tok = strtok(buf, "\t\n "); - int taxid = atoi(strtok(NULL, ":")); - // fprintf(stderr,"taxid: %d\n",taxid); - int nclass = -1; - int2int::iterator itit = lcadist.find(taxid); - if (itit != lcadist.end()) - nclass = itit->second; - double *valval = rawval; - std::map::iterator ititit = chris.find(taxid); - if (ititit != chris.end()) - valval = ititit->second; - double *dbl = getval(retmap, child, taxid, howmany); - double dbldbl[3 * howmany + 1]; // 3 because ct,ga,other - dbldbl[0] = dbl[0]; - for (int i = 0; i < 3 * howmany; i++) - dbldbl[i + 1] = dbl[1 + i] / dbl[0]; - - // orig[strlen(orig)-1] = '\0'; - // rollout results - tok = strtok(orig, "\t\n "); - fprintf(stdout, "%s\t%d:%.0f:%d:%f:%f\t", tok, taxid, dbldbl[0], nclass, valval[0], valval[1]); - for (int i = 0; i < 3 * howmany - 1; i++) - fprintf(stdout, "%f:", dbldbl[1 + i]); - fprintf(stdout, "%f", dbldbl[3 * howmany]); - while (((tok = strtok(NULL, "\t\n ")))) { - fprintf(stdout, "\t%s", tok); - } - fprintf(stdout, "\n"); - } - float postsize = retmap.size(); - fprintf(stderr, "\t-> pre: %f post:%f grownbyfactor: %f\n", presize, postsize, postsize / presize); - if (bgfp) - bgzf_close(bgfp); - if (fp != Z_NULL) - gzclose(fp); - return 0; -} - +//this very verbose function is translated by a thinking machine int main_print_all(int argc, char **argv) { - fprintf(stderr, "./metaDMG-cpp print_all file.bdamage.gz -names file.gz -nodes trestructure.gz -names fil.gz\n"); + fprintf(stderr, "./metaDMG-cpp print_all file.bdamage.gz -names file.gz -nodes trestructure.gz\n"); + if (argc <= 2) return 0; + char *infile_bdamage = NULL; char *infile_nodes = NULL; char *infile_names = NULL; - while (*(++argv)) { - if (strcasecmp("-names", *argv) == 0) + if (strcasecmp("-names", *argv) == 0) { + if (*(argv + 1) == NULL) { + fprintf(stderr, "\t-> Error: -names requires an argument, will exit\n"); + exit(1); + } infile_names = strdup(*(++argv)); - else if (strcasecmp("-nodes", *argv) == 0) + } else if (strcasecmp("-nodes", *argv) == 0) { + if (*(argv + 1) == NULL) { + fprintf(stderr, "\t-> Error: -nodes requires an argument, will exit\n"); + exit(1); + } infile_nodes = strdup(*(++argv)); - else + } else { infile_bdamage = strdup(*argv); + } + } + + if (infile_bdamage == NULL) { + fprintf(stderr, "\t-> Error: no input bdamage file provided, will exit\n"); + exit(1); + } + if (infile_names == NULL) { + fprintf(stderr, "\t-> Error: no names file provided, will exit\n"); + exit(1); } - fprintf(stderr, "infile_names: %s infile_bdamage: %s nodes: %s ", infile_names, infile_bdamage, infile_nodes); + fprintf(stderr, "infile_names: %s infile_bdamage: %s nodes: %s ", + infile_names ? infile_names : "NULL", + infile_bdamage ? infile_bdamage : "NULL", + infile_nodes ? infile_nodes : "NULL"); fprintf(stderr, "#VERSION:%s\n", METADAMAGE_VERSION); - // map of taxid -> taxid + int2int parent; - // map of taxid -> rank int2char rank; - // map of parent -> child taxids int2intvec child; if (infile_nodes != NULL) parse_nodes(infile_nodes, rank, parent, child, 1); + int howmany = 5; std::map retmap = load_bdamage3(infile_bdamage, howmany); fprintf(stderr, "\t-> number of entries in damage pattern file: %lu\n", retmap.size()); + int2char name = parse_names(infile_names); - BGZF *bgfp = NULL; - - float presize = retmap.size(); + size_t presize = retmap.size(); + + getval(retmap, child, 1, howmany); + + size_t postsize = retmap.size(); + fprintf(stderr, "\t-> pre: %lu post:%lu grownbyfactor: %f\n", + (unsigned long)presize, + (unsigned long)postsize, + presize ? (double)postsize / presize : 0.0); - getval(retmap, child, 1, howmany); // this will do everything - float postsize = retmap.size(); - fprintf(stderr, "\t-> pre: %f post:%f grownbyfactor: %f\n", presize, postsize, postsize / presize); for (std::map::iterator it = retmap.begin(); it != retmap.end(); it++) { int taxid = it->first; double *dbl = it->second; + char *myrank = NULL; char *myname = NULL; + int2char::iterator itc = rank.find(taxid); if (itc != rank.end()) myrank = itc->second; + itc = name.find(taxid); if (itc != name.end()) myname = itc->second; - double dbldbl[3 * howmany + 1]; // 3 because ct,ga,other - dbldbl[0] = dbl[0]; - for (int i = 0; i < 3 * howmany; i++) - dbldbl[i + 1] = dbl[1 + i] / dbl[0]; - if (dbldbl[0] > 0) { - fprintf(stdout, "%d:\"%s\":\"%s\":%.0f", taxid, myname, myrank, dbldbl[0]); + + if (dbl[0] > 0) { + double dbldbl[3 * 5 + 1]; + dbldbl[0] = dbl[0]; + for (int i = 0; i < 3 * howmany; i++) + dbldbl[i + 1] = dbl[1 + i] / dbl[0]; + + fprintf(stdout, "%d:\"%s\":\"%s\":%.0f", + taxid, + myname ? myname : "NULL", + myrank ? myrank : "NULL", + dbldbl[0]); for (int i = 0; i < 3 * howmany; i++) fprintf(stdout, "\t%f", dbldbl[1 + i]); fprintf(stdout, "\n"); - }; + } } - if (bgfp) - bgzf_close(bgfp); + free(infile_bdamage); + free(infile_nodes); + free(infile_names); return 0; } - int main_print_ugly(int argc, char **argv) { htsFormat *dingding2 = (htsFormat *)calloc(1, sizeof(htsFormat)); - if(1) - fprintf(stderr,"\t-> print_ugly functionality will be removed\n"); + fprintf(stderr,"\t-> print_ugly functionality will be removed\n"); - fprintf(stderr, "./metaDMG-cpp print_ugly file.bdamage.gz --names file.gz --nodes trestructure.gz --lcastat file.gz --out_prefix out\n"); + fprintf(stderr, "./metaDMG-cpp print_ugly file.bdamage.gz --names file.gz --nodes trestructure.gz --lcastat file.gz --out_prefix out\n"); if (argc <= 1) return 0; char *infile_bdamage = NULL; @@ -1201,6 +1117,7 @@ int main_print_ugly(int argc, char **argv) { char *infile_bam = NULL; char *out_prefix = NULL; int howmany; + int abe = 4;//<- monkey while (*(++argv)) { if (strcasecmp("--names", *argv) == 0) @@ -1324,9 +1241,9 @@ int main_print_ugly(int argc, char **argv) { std::map::iterator itold = retmap.find(it->first); size_t nalign = 0; if (itold == retmap.end()) { - int abe=2; - //fprintf(stderr, "[%s]\t-> Problem finding taxid: %d\n",__FUNCTION__, it->first); - // exit(0); + if(abe>0) + fprintf(stderr, "[%s]\t-> Problem finding taxid: %d will stop print this in: %d\n",__FUNCTION__, it->first,abe--); + } else nalign = itold->second.nal; char *myrank = NULL; @@ -1398,10 +1315,7 @@ int main(int argc, char **argv) { #endif fprintf(stderr, "./metaDMG-cpp pmd [other options]\n"); fprintf(stderr, "./metaDMG-cpp getdamage file.bam\n"); - fprintf(stderr, "./metaDMG-cpp mergedamage -b *bdamage.gz -r *rlens.gz\n"); fprintf(stderr, "./metaDMG-cpp index files.damage.gz\n"); - fprintf(stderr, "./metaDMG-cpp merge files.lca files.bdamage.gz\n"); - fprintf(stderr, "./metaDMG-cpp merge2 files.lca files.bdamage.gz christian.csv\n"); fprintf(stderr, "./metaDMG-cpp lca [many options]\n"); fprintf(stderr, "./metaDMG-cpp print bdamage.gz\n"); fprintf(stderr, "./metaDMG-cpp print2 [many options] bdamage.gz\n"); @@ -1439,10 +1353,6 @@ int main(int argc, char **argv) { main_aggregate(argc, argv); else if (!strcmp(argv[0], "print2")) main_print2(argc, argv); - else if (!strcmp(argv[0], "merge")) - main_merge(argc, argv); - else if (!strcmp(argv[0], "merge2")) - main_merge2(argc, argv); else if (!strcmp(argv[0], "lca")) main_lca(argc, argv); else{ diff --git a/misc/Likelihood.cpp b/misc/Likelihood.cpp index de2d438..b5a06bc 100644 --- a/misc/Likelihood.cpp +++ b/misc/Likelihood.cpp @@ -55,6 +55,7 @@ double log_exponential(double x, double loc, double scale){ } double compute_log_likelihood(const double DMGparam[], const void *dats){ + (void) dats; //Compute the log-likelihood across all positions and return the result /*double part1=0; double part2=0; @@ -132,7 +133,7 @@ int main() { int num_rows, num_cols; read_count_matrix("MycoBactBamSEOutSortMDSortN.mismatches.txt.gz", M3,tax_id,dir,&num_rows, &num_cols); - Alter_count_matrix(M3,tax_id,dir,num_rows,num_cols); + Alter_count_matrix(M3,tax_id,dir,num_rows); int n_priors=8; double* priors = (double*) malloc(n_priors*sizeof(double)); diff --git a/misc/M3Load.cpp b/misc/M3Load.cpp index b35f503..d632fe0 100644 --- a/misc/M3Load.cpp +++ b/misc/M3Load.cpp @@ -83,13 +83,13 @@ void read_count_matrix(const char* filename, double M3[MAX_ROWS][MAX_COLS],char gzclose(file); } -void Alter_count_matrix(double M3[MAX_ROWS][MAX_COLS],char tax_id[MAX_ROWS][MAX_COLS],char dir[MAX_ROWS][MAX_COLS],int num_rows,int num_cols){ +void Alter_count_matrix(double M3[MAX_ROWS][MAX_COLS],char tax_id[MAX_ROWS][MAX_COLS],char dir[MAX_ROWS][MAX_COLS],int num_rows){ + (void) tax_id; // Alters / adds additional columns to the count M^3 matrix, needed for log-likelihood calculations double C_total, G_total = 0; int k_sum_total = 0; - double N_in_group[MAX_ROWS]; int Max_N_in_group = 0; int Min_N_in_group = 0; @@ -151,7 +151,7 @@ int main() { int num_rows, num_cols; read_count_matrix("MycoBactBamSEOutSortMDSortN.mismatches.txt.gz", M3,tax_id,dir,&num_rows, &num_cols); printf("Read %d rows and %d columns\n", num_rows, num_cols); - Alter_count_matrix(M3,tax_id,dir,num_rows,num_cols); + Alter_count_matrix(M3,tax_id,dir,num_rows); return 0; } @@ -159,4 +159,4 @@ int main() { //g++ M3Load.cpp -std=c++11 -lm -lz -o M3 //g++ M3Load.cpp -c -std=c++11 -lm -lz -//g++ M3Load.cpp -std=c++11 -lm -lz -D __WITH_MAIN__ -o M3 \ No newline at end of file +//g++ M3Load.cpp -std=c++11 -lm -lz -D __WITH_MAIN__ -o M3 diff --git a/misc/M3Load.h b/misc/M3Load.h index fc7ed41..619943f 100644 --- a/misc/M3Load.h +++ b/misc/M3Load.h @@ -9,6 +9,6 @@ extern char dir[MAX_ROWS][MAX_COLS]; void read_count_matrix(const char* filename, double M3[MAX_ROWS][MAX_COLS],char tax_id[MAX_ROWS][MAX_COLS],char dir[MAX_ROWS][MAX_COLS],int* num_rows, int* num_cols); -void Alter_count_matrix(double M3[MAX_ROWS][MAX_COLS],char tax_id[MAX_ROWS][MAX_COLS],char dir[MAX_ROWS][MAX_COLS],int num_rows,int num_cols); +void Alter_count_matrix(double M3[MAX_ROWS][MAX_COLS],char tax_id[MAX_ROWS][MAX_COLS],char dir[MAX_ROWS][MAX_COLS],int num_rows); -#endif \ No newline at end of file +#endif diff --git a/misc/MAP.cpp b/misc/MAP.cpp index 96c0518..87b21ec 100644 --- a/misc/MAP.cpp +++ b/misc/MAP.cpp @@ -4,7 +4,6 @@ #include #include #include -#include #include "bfgs.h" #include "M3Load.h" @@ -35,7 +34,7 @@ int HelpPage(FILE *fp){ return 0; } -argStruct *getpars(int argc,char ** argv){ +argStruct *getpars(int /*argc*/,char ** argv){ argStruct *mypars = new argStruct; mypars->M3_matrix_print = NULL; mypars->OutputStat = NULL; @@ -57,7 +56,7 @@ argStruct *getpars(int argc,char ** argv){ return mypars; } -void MAP(double M3[MAX_ROWS][MAX_COLS],double* MAP){ +void MAP(double* MAP){ double lowbound[] = {0.00000001,0.00000001,0.00000001,2}; double upbound[] = {1-0.00000001,1-0.00000001,1-0.00000001,100000}; int nbd[] = {2,2,2,1}; @@ -226,7 +225,7 @@ const char** getColumnNames(int* colnumber) { columnNames[0] = "sample"; columnNames[1] = "tax_id"; - + /* int N_sum = 0; int N_sum_fwd = 0; int N_sum_rev = 0; @@ -244,7 +243,7 @@ const char** getColumnNames(int* colnumber) { K_sum_rev += M3[i][KCOL]; } } - + */ columnNames[2] = "N_x=1_forward"; columnNames[3] = "N_x=1_reverse"; columnNames[4] = "N_sum_total"; @@ -358,11 +357,11 @@ void M3Print_to_OutStat(int argc,char **argv){ const char* filename = mypars-> OutputStat; //"outputtest.txt.gz"; read_count_matrix(M3file, M3,tax_id,dir,&num_rows, &num_cols); - Alter_count_matrix(M3,tax_id,dir,num_rows,num_cols); + Alter_count_matrix(M3,tax_id,dir,num_rows); int numpars = 5; double* LlhRes = (double*) malloc(numpars*sizeof(double)); - MAP(M3,LlhRes); + MAP(LlhRes); //fprintf(stderr,"A: %f \t q: %f \t c: %f \t phi: %f \t llh: %f \n",LlhRes[0],LlhRes[1],LlhRes[2],LlhRes[3],LlhRes[4]); int colnumber = 0; @@ -370,7 +369,10 @@ void M3Print_to_OutStat(int argc,char **argv){ gzFile gz = Z_NULL; gz = gzopen(filename,"w"); - assert(gz!=Z_NULL); + if (gz == Z_NULL) { + fprintf(stderr, "\t-> Error: gz file handle is NULL, will exit\n"); + exit(1); + } if (columnNames != NULL) { for (int i = 0; i < colnumber; i++) { diff --git a/misc/Makefile b/misc/Makefile index 1f33f53..70466c0 100644 --- a/misc/Makefile +++ b/misc/Makefile @@ -84,8 +84,9 @@ $(LIBHTS): exit 1; \ fi +KMATH_CXXFLAGS := $(filter-out -Wall -Wextra,$(CXXFLAGS)) -w kmath.o: bfgs.cpp bfgs.h - $(CXX) -c $(CPPFLAGS) $(CXXFLAGS) bfgs.cpp -o kmath.o + $(CXX) -c $(CPPFLAGS) $(KMATH_CXXFLAGS) bfgs.cpp -o kmath.o Likelihood.o: Likelihood.cpp $(CXX) -c $(CPPFLAGS) $(CXXFLAGS) Likelihood.cpp -o Likelihood.o diff --git a/misc/compressbam.cpp b/misc/compressbam.cpp index 7afa029..2104c90 100644 --- a/misc/compressbam.cpp +++ b/misc/compressbam.cpp @@ -10,13 +10,10 @@ #include #include #include -#include #include #include #include -//#include "../shared.h" - #define VERSION "0.1" @@ -82,9 +79,16 @@ size_t getrefused(samFile *htsfp,bam_hdr_t *hdr,int *keeplist,int &nkeep,int min int VERB2[2]={5,5}; void writemod(const char *outfile ,bam_hdr_t *hdr,int *keeplist,samFile *htsfp,char *mycl,int minmatch){ BGZF *fp = NULL; - fp=bgzf_open(outfile,"w5");//write bgzf with compression level5. Maybe this should be changed. - bgzf_mt(fp,nthreads,256); - assert(fp!=NULL); + fp = bgzf_open(outfile, "w5"); + + if (fp == NULL) { + fprintf(stderr, "\t-> Error: failed to open BGZF file %s, will exit\n", outfile); + exit(1); + } + if (bgzf_mt(fp, nthreads, 256) != 0) { + fprintf(stderr, "\t-> Error: failed to initialize BGZF multithreading, will exit\n"); + exit(1); + } kstring_t kstmp ={0,0,NULL}; kstring_t newhdr_ks ={0,0,NULL}; @@ -99,11 +103,19 @@ void writemod(const char *outfile ,bam_hdr_t *hdr,int *keeplist,samFile *htsfp,c continue; // fprintf(stderr,"SQ: it->first: %d\n",it->first); kstmp.l =0; - assert(sam_hdr_find_line_pos(hdr,"SQ",i,&kstmp)==0); + + if (sam_hdr_find_line_pos(hdr, "SQ", i, &kstmp) != 0) { + fprintf(stderr, "\t-> Error: failed to find SQ line at position %d in header, will exit\n", i); + exit(1); + } + kputs(kstmp.s,&newhdr_ks); kputc('\n',&newhdr_ks); if(newhdr_ks.l>10000000){ - assert(bgzf_write(fp,newhdr_ks.s,newhdr_ks.l)==(ssize_t)newhdr_ks.l); + if (bgzf_write(fp, newhdr_ks.s, newhdr_ks.l) != (ssize_t)newhdr_ks.l) { + fprintf(stderr, "\t-> Error: failed to write expected number of bytes for new header, will exit\n"); + exit(1); + } newhdr_ks.l =0; } @@ -111,29 +123,44 @@ void writemod(const char *outfile ,bam_hdr_t *hdr,int *keeplist,samFile *htsfp,c for(int i=0;i Error: failed to find RG line at position %d in header, will exit\n", i); + exit(1); + } kputs(kstmp.s,&newhdr_ks); kputc('\n',&newhdr_ks); if(newhdr_ks.l>10000000){ - assert(bgzf_write(fp,newhdr_ks.s,newhdr_ks.l)==(ssize_t)newhdr_ks.l); + if (bgzf_write(fp, newhdr_ks.s, newhdr_ks.l) != (ssize_t)newhdr_ks.l) { + fprintf(stderr, "\t-> Error: failed to write expected number of bytes for new header, will exit\n"); + exit(1); + } newhdr_ks.l =0; } } for(int i=0;i Error: failed to find PG line at position %d in header, will exit\n", i); + exit(1); + } kputs(kstmp.s,&newhdr_ks); kputc('\n',&newhdr_ks); } for(int i=0;i Error: failed to find CO line at position %d in header, will exit\n", i); + exit(1); + } kputs(kstmp.s,&newhdr_ks); kputc('\n',&newhdr_ks); } - assert(bgzf_write(fp,newhdr_ks.s,newhdr_ks.l)==(ssize_t)newhdr_ks.l); + if (bgzf_write(fp, newhdr_ks.s, newhdr_ks.l) != (ssize_t)newhdr_ks.l) { + fprintf(stderr, "\t-> Error: failed to write expected number of bytes for new header, will exit\n"); + exit(1); + } free(newhdr_ks.s); free(kstmp.s); bgzf_close(fp); @@ -141,10 +168,16 @@ void writemod(const char *outfile ,bam_hdr_t *hdr,int *keeplist,samFile *htsfp,c samFile *tmpfp = NULL; tmpfp = hts_open(outfile,"r"); - assert(tmpfp!=NULL); + if (tmpfp == NULL) { + fprintf(stderr, "\t-> Error: failed to open file %s with hts_open, will exit\n", outfile); + exit(1); + } sam_hdr_t *newhdr = NULL; newhdr = sam_hdr_read(tmpfp); - assert(newhdr!=NULL); + if (newhdr == NULL) { + fprintf(stderr, "\t-> Error: failed to read SAM header from file, will exit\n"); + exit(1); + } sam_close(tmpfp); //fprintf(stderr,"newhdr: %s\n",sam_hdr_str(newhdr)); @@ -170,12 +203,18 @@ void writemod(const char *outfile ,bam_hdr_t *hdr,int *keeplist,samFile *htsfp,c continue; keeplist[i] = sam_hdr_name2tid(newhdr,sam_hdr_tid2name(hdr,i)); - assert(keeplist[i]>=0); + if (keeplist[i] < 0) { + fprintf(stderr, "\t-> Error: keeplist[%d] is negative (%d), will exit\n", i, keeplist[i]); + exit(1); + } } sam_hdr_add_pg(newhdr,"compressbam","VN",VERSION,"CL",mycl,NULL); int ret = sam_hdr_write(outhts,newhdr); - assert(ret>=0); + if (ret < 0) { + fprintf(stderr, "\t-> Error: return value ret is negative (%d), will exit\n", ret); + exit(1); + } fprintf(stderr,"\t-> Done writing new header as binary\n");fflush(stderr); //now mainloop bam1_t *aln = bam_init1(); @@ -201,7 +240,10 @@ void writemod(const char *outfile ,bam_hdr_t *hdr,int *keeplist,samFile *htsfp,c aln->core.tid = keeplist[aln->core.tid]; if(aln->core.mtid!=-1) aln->core.mtid = keeplist[aln->core.mtid]; - assert(sam_write1(outhts, newhdr,aln)>=0); + if (sam_write1(outhts, newhdr, aln) < 0) { + fprintf(stderr, "\t-> Error: failed to write alignment with sam_write1, will exit\n"); + exit(1); + } } sam_flush(htsfp); sam_hdr_destroy(newhdr); @@ -219,7 +261,7 @@ int main(int argc,char**argv){ if(argc==1){ fprintf(stderr,"./compressbam --threads --input --ref --output [--min-match ]\n"); - return -1; + return 0; } char *mycl = stringify_argv(argc,argv); fprintf(stderr,"\t-> compressbam: (%s;%s;%s): \'%s\'\n",__FILE__,__DATE__,__TIME__,mycl); diff --git a/misc/dfit_optim.cpp b/misc/dfit_optim.cpp index ac15483..2ad9adb 100644 --- a/misc/dfit_optim.cpp +++ b/misc/dfit_optim.cpp @@ -1,5 +1,4 @@ #include -#include #include #include #include @@ -10,7 +9,12 @@ double **read1_ugly_matrix(const char *fname){ fprintf(stderr,"\t-> Reading file: \'%s\'\n",fname); gzFile gz = Z_NULL; - assert((gz=gzopen(fname,"rb"))!=Z_NULL); + gz = gzopen(fname, "rb"); + if (gz == Z_NULL) { + fprintf(stderr, "\t-> Error: failed to open gz file %s, will exit\n", fname); + exit(1); + } + char buf[4096]; char *taxid = NULL; gzgets(gz,buf,4096); @@ -19,12 +23,25 @@ double **read1_ugly_matrix(const char *fname){ std::vector ncol; while(gzgets(gz,buf,4096)){ - char *tok = strtok(buf,"\t\n "); - if(taxid==NULL) + char *tok = strtok(buf, "\t\n "); + + if (tok == NULL) { + fprintf(stderr, "\t-> Error: failed to parse token from buffer, will exit\n"); + exit(1); + } + + if (taxid == NULL) { taxid = strdup(tok); - else - assert(strcmp(taxid,tok)==0); - + if (taxid == NULL) { + fprintf(stderr, "\t-> Error: failed to allocate memory for taxid, will exit\n"); + exit(1); + } + } else { + if (strcmp(taxid, tok) != 0) { + fprintf(stderr, "\t-> Error: taxid mismatch (%s vs %s), will exit\n", taxid, tok); + exit(1); + } + } tok = strtok(NULL,"\t\n "); int is5 = 1; if(strcmp(tok,"3'")==0) @@ -229,7 +246,14 @@ int main(int argc,char **argv){ double **dat = read1_ugly_matrix(fname); double pars[6] = {0.1,0.1,0.01,1000};//last one will contain the llh,and the ncall for the objective function optimoptim(pars,dat,nopt); - double stats[2+2*(int)dat[0][0]]; + + int n = 2 + 2 * (int)dat[0][0]; + + double *stats = (double *)malloc(n * sizeof(double)); + if (stats == NULL) { + fprintf(stderr, "\t-> Error: failed to allocate memory for stats, will exit\n"); + exit(1); + } getstat(dat,pars,stats); //printit @@ -249,6 +273,7 @@ int main(int argc,char **argv){ dx_conf = stats+2+nrows+ncycle; for(int i=0;i #include #include -#include #include #include "../shared.h" @@ -122,10 +121,14 @@ int2int getkeysint(const char *key,int value){ void doflush(queue *myq,int2int &keeplist,bam_hdr_t *hdr,samFile *outhts,int strict){ // fprintf(stderr,"flush: %lu strictk:%d\n",myq->l,strict); if(strict==1){//will only print specific match - for(size_t i=0;il;i++){ - int2int::iterator it=keeplist.find(myq->ary[i]->core.tid); - if(it!=keeplist.end()) - assert(sam_write1(outhts, hdr,myq->ary[i])>=0); + for (size_t i = 0; i < myq->l; i++) { + int2int::iterator it = keeplist.find(myq->ary[i]->core.tid); + if (it != keeplist.end()) { + if (sam_write1(outhts, hdr, myq->ary[i]) < 0) { + fprintf(stderr, "\t-> Error: failed to write alignment with sam_write1, will exit\n"); + exit(1); + } + } } } if(strict==0){ @@ -140,7 +143,10 @@ void doflush(queue *myq,int2int &keeplist,bam_hdr_t *hdr,samFile *outhts,int str } if(writedata>0){ for(size_t i=0;il;i++){ - assert(sam_write1(outhts, hdr,myq->ary[i])>=0); + if (sam_write1(outhts, hdr, myq->ary[i]) < 0) { + fprintf(stderr, "\t-> Error: failed to write alignment with sam_write1, will exit\n"); + exit(1); + } } } } @@ -180,7 +186,10 @@ void runextract_int2int(int2int &keeplist,samFile *htsfp,bam_hdr_t *hdr,const ch myq->l = 0; last=strdup(qname); } - assert(bam_copy1(myq->ary[myq->l],aln)!=NULL); + if (bam_copy1(myq->ary[myq->l], aln) == NULL) { + fprintf(stderr, "\t-> Error: bam_copy1 returned NULL, will exit\n"); + exit(1); + } myq->l++; if(myq->l==myq->m) expand_queue(myq); @@ -220,11 +229,17 @@ void runextract_readid(char2int &keeplist,samFile *htsfp,bam_hdr_t *hdr,const ch isthere=1; it->second = it->second + 1; } - if(complement==0&&isthere==1) - assert(sam_write1(outhts, hdr,aln)>=0); - else if(complement==1&&isthere==0) - assert(sam_write1(outhts, hdr,aln)>=0); - + if(complement==0&&isthere==1){ + if (sam_write1(outhts, hdr, aln) < 0) { + fprintf(stderr, "\t-> Error: failed to write alignment with sam_write1, will exit\n"); + exit(1); + } + }else if(complement==1&&isthere==0){ + if (sam_write1(outhts, hdr, aln) < 0) { + fprintf(stderr, "\t-> Error: failed to write alignment with sam_write1, will exit\n"); + exit(1); + } + } } for(char2int::iterator it=keeplist.begin();it!=keeplist.end();it++){ if(it->second==0) @@ -290,7 +305,14 @@ int main_byrefid(int argc,char**argv){ ++argv; } - fprintf(stderr,"\t-> key: %s \n\t-> hts: %s \n\t-> type: %s \n\t-> outfile: %s\n\t-> strict: %d\n",keyfile,hts,type,outfile,strict); + fprintf(stderr, + "\t-> key: %s \n\t-> hts: %s \n\t-> type: %s \n\t-> outfile: %s\n\t-> strict: %d\n", + keyfile ? keyfile : "NULL", + hts ? hts : "NULL", + type ? type : "NULL", + outfile ? outfile : "NULL", + strict + ); //open inputfile and parse header samFile *htsfp = hts_open(hts,"r"); @@ -363,7 +385,7 @@ int main_bytaxid(int argc,char**argv){ char *hts = NULL; char *taxid = NULL; char *nodefile = NULL; - char out_mode[5] = "wb"; + // char out_mode[5] = "wb"; char *acc2tax = NULL; int strict = 0; char *type = NULL; @@ -396,7 +418,21 @@ int main_bytaxid(int argc,char**argv){ ++argv; } - fprintf(stderr,"\t-> key: %s \n\t-> hts: %s \n\t-> nodefile: %s \n\t-> acc2tax: %s \n\t-> taxid: %s\n\t-> taxnames: %s \n\t-> strict: %d \n\t-> type: %s \n\t-> outfile: %s\n\t-> forcedump:%d\n\t-> accout:%s\n\t-> names: %s\n",keyfile,hts,nodefile, acc2tax, taxid,taxnames,strict,type,outfile,forcedump,accout,names); + fprintf(stderr, + "\t-> key: %s \n\t-> hts: %s \n\t-> nodefile: %s \n\t-> acc2tax: %s \n\t-> taxid: %s\n\t-> taxnames: %s \n\t-> strict: %d \n\t-> type: %s \n\t-> outfile: %s\n\t-> forcedump:%d\n\t-> accout:%s\n\t-> names: %s\n", + keyfile ? keyfile : "NULL", + hts ? hts : "NULL", + nodefile ? nodefile : "NULL", + acc2tax ? acc2tax : "NULL", + taxid ? taxid : "NULL", + taxnames ? taxnames : "NULL", + strict, + type ? type : "NULL", + outfile ? outfile : "NULL", + forcedump, + accout ? accout : "NULL", + names ? names : "NULL" + ); if(taxid ==NULL&&taxnames==NULL){ fprintf(stderr,"\t-> Need to supply -taxid and/or -taxnames\n"); @@ -509,7 +545,10 @@ int main_bytaxid(int argc,char**argv){ int2int keeplist; for(char2int::iterator it=cmap.begin();it!=cmap.end();it++){ int tokeep = sam_hdr_name2tid(hdr,it->first); - assert(tokeep>=0); + if (tokeep < 0) { + fprintf(stderr, "\t-> Error: tokeep is negative (%d), will exit\n", tokeep); + exit(1); + } keeplist[tokeep] =1; } diff --git a/mrand.cpp b/mrand.cpp index 2b4a069..2558613 100644 --- a/mrand.cpp +++ b/mrand.cpp @@ -3,7 +3,6 @@ #include #include #include "mrand.h" -#include mrand_t *mrand_alloc(int type_a,long int seedval){ mrand_t *ret = (mrand_t *) malloc(sizeof(mrand_t)); @@ -46,10 +45,13 @@ mrand_t *mrand_alloc(int type_a,long int seedval){ } double mrand_pop(mrand_t *mr){ - double res; + double res = -1; if(mr->type==0){ #if defined(__linux__) || defined(__unix__) drand48_r((struct drand48_data*)&mr->buf0,&res); +#else + fprintf(stderr,"\t-> Error: mr->type==0 is not supported on this platform, will exit\n"); + exit(1); #endif } else if(mr->type==1){ @@ -72,17 +74,22 @@ double mrand_pop(mrand_t *mr){ fprintf(stderr,"Random parameter %d is not supported\n",mr->type); exit(0); } - assert(res!=0.0); - assert(res!=1.0); + if (!(res > 0.0 && res < 1.0)) { + fprintf(stderr,"\t-> Unclear and unobvious subtle strange error message, will exit mrand_pop\n"); + exit(1); + } return res; } long mrand_pop_long(mrand_t *mr){ - long res; + long res = 0; if(mr->type==0){ - #if defined(__linux__) || defined(__unix__) +#if defined(__linux__) || defined(__unix__) lrand48_r((struct drand48_data*)&mr->buf0,&res); res >>= 4; - #endif +#else + fprintf(stderr, "\t-> Error: mr->type==0 is not supported on this platform, will exit\n"); + exit(1); +#endif } else if(mr->type==1){ res = mr->distrInt(mr->eng); diff --git a/ngsLCA.cpp b/ngsLCA.cpp index 1dadda4..1a77f15 100644 --- a/ngsLCA.cpp +++ b/ngsLCA.cpp @@ -11,7 +11,6 @@ int mod_out[] = {1333996, 1333996, 1582270, 1914213, 1917265, 1915309, 263865, 2 #include // for gzprintf, gzFile #include -#include // for assert #include // for pow #include // for fprintf, stderr, NULL, FILE #include // for free, exit, calloc @@ -63,7 +62,14 @@ char2int setlevels(int norank2species, char *key, int &value) { } if (norank2species) { char2int::iterator it = c2i.find(spec); - assert(it != c2i.end()); + if (it == c2i.end()) { + fprintf(stderr, "\t-> Error: iterator reached end in c2i, will exit\n"); + exit(1); + } + if (it == c2i.end()) { + fprintf(stderr, "\t-> Error: iterator reached end in c2i, will exit\n"); + exit(1); + } c2i[strdup("no rank")] = it->second; } fprintf(stderr, "\t-> Number of entries with level information: %lu \n", c2i.size()); @@ -189,7 +195,10 @@ void adder(int taxid, int readlengths, float gccontent) { } int do_lca(std::vector &taxids, int2int &parent) { // fprintf(stderr,"\t-> [%s] with number of taxids: %lu\n",__func__,taxids.size()); - assert(taxids.size() > 0); + if (taxids.size() == 0) { + fprintf(stderr, "\t-> Error: taxids vector is empty, will exit\n"); + exit(1); + } if (taxids.size() == 1) { int taxa = taxids[0]; if (parent.count(taxa) == 0) { @@ -236,7 +245,7 @@ int do_lca(std::vector &taxids, int2int &parent) { // now counts contain how many time a node is traversed to the root int2int dist2root; for (int2int::iterator it = counter.begin(); it != counter.end(); it++) - if (it->second == taxids.size()) + if (it->second ==(int) taxids.size()) dist2root[nodes2root(it->first, parent)] = it->first; for (int2int::iterator it = dist2root.begin(); 0 && it != dist2root.end(); it++) fprintf(stderr, "%d\t->%d\n", it->first, it->second); @@ -255,7 +264,10 @@ void print_chain1(kstring_t *kstr, int taxa, int2char &rank, int2char &name_map, if (it1 == name_map.end()) { fprintf(stderr, "\t-> Problem finding taxaid:%d\n", taxa); } - assert(it2 != rank.end()); + if (it2 == rank.end()) { + fprintf(stderr, "\t-> Error: iterator reached end in rank, will exit\n"); + exit(1); + } if (it1 == name_map.end() || it2 == rank.end()) { fprintf(stderr, "taxa: %d %s doesnt exists will exit\n", taxa, it1->second); exit(1); @@ -265,7 +277,10 @@ void print_chain1(kstring_t *kstr, int taxa, int2char &rank, int2char &name_map, void print_rank(FILE *fp, int taxa, int2char &rank) { int2char::iterator it2 = rank.find(taxa); - assert(it2 != rank.end()); + if (it2 == rank.end()) { + fprintf(stderr, "\t-> Error: iterator reached end in rank, will exit\n"); + exit(1); + } fprintf(fp, "taxa: %d rank %s\n", taxa, it2->second); } @@ -276,7 +291,10 @@ void print_chain(kstring_t *kstr, int taxa, int2int &parent, int2char &rank, int while (1) { print_chain1(kstr, taxa, rank, name_map,sep); int2int::iterator it = parent.find(taxa); - assert(it != parent.end()); + if (it == parent.end()) { + fprintf(stderr, "\t-> Error: iterator reached end in parent, will exit\n"); + exit(1); + } if (taxa == it->second) //<- catch root break; taxa = it->second; @@ -306,9 +324,10 @@ int get_species1(int taxa, int2int &parent, int2char &rank) { if (it == rank.end()) return -1; char *val = it->second; - if (val == NULL) - fprintf(stderr, "taxa:%d\n", taxa); - assert(val); + if (val == NULL) { + fprintf(stderr, "\t-> Error: rank value is NULL for taxa:%d, will exit\n", taxa); + exit(1); + } if (!strcmp(val, "species")) return taxa; int next = parent[taxa]; @@ -370,7 +389,7 @@ int refToInt[256] = { float gccontent(char *seq) { int counts[5] = {0, 0, 0, 0, 0}; for (size_t i = 0; i < strlen(seq); i++) - counts[refToInt[seq[i]]]++; + counts[refToInt[(unsigned char)seq[i]]]++; float gcs = counts[1] + counts[2]; float tot = counts[0] + counts[1] + counts[2] + counts[3]; @@ -383,7 +402,7 @@ float gccontent(bam1_t *aln) { int len = aln->core.l_qseq; uint8_t *q = bam_get_seq(aln); for (int i = 0; i < len; i++) - counts[refToInt[seq_nt16_str[bam_seqi(q, i)]]]++; + counts[refToInt[(unsigned char)seq_nt16_str[bam_seqi(q, i)]]]++; float gcs = counts[1] + counts[2]; float tot = counts[0] + counts[1] + counts[2] + counts[3]; @@ -395,7 +414,10 @@ int printonce = 1; std::vector purge(std::vector &taxids, std::vector &editdist) { if (printonce == 1) fprintf(stderr, "\t-> purging taxids oldsize:%lu\n", taxids.size()); - assert(taxids.size() == editdist.size()); + if (taxids.size() != editdist.size()) { + fprintf(stderr, "\t-> Error: taxids and editdist sizes do not match, will exit\n"); + exit(1); + } std::vector tmpnewvec; int mylow = *std::min_element(editdist.begin(), editdist.end()); for (size_t i = 0; i < taxids.size(); i++) @@ -408,7 +430,10 @@ std::vector purge(std::vector &taxids, std::vector &editdist) { void hts(gzFile fp, samFile *fp_in, int2int &i2i, int2int &parent, bam_hdr_t *hdr, int2char &rank, int2char &name_map, int minmapq, int discard, int editMin, int editMax, double scoreLow, double scoreHigh, int minlength, int lca_rank, char *prefix, int howmany, samFile *fp_usedreads, int skipnorank, int2int &rank2level, int nthreads, int weighttype,long maxreads,samFile *fp_famout,int rlens_flat_out) { fprintf(stderr, "[%s] \t-> editMin:%d editmMax:%d scoreLow:%f scoreHigh:%f minlength:%d discard: %d prefix: %s howmany: %d skipnorank: %d weighttype: %d maxreads: %ld rlens_flat_out: %d\n", __FUNCTION__, editMin, editMax, scoreLow, scoreHigh, minlength, discard, prefix, howmany, skipnorank, weighttype,maxreads,rlens_flat_out); - assert(fp_in != NULL); + if (fp_in == NULL) { + fprintf(stderr, "\t-> Error: fp_in is NULL (failed to open input file), will exit\n"); + exit(1); + } damage *dmg = new damage(howmany, nthreads, 13); bam1_t *aln = bam_init1(); // initialize an alignment @@ -463,7 +488,10 @@ void hts(gzFile fp, samFile *fp_in, int2int &i2i, int2int &parent, bam_hdr_t *hd // fprintf(stderr,"length of taxids:%lu and other:%lu minedit:%d\n",taxids.size(),editdist.size(),*std::min_element(editdist.begin(),editdist.end())); size_t size = taxids.size(); - assert(size == myq->l); + if (size != myq->l) { + fprintf(stderr, "\t-> Error: size does not match myq->l, will exit\n"); + exit(1); + } if (editMin == -1 && editMax == -1) taxids = purge(taxids, editdist); nreads++; @@ -487,13 +515,20 @@ void hts(gzFile fp, samFile *fp_in, int2int &i2i, int2int &parent, bam_hdr_t *hd // lca_rank is integer and is different from minus one int2int::iterator myit = rank2level.find(lca); - assert(myit != rank2level.end()); + if (myit == rank2level.end()) { + fprintf(stderr, "\t-> Error: iterator reached end in rank2level, will exit\n"); + exit(1); + } //write to family out if (myit->second != -1 && (myit->second <= 16)) { //family is 16, see rank2level.txt for (size_t i = 0; i < myq->l; i++) - if(fp_famout) - assert(sam_write1(fp_famout, hdr, myq->ary[i]) >= 0); + if(fp_famout){ + if (sam_write1(fp_famout, hdr, myq->ary[i]) < 0) { + fprintf(stderr, "\t-> Error: failed to write record with sam_write1, will exit\n"); + exit(1); + } + } } //standard analyses if (myit->second != -1 && (myit->second <= lca_rank)) { @@ -502,7 +537,10 @@ void hts(gzFile fp, samFile *fp_in, int2int &i2i, int2int &parent, bam_hdr_t *hd for (size_t i = 0; i < myq->l; i++) { int2int::iterator it2k = i2i.find(myq->ary[i]->core.tid); - assert(it2k != i2i.end()); + if (it2k == i2i.end()) { + fprintf(stderr, "\t-> Error: iterator reached end in i2i, will exit\n"); + exit(1); + } int2char::iterator ititit = rank.find(it2k->second); if (ititit == rank.end()) { fprintf(stderr, "\t-> Potential problem no rank for taxid: %d\n", it2k->second); @@ -511,8 +549,12 @@ void hts(gzFile fp, samFile *fp_in, int2int &i2i, int2int &parent, bam_hdr_t *hd if (skipnorank == 1 && strcasecmp(ititit->second, "no rank") == 0) continue; dmg->damage_analysis(myq->ary[i], it2k->second, weighttype == 0 ? 1 : (1.0 / (float)myq->l)); - if (fp_usedreads) - assert(sam_write1(fp_usedreads, hdr, myq->ary[i]) >= 0); + if (fp_usedreads){ + if (sam_write1(fp_usedreads, hdr, myq->ary[i]) < 0) { + fprintf(stderr, "\t-> Error: failed to write record with sam_write1 (fp_usedreads), will exit\n"); + exit(1); + } + } } } } @@ -576,7 +618,10 @@ void hts(gzFile fp, samFile *fp_in, int2int &i2i, int2int &parent, bam_hdr_t *hd } else it2missing->second = it2missing->second + 1; } else { - assert(bam_copy1(myq->ary[myq->l], aln) != NULL); + if (bam_copy1(myq->ary[myq->l], aln) == NULL) { + fprintf(stderr, "\t-> Error: bam_copy1 returned NULL, will exit\n"); + exit(1); + } myq->l++; if (myq->l == myq->m) expand_queue(myq); @@ -590,7 +635,10 @@ void hts(gzFile fp, samFile *fp_in, int2int &i2i, int2int &parent, bam_hdr_t *hd } if (taxids.size() > 0 && skip == 0) { size_t size = taxids.size(); - assert(size == myq->l); + if (size != myq->l) { + fprintf(stderr, "\t-> Error: size does not match myq->l, will exit\n"); + exit(1); + } if (editMin == -1 && editMax == -1) taxids = purge(taxids, editdist); @@ -612,13 +660,20 @@ void hts(gzFile fp, samFile *fp_in, int2int &i2i, int2int &parent, bam_hdr_t *hd } // lca_rank is integer and is different from minus one int2int::iterator myit = rank2level.find(lca); - assert(myit != rank2level.end()); + if (myit == rank2level.end()) { + fprintf(stderr, "\t-> Error: iterator reached end in rank2level, will exit\n"); + exit(1); + } //write to family out if (myit->second != -1 && (myit->second <= 16)) { //family is 16, see rank2level.txt for (size_t i = 0; i < myq->l; i++) - if(fp_famout) - assert(sam_write1(fp_famout, hdr, myq->ary[i]) >= 0); + if(fp_famout){ + if (sam_write1(fp_famout, hdr, myq->ary[i]) < 0) { + fprintf(stderr, "\t-> Error: failed to write record with sam_write1 (fp_famout), will exit\n"); + exit(1); + } + } } if (myit->second != -1 && (myit->second <= lca_rank)) { adder(lca, strlen(seq), gccontent(seq)); @@ -627,7 +682,10 @@ void hts(gzFile fp, samFile *fp_in, int2int &i2i, int2int &parent, bam_hdr_t *hd for (size_t i = 0; i < myq->l; i++) { // dmg->damage_analysis(myq->ary[i],myq->ary[i]->core.tid); int2int::iterator ittt = i2i.find(myq->ary[i]->core.tid); - assert(ittt != i2i.end()); + if (ittt == i2i.end()) { + fprintf(stderr, "\t-> Error: iterator reached end in i2i (ittt), will exit\n"); + exit(1); + } int2char::iterator ititit = rank.find(ittt->second); if (ititit == rank.end()) { fprintf(stderr, "\t-> Potential problem no rank for taxid: %d\n", ititit->first); @@ -638,8 +696,12 @@ void hts(gzFile fp, samFile *fp_in, int2int &i2i, int2int &parent, bam_hdr_t *hd continue; // fprintf(stderr,"uaua\n"); dmg->damage_analysis(myq->ary[i], ittt->second, weighttype == 0 ? 1 : (1.0 / (float)myq->l)); - if (fp_usedreads) - assert(sam_write1(fp_usedreads, hdr, myq->ary[i]) >= 0); + if (fp_usedreads){ + if (sam_write1(fp_usedreads, hdr, myq->ary[i]) < 0) { + fprintf(stderr, "\t-> Error: failed to write record with sam_write1 (fp_usedreads), will exit\n"); + exit(1); + } + } } } } diff --git a/ngsLCA_cli.cpp b/ngsLCA_cli.cpp index 3524c82..eda3391 100644 --- a/ngsLCA_cli.cpp +++ b/ngsLCA_cli.cpp @@ -1,6 +1,5 @@ #include "ngsLCA_cli.h" -#include // for assert #include // for hts_open #include // for pthread_mutex_lock, pthre... #include @@ -125,9 +124,17 @@ char2int *ass2bin(const char *fname, int redo) { while (sizeof(int) == gzread(FP, &key_l, sizeof(int))) { char *key = (char *)calloc(key_l + 1, sizeof(char)); - assert((key_l = gzread(FP, key, key_l))); - int val; - assert(sizeof(int) == gzread(FP, &val, sizeof(int))); + key_l = gzread(FP, key, key_l); + if (key_l <= 0) { + fprintf(stderr, "\t-> Error: failed to read key with gzread, will exit\n"); + exit(1); + } + + int val; + if (gzread(FP, &val, sizeof(int)) != sizeof(int)) { + fprintf(stderr, "\t-> Error: failed to read expected number of bytes for int value with gzread, will exit\n"); + exit(1); + } if (cm->find(key) != cm->end()) { fprintf(stderr, "\t-> Duplicate entries found \'%s\'\n", key); } else @@ -275,12 +282,18 @@ pars *get_pars(int argc, char **argv) { snprintf(buf, 1024, "%s.lca.gz", p->outnames); fprintf(stderr, "\t-> Will output lca results in file:\t\t\'%s\'\n", buf); p->fp1 = gzopen(buf, "wb"); - assert(p->fp1); + if (!p->fp1) { + fprintf(stderr, "\t-> Error: p->fp1 is NULL (file pointer not initialized), will exit\n"); + exit(1); + } snprintf(buf, 1024, "%s.stat.gz", p->outnames); fprintf(stderr, "\t-> Will output lca distribution in file:\t\t\'%s\'\n", buf); p->fp_lcadist = NULL; p->fp_lcadist = gzopen(buf, "wb"); - assert(p->fp_lcadist); + if (!p->fp_lcadist) { + fprintf(stderr, "\t-> Error: p->fp_lcadist is NULL (file pointer not initialized), will exit\n"); + exit(1); + } snprintf(buf, 1024, "%s.wlca.gz", p->outnames); fprintf(stderr, "\t-> Will output lca weight in file:\t\t\'%s\'\n", buf); // p->fp2 = gzopen(buf,"wb"); @@ -288,7 +301,10 @@ pars *get_pars(int argc, char **argv) { snprintf(buf, 1024, "%s.log", p->outnames); fprintf(stderr, "\t-> Will output log info (problems) in file:\t\'%s\'\n", buf); p->fp3 = fopen(buf, "wb"); - assert(p->fp3); + if (!p->fp3) { + fprintf(stderr, "\t-> Error: p->fp3 is NULL (file pointer not initialized), will exit\n"); + exit(1); + } #endif if (make_used_reads) { snprintf(buf, 1024, "%s.usedreads.bam", p->outnames); @@ -334,8 +350,10 @@ void print_pars(FILE *fp, pars *p) { #ifdef __WITH_MAIN__ int main(int argc, char **argv) { pars *p = get_pars(--argc, ++argv); - assert(p); - + if (!p) { + fprintf(stderr, "\t-> Error: get_pars returned NULL, will exit\n"); + exit(1); + } print_pars(stdout, p); } #endif diff --git a/profile.cpp b/profile.cpp index 300eb2e..eb5f47f 100644 --- a/profile.cpp +++ b/profile.cpp @@ -8,7 +8,6 @@ #include // for int32_t, uint32_t, uint8_t #include // for max -#include // for assert #include // for log10 #include // for strlen, strtok, memset, strdup #include // for vector @@ -59,7 +58,10 @@ void destroy_damage(damage *dmg) { BGZF *my_bgzf_open(const char *name, int nthreads) { BGZF *ret = NULL; ret = bgzf_open(name, "wb"); - assert(ret != NULL); + if (ret == NULL) { + fprintf(stderr, "\t-> Error: ret is NULL, will exit\n"); + exit(1); + } if (nthreads > 1) { fprintf(stderr, "\t-> Setting threads to: %d \n", nthreads); bgzf_mt(ret, nthreads, 64); @@ -300,8 +302,14 @@ inline void increaseCounters(const bam1_t *b, const char *reconstructedReference //this function will not get a commen void mysuperduper_function(triple &t, size_t old_m, size_t new_m) { - assert(t.rlens!=NULL&&old_m Error: invalid state (t.rlens NULL or old_m >= new_m), will exit\n"); + exit(1); + } + if ((size_t)t.rlens_m != old_m) { + fprintf(stderr, "\t-> Error: t.rlens_m does not match old_m, will exit\n"); + exit(1); + } size_t *tmp = new size_t[new_m]; memset(tmp, 0, new_m * sizeof(size_t)); memcpy(tmp, t.rlens, t.rlens_m * sizeof(size_t)); diff --git a/shared.cpp b/shared.cpp index df09803..05d2980 100644 --- a/shared.cpp +++ b/shared.cpp @@ -8,7 +8,6 @@ #include // for isatty #include // for gzclose, gzgets, gzopen, Z_NULL, gzFile -#include // for assert #include // for fprintf, stderr, snprintf, NULL, size_t #include // for atoi, exit, malloc #include // for time @@ -39,19 +38,33 @@ BGZF *getbgzf(const char *str1, const char *mode, int nthreads) { } BGZF *getbgzf2(const char *str1, const char *str2, const char *mode, int nthreads) { - unsigned tmp_l = strlen(str1) + strlen(str2) + 5; - char tmp[tmp_l]; + unsigned tmp_l = strlen(str1) + strlen(str2) + 1; // +1 for '\0' + + char *tmp = (char*)malloc(tmp_l); + if (!tmp) { + fprintf(stderr, "\t-> Error: failed to allocate memory for tmp, will exit\n"); + exit(1); + } + snprintf(tmp, tmp_l, "%s%s", str1, str2); - return getbgzf(tmp, mode, nthreads); + + BGZF *ret = getbgzf(tmp, mode, nthreads); + free(tmp); + return ret; } BGZF *getbgzf3(const char *str1, const char *str2, const char *str3, const char *mode, int nthreads) { - unsigned tmp_l = strlen(str1) + strlen(str2) + strlen(str3) + 5; - char tmp[tmp_l]; + size_t tmp_l = strlen(str1) + strlen(str2) + strlen(str3) + 1; + char *tmp = (char *)malloc(tmp_l); + if (tmp == NULL) { + fprintf(stderr, "\t-> Error: failed to allocate memory for tmp, will exit\n"); + exit(1); + } snprintf(tmp, tmp_l, "%s%s%s", str1, str2, str3); - return getbgzf(tmp, mode, nthreads); + BGZF *ret = getbgzf(tmp, mode, nthreads); + free(tmp); + return ret; } - char *getfilename4(const char *str1, const char *str2, const char *str3, const char *str4) { unsigned tmp_l = strlen(str1) + strlen(str2) + strlen(str3) + strlen(str4) + 5; char *tmp=(char*)calloc(tmp_l,1); @@ -60,12 +73,20 @@ char *getfilename4(const char *str1, const char *str2, const char *str3, const c } - BGZF *getbgzf4(const char *str1, const char *str2, const char *str3, const char *str4, const char *mode, int nthreads) { - unsigned tmp_l = strlen(str1) + strlen(str2) + strlen(str3) + strlen(str4) + 5; - char tmp[tmp_l]; + size_t tmp_l = strlen(str1) + strlen(str2) + strlen(str3) + strlen(str4) + 1; + + char *tmp = (char *)malloc(tmp_l); + if (tmp == NULL) { + fprintf(stderr, "\t-> Error: failed to allocate memory for tmp, will exit\n"); + exit(1); + } + snprintf(tmp, tmp_l, "%s%s%s%s", str1, str2, str3, str4); - return getbgzf(tmp, mode, nthreads); + + BGZF *ret = getbgzf(tmp, mode, nthreads); + free(tmp); + return ret; } int fexists(const char *str) { ///@param str Filename given as a string. @@ -75,12 +96,20 @@ int fexists(const char *str) { ///@param str Filename given as a string. } int fexists2(const char *str1, const char *str2) { - unsigned tmp_l = strlen(str1) + strlen(str2) + 5; - char tmp[tmp_l]; + size_t tmp_l = strlen(str1) + strlen(str2) + 1; + + char *tmp = (char *)malloc(tmp_l); + if (tmp == NULL) { + fprintf(stderr, "\t-> Error: failed to allocate memory for tmp, will exit\n"); + exit(1); + } + snprintf(tmp, tmp_l, "%s%s", str1, str2); - return fexists(tmp); -} + int ret = fexists(tmp); + free(tmp); + return ret; +} size_t fsize(const char *fname) { struct stat st; stat(fname, &st); @@ -88,22 +117,40 @@ size_t fsize(const char *fname) { } int fexists3(const char *str1, const char *str2, const char *str3) { - unsigned tmp_l = strlen(str1) + strlen(str2) + strlen(str3) + 5; - char tmp[tmp_l]; + size_t tmp_l = strlen(str1) + strlen(str2) + strlen(str3) + 1; + + char *tmp = (char *)malloc(tmp_l); + if (tmp == NULL) { + fprintf(stderr, "\t-> Error: failed to allocate memory for tmp, will exit\n"); + exit(1); + } + snprintf(tmp, tmp_l, "%s%s%s", str1, str2, str3); + size_t fs = fsize(tmp); - return fexists(tmp) && fs > 0; + int ret = fexists(tmp) && fs > 0; + + free(tmp); + return ret; } int fexists4(const char *str1, const char *str2, const char *str3, const char *str4) { - unsigned tmp_l = strlen(str1) + strlen(str2) + strlen(str3) + strlen(str4) + 5; - char tmp[tmp_l]; + size_t tmp_l = strlen(str1) + strlen(str2) + strlen(str3) + strlen(str4) + 1; + + char *tmp = (char *)malloc(tmp_l); + if (tmp == NULL) { + fprintf(stderr, "\t-> Error: failed to allocate memory for tmp, will exit\n"); + exit(1); + } + snprintf(tmp, tmp_l, "%s%s%s%s", str1, str2, str3, str4); - // fprintf(stderr,"\t-> checking if : %s exists\n",tmp ); + size_t fs = fsize(tmp); - return fexists(tmp) && fs > 0; -} + int ret = fexists(tmp) && fs > 0; + free(tmp); + return ret; +} // usefull little function to split char *strpop(char **str, char split) { char *tok = *str; @@ -205,7 +252,6 @@ void parse_nodes(const char *fname, int2char &rank, int2int &parent, int2intvec exit(1); } char buf[4096]; - int at = 0; char **toks = new char *[5]; while (gzgets(gz, buf, 4096)) { @@ -282,6 +328,7 @@ int SIG_COND = 1; int2int *bamRefId2tax(bam_hdr_t *hdr, char *acc2taxfile, char *bamfile, int2int &errmap, char *tempfolder, int usedump, char *filteredAcc2taxfile, char2int *acc2taxidmap) { + (void) errmap; fprintf(stderr, "\t-> Starting to extract (acc->taxid) from binary file: '%s'\n", acc2taxfile); fflush(stderr); @@ -365,8 +412,11 @@ int2int *bamRefId2tax(bam_hdr_t *hdr, char *acc2taxfile, char *bamfile, int2int nprocs++; if (fp != NULL) - assert(bgzf_write(fp, &valinbam, sizeof(int)) == sizeof(int) && - bgzf_write(fp, &val, sizeof(int)) == sizeof(int)); + if (bgzf_write(fp, &valinbam, sizeof(int)) != sizeof(int) || + bgzf_write(fp, &val, sizeof(int)) != sizeof(int)) { + fprintf(stderr, "\t-> Error: failed to write expected number of bytes with bgzf_write, will exit\n"); + exit(1); + } if (am->find(valinbam) != am->end()) fprintf(stderr, "\t-> Duplicate entries found '%s'\n", key); @@ -385,8 +435,11 @@ int2int *bamRefId2tax(bam_hdr_t *hdr, char *acc2taxfile, char *bamfile, int2int } else { int valinbam, val; while (bgzf_read(fp, &valinbam, sizeof(int))) { - assert(bgzf_read(fp, &val, sizeof(int)) == sizeof(int)); - (*am)[valinbam] = val; + if (bgzf_read(fp, &val, sizeof(int)) != sizeof(int)) { + fprintf(stderr, "\t-> Error: failed to read expected number of bytes for int value with bgzf_read, will exit\n"); + exit(1); + } + (*am)[valinbam] = val; } }