Added logic to keep dups and failed QC for multiBamCov. Thanks to Chip Stewart.

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;

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;
     

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;