diff --git a/base.cpp b/base.cpp index 874ee71..c66ca56 100644 --- a/base.cpp +++ b/base.cpp @@ -1,49 +1,48 @@ - /* - * TODO: - * check that input bams are sorted - * wiggle output - */ +/* + * TODO: + * check that input bams are sorted + * wiggle output + */ +#include +#include +#include +#include +#include +#include +#include +#include +#include "OptionParser.h" +#include "interval.h" - #include - #include - #include - #include - #include - #include - #include - #include - #include "OptionParser.h" - #include "interval.h" - - using namespace BamTools; - using namespace std; +using namespace BamTools; +using namespace std; - // types - typedef uint16_t DepthType; - typedef size_t CountType; +// types +typedef uint16_t DepthType; +typedef size_t CountType; - const char ref_char = '>'; +const char ref_char = '>'; - const string VERSION = "%prog 0.5" - "\nCopyright (C) 2014-2019 Giovanni Birolo and Andrea Telatin\n" +const string VERSION = "%prog 0.4" + "\nCopyright (C) 2014-2019 Giovanni Birolo and Andrea Telatin\n" "License MIT" ".\n" "This is free software: you are free to change and redistribute it.\n" "There is NO WARRANTY, to the extent permitted by law."; - #define debug if(true) +#define debug if(false) - struct CovEnd { +struct CovEnd { PositionType end; bool rev; // order for queuing bool operator<(const CovEnd &o) const { return end > o.end; } - }; +}; - struct Alignments { +struct Alignments { BamMultiReader input_bams; //map mapq_distr; //map > cigar_distr; @@ -68,6 +67,7 @@ bool get_next_alignment(BamAlignment & alignment) { bool more_alignments; do { + debug cerr << "[M] " << alignment.Name << ":" << alignment.Position << " | Is mapped? " << alignment.IsMapped() << endl; more_alignments = input_bams.GetNextAlignmentCore(alignment); //++mapq_distr[alignment.MapQuality]; } while (more_alignments && alignment.MapQuality < min_mapq); @@ -77,9 +77,9 @@ } vector get_ref_data() { return input_bams.GetReferenceData(); } int get_ref_id(const string &ref) { return input_bams.GetReferenceID(ref); } - }; +}; - struct Coverage { +struct Coverage { DepthType f = 0, r = 0; operator DepthType() const { return f + r; } @@ -96,9 +96,9 @@ else --f; } - }; +}; - class Output { +class Output { public: enum Format {BED, COUNTS}; @@ -163,9 +163,9 @@ const int minlen; Interval last_interval; Coverage last_coverage; - }; +}; - int main(int argc, char *argv[]) { +int main(int argc, char *argv[]) { // general options optparse::OptionParser parser = optparse::OptionParser().description("Computes coverage from alignments").usage("%prog [options] [BAM]...").version(VERSION); @@ -176,7 +176,7 @@ parser.add_option("-x", "--max-cov").metavar("MAXCOV").type("int").set_default("100000").help("print BED feature only if the coverage is lower than MAXCOV (default: %default)"); parser.add_option("-l", "--min-len").metavar("MINLEN").type("int").set_default("1").help("print BED feature only if its length is bigger (or equal to) than MINLELN (default: %default)"); parser.add_option("-f", "--flatten").action("store_true").set_default("0").help("Flatten adjacent BED features (usually when specifying -m/-x)"); - + parser.add_option("-v") .action("version") .help("prints program version"); // output options parser.add_option("--output-strands").action("store_true").set_default("0").help("outputs coverage and stats separately for each strand"); @@ -210,13 +210,14 @@ Output output(&cout, f, options.get("output_strands"), minimum_coverage, maximum_coverage, minimum_length); Alignments input(parser.args(), options.get("min_mapq")); - + // main alignment parsing loop BamAlignment alignment; bool more_alignments = input.get_next_alignment(alignment); for (const auto &ref : input.get_ref_data()) { // loop on reference // init new reference data const int ref_id = input.get_ref_id(ref.RefName); + debug cerr << "[R] Reference: " << ref_id << endl; PositionType last_pos = 0; priority_queue coverage_ends; Coverage coverage; @@ -254,13 +255,17 @@ debug cerr << "[<] Coverage is " << coverage << " from " << next_change << endl; last_pos = next_change; } while (last_pos != ref.RefLength); - + debug cerr << "[_] Completed at " << ref.RefName << ":" << last_pos << endl; // reference ended assert(coverage_ends.empty()); } + while (more_alignments && !(alignment.IsMapped()) ) { + debug cerr << "[8] Emtpying trash" << endl; + more_alignments = input.get_next_alignment(alignment); + } assert(!more_alignments); } catch (string msg) { parser.error(msg); } return 0; - } +}