From cc031ba9de39d8a6ba67887bcae82973ea05e4cb Mon Sep 17 00:00:00 2001
From: Aaron <aaronquinlan@gmail.com>
Date: Wed, 22 Jun 2011 21:24:23 -0400
Subject: [PATCH] Added minQual (-q) and force proper pairs (-p) to
multiBamCov.
---
src/multiBamCov/multiBamCov.cpp | 16 +++++++++++-----
src/multiBamCov/multiBamCov.h | 6 +++++-
src/multiBamCov/multiBamCovMain.cpp | 19 ++++++++++++++++---
3 files changed, 32 insertions(+), 9 deletions(-)
diff --git a/src/multiBamCov/multiBamCov.cpp b/src/multiBamCov/multiBamCov.cpp
index e19b3c9f..57006775 100644
--- a/src/multiBamCov/multiBamCov.cpp
+++ b/src/multiBamCov/multiBamCov.cpp
@@ -17,10 +17,12 @@
/*
Constructor
*/
-MultiCovBam::MultiCovBam(const vector<string> &bam_files, const string bed_file)
+MultiCovBam::MultiCovBam(const vector<string> &bam_files, const string bed_file, int minQual, bool properOnly)
:
_bam_files(bam_files),
-_bed_file(bed_file)
+_bed_file(bed_file),
+_minQual(minQual),
+_properOnly(properOnly)
{
_bed = new BedFile(_bed_file);
LoadBamFileMap();
@@ -72,9 +74,13 @@ void MultiCovBam::CollectCoverage()
BamAlignment al;
while ( reader.GetNextAlignment(al) )
{
- // lookup the offset of the file name and tabulate coverage
- // for the appropriate file
- counts[bamFileMap[al.Filename]]++;
+ if (al.MapQuality >= _minQual) {
+ if (_properOnly && !IsProperPair())
+ continue;
+ // lookup the offset of the file name and tabulate coverage
+ // for the appropriate file
+ counts[bamFileMap[al.Filename]]++;
+ }
}
// report the cov at this interval for each file and reset
_bed->reportBedTab(bed);
diff --git a/src/multiBamCov/multiBamCov.h b/src/multiBamCov/multiBamCov.h
index 2e8f0fb6..95970545 100644
--- a/src/multiBamCov/multiBamCov.h
+++ b/src/multiBamCov/multiBamCov.h
@@ -30,7 +30,7 @@ class MultiCovBam {
public:
// constructor
- MultiCovBam(const vector<string> &bam_files, const string bed_file);
+ MultiCovBam(const vector<string> &bam_files, const string bed_file, int minQual, bool properOnly);
// destructor
~MultiCovBam(void);
@@ -45,6 +45,10 @@ private:
vector<string> _bam_files;
string _bed_file;
BedFile *_bed;
+
+ // attributes to control what is counted
+ int _minQual;
+ bool _properOnly;
map<string, int> bamFileMap;
diff --git a/src/multiBamCov/multiBamCovMain.cpp b/src/multiBamCov/multiBamCovMain.cpp
index 66552e69..4e00ffdf 100644
--- a/src/multiBamCov/multiBamCovMain.cpp
+++ b/src/multiBamCov/multiBamCovMain.cpp
@@ -32,11 +32,13 @@ int main(int argc, char* argv[]) {
// input files
string bedFile;
vector<string> bamFiles;
+ minQual = 0;
// input arguments
bool haveBed = false;
bool haveBams = false;
-
+ bool properOnly = false;
+
// check to see if we should print out some help
if(argc <= 1) showHelp = true;
@@ -77,10 +79,19 @@ int main(int argc, char* argv[]) {
i--;
}
}
+ else if(PARAMETER_CHECK("-q", 2, parameterLength)) {
+ if ((i+1) < argc) {
+ minQual = atoi(argv[i + 1]);
+ i++;
+ }
+ }
+ else if(PARAMETER_CHECK("-p", 2, parameterLength)) {
+ properOnly = true;
+ }
}
if (!showHelp) {
- MultiCovBam *mc = new MultiCovBam(bamFiles, bedFile);
+ MultiCovBam *mc = new MultiCovBam(bamFiles, bedFile, minQual, properOnly);
mc->CollectCoverage();
delete mc;
return 0;
@@ -102,11 +113,13 @@ void ShowHelp(void) {
cerr << "Options: " << endl;
- cerr << "\t-bams\t" << "The bam files." << endl << endl;
+ cerr << "\t-bams\t" << "The bam files." << endl << endl;
cerr << "\t-bed\t" << "The bed file." << endl << endl;
+ cerr << "\t-q [INT]\t" << "Minimum mapping quality allowed. Default is 0." << endl << endl;
+ cerr << "\t-p\t" << "Omly count proper pairs. Default is to count all alignments >= -q" << endl << endl;
// end the program here
exit(1);
--
GitLab