|
0
|
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 }
|
