Skip to content

Commit

Permalink
keep bitset in order to only write mapped positions. no more writing …
Browse files Browse the repository at this point in the history
…default values in gaps
  • Loading branch information
glennhickey committed Jul 24, 2013
1 parent ed7d92c commit 2c1120a
Show file tree
Hide file tree
Showing 2 changed files with 52 additions and 10 deletions.
35 changes: 25 additions & 10 deletions liftover/impl/halWiggleLiftover.cpp
Expand Up @@ -208,6 +208,8 @@ void WiggleLiftover::write()
{
const Sequence* outSequence = NULL;
hal_size_t ogSize = _tgtGenome->getSequenceLength();
bool needHeader = true;
hal_index_t prevPos = NULL_INDEX;
for (hal_size_t i = 0; i < _outVals.getNumTiles(); ++i)
{
if (_outVals.isTileEmpty(i) == false)
Expand All @@ -216,18 +218,31 @@ void WiggleLiftover::write()
for (hal_size_t j = 0; pos < ogSize && j < _outVals.getTileSize(); ++j,
++pos)
{
if (outSequence == NULL || pos < outSequence->getStartPosition() ||
pos > outSequence->getEndPosition())
if (_outVals.exists(pos) == true)
{
outSequence = _tgtGenome->getSequenceBySite(pos);
assert(outSequence != NULL);
*_outStream << "fixedStep"
<< "\tchrom=" << outSequence->getName()
<< "\tstart="
<< (1 + pos - outSequence->getStartPosition())
<< "\tstep=1\n";
if (outSequence == NULL || pos < outSequence->getStartPosition() ||
pos > outSequence->getEndPosition())
{
outSequence = _tgtGenome->getSequenceBySite(pos);
assert(outSequence != NULL);
needHeader = true;
}
else if (pos != prevPos + 1)
{
needHeader = true;
}
if (needHeader == true)
{
*_outStream << "fixedStep"
<< "\tchrom=" << outSequence->getName()
<< "\tstart="
<< (1 + pos - outSequence->getStartPosition())
<< "\tstep=1\n";
needHeader = false;
}
*_outStream << _outVals.get(pos) << '\n';
prevPos = pos;
}
*_outStream << _outVals.get(pos) << '\n';
}
}
}
Expand Down
27 changes: 27 additions & 0 deletions liftover/inc/halWiggleTiles.h
Expand Up @@ -42,6 +42,11 @@ class WiggleTiles
* pos which is set to val*/
void set(hal_index_t pos, T val);

/** Test if a value was written to using set(). Doing a get() isn't
* really sufficient because it will return the default value in the case
* where it was not set, or if it was set with the default value */
bool exists(hal_index_t pos) const;

/** Methods to get basic structure info */
hal_size_t getGenomeSize() const;
hal_size_t getTileSize() const;
Expand All @@ -52,6 +57,7 @@ class WiggleTiles
protected:

std::vector<std::vector<T> > _tiles;
std::vector<std::vector<bool> > _bits;
hal_size_t _tileSize;
hal_size_t _genomeSize;
hal_size_t _lastTileSize;
Expand Down Expand Up @@ -91,14 +97,17 @@ inline void WiggleTiles<T>::init(hal_size_t genomeSize, T defaultValue,
for (size_t i = 0; i < _tiles.size(); ++i)
{
_tiles[i].clear();
_bits[i].clear();
}
_tiles.resize(numTiles);
_bits.resize(numTiles);
}

template <class T>
inline void WiggleTiles<T>::clear()
{
_tiles.clear();
_bits.clear();
}

template <class T>
Expand Down Expand Up @@ -129,9 +138,27 @@ inline void WiggleTiles<T>::set(hal_index_t pos, T val)
{
hal_size_t len = tile == _tiles.size() - 1 ? _lastTileSize : _tileSize;
_tiles[tile].assign(len, _defaultValue);
_bits[tile].assign(len, false);
}
hal_size_t offset = pos % _tileSize;
_tiles[tile][offset] = val;
_bits[tile][offset] = true;
}

template <class T>
inline bool WiggleTiles<T>::exists(hal_index_t pos) const
{
assert(pos < _genomeSize);
hal_size_t tile = pos / _tileSize;
assert(tile < _tiles.size());
assert(_tiles[tile].size() == 0 || _tiles[tile].size() == _tileSize ||
_tiles[tile].size() == _lastTileSize);
if (_tiles[tile].size() == 0)
{
return false;
}
hal_size_t offset = pos % _tileSize;
return _bits[tile][offset];
}

template <class T>
Expand Down

0 comments on commit 2c1120a

Please sign in to comment.