Mercurial > repos > aaronquinlan > multi_intersect

comparison BEDTools-Version-2.14.3/src/utils/BamTools-Ancillary/BamAncillary.cpp @ 0:dfcd8b6c1bda

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author aaronquinlan
date Thu, 03 Nov 2011 10:25:04 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:dfcd8b6c1bda
1 /*****************************************************************************
2 bamAncillary.cpp
3
4 (c) 2009 - Aaron Quinlan
5 Hall Laboratory
6 Department of Biochemistry and Molecular Genetics
7 University of Virginia
8 aaronquinlan@gmail.com
9
10 Licensed under the GNU General Public License 2.0 license.
11 ******************************************************************************/
12 #include "BamAncillary.h"
13 using namespace std;
14
15 // 10 15 20 25 30 4000
16 // acccctttggacct---ataggga.................aaaa
17 // acccc---ggaccttttataggga.................aaaa
18 // 5M 3D 6M 2I 7M 20N 4M
19
20 namespace BamTools {
21 void getBamBlocks(const BamAlignment &bam, const RefVector &refs,
22 vector<BED> &blocks, bool breakOnDeletionOps) {
23
24 CHRPOS currPosition = bam.Position;
25 CHRPOS blockStart = bam.Position;
26 string chrom = refs.at(bam.RefID).RefName;
27 string name = bam.Name;
28 string strand = "+";
29 string score = ToString(bam.MapQuality);
30 char prevOp = '\0';
31 if (bam.IsReverseStrand()) strand = "-";
32 bool blocksFound = false;
33
34 vector<CigarOp>::const_iterator cigItr = bam.CigarData.begin();
35 vector<CigarOp>::const_iterator cigEnd = bam.CigarData.end();
36 for ( ; cigItr != cigEnd; ++cigItr ) {
37 if (cigItr->Type == 'M') {
38 currPosition += cigItr->Length;
39 // we only want to create a new block if the current M op
40 // was preceded by an N op or a D op (and we are breaking on D ops)
41 if ((prevOp == 'D' && breakOnDeletionOps == true) || (prevOp == 'N')) {
42 blocks.push_back( BED(chrom, blockStart, currPosition, name, score, strand) );
43 blockStart = currPosition;
44 }
45 }
46 else if (cigItr->Type == 'D') {
47 if (breakOnDeletionOps == false)
48 currPosition += cigItr->Length;
49 else {
50 currPosition += cigItr->Length;
51 blockStart = currPosition;
52 }
53 }
54 else if (cigItr->Type == 'N') {
55 currPosition += cigItr->Length;
56 blockStart = currPosition;
57 }
58 else if (cigItr->Type == 'S' || cigItr->Type == 'H' || cigItr->Type == 'P' || cigItr->Type == 'I') {
59 // do nothing
60 }
61 else {
62 cerr << "Input error: invalid CIGAR type (" << cigItr->Type
63 << ") for: " << bam.Name << endl;
64 exit(1);
65 }
66 prevOp = cigItr->Type;
67 }
68 // if there were no splits, we just create a block representing the contiguous alignment.
69 if (blocksFound == false) {
70 blocks.push_back( BED(chrom, bam.Position, currPosition, name, score, strand) );
71 }
72 }
73 }