From 5e198ba9cd8cdf4b4cdce7f1011ff00aca4b2a04 Mon Sep 17 00:00:00 2001 From: Your Name Date: Sat, 25 Apr 2026 22:13:37 +0200 Subject: [PATCH 1/4] Drop getdamage .res.gz output --- R/plotdamage.R | 2 +- main_getdamage.cpp | 2 +- test/output.md5 | 2 -- test/output.md5.macos | 2 -- 4 files changed, 2 insertions(+), 6 deletions(-) diff --git a/R/plotdamage.R b/R/plotdamage.R index f58b110..b2a9683 100644 --- a/R/plotdamage.R +++ b/R/plotdamage.R @@ -38,7 +38,7 @@ print.args<-function(args,des){ # NULL is an non-optional argument, NA is an optional argument with no default, others are the default arguments args<-list(file=NULL,outfile=NA,type=1) #if no argument are given prints the need arguments and the optional ones with default -des<-list(file=" The damage file, typicaly called meta.res.gz",outfile="name of output file",type="type=1 -r 1,type=1 -r 0 in ./metaDMG-cpp getdamage [maybe not used here yet]") +des<-list(file=" The damage file, typicaly called meta.bdamage.gz",outfile="name of output file",type="type=1 -r 1,type=1 -r 0 in ./metaDMG-cpp getdamage [maybe not used here yet]") ###################################### #######get arguments and add to workspace diff --git a/main_getdamage.cpp b/main_getdamage.cpp index f67a462..0643cfb 100644 --- a/main_getdamage.cpp +++ b/main_getdamage.cpp @@ -181,7 +181,7 @@ static void write_getdamage_outputs(const getdamage_args &args, bam_hdr_t *hdr, damage *dmg) { dmg->printit(stdout, args.printLength); - dmg->write(args.onam, args.runmode == 1 ? hdr : NULL); + (void)hdr; dmg->bwrite(args.onam, args.rlens_flat_out); } diff --git a/test/output.md5 b/test/output.md5 index 52ff38c..4b3a04f 100644 --- a/test/output.md5 +++ b/test/output.md5 @@ -5,11 +5,9 @@ # 7a7363a09607cacbad8c432fdd46b26e output/test_dfit_local.dfit.fix 783c885b2ee8ac6469ab070b972e73d8 output/test_getdamage_global.bdamage 3398e353af15d1646228cb8bbaf4a5a8 output/test_getdamage_global.bdamage.tsv -f0882ca11587cea63aa120d1a64c8d14 output/test_getdamage_global.res ba1415f0e1889f819b430651dc56b4ff output/test_getdamage_global.stat 71308047bbc8f703fdd65c798de91cfa output/test_getdamage_local.bdamage 93ec5114d069f984ccaaa5079c15db1e output/test_getdamage_local.bdamage.tsv -d95e54a3886b9f411bde8d96de4af32e output/test_getdamage_local.res 5f51536e977e54529cb2971462561315 output/test_getdamage_local.stat 443b35425bc868259fde18dbe86011e3 output/test_lca.bdamage # d4532cede0b424946f92559b2e737899 output/test_lca.lca##this has been commented out because eithput is random between reqad and the reverse complement diff --git a/test/output.md5.macos b/test/output.md5.macos index 82fdce6..f2d2253 100644 --- a/test/output.md5.macos +++ b/test/output.md5.macos @@ -3,11 +3,9 @@ 197881fb3261dbc475146a8b896b2118 output/test_dfit_global.dfit.fix 783c885b2ee8ac6469ab070b972e73d8 output/test_getdamage_global.bdamage 3398e353af15d1646228cb8bbaf4a5a8 output/test_getdamage_global.bdamage.tsv -f0882ca11587cea63aa120d1a64c8d14 output/test_getdamage_global.res ba1415f0e1889f819b430651dc56b4ff output/test_getdamage_global.stat 71308047bbc8f703fdd65c798de91cfa output/test_getdamage_local.bdamage 93ec5114d069f984ccaaa5079c15db1e output/test_getdamage_local.bdamage.tsv -d95e54a3886b9f411bde8d96de4af32e output/test_getdamage_local.res 5f51536e977e54529cb2971462561315 output/test_getdamage_local.stat 443b35425bc868259fde18dbe86011e3 output/test_lca.bdamage # d4532cede0b424946f92559b2e737899 output/test_lca.lca##this has been commented out because eithput is random between reqad and the reverse complement From 2414a5d0ac96e740f16a871488dfd381d36672d7 Mon Sep 17 00:00:00 2001 From: Your Name Date: Sat, 25 Apr 2026 23:31:29 +0200 Subject: [PATCH 2/4] Add filter_bdamage command with resolve and companion filtering --- README.md | 28 ++ main_filter_bdamage.cpp | 826 ++++++++++++++++++++++++++++++++++++++++ main_filter_bdamage.h | 3 + metaDMG.cpp | 4 + 4 files changed, 861 insertions(+) create mode 100644 main_filter_bdamage.cpp create mode 100644 main_filter_bdamage.h diff --git a/README.md b/README.md index 5b269d0..9d020de 100644 --- a/README.md +++ b/README.md @@ -70,6 +70,34 @@ Top-level and subcommand help: ./metaDMG-cpp -h ``` +## Filter existing bdamage output +`filter_bdamage` creates a subset of an existing `.bdamage.gz` file and (when available) matching `.stat.gz` and `.rlens.gz` files. + +Important ID semantics: +- For `lca` output, IDs are taxids. +- For `getdamage --run_mode 1`, IDs are BAM reference offsets (tid), not taxids. + +Basic usage: +``` +./metaDMG-cpp filter_bdamage in.bdamage.gz --id 11 --out_prefix subset +``` + +Taxonomic expansion with descendants (from nodes): +``` +./metaDMG-cpp filter_bdamage in.bdamage.gz --id 11 --out_prefix subset --nodes nodes.dmp.gz +``` +This includes taxid `11` and all descendants of `11` from `nodes.dmp.gz`. + +Resolve by reference/accession: +``` +./metaDMG-cpp filter_bdamage in.bdamage.gz --resolve ref2 --bam reads.bam --out_prefix subset +./metaDMG-cpp filter_bdamage in.bdamage.gz --resolve ACC123 --acc2tax acc2taxid.map.gz --out_prefix subset +``` + +Companion files (`stat` and `rlens`): +- Auto-detected by prefix when omitted: `X.bdamage.gz` -> `X.stat.gz` and `X.rlens.gz`. +- Can be set explicitly: `--stat file.stat.gz --rlens file.rlens.gz`. + ## Damage analysis (non-taxonomically assigned) metaDMG-cpp calculates substitutions between reads and reference sequences. It can operate in two modes: diff --git a/main_filter_bdamage.cpp b/main_filter_bdamage.cpp new file mode 100644 index 0000000..ec13904 --- /dev/null +++ b/main_filter_bdamage.cpp @@ -0,0 +1,826 @@ +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "main_filter_bdamage.h" +#include "shared.h" + +struct filter_bdamage_args { + char *infile; + char *out_prefix; + char *id_file; + char *resolve_file; + char *bam_file; + char *acc2tax_file; + char *nodes_file; + char *stat_infile; + char *rlens_infile; + int min_reads; + int max_reads; + int nthreads; + int exclude; + int strict_resolve; + std::set ids; + std::vector selectors; +}; + +static int usage_filter_bdamage(FILE *fp) { + fprintf(fp, "Usage: ./metaDMG-cpp filter_bdamage [options]\n"); + fprintf(fp, "\nOptions:\n"); + fprintf(fp, " -o/--out_prefix STR output prefix (default: filtered)\n"); + fprintf(fp, " -n/--threads INT BGZF read/write threads (default: 4)\n"); + fprintf(fp, " --id INT include this numeric id (repeatable)\n"); + fprintf(fp, " --id_list STR comma-separated numeric ids/selectors\n"); + fprintf(fp, " --id_file FILE newline separated ids/selectors (first column used)\n"); + fprintf(fp, " --resolve STR selector to resolve (repeatable)\n"); + fprintf(fp, " --resolve_list STR comma-separated selectors to resolve\n"); + fprintf(fp, " --resolve_file FILE newline separated selectors (first column used)\n"); + fprintf(fp, " --bam FILE BAM/CRAM/SAM for refname -> bam-offset (tid) resolve\n"); + fprintf(fp, " --acc2tax FILE accession2taxid map for accession -> taxid resolve\n"); + fprintf(fp, " --nodes FILE expand selected taxid to all descendants from nodes file\n"); + fprintf(fp, " --stat FILE companion stat.gz to filter (auto-inferred if omitted)\n"); + fprintf(fp, " --rlens FILE companion rlens.gz to filter (auto-inferred if omitted)\n"); + fprintf(fp, " --strict_resolve INT 1: fail on unresolved selectors (default), 0: ignore\n"); + fprintf(fp, " --exclude exclude listed ids/selectors instead of including\n"); + fprintf(fp, " --min_reads INT keep entries with nreads >= INT\n"); + fprintf(fp, " --max_reads INT keep entries with nreads <= INT\n"); + fprintf(fp, " -h/--help show this help\n"); + fprintf(fp, "\nOutput:\n"); + fprintf(fp, " .bdamage.gz\n"); + fprintf(fp, " .stat.gz (if companion found/specified)\n"); + fprintf(fp, " .rlens.gz (if companion found/specified)\n"); + fprintf(fp, "\nExamples:\n"); + fprintf(fp, " ./metaDMG-cpp filter_bdamage in.bdamage.gz --id 11 --out_prefix subset\n"); + fprintf(fp, " ./metaDMG-cpp filter_bdamage in.bdamage.gz --resolve ref2 --bam reads.bam --out_prefix subset\n"); + fprintf(fp, " ./metaDMG-cpp filter_bdamage in.bdamage.gz --resolve ACC123 --acc2tax acc2taxid.map.gz --nodes nodes.dmp.gz --out_prefix subset\n"); + fprintf(fp, " ./metaDMG-cpp filter_bdamage in.bdamage.gz --id_list 2,10,42 --exclude --min_reads 100 --out_prefix subset\n"); + return 0; +} + +static filter_bdamage_args init_filter_bdamage_args() { + filter_bdamage_args args; + args.infile = NULL; + args.out_prefix = strdup("filtered"); + args.id_file = NULL; + args.resolve_file = NULL; + args.bam_file = NULL; + args.acc2tax_file = NULL; + args.nodes_file = NULL; + args.stat_infile = NULL; + args.rlens_infile = NULL; + args.min_reads = -1; + args.max_reads = -1; + args.nthreads = 4; + args.exclude = 0; + args.strict_resolve = 1; + return args; +} + +static int parse_int(const char *txt, int *out) { + if (txt == NULL || out == NULL) + return 1; + + errno = 0; + char *end = NULL; + long val = strtol(txt, &end, 10); + if (end == txt || *end != '\0' || errno != 0 || val < INT_MIN || val > INT_MAX) + return 1; + *out = (int)val; + return 0; +} + +static int has_suffix(const std::string &txt, const std::string &suffix) { + if (txt.size() < suffix.size()) + return 0; + return txt.compare(txt.size() - suffix.size(), suffix.size(), suffix) == 0; +} + +static void add_selector_token(const char *tok, filter_bdamage_args &args) { + if (tok == NULL || tok[0] == '\0') + return; + + int id = 0; + if (parse_int(tok, &id) == 0) + args.ids.insert(id); + else + args.selectors.push_back(std::string(tok)); +} + +static int add_csv_tokens(const char *csv, filter_bdamage_args &args) { + if (csv == NULL) + return 1; + + char *tmp = strdup(csv); + char *token = strtok(tmp, ","); + while (token != NULL) { + add_selector_token(token, args); + token = strtok(NULL, ","); + } + free(tmp); + return 0; +} + +static int load_tokens_from_file(const char *path, filter_bdamage_args &args) { + if (path == NULL) + return 0; + + FILE *fp = fopen(path, "r"); + if (fp == NULL) { + fprintf(stderr, "\t-> Error: could not open file: %s\n", path); + return 1; + } + + char line[4096]; + while (fgets(line, sizeof(line), fp) != NULL) { + char *ptr = line; + while (*ptr == ' ' || *ptr == '\t') + ptr++; + if (*ptr == '\0' || *ptr == '\n' || *ptr == '#') + continue; + + char *tok = strtok(ptr, "\t\n "); + if (tok == NULL) + continue; + + add_selector_token(tok, args); + } + + fclose(fp); + return 0; +} + +static int parse_filter_bdamage_args(int argc, char **argv, filter_bdamage_args &args) { + for (int i = 1; i < argc; i++) { + if (!strcasecmp(argv[i], "-h") || !strcasecmp(argv[i], "--help")) { + usage_filter_bdamage(stdout); + return 2; + } else if (!strcasecmp(argv[i], "-o") || !strcasecmp(argv[i], "--out_prefix")) { + if (i + 1 >= argc) { + fprintf(stderr, "\t-> Error: missing value for %s\n", argv[i]); + return 1; + } + free(args.out_prefix); + args.out_prefix = strdup(argv[++i]); + } else if (!strcasecmp(argv[i], "-n") || !strcasecmp(argv[i], "--threads")) { + if (i + 1 >= argc) { + fprintf(stderr, "\t-> Error: missing value for %s\n", argv[i]); + return 1; + } + args.nthreads = atoi(argv[++i]); + } else if (!strcasecmp(argv[i], "--id")) { + if (i + 1 >= argc) { + fprintf(stderr, "\t-> Error: missing value for --id\n"); + return 1; + } + int id = 0; + if (parse_int(argv[++i], &id) != 0) { + fprintf(stderr, "\t-> Error: invalid value for --id: %s\n", argv[i]); + return 1; + } + args.ids.insert(id); + } else if (!strcasecmp(argv[i], "--id_list")) { + if (i + 1 >= argc) { + fprintf(stderr, "\t-> Error: missing value for --id_list\n"); + return 1; + } + if (add_csv_tokens(argv[++i], args) != 0) + return 1; + } else if (!strcasecmp(argv[i], "--id_file")) { + if (i + 1 >= argc) { + fprintf(stderr, "\t-> Error: missing value for --id_file\n"); + return 1; + } + free(args.id_file); + args.id_file = strdup(argv[++i]); + } else if (!strcasecmp(argv[i], "--resolve")) { + if (i + 1 >= argc) { + fprintf(stderr, "\t-> Error: missing value for --resolve\n"); + return 1; + } + add_selector_token(argv[++i], args); + } else if (!strcasecmp(argv[i], "--resolve_list")) { + if (i + 1 >= argc) { + fprintf(stderr, "\t-> Error: missing value for --resolve_list\n"); + return 1; + } + if (add_csv_tokens(argv[++i], args) != 0) + return 1; + } else if (!strcasecmp(argv[i], "--resolve_file")) { + if (i + 1 >= argc) { + fprintf(stderr, "\t-> Error: missing value for --resolve_file\n"); + return 1; + } + free(args.resolve_file); + args.resolve_file = strdup(argv[++i]); + } else if (!strcasecmp(argv[i], "--bam")) { + if (i + 1 >= argc) { + fprintf(stderr, "\t-> Error: missing value for --bam\n"); + return 1; + } + free(args.bam_file); + args.bam_file = strdup(argv[++i]); + } else if (!strcasecmp(argv[i], "--acc2tax")) { + if (i + 1 >= argc) { + fprintf(stderr, "\t-> Error: missing value for --acc2tax\n"); + return 1; + } + free(args.acc2tax_file); + args.acc2tax_file = strdup(argv[++i]); + } else if (!strcasecmp(argv[i], "--nodes")) { + if (i + 1 >= argc) { + fprintf(stderr, "\t-> Error: missing value for --nodes\n"); + return 1; + } + free(args.nodes_file); + args.nodes_file = strdup(argv[++i]); + } else if (!strcasecmp(argv[i], "--stat")) { + if (i + 1 >= argc) { + fprintf(stderr, "\t-> Error: missing value for --stat\n"); + return 1; + } + free(args.stat_infile); + args.stat_infile = strdup(argv[++i]); + } else if (!strcasecmp(argv[i], "--rlens")) { + if (i + 1 >= argc) { + fprintf(stderr, "\t-> Error: missing value for --rlens\n"); + return 1; + } + free(args.rlens_infile); + args.rlens_infile = strdup(argv[++i]); + } else if (!strcasecmp(argv[i], "--strict_resolve")) { + if (i + 1 >= argc) { + fprintf(stderr, "\t-> Error: missing value for --strict_resolve\n"); + return 1; + } + args.strict_resolve = atoi(argv[++i]); + } else if (!strcasecmp(argv[i], "--exclude")) { + args.exclude = 1; + } else if (!strcasecmp(argv[i], "--min_reads")) { + if (i + 1 >= argc) { + fprintf(stderr, "\t-> Error: missing value for --min_reads\n"); + return 1; + } + args.min_reads = atoi(argv[++i]); + } else if (!strcasecmp(argv[i], "--max_reads")) { + if (i + 1 >= argc) { + fprintf(stderr, "\t-> Error: missing value for --max_reads\n"); + return 1; + } + args.max_reads = atoi(argv[++i]); + } else if (argv[i][0] == '-') { + fprintf(stderr, "\t-> Error: unknown option: %s\n", argv[i]); + return 1; + } else { + if (args.infile == NULL) + args.infile = strdup(argv[i]); + else { + fprintf(stderr, "\t-> Error: multiple input files provided: %s and %s\n", args.infile, argv[i]); + return 1; + } + } + } + + if (args.infile == NULL) { + fprintf(stderr, "\t-> Error: missing input .bdamage.gz file\n"); + return 1; + } + + if (args.nthreads < 1) { + fprintf(stderr, "\t-> Error: --threads must be >= 1\n"); + return 1; + } + if (args.min_reads < -1) { + fprintf(stderr, "\t-> Error: --min_reads must be -1 or greater\n"); + return 1; + } + if (args.max_reads < -1) { + fprintf(stderr, "\t-> Error: --max_reads must be -1 or greater\n"); + return 1; + } + if (args.min_reads != -1 && args.max_reads != -1 && args.min_reads > args.max_reads) { + fprintf(stderr, "\t-> Error: --min_reads cannot be greater than --max_reads\n"); + return 1; + } + if (!(args.strict_resolve == 0 || args.strict_resolve == 1)) { + fprintf(stderr, "\t-> Error: --strict_resolve must be 0 or 1\n"); + return 1; + } + + if (load_tokens_from_file(args.id_file, args) != 0) + return 1; + if (load_tokens_from_file(args.resolve_file, args) != 0) + return 1; + + return 0; +} + +static int maybe_infer_companion_inputs(filter_bdamage_args &args) { + if (args.infile == NULL) + return 0; + + std::string in(args.infile); + std::string base; + if (has_suffix(in, ".bdamage.gz")) + base = in.substr(0, in.size() - strlen(".bdamage.gz")); + else if (has_suffix(in, ".bdamage")) + base = in.substr(0, in.size() - strlen(".bdamage")); + else + return 0; + + if (args.stat_infile == NULL) { + std::string s = base + ".stat.gz"; + if (fexists(s.c_str())) + args.stat_infile = strdup(s.c_str()); + } + if (args.rlens_infile == NULL) { + std::string s = base + ".rlens.gz"; + if (fexists(s.c_str())) + args.rlens_infile = strdup(s.c_str()); + } + return 0; +} + +static int load_bam_lookup(const char *bam_file, std::map &lookup) { + if (bam_file == NULL) + return 0; + + samFile *samfp = sam_open_format(bam_file, "r", NULL); + if (samfp == NULL) { + fprintf(stderr, "\t-> Error: could not open BAM/SAM/CRAM for resolve: %s\n", bam_file); + return 1; + } + bam_hdr_t *hdr = sam_hdr_read(samfp); + if (hdr == NULL) { + fprintf(stderr, "\t-> Error: could not read header for resolve: %s\n", bam_file); + sam_close(samfp); + return 1; + } + + for (int i = 0; i < hdr->n_targets; i++) + lookup[std::string(hdr->target_name[i])] = i; + + bam_hdr_destroy(hdr); + sam_close(samfp); + return 0; +} + +static int load_acc2tax_lookup(const char *acc2tax_file, std::map &lookup) { + if (acc2tax_file == NULL) + return 0; + + gzFile fp = gzopen(acc2tax_file, "rb"); + if (fp == Z_NULL) { + fprintf(stderr, "\t-> Error: could not open acc2tax for resolve: %s\n", acc2tax_file); + return 1; + } + + char line[8192]; + int at = 0; + while (gzgets(fp, line, sizeof(line)) != NULL) { + at++; + if (at == 1) + continue; + + char *tok = strtok(line, "\t\n "); + if (tok == NULL) + continue; + char *key = strtok(NULL, "\t\n "); + if (key == NULL) + continue; + char *val = strtok(NULL, "\t\n "); + if (val == NULL) + continue; + + lookup[std::string(key)] = atoi(val); + } + + gzclose(fp); + return 0; +} + +static int resolve_selectors(filter_bdamage_args &args) { + if (args.selectors.empty()) + return 0; + + std::map lookup; + if (load_bam_lookup(args.bam_file, lookup) != 0) + return 1; + if (load_acc2tax_lookup(args.acc2tax_file, lookup) != 0) + return 1; + + int unresolved = 0; + for (size_t i = 0; i < args.selectors.size(); i++) { + const std::string &sel = args.selectors[i]; + std::map::iterator it = lookup.find(sel); + if (it == lookup.end()) { + unresolved++; + fprintf(stderr, "\t-> resolve: selector not found: %s\n", sel.c_str()); + continue; + } + args.ids.insert(it->second); + } + + if (unresolved > 0 && args.strict_resolve == 1) { + fprintf(stderr, "\t-> Error: %d selector(s) could not be resolved\n", unresolved); + return 1; + } + + if (unresolved > 0) + fprintf(stderr, "\t-> resolve: ignored unresolved selectors because --strict_resolve 0\n"); + + return 0; +} + +static void expand_descendants_from_seed(int seed, int2intvec &child, std::set &ids) { + std::vector stack; + stack.push_back(seed); + + while (!stack.empty()) { + int cur = stack.back(); + stack.pop_back(); + + int2intvec::iterator it = child.find(cur); + if (it == child.end()) + continue; + + for (size_t i = 0; i < it->second.size(); i++) { + const int kid = it->second[i]; + if (ids.insert(kid).second) + stack.push_back(kid); + } + } +} + +static int maybe_expand_ids_from_nodes(filter_bdamage_args &args) { + if (args.nodes_file == NULL || args.ids.empty()) + return 0; + + int2char rank; + int2int parent; + int2intvec child; + parse_nodes(args.nodes_file, rank, parent, child, 1); + + const size_t before = args.ids.size(); + std::vector seeds(args.ids.begin(), args.ids.end()); + for (size_t i = 0; i < seeds.size(); i++) + expand_descendants_from_seed(seeds[i], child, args.ids); + + for (int2char::iterator it = rank.begin(); it != rank.end(); it++) + free(it->second); + + fprintf(stderr, + "\t-> nodes expansion with %s: %zu seed ids -> %zu ids after descendants\n", + args.nodes_file, + before, + args.ids.size()); + return 0; +} + +static int keep_entry(int id, int nreads, const filter_bdamage_args &args) { + int keep = 1; + + if (!args.ids.empty()) { + const int in_set = args.ids.find(id) != args.ids.end(); + keep = args.exclude ? !in_set : in_set; + } + + if (keep && args.min_reads != -1 && nreads < args.min_reads) + keep = 0; + if (keep && args.max_reads != -1 && nreads > args.max_reads) + keep = 0; + + return keep; +} + +static int filter_rlens_companion(const char *infile, const char *outfile, const std::set &kept_ids) { + if (infile == NULL || outfile == NULL) + return 0; + + gzFile in = gzopen(infile, "rb"); + if (in == Z_NULL) { + fprintf(stderr, "\t-> Error: could not open rlens companion: %s\n", infile); + return 1; + } + gzFile out = gzopen(outfile, "wb"); + if (out == Z_NULL) { + fprintf(stderr, "\t-> Error: could not open rlens output: %s\n", outfile); + gzclose(in); + return 1; + } + + char line[8192]; + size_t at = 0; + size_t kept = 0; + while (gzgets(in, line, sizeof(line)) != NULL) { + at++; + if (at == 1) { + gzputs(out, line); + continue; + } + + char copy[8192]; + strncpy(copy, line, sizeof(copy)); + copy[sizeof(copy) - 1] = '\0'; + char *tok = strtok(copy, "\t\n "); + if (tok == NULL) + continue; + + int id = 0; + if (parse_int(tok, &id) != 0) + continue; + + if (kept_ids.find(id) != kept_ids.end()) { + gzputs(out, line); + kept++; + } + } + + gzclose(in); + gzclose(out); + + fprintf(stderr, "\t-> Companion rlens filtered: %s -> %s (kept lines: %zu)\n", infile, outfile, kept); + return 0; +} + +static int resolve_stat_id(const char *tok, + const std::set &kept_ids, + const std::map &bam_lookup, + int has_bam_lookup, + int &is_keep, + int &is_resolved) { + is_keep = 0; + is_resolved = 0; + + if (tok == NULL) + return 0; + + int id = 0; + if (parse_int(tok, &id) == 0) { + is_resolved = 1; + is_keep = kept_ids.find(id) != kept_ids.end(); + return 0; + } + + if (strcmp(tok, "global") == 0) { + is_resolved = 1; + is_keep = kept_ids.find(0) != kept_ids.end(); + return 0; + } + + if (has_bam_lookup) { + std::map::const_iterator it = bam_lookup.find(tok); + if (it != bam_lookup.end()) { + is_resolved = 1; + is_keep = kept_ids.find(it->second) != kept_ids.end(); + return 0; + } + } + + return 0; +} + +static int filter_stat_companion(const char *infile, + const char *outfile, + const std::set &kept_ids, + const std::map &bam_lookup, + int has_bam_lookup) { + if (infile == NULL || outfile == NULL) + return 0; + + gzFile in = gzopen(infile, "rb"); + if (in == Z_NULL) { + fprintf(stderr, "\t-> Error: could not open stat companion: %s\n", infile); + return 1; + } + gzFile out = gzopen(outfile, "wb"); + if (out == Z_NULL) { + fprintf(stderr, "\t-> Error: could not open stat output: %s\n", outfile); + gzclose(in); + return 1; + } + + char line[8192]; + size_t at = 0; + size_t kept = 0; + size_t unresolved = 0; + + while (gzgets(in, line, sizeof(line)) != NULL) { + at++; + if (at == 1) { + gzputs(out, line); + continue; + } + + char copy[8192]; + strncpy(copy, line, sizeof(copy)); + copy[sizeof(copy) - 1] = '\0'; + char *tok = strtok(copy, "\t\n "); + if (tok == NULL) + continue; + + int is_keep = 0; + int is_resolved = 0; + resolve_stat_id(tok, kept_ids, bam_lookup, has_bam_lookup, is_keep, is_resolved); + if (!is_resolved) { + unresolved++; + continue; + } + if (is_keep) { + gzputs(out, line); + kept++; + } + } + + gzclose(in); + gzclose(out); + + fprintf(stderr, "\t-> Companion stat filtered: %s -> %s (kept lines: %zu)\n", infile, outfile, kept); + if (unresolved > 0 && !has_bam_lookup) + fprintf(stderr, "\t-> Companion stat: %zu unresolved id labels (tip: provide --bam for refname->tid mapping)\n", unresolved); + else if (unresolved > 0) + fprintf(stderr, "\t-> Companion stat: %zu unresolved lines were skipped\n", unresolved); + + return 0; +} + +int main_filter_bdamage(int argc, char **argv) { + if (argc == 1 || (argc == 2 && (!strcasecmp(argv[1], "-h") || !strcasecmp(argv[1], "--help")))) + return usage_filter_bdamage(stderr); + + filter_bdamage_args args = init_filter_bdamage_args(); + const int parse_rc = parse_filter_bdamage_args(argc, argv, args); + if (parse_rc == 2) + return 0; + if (parse_rc != 0) + return 1; + + if (maybe_infer_companion_inputs(args) != 0) + return 1; + if (resolve_selectors(args) != 0) + return 1; + if (maybe_expand_ids_from_nodes(args) != 0) + return 1; + + char out_bdamage[4096]; + char out_stat[4096]; + char out_rlens[4096]; + snprintf(out_bdamage, sizeof(out_bdamage), "%s.bdamage.gz", args.out_prefix); + snprintf(out_stat, sizeof(out_stat), "%s.stat.gz", args.out_prefix); + snprintf(out_rlens, sizeof(out_rlens), "%s.rlens.gz", args.out_prefix); + + fprintf(stderr, + "\t-> filter_bdamage infile: %s out_prefix: %s ids: %zu selectors: %zu exclude: %d min_reads: %d max_reads: %d nthreads: %d\n", + args.infile, + args.out_prefix, + args.ids.size(), + args.selectors.size(), + args.exclude, + args.min_reads, + args.max_reads, + args.nthreads); + fprintf(stderr, + "\t-> companions stat:%s rlens:%s\n", + args.stat_infile ? args.stat_infile : "(none)", + args.rlens_infile ? args.rlens_infile : "(none)"); + + BGZF *in = bgzf_open(args.infile, "r"); + if (in == NULL) { + fprintf(stderr, "\t-> Error: failed to open input file: %s\n", args.infile); + return 1; + } + BGZF *out = bgzf_open(out_bdamage, "w"); + if (out == NULL) { + fprintf(stderr, "\t-> Error: failed to open output file: %s\n", out_bdamage); + bgzf_close(in); + return 1; + } + + if (args.nthreads > 1) { + bgzf_mt(in, args.nthreads, 256); + bgzf_mt(out, args.nthreads, 256); + } + + int printlength = 0; + if (bgzf_read(in, &printlength, sizeof(int)) != sizeof(int)) { + fprintf(stderr, "\t-> Error: failed to read printlength from input bdamage file\n"); + bgzf_close(in); + bgzf_close(out); + return 1; + } + if (printlength <= 0) { + fprintf(stderr, "\t-> Error: invalid printlength in input bdamage file: %d\n", printlength); + bgzf_close(in); + bgzf_close(out); + return 1; + } + + if (bgzf_write(out, &printlength, sizeof(int)) != sizeof(int)) { + fprintf(stderr, "\t-> Error: failed to write printlength to output bdamage file\n"); + bgzf_close(in); + bgzf_close(out); + return 1; + } + + const size_t nvals = (size_t)printlength * 16 * 2; + const size_t datab = nvals * sizeof(float); + std::vector data(nvals); + + std::set kept_ids; + size_t total = 0; + size_t kept = 0; + while (1) { + int ref_nreads[2]; + const int nread = bgzf_read(in, ref_nreads, 2 * sizeof(int)); + if (nread == 0) + break; + if (nread != 2 * (int)sizeof(int)) { + fprintf(stderr, "\t-> Error: truncated/corrupt bdamage file (record header)\n"); + bgzf_close(in); + bgzf_close(out); + return 1; + } + + if (bgzf_read(in, data.data(), datab) != (ssize_t)datab) { + fprintf(stderr, "\t-> Error: truncated/corrupt bdamage file (record payload)\n"); + bgzf_close(in); + bgzf_close(out); + return 1; + } + + total++; + if (!keep_entry(ref_nreads[0], ref_nreads[1], args)) + continue; + + if (bgzf_write(out, ref_nreads, 2 * sizeof(int)) != 2 * sizeof(int)) { + fprintf(stderr, "\t-> Error: failed writing record header to output bdamage file\n"); + bgzf_close(in); + bgzf_close(out); + return 1; + } + if (bgzf_write(out, data.data(), datab) != (ssize_t)datab) { + fprintf(stderr, "\t-> Error: failed writing record payload to output bdamage file\n"); + bgzf_close(in); + bgzf_close(out); + return 1; + } + kept_ids.insert(ref_nreads[0]); + kept++; + } + + bgzf_close(in); + bgzf_close(out); + + std::map bam_lookup; + int has_bam_lookup = 0; + if (args.bam_file != NULL) { + if (load_bam_lookup(args.bam_file, bam_lookup) != 0) + return 1; + has_bam_lookup = 1; + } + + if (args.stat_infile != NULL) { + if (filter_stat_companion(args.stat_infile, out_stat, kept_ids, bam_lookup, has_bam_lookup) != 0) + return 1; + } + if (args.rlens_infile != NULL) { + if (filter_rlens_companion(args.rlens_infile, out_rlens, kept_ids) != 0) + return 1; + } + + fprintf(stderr, + "\t-> filter_bdamage done. Input entries: %zu Kept entries: %zu Removed entries: %zu Output: %s\n", + total, + kept, + total - kept, + out_bdamage); + + free(args.infile); + free(args.out_prefix); + if (args.id_file) + free(args.id_file); + if (args.resolve_file) + free(args.resolve_file); + if (args.bam_file) + free(args.bam_file); + if (args.acc2tax_file) + free(args.acc2tax_file); + if (args.nodes_file) + free(args.nodes_file); + if (args.stat_infile) + free(args.stat_infile); + if (args.rlens_infile) + free(args.rlens_infile); + + return 0; +} diff --git a/main_filter_bdamage.h b/main_filter_bdamage.h new file mode 100644 index 0000000..bbf5320 --- /dev/null +++ b/main_filter_bdamage.h @@ -0,0 +1,3 @@ +#pragma once + +int main_filter_bdamage(int argc, char **argv); diff --git a/metaDMG.cpp b/metaDMG.cpp index b416295..11700eb 100644 --- a/metaDMG.cpp +++ b/metaDMG.cpp @@ -6,6 +6,7 @@ #include "Aggregate_stat.h" // for main_aggregate #include "main_dfit.h" // for main_dfit +#include "main_filter_bdamage.h" #include "main_getdamage.h" #include "main_pmd.h" // for main_pmd #include "main_print.h" @@ -32,6 +33,7 @@ static void print_top_level_usage(FILE *fp) { fprintf(fp, "./metaDMG-cpp print_ugly [many options] bdamage.gz\n"); fprintf(fp, "./metaDMG-cpp dfit [many options] bdamage.gz\n"); fprintf(fp, "./metaDMG-cpp aggregate [many options] bdamage.gz\n"); + fprintf(fp, "./metaDMG-cpp filter_bdamage [options] in.bdamage.gz\n"); fprintf(fp, "\nUse './metaDMG-cpp --help' for command-specific help.\n"); } @@ -75,6 +77,8 @@ int main(int argc, char **argv) { rc = main_dfit(argc, argv); else if (!strcmp(argv[0], "aggregate")) rc = main_aggregate(argc, argv); + else if (!strcmp(argv[0], "filter_bdamage")) + rc = main_filter_bdamage(argc, argv); else if (!strcmp(argv[0], "print2")) rc = main_print2(argc, argv); else if (!strcmp(argv[0], "lca")) From ee4923f6b8fb4653fc4754466180706bbb97f4a5 Mon Sep 17 00:00:00 2001 From: Your Name Date: Sun, 26 Apr 2026 00:29:44 +0200 Subject: [PATCH 3/4] Fix companion autodetect and add filter_bdamage tests --- main_filter_bdamage.cpp | 26 ++++++++++++++++----- test/testAll2.sh | 51 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 71 insertions(+), 6 deletions(-) diff --git a/main_filter_bdamage.cpp b/main_filter_bdamage.cpp index ec13904..047034c 100644 --- a/main_filter_bdamage.cpp +++ b/main_filter_bdamage.cpp @@ -339,15 +339,29 @@ static int maybe_infer_companion_inputs(filter_bdamage_args &args) { return 0; if (args.stat_infile == NULL) { - std::string s = base + ".stat.gz"; - if (fexists(s.c_str())) - args.stat_infile = strdup(s.c_str()); + std::vector candidates; + candidates.push_back(base + ".stat.gz"); + candidates.push_back(base + ".stat"); + for (size_t i = 0; i < candidates.size(); i++) { + if (fexists(candidates[i].c_str())) { + args.stat_infile = strdup(candidates[i].c_str()); + break; + } + } } + if (args.rlens_infile == NULL) { - std::string s = base + ".rlens.gz"; - if (fexists(s.c_str())) - args.rlens_infile = strdup(s.c_str()); + std::vector candidates; + candidates.push_back(base + ".rlens.gz"); + candidates.push_back(base + ".rlens"); + for (size_t i = 0; i < candidates.size(); i++) { + if (fexists(candidates[i].c_str())) { + args.rlens_infile = strdup(candidates[i].c_str()); + break; + } + } } + return 0; } diff --git a/test/testAll2.sh b/test/testAll2.sh index 9cbc4d2..7cd2d15 100755 --- a/test/testAll2.sh +++ b/test/testAll2.sh @@ -103,6 +103,14 @@ assert_gzip_contains() { fi } +assert_gzip_not_contains() { + local file="$1" + local pattern="$2" + + if gunzip -c "${file}" | grep -F -- "${pattern}" >/dev/null; then + mark_fail "Unexpected pattern found in ${file}: ${pattern}" + fi +} test_getdamage() { run_logged "Running getdamage global" \ "${PRG}" getdamage --run_mode 0 --min_length 35 --print_length 5 \ @@ -262,6 +270,48 @@ test_data2_getdamage() { assert_gzip_contains output_data2_gd/sam5.rlens.gz $'0\t607:1' } + +test_filter_bdamage() { + run_logged "Running filter_bdamage nodes expansion (lca taxid)" \ + "${PRG}" filter_bdamage output_data2/sam3.bdamage.gz \ + --id 10 --nodes data2/nodes.dmp --out_prefix output_data2/filter_nodes + + require_file output_data2/filter_nodes.bdamage.gz || return 0 + require_file output_data2/filter_nodes.stat.gz || return 0 + require_file output_data2/filter_nodes.rlens.gz || return 0 + + assert_gzip_contains output_data2/filter_nodes.stat.gz $'11\t1\t30.000000\t0.000000\t0.000000\t0.000000\t"l__Bacteria"\t"subspecies"' + assert_gzip_not_contains output_data2/filter_nodes.stat.gz $'10\t1\t30.000000\t0.000000\t0.000000\t0.000000\t"l__Bacteria"\t"species"' + + if ! "${PRG}" print output_data2/filter_nodes.bdamage.gz \ + 1>output_data2/filter_nodes.bdamage.tsv 2>>"${LOG}"; then + mark_fail "Problem running print on output_data2/filter_nodes.bdamage.gz" + fi + assert_file_line_count output_data2/filter_nodes.bdamage.tsv 61 + + run_logged "Running filter_bdamage resolve with BAM (getdamage local)" \ + "${PRG}" filter_bdamage output_data2_gd/sam2.bdamage.gz \ + --resolve ref2 --bam data2/sam2.sam --out_prefix output_data2_gd/filter_ref2 + + require_file output_data2_gd/filter_ref2.bdamage.gz || return 0 + require_file output_data2_gd/filter_ref2.stat.gz || return 0 + require_file output_data2_gd/filter_ref2.rlens.gz || return 0 + + assert_gzip_contains output_data2_gd/filter_ref2.stat.gz $'ref2\t1\t30.000000\t0.000000\t0.000000\t0.000000\tNA\tNA' + assert_gzip_not_contains output_data2_gd/filter_ref2.stat.gz $'ref1\t1\t30.000000\t0.000000\t0.000000\t0.000000\tNA\tNA' + assert_gzip_not_contains output_data2_gd/filter_ref2.stat.gz $'ref3\t1\t30.000000\t0.000000\t0.000000\t0.000000\tNA\tNA' + + assert_gzip_contains output_data2_gd/filter_ref2.rlens.gz $'1\t30:1' + assert_gzip_not_contains output_data2_gd/filter_ref2.rlens.gz $'0\t30:1' + assert_gzip_not_contains output_data2_gd/filter_ref2.rlens.gz $'2\t30:1' + + if ! "${PRG}" print output_data2_gd/filter_ref2.bdamage.gz \ + 1>output_data2_gd/filter_ref2.bdamage.tsv 2>>"${LOG}"; then + mark_fail "Problem running print on output_data2_gd/filter_ref2.bdamage.gz" + fi + assert_file_line_count output_data2_gd/filter_ref2.bdamage.tsv 11 +} + test_compressbam() { local compressbam="../misc/compressbam" local header_all="output_compressbam/basic.all.header.sam" @@ -455,6 +505,7 @@ main() { test_prints test_data2 test_data2_getdamage + test_filter_bdamage test_compressbam test_extract_reads validate_checksums From cf30da14565de715b3f1d7b2b2be2eae7a0fdeed Mon Sep 17 00:00:00 2001 From: Your Name Date: Sun, 26 Apr 2026 00:50:45 +0200 Subject: [PATCH 4/4] Remove print2 command and improve print help text --- main_print.cpp | 308 +++---------------------------------------------- main_print.h | 1 - metaDMG.cpp | 3 - 3 files changed, 18 insertions(+), 294 deletions(-) diff --git a/main_print.cpp b/main_print.cpp index 68001e8..7b07202 100644 --- a/main_print.cpp +++ b/main_print.cpp @@ -160,14 +160,23 @@ mydata2 getval_stats(std::map &retmap, int2intvec &child, int taxi } static int usage_print(FILE *fp) { - fprintf(fp, "Usage: ./metaDMG-cpp print file.bdamage.gz [options]\n"); - fprintf(fp, "Options: -names FILE -bam FILE -nodes FILE -r INT -howmany INT -ctga -countout -doOld\n"); - return 0; -} - -static int usage_print2(FILE *fp) { - fprintf(fp, "Usage: ./metaDMG-cpp print2 file.bdamage.gz [options]\n"); - fprintf(fp, "Options: -acc2tax FILE -bam FILE -nodes FILE -r INT -howmany INT -ctga -countout -doOld\n"); + fprintf(fp, "Usage: ./metaDMG-cpp print file.bdamage.gz [options]\n\n"); + fprintf(fp, "Options:\n"); + fprintf(fp, " -names FILE | -acc2tax FILE Mapping from id to display name\n"); + fprintf(fp, " -bam FILE Use BAM header reference names\n"); + fprintf(fp, " -nodes FILE nodes.dmp(.gz) for taxonomy roll-up with -r\n"); + fprintf(fp, " -r INT Select one taxid/reference id (with descendants if -nodes)\n"); + fprintf(fp, " -howmany INT Number of positions to print (default: 15)\n"); + fprintf(fp, " -ctga Print compact CT/GA-only output\n"); + fprintf(fp, " -countout Print raw counts instead of normalized frequencies\n"); + fprintf(fp, " -doOld Use legacy mode (default behavior)\n"); + fprintf(fp, " -h, --help Show this help message\n\n"); + fprintf(fp, "Examples:\n"); + fprintf(fp, " ./metaDMG-cpp print sample.bdamage.gz\n"); + fprintf(fp, " ./metaDMG-cpp print sample.bdamage.gz -names names.dmp.gz\n"); + fprintf(fp, " ./metaDMG-cpp print sample.bdamage.gz -acc2tax acc2tax.tsv.gz -r 11 -nodes nodes.dmp.gz\n"); + fprintf(fp, " ./metaDMG-cpp print sample.bdamage.gz -bam reads.bam -ctga -howmany 25\n"); + fprintf(fp, " ./metaDMG-cpp print sample.bdamage.gz -countout\n"); return 0; } @@ -199,7 +208,7 @@ int main_print(int argc, char **argv) { while (*(++argv)) { if (!strcasecmp("-h", *argv) || !strcasecmp("--help", *argv)) return usage_print(stderr); - else if (strcasecmp("-names", *argv) == 0) + else if (strcasecmp("-names", *argv) == 0 || strcasecmp("-acc2tax", *argv) == 0) infile_names = strdup(*(++argv)); else if (strcasecmp("-bam", *argv) == 0) inbam = strdup(*(++argv)); @@ -455,287 +464,6 @@ int main_print(int argc, char **argv) { return 0; } -int main_print2(int argc, char **argv) { - if (argc == 1 || (argc == 2 && (!strcasecmp(argv[1], "-h") || !strcasecmp(argv[1], "--help")))) - return usage_print2(stderr); - char *infile = NULL; - char *inbam = NULL; - char *acc2tax = NULL; - int ctga = 0; // only print ctga errors - int search = -1; - int countout = 0; - char *infile_nodes = NULL; - int howmany = 15; - int doold = 0; - while (*(++argv)) { - if (!strcasecmp("-h", *argv) || !strcasecmp("--help", *argv)) - return usage_print2(stderr); - else if (strcasecmp("-acc2tax", *argv) == 0) - acc2tax = strdup(*(++argv)); - else if (strcasecmp("-bam", *argv) == 0) - inbam = strdup(*(++argv)); - else if (strcasecmp("-r", *argv) == 0) - search = atoi(*(++argv)); - else if (strcasecmp("-howmany", *argv) == 0) - howmany = atoi(*(++argv)); - else if (strcasecmp("-ctga", *argv) == 0) - ctga = 1; - else if (strcasecmp("-doOld", *argv) == 0) - doold = 1; - else if (strcasecmp("-countout", *argv) == 0) - countout = 1; - else if (strcasecmp("-nodes", *argv) == 0) - infile_nodes = strdup(*(++argv)); - else - infile = strdup(*argv); - } - - 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); - - // 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); - if (search != -1 && doold == 0) { - std::map retmap = load_bdamage3(infile, howmany); - double *dbl = getval(retmap, child, search, howmany); - double *dbldbl = new double [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]; - - fprintf(stdout, "%d\t%.0f", search, dbldbl[0]); - for (int i = 0; i < 3 * howmany; i++) - fprintf(stdout, "\t%f", dbldbl[1 + i]); - fprintf(stdout, "\n"); - delete [] dbldbl; - return 0; - } - - BGZF *bgfp = NULL; - samFile *samfp = NULL; - bam_hdr_t *hdr = NULL; - - if (((bgfp = bgzf_open(infile, "r"))) == NULL) { - fprintf(stderr, "Could not open input bdamage.gz file: %s\n", infile); - 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; - } - } - - int printlength; - 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]; - char *type_name = NULL; - if (hdr != NULL) - type_name = strdup("Reference"); - else if (acc2tax != NULL) - type_name = strdup("FunkyName"); - else - type_name = strdup("taxid"); - - if (ctga == 0) { - fprintf(stdout, "%s\tNalignments\tDirection\tPos\tAA\tAC\tAG\tAT\tCA\tCC\tCG\tCT\tGA\tGC\tGG\tGT\tTA\tTC\tTG\tTT\n", type_name); - } else { - fprintf(stdout, "%s\tNalignment", type_name); - for (int i = 0; i < howmany; i++) - fprintf(stdout, "\tCT_%d", i); - for (int i = 0; i < howmany; i++) - fprintf(stdout, "\tGA_%d", i); - fprintf(stdout, "\n"); - } - - float data[16]; - - while (1) { - int nread = bgzf_read(bgfp, ref_nreads, 2 * sizeof(int)); - 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]); - 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]) { - if (ctga == 0) { - if (hdr != NULL) - fprintf(stdout, "%s\t%d\t5\'\t%d", hdr->target_name[ref_nreads[0]], ref_nreads[1], at); - else if (acc2tax != NULL) { - int2char::iterator itt = name_map.find(ref_nreads[0]); - if (itt == name_map.end()) { - fprintf(stderr, "\t-> Problem finding taxid: \'%d' in namedatabase: \'%s\'\n", ref_nreads[0], acc2tax); - exit(0); - } - fprintf(stdout, "\"%s\"\t%d\t5\'\t%d", itt->second, ref_nreads[1], at); - } else - fprintf(stdout, "%d\t%d\t5\'\t%d", ref_nreads[0], ref_nreads[1], at); - } else { - if (at == 0) - fprintf(stdout, "%d\t%d", ref_nreads[0], ref_nreads[1]); - } - if (countout == 1) { - for (int i = 0; i < 16; i++) - fprintf(stdout, "\t%f", data[i]); - fprintf(stdout, "\n"); - } else { - float flt[16]; - - for (int i = 0; i < 4; i++) { - double tsum = 0; - for (int j = 0; j < 4; j++) { - tsum += data[i * 4 + j]; - flt[i * 4 + j] = data[i * 4 + j]; - } - if (tsum == 0) tsum = 1; - for (int j = 0; j < 4; j++) - flt[i * 4 + j] /= tsum; - } - - if (ctga == 0) { - for (int j = 0; j < 16; j++) - fprintf(stdout, "\t%f", flt[j]); - fprintf(stdout, "\n"); - } else - ctgas[at] = flt[7]; - } - } - } - if (search == -1 || search == ref_nreads[0]) { - if (ctga == 1) { - for (int i = 0; i < howmany; i++) - fprintf(stdout, "\t%f", ctgas[i]); - } - } - - 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]) { - if (ctga == 0) { - if (hdr != NULL) - fprintf(stdout, "%s\t%d\t3\'\t%d", hdr->target_name[ref_nreads[0]], ref_nreads[1], at); - else if (acc2tax != NULL) { - int2char::iterator itt = name_map.find(ref_nreads[0]); - if (itt == name_map.end()) { - fprintf(stderr, "\t-> Problem finding taxid: \'%d' in namedatabase: \'%s\'\n", ref_nreads[0], acc2tax); - exit(0); - } - fprintf(stdout, "\"%s\"\t%d\t3\'\t%d", itt->second, ref_nreads[1], at); - } else - fprintf(stdout, "%d\t%d\t3\'\t%d", ref_nreads[0], ref_nreads[1], at); - } - if (countout == 1) { - for (int i = 0; i < 16; i++) - fprintf(stdout, "\t%f", data[i]); - fprintf(stdout, "\n"); - } else { - float flt[16]; - - for (int i = 0; i < 4; i++) { - double tsum = 0; - for (int j = 0; j < 4; j++) { - tsum += data[i * 4 + j]; - flt[i * 4 + j] = data[i * 4 + j]; - } - if (tsum == 0) tsum = 1; - for (int j = 0; j < 4; j++) - flt[i * 4 + j] /= tsum; - } - if (ctga == 0) { - for (int j = 0; j < 16; j++) - fprintf(stdout, "\t%f", flt[j]); - fprintf(stdout, "\n"); - } else - ctgas[at + printlength] = flt[8]; - } - } - } - - if (search == -1 || search == ref_nreads[0]) { - if (ctga == 1) { - for (int i = 0; i < howmany; i++) - fprintf(stdout, "\t%f", ctgas[printlength + i]); - fprintf(stdout, "\n"); - } - } - free(ctgas); - } - //clean up - for(int2char::iterator it=name_map.begin();it!=name_map.end();it++) - free(it->second); - for(int2char::iterator it=rank.begin();it!=rank.end();it++) - free(it->second); - - if (bgfp) - bgzf_close(bgfp); - if (hdr) - bam_hdr_destroy(hdr); - if (samfp) - sam_close(samfp); - if(type_name) - free(type_name); - if(infile_nodes) - free(infile_nodes); - if(infile) - free(infile); - if(inbam) - free(inbam); - if(acc2tax) - free(acc2tax); - return 0; -} int2int getlcadist(char *fname) { // fprintf(stderr,"fname: %s\n",fname); diff --git a/main_print.h b/main_print.h index f613025..425f2da 100644 --- a/main_print.h +++ b/main_print.h @@ -11,6 +11,5 @@ std::map getval_full_norec(std::map &retmap, int2int mydata2 getval_stats(std::map &retmap, int2intvec &child, int taxid); int main_print(int argc, char **argv); -int main_print2(int argc, char **argv); int main_print_all(int argc, char **argv); int main_print_ugly(int argc, char **argv); diff --git a/metaDMG.cpp b/metaDMG.cpp index 11700eb..9e20463 100644 --- a/metaDMG.cpp +++ b/metaDMG.cpp @@ -28,7 +28,6 @@ static void print_top_level_usage(FILE *fp) { fprintf(fp, "./metaDMG-cpp index files.damage.gz\n"); fprintf(fp, "./metaDMG-cpp lca [many options]\n"); fprintf(fp, "./metaDMG-cpp print bdamage.gz\n"); - fprintf(fp, "./metaDMG-cpp print2 [many options] bdamage.gz\n"); fprintf(fp, "./metaDMG-cpp print_all [many options] bdamage.gz\n"); fprintf(fp, "./metaDMG-cpp print_ugly [many options] bdamage.gz\n"); fprintf(fp, "./metaDMG-cpp dfit [many options] bdamage.gz\n"); @@ -79,8 +78,6 @@ int main(int argc, char **argv) { rc = main_aggregate(argc, argv); else if (!strcmp(argv[0], "filter_bdamage")) rc = main_filter_bdamage(argc, argv); - else if (!strcmp(argv[0], "print2")) - rc = main_print2(argc, argv); else if (!strcmp(argv[0], "lca")) rc = main_lca(argc, argv); else{