diff BEDTools-Version-2.14.3/src/pairToBed/pairToBed.h @ 0:dfcd8b6c1bda

Uploaded
author aaronquinlan
date Thu, 03 Nov 2011 10:25:04 -0400
parents
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/BEDTools-Version-2.14.3/src/pairToBed/pairToBed.h	Thu Nov 03 10:25:04 2011 -0400
@@ -0,0 +1,161 @@
+/*****************************************************************************
+  pairToBed.h
+
+  (c) 2009 - Aaron Quinlan
+  Hall Laboratory
+  Department of Biochemistry and Molecular Genetics
+  University of Virginia
+  aaronquinlan@gmail.com
+
+  Licenced under the GNU General Public License 2.0 license.
+******************************************************************************/
+#ifndef INTERSECTBED_H
+#define INTERSECTBED_H
+
+#include "api/BamReader.h"
+#include "api/BamWriter.h"
+#include "api/BamAux.h"
+using namespace BamTools;
+
+#include "bedFile.h"
+#include "bedFilePE.h"
+#include <vector>
+#include <iostream>
+#include <fstream>
+
+using namespace std;
+
+
+
+/**************************************************
+Helper function protoypes
+**************************************************/
+void IsCorrectMappingForBEDPE (const BamAlignment &bam, const RefVector &refs, BEDPE &a);
+
+
+
+//************************************************
+// Class methods and elements
+//************************************************
+class BedIntersectPE {
+
+public:
+
+    // constructor
+    BedIntersectPE(string bedAFilePE, string bedBFile, float overlapFraction,
+        string searchType, bool sameStrand, bool diffStrand, bool bamInput, bool bamOutput, bool uncompressedBam, bool useEditDistance);
+    // destructor
+    ~BedIntersectPE(void);
+
+    void FindOverlaps(const BEDPE &, vector<BED> &hits1, vector<BED> &hits2, const string &type);
+
+    bool FindOneOrMoreOverlaps(const BEDPE &, const string &type);
+
+    void FindSpanningOverlaps(const BEDPE &a, vector<BED> &hits, const string &type);
+    bool FindOneOrMoreSpanningOverlaps(const BEDPE &a, const string &type);
+
+    void IntersectBedPE();
+    void IntersectBamPE(string bamFile);
+
+    void DetermineBedPEInput();
+
+private:
+
+    string _bedAFilePE;
+    string _bedBFile;
+    float _overlapFraction;
+    string _searchType;
+    bool _sameStrand;
+    bool _diffStrand;
+    bool _useEditDistance;
+    bool _bamInput;
+    bool _bamOutput;
+    bool  _isUncompressedBam;
+
+    // instance of a paired-end bed file class.
+    BedFilePE *_bedA;
+
+    // instance of a bed file class.
+    BedFile *_bedB;
+
+    inline
+    void ConvertBamToBedPE(const BamAlignment &bam1, const BamAlignment &bam2, const RefVector &refs, BEDPE &a) {
+
+        // initialize BEDPE variables
+        a.start1 = a.start2 = a.end1 = a.end2 = -1;
+        a.chrom1 = a.chrom2 = ".";
+        a.strand1 = a.strand2 = '.';
+        uint32_t editDistance1, editDistance2;
+        editDistance1 = editDistance2 = 0;
+
+        // take the qname from end 1.
+        a.name = bam1.Name;
+
+        // end 1
+        if (bam1.IsMapped()) {
+            a.chrom1  = refs.at(bam1.RefID).RefName;
+            a.start1  = bam1.Position;
+            a.end1    = bam1.GetEndPosition(false, false);
+            a.strand1 = "+";
+            if (bam1.IsReverseStrand()) a.strand1 = "-";
+
+            // extract the edit distance from the NM tag
+            // if possible. otherwise, complain.
+            if (_useEditDistance == true) {
+                if (bam1.GetTag("NM", editDistance1) == false) {
+                    cerr << "The edit distance tag (NM) was not found in the BAM file.  Please disable -ed.  Exiting\n";
+                    exit(1);
+                }
+            }
+        }
+
+        // end 2
+        if (bam2.IsMapped()) {
+            a.chrom2  = refs.at(bam2.RefID).RefName;
+            a.start2  = bam2.Position;
+            a.end2    = bam2.GetEndPosition(false, false);
+            a.strand2 = "+";
+            if (bam2.IsReverseStrand()) a.strand2 = "-";
+
+            // extract the edit distance from the NM tag
+            // if possible. otherwise, complain.
+            if (_useEditDistance == true) {
+                if (bam2.GetTag("NM", editDistance2) == false) {
+                    cerr << "The edit distance tag (NM) was not found in the BAM file.  Please disable -ed.  Exiting\n";
+                    exit(1);
+                }
+            }
+        }
+
+        // swap the ends if necessary
+        if ( a.chrom1 > a.chrom2 || ((a.chrom1 == a.chrom2) && (a.start1 > a.start2)) ) {
+            swap(a.chrom1, a.chrom2);
+            swap(a.start1, a.start2);
+            swap(a.end1, a.end2);
+            swap(a.strand1, a.strand2);
+        }
+
+        // compute the minimum mapping quality b/w the two ends of the pair.
+        a.score = "0";
+        if (_useEditDistance == false) {
+            if (bam1.IsMapped() == true && bam2.IsMapped() == true)
+                a.score = ToString(min(bam1.MapQuality, bam2.MapQuality));
+        }
+        // BEDPE using edit distance
+        else {
+            if (bam1.IsMapped() == true && bam2.IsMapped() == true)
+                a.score = ToString((int) (editDistance1 + editDistance2));
+            else if (bam1.IsMapped() == true)
+                a.score = ToString((int) editDistance1);
+            else if (bam2.IsMapped() == true)
+                a.score = ToString((int) editDistance2);
+        }
+    };
+
+    inline
+    void ProcessBamBlock (const BamAlignment &bam1, const BamAlignment &bam2,
+                                          const RefVector &refs,
+                                          BamWriter &writer);
+};
+
+#endif /* PEINTERSECTBED_H */