Skip to content
Snippets Groups Projects
Commit 58b1ce0b authored by Aaron's avatar Aaron
Browse files

80 char width

parent baa7e8b5
No related branches found
No related tags found
No related merge requests found
...@@ -102,7 +102,9 @@ void BedFile::Open(void) { ...@@ -102,7 +102,9 @@ void BedFile::Open(void) {
_bedStream = new igzstream(bedFile.c_str(), ios::in); _bedStream = new igzstream(bedFile.c_str(), ios::in);
} }
if ( !(_bedStream->good()) ) { if ( !(_bedStream->good()) ) {
cerr << "Error: The requested bed file (" << bedFile << ") could not be opened. Exiting!" << endl; cerr << "Error: The requested bed file ("
<< bedFile
<< ") could not be opened. Exiting!" << endl;
exit (1); exit (1);
} }
} }
...@@ -270,9 +272,12 @@ bool BedFile::GetNextMergedBed(BED &merged_bed) { ...@@ -270,9 +272,12 @@ bool BedFile::GetNextMergedBed(BED &merged_bed) {
} }
void BedFile::allHits(string chrom, CHRPOS start, CHRPOS end, string strand, void BedFile::allHits(string chrom, CHRPOS start,
vector<BED> &hits, bool sameStrand, bool diffStrand, CHRPOS end, string strand,
float overlapFraction, bool reciprocal) { vector<BED> &hits, bool sameStrand,
bool diffStrand, float overlapFraction,
bool reciprocal)
{
BIN startBin, endBin; BIN startBin, endBin;
startBin = (start >> _binFirstShift); startBin = (start >> _binFirstShift);
...@@ -297,7 +302,10 @@ void BedFile::allHits(string chrom, CHRPOS start, CHRPOS end, string strand, ...@@ -297,7 +302,10 @@ void BedFile::allHits(string chrom, CHRPOS start, CHRPOS end, string strand,
CHRPOS e = min(end, bedItr->end); CHRPOS e = min(end, bedItr->end);
int overlapBases = (e - s); int overlapBases = (e - s);
// 1. is there sufficient overlap w.r.t A? // 1. is there sufficient overlap w.r.t A?
if ( (float) overlapBases / (float) aLength >= overlapFraction) { if ( (float) overlapBases
/
(float) aLength >= overlapFraction)
{
CHRPOS bLength = (bedItr->end - bedItr->start); CHRPOS bLength = (bedItr->end - bedItr->start);
float bOverlap = ( (float) overlapBases / (float) bLength ); float bOverlap = ( (float) overlapBases / (float) bLength );
bool strands_are_same = (strand == bedItr->strand); bool strands_are_same = (strand == bedItr->strand);
...@@ -352,7 +360,10 @@ bool BedFile::anyHits(string chrom, CHRPOS start, CHRPOS end, string strand, ...@@ -352,7 +360,10 @@ bool BedFile::anyHits(string chrom, CHRPOS start, CHRPOS end, string strand,
CHRPOS e = min(end, bedItr->end); CHRPOS e = min(end, bedItr->end);
int overlapBases = (e - s); int overlapBases = (e - s);
// 1. is there sufficient overlap w.r.t A? // 1. is there sufficient overlap w.r.t A?
if ( (float) overlapBases / (float) aLength >= overlapFraction) { if ( (float) overlapBases
/
(float) aLength >= overlapFraction)
{
CHRPOS bLength = (bedItr->end - bedItr->start); CHRPOS bLength = (bedItr->end - bedItr->start);
float bOverlap = ( (float) overlapBases / (float) bLength ); float bOverlap = ( (float) overlapBases / (float) bLength );
bool strands_are_same = (strand == bedItr->strand); bool strands_are_same = (strand == bedItr->strand);
...@@ -393,9 +404,9 @@ void BedFile::countHits(const BED &a, bool sameStrand, bool diffStrand, bool cou ...@@ -393,9 +404,9 @@ void BedFile::countHits(const BED &a, bool sameStrand, bool diffStrand, bool cou
// loop through each bin at this level of the hierarchy // loop through each bin at this level of the hierarchy
BIN offset = _binOffsetsExtended[i]; BIN offset = _binOffsetsExtended[i];
for (BIN j = (startBin+offset); j <= (endBin+offset); ++j) { for (BIN j = (startBin+offset); j <= (endBin+offset); ++j) {
// loop through each feature in this chrom/bin and see if it overlaps // loop through each feature in this chrom/bin and
// with the feature that was passed in. if so, add the feature to // see if it overlaps with the feature that was passed in.
// the list of hits. // if so, add the feature to the list of hits.
vector<BEDCOV>::iterator bedItr = bedCovMap[a.chrom][j].begin(); vector<BEDCOV>::iterator bedItr = bedCovMap[a.chrom][j].begin();
vector<BEDCOV>::iterator bedEnd = bedCovMap[a.chrom][j].end(); vector<BEDCOV>::iterator bedEnd = bedCovMap[a.chrom][j].end();
for (; bedItr != bedEnd; ++bedItr) { for (; bedItr != bedEnd; ++bedItr) {
...@@ -408,8 +419,9 @@ void BedFile::countHits(const BED &a, bool sameStrand, bool diffStrand, bool cou ...@@ -408,8 +419,9 @@ void BedFile::countHits(const BED &a, bool sameStrand, bool diffStrand, bool cou
{ {
continue; continue;
} }
else if (overlaps(bedItr->start, bedItr->end, a.start, a.end) > 0) { else if (overlaps(bedItr->start, bedItr->end, a.start, a.end)
> 0)
{
bedItr->count++; bedItr->count++;
if (countsOnly == false) { if (countsOnly == false) {
if (a.zeroLength == false) { if (a.zeroLength == false) {
...@@ -417,9 +429,10 @@ void BedFile::countHits(const BED &a, bool sameStrand, bool diffStrand, bool cou ...@@ -417,9 +429,10 @@ void BedFile::countHits(const BED &a, bool sameStrand, bool diffStrand, bool cou
bedItr->depthMap[a.end].ends++; bedItr->depthMap[a.end].ends++;
} }
else { else {
// correct for the fact that we artificially expanded the zeroLength feature // correct for the fact that we artificially
// expanded the zeroLength feature
bedItr->depthMap[a.start+2].starts++; bedItr->depthMap[a.start+2].starts++;
bedItr->depthMap[a.end-1].ends++; bedItr->depthMap[a.end-1].ends++;
} }
if (a.start < bedItr->minOverlapStart) { if (a.start < bedItr->minOverlapStart) {
...@@ -456,15 +469,16 @@ void BedFile::countSplitHits(const vector<BED> &bedBlocks, bool sameStrand, bool ...@@ -456,15 +469,16 @@ void BedFile::countSplitHits(const vector<BED> &bedBlocks, bool sameStrand, bool
// loop through each bin at this level of the hierarchy // loop through each bin at this level of the hierarchy
BIN offset = _binOffsetsExtended[i]; BIN offset = _binOffsetsExtended[i];
for (BIN j = (startBin+offset); j <= (endBin+offset); ++j) { for (BIN j = (startBin+offset); j <= (endBin+offset); ++j) {
// loop through each feature in this chrom/bin and see if it
// loop through each feature in this chrom/bin and see if it overlaps // overlaps with the feature that was passed in.
// with the feature that was passed in. if so, add the feature to // if so, add the feature to the list of hits.
// the list of hits. vector<BEDCOV>::iterator
vector<BEDCOV>::iterator bedItr = bedCovMap[blockItr->chrom][j].begin(); bedItr = bedCovMap[blockItr->chrom][j].begin();
vector<BEDCOV>::iterator bedEnd = bedCovMap[blockItr->chrom][j].end(); vector<BEDCOV>::iterator
bedEnd = bedCovMap[blockItr->chrom][j].end();
for (; bedItr != bedEnd; ++bedItr) { for (; bedItr != bedEnd; ++bedItr) {
bool strands_are_same =
bool strands_are_same = (blockItr->strand == bedItr->strand); (blockItr->strand == bedItr->strand);
// skip the hit if not on the same strand (and we care) // skip the hit if not on the same strand (and we care)
if ((sameStrand == true && strands_are_same == false) || if ((sameStrand == true && strands_are_same == false) ||
(diffStrand == true && strands_are_same == true) (diffStrand == true && strands_are_same == true)
...@@ -472,14 +486,17 @@ void BedFile::countSplitHits(const vector<BED> &bedBlocks, bool sameStrand, bool ...@@ -472,14 +486,17 @@ void BedFile::countSplitHits(const vector<BED> &bedBlocks, bool sameStrand, bool
{ {
continue; continue;
} }
else if (overlaps(bedItr->start, bedItr->end, blockItr->start, blockItr->end) > 0) { else if (overlaps(bedItr->start, bedItr->end,
blockItr->start, blockItr->end) > 0)
{
if (countsOnly == false) { if (countsOnly == false) {
if (blockItr->zeroLength == false) { if (blockItr->zeroLength == false) {
bedItr->depthMap[blockItr->start+1].starts++; bedItr->depthMap[blockItr->start+1].starts++;
bedItr->depthMap[blockItr->end].ends++; bedItr->depthMap[blockItr->end].ends++;
} }
else { else {
// correct for the fact that we artificially expanded the zeroLength feature // correct for the fact that we artificially
// expanded the zeroLength feature
bedItr->depthMap[blockItr->start+2].starts++; bedItr->depthMap[blockItr->start+2].starts++;
bedItr->depthMap[blockItr->end-1].ends++; bedItr->depthMap[blockItr->end-1].ends++;
} }
...@@ -501,8 +518,8 @@ void BedFile::countSplitHits(const vector<BED> &bedBlocks, bool sameStrand, bool ...@@ -501,8 +518,8 @@ void BedFile::countSplitHits(const vector<BED> &bedBlocks, bool sameStrand, bool
set< vector<BEDCOV>::iterator >::iterator validHitsItr = validHits.begin(); set< vector<BEDCOV>::iterator >::iterator validHitsItr = validHits.begin();
set< vector<BEDCOV>::iterator >::iterator validHitsEnd = validHits.end(); set< vector<BEDCOV>::iterator >::iterator validHitsEnd = validHits.end();
for (; validHitsItr != validHitsEnd; ++validHitsItr) for (; validHitsItr != validHitsEnd; ++validHitsItr)
// the validHitsItr points to another itr, hence the (*itr)-> dereferencing. // the validHitsItr points to another itr, hence
// ugly, but that's C++. // the (*itr)-> dereferencing. ugly, but that's C++.
(*validHitsItr)->count++; (*validHitsItr)->count++;
} }
...@@ -520,11 +537,13 @@ void BedFile::countListHits(const BED &a, int index, bool sameStrand, bool diffS ...@@ -520,11 +537,13 @@ void BedFile::countListHits(const BED &a, int index, bool sameStrand, bool diffS
BIN offset = _binOffsetsExtended[i]; BIN offset = _binOffsetsExtended[i];
for (BIN j = (startBin+offset); j <= (endBin+offset); ++j) { for (BIN j = (startBin+offset); j <= (endBin+offset); ++j) {
// loop through each feature in this chrom/bin and see if it overlaps // loop through each feature in this chrom/bin and see if it
// with the feature that was passed in. if so, add the feature to // overlaps with the feature that was passed in. if so,
// the list of hits. // add the feature tothe list of hits.
vector<BEDCOVLIST>::iterator bedItr = bedCovListMap[a.chrom][j].begin(); vector<BEDCOVLIST>::iterator
vector<BEDCOVLIST>::iterator bedEnd = bedCovListMap[a.chrom][j].end(); bedItr = bedCovListMap[a.chrom][j].begin();
vector<BEDCOVLIST>::iterator
bedEnd = bedCovListMap[a.chrom][j].end();
for (; bedItr != bedEnd; ++bedItr) { for (; bedItr != bedEnd; ++bedItr) {
bool strands_are_same = (a.strand == bedItr->strand); bool strands_are_same = (a.strand == bedItr->strand);
...@@ -535,16 +554,19 @@ void BedFile::countListHits(const BED &a, int index, bool sameStrand, bool diffS ...@@ -535,16 +554,19 @@ void BedFile::countListHits(const BED &a, int index, bool sameStrand, bool diffS
{ {
continue; continue;
} }
else if (overlaps(bedItr->start, bedItr->end, a.start, a.end) > 0) { else if (overlaps(bedItr->start, bedItr->end,
a.start, a.end) > 0)
{
bedItr->counts[index]++; bedItr->counts[index]++;
if (a.zeroLength == false) { if (a.zeroLength == false) {
bedItr->depthMapList[index][a.start+1].starts++; bedItr->depthMapList[index][a.start+1].starts++;
bedItr->depthMapList[index][a.end].ends++; bedItr->depthMapList[index][a.end].ends++;
} }
else { else {
// correct for the fact that we artificially expanded the zeroLength feature // correct for the fact that we artificially expanded
// the zeroLength feature
bedItr->depthMapList[index][a.start+2].starts++; bedItr->depthMapList[index][a.start+2].starts++;
bedItr->depthMapList[index][a.end-1].ends++; bedItr->depthMapList[index][a.end-1].ends++;
} }
if (a.start < bedItr->minOverlapStarts[index]) { if (a.start < bedItr->minOverlapStarts[index]) {
...@@ -656,7 +678,10 @@ void BedFile::loadBedFileIntoMapNoBin() { ...@@ -656,7 +678,10 @@ void BedFile::loadBedFileIntoMapNoBin() {
// sort the BED entries for each chromosome // sort the BED entries for each chromosome
// in ascending order of start position // in ascending order of start position
for (masterBedMapNoBin::iterator m = this->bedMapNoBin.begin(); m != this->bedMapNoBin.end(); ++m) { for (masterBedMapNoBin::iterator m = this->bedMapNoBin.begin();
m != this->bedMapNoBin.end();
++m)
{
sort(m->second.begin(), m->second.end(), sortByStart); sort(m->second.begin(), m->second.end(), sortByStart);
} }
} }
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment