Mercurial > repos > aaronquinlan > multi_intersect

comparison BEDTools-Version-2.14.3/src/flankBed/flankBed.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 flankBed.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 "lineFileUtilities.h"
13 #include "flankBed.h"
14
15
16 BedFlank::BedFlank(string &bedFile, string &genomeFile, bool forceStrand, float leftFlank, float rightFlank, bool fractional) {
17
18 _bedFile = bedFile;
19 _genomeFile = genomeFile;
20 _forceStrand = forceStrand;
21 _leftFlank = leftFlank;
22 _rightFlank = rightFlank;
23 _fractional = fractional;
24
25 _bed = new BedFile(bedFile);
26 _genome = new GenomeFile(genomeFile);
27
28 // get going, slop it up.
29 FlankBed();
30 }
31
32
33 BedFlank::~BedFlank(void) {
34
35 }
36
37
38 void BedFlank::FlankBed() {
39
40 int lineNum = 0;
41 BED bedEntry, nullBed; // used to store the current BED line from the BED file.
42 BedLineStatus bedStatus;
43
44 _bed->Open();
45 bedStatus = _bed->GetNextBed(bedEntry, lineNum);
46 while (bedStatus != BED_INVALID) {
47 if (bedStatus == BED_VALID) {
48
49 int leftFlank = _leftFlank;
50 int rightFlank = _rightFlank;
51 if (_fractional == true) {
52 leftFlank = (int) (_leftFlank * bedEntry.size());
53 rightFlank = (int) (_rightFlank * bedEntry.size());
54 }
55
56 if ((_forceStrand == false) || (bedEntry.strand == "+"))
57 {
58 AddFlank(bedEntry, leftFlank, rightFlank);
59 }
60 else if ((_forceStrand == true) && (bedEntry.strand == "-" ))
61 {
62 AddStrandedFlank(bedEntry, leftFlank, rightFlank);
63 }
64 bedEntry = nullBed;
65 }
66 bedStatus = _bed->GetNextBed(bedEntry, lineNum);
67 }
68 _bed->Close();
69 }
70
71
72 void BedFlank::AddFlank(BED &bed, int leftFlank, int rightFlank) {
73
74 int chromSize = _genome->getChromSize(bed.chrom);
75 if (chromSize == -1) {
76 cerr << "ERROR: chrom \"" << bed.chrom << "\" not found in genome file. Exiting." << endl;
77 exit(1);
78 }
79
80 // init. our left and right flanks to the original BED entry.
81 // we'll create the flanks from these coordinates.
82 BED left = bed;
83 BED right = bed;
84
85 // make the left flank (if necessary)
86 if (leftFlank > 0) {
87 if ( (static_cast<int>(left.start) - leftFlank) > 0)
88 {
89 left.end = left.start;
90 left.start -= leftFlank;
91 }
92 else
93 {
94 left.end = left.start;
95 left.start = 0;
96 }
97 // report the left flank
98 _bed->reportBedNewLine(left);
99 }
100
101 // make the left flank (if necessary)
102 if (rightFlank > 0) {
103 if ( (static_cast<int>(right.end) + (rightFlank+1)) <= static_cast<int>(chromSize))
104 {
105 right.start = right.end;
106 right.end += (rightFlank);
107 }
108 else {
109 right.start = right.end;
110 right.end += chromSize;
111 }
112 // report the right flank
113 _bed->reportBedNewLine(right);
114 }
115 }
116
117
118 void BedFlank::AddStrandedFlank(BED &bed, int leftFlank, int rightFlank) {
119
120 int chromSize = _genome->getChromSize(bed.chrom);
121 if (chromSize == -1) {
122 cerr << "ERROR: chrom \"" << bed.chrom << "\" not found in genome file. Exiting." << endl;
123 exit(1);
124 }
125
126 // init. our left and right flanks to the original BED entry.
127 // we'll create the flanks from these coordinates.
128 BED left = bed;
129 BED right = bed;
130
131 // make the left flank (if necessary)
132 if (rightFlank > 0) {
133 if ( (static_cast<int>(left.start) - rightFlank) > 0)
134 {
135 left.end = left.start;
136 left.start -= rightFlank;
137 }
138 else
139 {
140 left.end = left.start;
141 left.start = 0;
142 }
143 // report the left flank
144 _bed->reportBedNewLine(left);
145 }
146
147 // make the left flank (if necessary)
148 if (leftFlank > 0) {
149 if ( (static_cast<int>(right.end) + leftFlank) <= static_cast<int>(chromSize))
150 {
151 right.start = right.end;
152 right.end += leftFlank;
153 }
154 else {
155 right.start = right.end;
156 right.end = chromSize;
157 }
158 // report the right flank
159 _bed->reportBedNewLine(right);
160 }
161 }
162
163