Efficient MAF revamp dump

Passes Dent's old mafValidator, but thorough tests are showing a few discrepancies.
ComparativeGenomicsToolkit · Aug 23, 2018 · 0b490aa · 0b490aa
1 parent 276eb8c
commit 0b490aa
Show file tree

Hide file tree

Showing 4 changed files with 1,758 additions and 5 deletions.
diff --git a/.travis.yml b/.travis.yml
@@ -6,6 +6,7 @@ before_install:
   - if [[ "$TRAVIS_OS_NAME" == "linux" ]]; then sudo apt-get install -y libhdf5-7 libhdf5-dev; fi
   - if [[ "$TRAVIS_OS_NAME" == "osx" ]]; then brew install hdf5 || echo "a brew error code when installing gcc is expected"; fi
   - git clone https://github.com/ComparativeGenomicsToolkit/sonLib.git
+  - pip install newick attrs pytest
 install:
   - sh -c 'cd sonLib && make'
 script:

diff --git a/Makefile b/Makefile
@@ -15,6 +15,7 @@ clean.%:
 	cd $* && ${MAKE} clean
 
 test: all
+	pytest maf/impl/naiveLiftUp.py
 	python allTests.py
 
 doxy : ${modules:%=doxy.%}

diff --git a/extract/impl/halAlignedExtract.cpp b/extract/impl/halAlignedExtract.cpp
@@ -14,7 +14,8 @@ using namespace hal;
 
 static void extractAlignedRegions(const Genome* genome,
                                   ostream* bedStream,
-                                  bool complement);
+                                  bool complement,
+                                  bool viewParentCoords);
 
 static CLParserPtr initParser()
 {
@@ -30,6 +31,8 @@ static CLParserPtr initParser()
                                "genome that are *unaligned* to the parent. "
                                "ie all intervals that are not returned with"
                                " the default setting.", false);
+  optionsParser->addOptionFlag("viewParentCoords", "view the corresponding parent "
+                               "coordinates for aligned regions", false);
 
   return optionsParser;
 }
@@ -42,13 +45,15 @@ int main(int argc, char** argv)
   string genomeName;
   string outBedPath;
   bool complement;
+  bool viewParentCoords;
   try
   {
     optionsParser->parseOptions(argc, argv);
     halPath = optionsParser->getArgument<string>("halPath");
     genomeName = optionsParser->getArgument<string>("genome");
     outBedPath = optionsParser->getOption<string>("alignedFile");
     complement = optionsParser->getFlag("complement");
+    viewParentCoords = optionsParser->getFlag("viewParentCoords");
   }
   catch(exception& e)
   {
@@ -63,7 +68,7 @@ int main(int argc, char** argv)
                                                              optionsParser);
     if (inAlignment->getNumGenomes() == 0)
     {
-      throw hal_exception("input hal alignmenet is empty");
+      throw hal_exception("input hal alignment is empty");
     }
 
     const Genome* genome = inAlignment->openGenome(genomeName);
@@ -87,7 +92,7 @@ int main(int argc, char** argv)
         throw hal_exception(string("Error opening ") + outBedPath + 
                             " for writing");
       }
-      extractAlignedRegions(genome, bedStream, complement);
+      extractAlignedRegions(genome, bedStream, complement, viewParentCoords);
     }
   }
   catch(hal_exception& e)
@@ -105,7 +110,7 @@ int main(int argc, char** argv)
 }
 
 void extractAlignedRegions(const Genome* genome, ostream* bedStream,
-                           bool complement)
+                           bool complement, bool viewParentCoords)
 {
   assert(genome && bedStream);
   if (genome->getNumTopSegments() > 0)
@@ -122,7 +127,22 @@ void extractAlignedRegions(const Genome* genome, ostream* bedStream,
                    << topSeg->getStartPosition() - sequence->getStartPosition()
                    << '\t' 
                    << (1 + topSeg->getEndPosition() -
-                       sequence->getStartPosition()) << '\n';
+                       sequence->getStartPosition());
+        if (!complement && viewParentCoords) {
+            BottomSegmentIteratorConstPtr botSeg = genome->getParent()->getBottomSegmentIterator();
+            botSeg->toParent(topSeg);
+            if (botSeg->getReversed()) {
+                // Ensure the parent segment coordinates always have start < end.
+                botSeg->toReverse();
+            }
+            const Sequence *parentSequence = botSeg->getSequence();
+            *bedStream << '\t'
+                       << parentSequence->getName() << '\t'
+                       << botSeg->getStartPosition() - parentSequence->getStartPosition() << '\t'
+                       << 1 + botSeg->getEndPosition() - parentSequence->getStartPosition() << '\t'
+                       << (topSeg->getParentReversed() ? "-" : "+");
+        }
+        *bedStream << '\n';
       }
     }
   }