diff BEDTools-Version-2.14.3/src/utils/BamTools-Ancillary/BamAncillary.cpp @ 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/utils/BamTools-Ancillary/BamAncillary.cpp	Thu Nov 03 10:25:04 2011 -0400
@@ -0,0 +1,73 @@
+/*****************************************************************************
+  bamAncillary.cpp
+
+  (c) 2009 - Aaron Quinlan
+  Hall Laboratory
+  Department of Biochemistry and Molecular Genetics
+  University of Virginia
+  aaronquinlan@gmail.com
+
+  Licensed under the GNU General Public License 2.0 license.
+******************************************************************************/
+#include "BamAncillary.h"
+using namespace std;
+
+// 10   15   20      25    30               4000
+// acccctttggacct---ataggga.................aaaa
+// acccc---ggaccttttataggga.................aaaa
+// 5M   3D 6M    2I 7M      20N             4M
+
+namespace BamTools {
+    void getBamBlocks(const BamAlignment &bam, const RefVector &refs,
+                      vector<BED> &blocks, bool breakOnDeletionOps) {
+
+        CHRPOS currPosition = bam.Position;
+        CHRPOS blockStart   = bam.Position;
+        string chrom        = refs.at(bam.RefID).RefName;
+        string name         = bam.Name;
+        string strand       = "+";
+        string score        = ToString(bam.MapQuality);
+        char  prevOp        = '\0';
+        if (bam.IsReverseStrand()) strand = "-";
+        bool blocksFound = false;
+
+        vector<CigarOp>::const_iterator cigItr = bam.CigarData.begin();
+        vector<CigarOp>::const_iterator cigEnd = bam.CigarData.end();
+        for ( ; cigItr != cigEnd; ++cigItr ) {
+            if (cigItr->Type == 'M') {
+                currPosition += cigItr->Length;
+                // we only want to create a new block if the current M op
+                // was preceded by an N op or a D op (and we are breaking on D ops)
+                if ((prevOp == 'D' && breakOnDeletionOps == true) || (prevOp == 'N')) {
+                    blocks.push_back( BED(chrom, blockStart, currPosition, name, score, strand) );
+                    blockStart = currPosition;
+                }
+            }
+            else if (cigItr->Type == 'D') {
+                if (breakOnDeletionOps == false)
+                    currPosition += cigItr->Length;
+                else {
+                    currPosition += cigItr->Length;
+                    blockStart    = currPosition;
+                }
+            }
+            else if (cigItr->Type == 'N') {
+                currPosition += cigItr->Length;
+                blockStart    = currPosition;
+            }
+            else if (cigItr->Type == 'S' || cigItr->Type == 'H' || cigItr->Type == 'P' || cigItr->Type == 'I') {
+                // do nothing
+            }
+            else {
+                cerr << "Input error: invalid CIGAR type (" << cigItr->Type
+                    << ") for: " << bam.Name << endl;
+                exit(1);
+            }
+            prevOp = cigItr->Type;
+        }
+        // if there were no splits, we just create a block representing the contiguous alignment.
+        if (blocksFound == false) {
+            blocks.push_back( BED(chrom, bam.Position, currPosition, name, score, strand) );
+        }
+    }
+}