Permalink
Browse files

change halWiggleLiftovers -append option to merge wiggle data in order

  • Loading branch information...
glennhickey committed Nov 2, 2013
1 parent 6bfd770 commit 560f1a790744502a5251aa748a96179b7ccb5c9d
@@ -8,6 +8,7 @@
#include <cassert>
#include "halWiggleLiftover.h"
#include "halBlockMapper.h"
+#include "halWiggleLoader.h"
using namespace std;
using namespace hal;
@@ -25,6 +26,15 @@ WiggleLiftover::~WiggleLiftover()
}
+void WiggleLiftover::preloadOutput(AlignmentConstPtr alignment,
+ const Genome* tgtGenome,
+ istream* inputFile)
+{
+ WiggleLoader loader;
+ _outVals.init(tgtGenome->getSequenceLength(), DefaultValue, DefaultTileSize);
+ loader.load(alignment, tgtGenome, inputFile, &_outVals);
+}
+
void WiggleLiftover::convert(AlignmentConstPtr alignment,
const Genome* srcGenome,
istream* inputFile,
@@ -56,7 +66,12 @@ void WiggleLiftover::convert(AlignmentConstPtr alignment,
inputSet.insert(_srcGenome);
inputSet.insert(_tgtGenome);
getGenomesInSpanningTree(inputSet, _tgtSet);
- _outVals.init(tgtGenome->getSequenceLength(), DefaultValue, DefaultTileSize);
+ // if not init'd by preload()...
+ if (_outVals.getGenomeSize() == 0)
+ {
+ _outVals.init(tgtGenome->getSequenceLength(), DefaultValue,
+ DefaultTileSize);
+ }
scan(inputFile);
write();
_outVals.clear();
@@ -24,7 +24,10 @@ static CLParserPtr initParser()
" to stream to standard output.");
optionsParser->addOptionFlag("noDupes", "do not map between duplications in"
" graph.", false);
- optionsParser->addOptionFlag("append", "append results to tgtWig", false);
+ optionsParser->addOptionFlag("append", "append/merge results into tgtWig. "
+ "Note that the entire tgtWig file will be loaded into"
+ " memory then overwritten, so this data can be lost "
+ "in event of a crash", false);
/* optionsParser->addOptionFlag("unique",
"only map block if its left-most paralog is in"
"the input. this "
@@ -107,8 +110,19 @@ int main(int argc, char** argv)
throw hal_exception("Error opening srcWig, " + srcWigPath);
}
}
-
- ios_base::openmode mode = append ? ios::out | ios::app : ios_base::out;
+
+ WiggleLiftover liftover;
+ if (append == true && tgtWigPath != "stdout")
+ {
+ // load the wig data into memory so that it can be properly merged
+ // with the new data from the liftover.
+ ifstream tgtWig(tgtWigPath.c_str());
+ if (tgtWig)
+ {
+ liftover.preloadOutput(alignment, tgtGenome, &tgtWig);
+ }
+ }
+
ofstream tgtWig;
ostream* tgtWigPtr;
if (tgtWigPath == "stdout")
@@ -117,15 +131,14 @@ int main(int argc, char** argv)
}
else
{
- tgtWig.open(tgtWigPath.c_str(), mode);
+ tgtWig.open(tgtWigPath.c_str());
tgtWigPtr = &tgtWig;
if (!tgtWig)
{
throw hal_exception("Error opening tgtWig, " + tgtWigPath);
}
}
-
- WiggleLiftover liftover;
+
liftover.convert(alignment, srcGenome, srcWigPtr, tgtGenome, tgtWigPtr,
!noDupes, unique);
}
@@ -0,0 +1,64 @@
+/*
+ * Copyright (C) 2013 by Glenn Hickey (hickey@soe.ucsc.edu)
+ *
+ * Released under the MIT license, see LICENSE.txt
+ */
+
+#include <deque>
+#include <cassert>
+#include "halWiggleLoader.h"
+#include "halWiggleLiftover.h"
+
+using namespace std;
+using namespace hal;
+
+
+WiggleLoader::WiggleLoader()
+{
+
+}
+
+WiggleLoader::~WiggleLoader()
+{
+
+}
+
+void WiggleLoader::load(AlignmentConstPtr alignment,
+ const Genome* genome,
+ istream* inputFile,
+ WiggleTiles<double>* vals)
+{
+ _alignment = alignment;
+ _srcGenome = genome;
+ _srcSequence = NULL;
+ _vals = vals;
+ scan(inputFile);
+}
+
+void WiggleLoader::visitHeader()
+{
+ _srcSequence = _srcGenome->getSequence(_sequenceName);
+ if (_srcSequence == NULL)
+ {
+ stringstream ss;
+ ss << "Sequence " << _sequenceName << " not found in genome "
+ << _srcGenome->getName();
+ throw hal_exception(ss.str());
+ }
+}
+
+void WiggleLoader::visitLine()
+{
+ if (_srcSequence == NULL)
+ {
+ throw hal_exception("Missing Wig header");
+ }
+
+ hal_index_t absFirst = _first + _srcSequence->getStartPosition();
+ hal_index_t absLast = _last + _srcSequence->getStartPosition();
+
+ for (hal_index_t absPos = absFirst; absPos <= absLast; ++absPos)
+ {
+ _vals->set(absPos, _value);
+ }
+}
@@ -22,7 +22,11 @@ class WiggleLiftover : public WiggleScanner
WiggleLiftover();
virtual ~WiggleLiftover();
-
+
+ void preloadOutput(AlignmentConstPtr alignment,
+ const Genome* tgtGenome,
+ std::istream* inputFile);
+
void convert(AlignmentConstPtr alignment,
const Genome* srcGenome,
std::istream* inputFile,
@@ -0,0 +1,47 @@
+/*
+ * Copyright (C) 2013 by Glenn Hickey (hickey@soe.ucsc.edu)
+ *
+ * Released under the MIT license, see LICENSE.txt
+ */
+
+#ifndef _HALWIGGLELOADER_H
+#define _HALWIGGLELOADER_H
+
+#include <vector>
+#include <string>
+#include <fstream>
+#include <iostream>
+#include "halWiggleScanner.h"
+#include "halWiggleTiles.h"
+
+namespace hal {
+
+/** Quick hack to load a wiggle into memory, in order to have wiggleLiftover
+ * --append option work better. Ideally would be base class of WiggleLiftover
+ * but don't have time to refactor right now */
+class WiggleLoader : public WiggleScanner
+{
+public:
+
+ WiggleLoader();
+ virtual ~WiggleLoader();
+
+ void load(AlignmentConstPtr alignment,
+ const Genome* genome,
+ std::istream* inputFile,
+ WiggleTiles<double>* vals);
+
+
+protected:
+
+ virtual void visitLine();
+ virtual void visitHeader();
+
+ AlignmentConstPtr _alignment;
+ const Genome* _srcGenome;
+ const Sequence* _srcSequence;
+ WiggleTiles<double>* _vals;
+};
+
+}
+#endif

0 comments on commit 560f1a7

Please sign in to comment.