|
0
|
1 /*****************************************************************************
|
|
|
2 complementBed.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 Licenced under the GNU General Public License 2.0 license.
|
|
|
11 ******************************************************************************/
|
|
|
12 #include "lineFileUtilities.h"
|
|
|
13 #include "complementBed.h"
|
|
|
14
|
|
|
15 BedComplement::BedComplement(string &bedFile, string &genomeFile) {
|
|
|
16
|
|
|
17 _bedFile = bedFile;
|
|
|
18 _genomeFile = genomeFile;
|
|
|
19
|
|
|
20 _bed = new BedFile(bedFile);
|
|
|
21 _genome = new GenomeFile(genomeFile);
|
|
|
22
|
|
|
23 }
|
|
|
24
|
|
|
25
|
|
|
26 BedComplement::~BedComplement(void) {
|
|
|
27 }
|
|
|
28
|
|
|
29
|
|
|
30 //
|
|
|
31 // Merge overlapping BED entries into a single entry
|
|
|
32 //
|
|
|
33 void BedComplement::ComplementBed() {
|
|
|
34
|
|
|
35 // load the "B" bed file into a map so
|
|
|
36 // that we can easily compare "A" to it for overlaps
|
|
|
37 _bed->loadBedFileIntoMapNoBin();
|
|
|
38
|
|
|
39 // get a list of the chroms in the user's genome
|
|
|
40 vector<string> chromList = _genome->getChromList();
|
|
|
41
|
|
|
42 // process each chrom in the genome
|
|
|
43 for (size_t c = 0; c < chromList.size(); ++c) {
|
|
|
44 string currChrom = chromList[c];
|
|
|
45
|
|
|
46 // create a "bit vector" for the chrom
|
|
|
47 CHRPOS currChromSize = _genome->getChromSize(currChrom);
|
|
|
48 vector<bool> chromMasks(currChromSize, 0);
|
|
|
49
|
|
|
50 // mask the chrom for every feature in the BED file
|
|
|
51 bedVector::const_iterator bItr = _bed->bedMapNoBin[currChrom].begin();
|
|
|
52 bedVector::const_iterator bEnd = _bed->bedMapNoBin[currChrom].end();
|
|
|
53 for (; bItr != bEnd; ++bItr) {
|
|
|
54 if (bItr->end > currChromSize) {
|
|
|
55 cout << "Warninge: end of BED entry exceeds chromosome length. Please correct." << endl;
|
|
|
56 _bed->reportBedNewLine(*bItr);
|
|
|
57 exit(1);
|
|
|
58 }
|
|
|
59
|
|
|
60 // mask all of the positions spanned by this BED entry.
|
|
|
61 for (CHRPOS b = bItr->start; b < bItr->end; b++)
|
|
|
62 chromMasks[b] = 1;
|
|
|
63 }
|
|
|
64
|
|
|
65 // report the unmasked, that is, complemented parts of the chrom
|
|
|
66 CHRPOS i = 0;
|
|
|
67 CHRPOS start;
|
|
|
68 while (i < chromMasks.size()) {
|
|
|
69 if (chromMasks[i] == 0) {
|
|
|
70 start = i;
|
|
|
71 while ((chromMasks[i] == 0) && (i < chromMasks.size()))
|
|
|
72 i++;
|
|
|
73
|
|
|
74 if (start > 0)
|
|
|
75 cout << currChrom << "\t" << start << "\t" << i << endl;
|
|
|
76 else
|
|
|
77 cout << currChrom << "\t" << 0 << "\t" << i << endl;
|
|
|
78 }
|
|
|
79 i++;
|
|
|
80 }
|
|
|
81 }
|
|
|
82 }
|
|
|
83
|
