diff --git a/src/multiBamCov/multiBamCov.cpp b/src/multiBamCov/multiBamCov.cpp index 7b287d4db6b941fd02c12647b9496214a2a9b42f..4c3c7396d274915d781347a02e556b9be91de73d 100644 --- a/src/multiBamCov/multiBamCov.cpp +++ b/src/multiBamCov/multiBamCov.cpp @@ -17,12 +17,16 @@ /* Constructor */ -MultiCovBam::MultiCovBam(const vector<string> &bam_files, const string bed_file, int minQual, bool properOnly) +MultiCovBam::MultiCovBam(const vector<string> &bam_files, const string bed_file, + int minQual, bool properOnly, + bool keepDuplicates, bool keepFailedQC) : _bam_files(bam_files), _bed_file(bed_file), _minQual(minQual), -_properOnly(properOnly) +_properOnly(properOnly), +_keepDuplicates(keepDuplicates), +_keepFailedQC(keepFailedQC) { _bed = new BedFile(_bed_file); LoadBamFileMap(); @@ -76,8 +80,12 @@ void MultiCovBam::CollectCoverage() BamAlignment al; while ( reader.GetNextAlignment(al) ) { + bool duplicate = al.IsDuplicate(); + bool failedQC = al.IsFailedQC(); + if (_keepDuplicates) duplicate = false; + if (_keepFailedQC) failedQC = false; // map qual must exceed minimum - if (al.MapQuality >= _minQual) { + if ((al.MapQuality >= _minQual) && (!duplicate) && (!failedQC)) { // ignore if not properly paired and we actually care. if (_properOnly && !al.IsProperPair()) continue; diff --git a/src/multiBamCov/multiBamCov.h b/src/multiBamCov/multiBamCov.h index 95970545e9305e93c87c5537524966580b99a445..12a3bf9f9fa7332fd1af1f27395044471c06e9dd 100644 --- a/src/multiBamCov/multiBamCov.h +++ b/src/multiBamCov/multiBamCov.h @@ -30,7 +30,9 @@ class MultiCovBam { public: // constructor - MultiCovBam(const vector<string> &bam_files, const string bed_file, int minQual, bool properOnly); + MultiCovBam(const vector<string> &bam_files, const string bed_file, + int minQual, bool properOnly, + bool keepDuplicates, bool keepFailedQC); // destructor ~MultiCovBam(void); @@ -49,6 +51,9 @@ private: // attributes to control what is counted int _minQual; bool _properOnly; + bool _keepDuplicates; + bool _keepFailedQC; + map<string, int> bamFileMap; diff --git a/src/multiBamCov/multiBamCovMain.cpp b/src/multiBamCov/multiBamCovMain.cpp index 55d93cbd2388b84a53f8844e542c0cdda20bb706..69b340dc975ba4c9cc7d8215bb03ea14e5e881a8 100644 --- a/src/multiBamCov/multiBamCovMain.cpp +++ b/src/multiBamCov/multiBamCovMain.cpp @@ -38,6 +38,8 @@ int main(int argc, char* argv[]) { bool haveBed = false; bool haveBams = false; bool properOnly = false; + bool keepDuplicates = false; + bool keepFailedQC = false; // check to see if we should print out some help if(argc <= 1) showHelp = true; @@ -88,10 +90,21 @@ int main(int argc, char* argv[]) { else if(PARAMETER_CHECK("-p", 2, parameterLength)) { properOnly = true; } + else if(PARAMETER_CHECK("-D", 2, parameterLength)) { + keepDuplicates = true; + } + + else if(PARAMETER_CHECK("-F", 2, parameterLength)) { + keepFailedQC = true; + } + else { + cerr << endl << "*****ERROR: Unrecognized parameter: " << argv[i] << " *****" << endl << endl; + showHelp = true; + } } if (!showHelp) { - MultiCovBam *mc = new MultiCovBam(bamFiles, bedFile, minQual, properOnly); + MultiCovBam *mc = new MultiCovBam(bamFiles, bedFile, minQual, properOnly, keepDuplicates, keepFailedQC); mc->CollectCoverage(); delete mc; return 0; @@ -119,6 +132,10 @@ void ShowHelp(void) { cerr << "\t-q\t" << "Minimum mapping quality allowed. Default is 0." << endl << endl; + cerr << "\t-D\t" << "Include duplicate-marked reads. Default is to count non-duplicates only" << endl << endl; + + cerr << "\t-F\t" << "Include failed-QC reads. Default is to count pass-QC reads only" << endl << endl; + cerr << "\t-p\t" << "Only count proper pairs. Default is to count all alignments with MAPQ" << endl; cerr << "\t\t" << "greater than the -q argument, regardless of the BAM FLAG field." << endl << endl;