hack in halGetMaf to the C interface in chains

ComparativeGenomicsToolkit · Aug 13, 2013 · 109b87b · 109b87b
1 parent f17ac47
commit 109b87b
Show file tree

Hide file tree

Showing 4 changed files with 124 additions and 15 deletions.
diff --git a/Makefile b/Makefile
@@ -1,5 +1,5 @@
 # order is important, libraries first
-modules = api stats randgen validate mutations fasta alignability lod chain liftover maf extract analysis phyloP
+modules = api stats randgen validate mutations fasta maf alignability lod chain liftover extract analysis phyloP
 
 .PHONY: all %.all clean %.clean doxy %.doxy
 

diff --git a/chain/Makefile b/chain/Makefile
@@ -12,12 +12,12 @@ libTestsHeaders = $(wildcard test/*.h)
 libHalTestsAll := $(wildcard ../api/tests/*.cpp)
 libHalTests = $(subst ../api/tests/allTests.cpp,,${libHalTestsAll})
 
-all : ${libPath}/halChain.a ${binPath}/hal2chain test/blockVizTest test/blockVizBed test/blockInterpolateTest ${binPath}/halChainTests ${binPath}/hal2assemblyHub.py
+all : ${libPath}/halChain.a ${binPath}/hal2chain test/blockVizTest test/blockVizBed test/blockVizMaf test/blockInterpolateTest ${binPath}/halChainTests ${binPath}/hal2assemblyHub.py
 
 clean : 
-	rm -f ${libPath}/halChain.a ${libPath}/*.h ${binPath}/hal2chain test/blockVizTest test/BlockVizBed test/blockInterpolateTest ${binPath}/hal2assemblyHub.py
+	rm -f ${libPath}/halChain.a ${libPath}/*.h ${binPath}/hal2chain test/blockVizTest test/BlockVizBed test/blockVizMaf test/blockInterpolateTest ${binPath}/hal2assemblyHub.py
 
-${libPath}/halChain.a : ${libSources} ${libHeaders} ${libPath}/halLib.a ${basicLibsDependencies} 
+${libPath}/halChain.a : ${libSources} ${libHeaders} ${libPath}/halLib.a ${libPath}/halMaf.a ${basicLibsDependencies} 
 	cp ${libHeaders} ${libPath}/
 	rm -f *.o
 	${cpp} ${cppflags} -I inc -I hdf5_impl -I impl -I ${libPath}/ -c ${libSources}
@@ -26,23 +26,27 @@ ${libPath}/halChain.a : ${libSources} ${libHeaders} ${libPath}/halLib.a ${basicL
 	rm *.o
 	mv halChain.a ${libPath}/
 
-${binPath}/hal2chain : impl/hal2chain.cpp ${libPath}/halChain.a ${libPath}/halLib.a ${basicLibsDependencies}
-	${cpp} ${cppflags} -I inc -I impl -I ${libPath} -I impl -I tests -o ${binPath}/hal2chain impl/hal2chain.cpp ${libPath}/halChain.a ${libPath}/halLib.a ${basicLibs}
+${binPath}/hal2chain : impl/hal2chain.cpp ${libPath}/halMaf.a ${libPath}/halChain.a ${libPath}/halLib.a ${basicLibsDependencies}
+	${cpp} ${cppflags} -I inc -I impl -I ${libPath} -I impl -I tests -o ${binPath}/hal2chain impl/hal2chain.cpp ${libPath}/halMaf.a ${libPath}/halChain.a ${libPath}/halLib.a ${basicLibs}
 
-test/blockVizTest: test/blockVizTest.c ${libPath}/halChain.a ${libPath}/halLod.a ${libPath}/halLib.a ${basicLib\
+test/blockVizTest: test/blockVizTest.c ${libPath}/halMaf.a ${libPath}/halChain.a ${libPath}/halLod.a ${libPath}/halLib.a ${basicLib\
 sDependencies}
-	${cpp} ${cppflags} -I inc -I impl -I ${libPath} -I impl -I tests -o test/blockVizTest test/blockVizTest.c ${libPath}/halChain.a ${libPath}/halLod.a ${libPath}/halLib.a ${basicLibs}
+	${cpp} ${cppflags} -I inc -I impl -I ${libPath} -I impl -I tests -o test/blockVizTest test/blockVizTest.c ${libPath}/halMaf.a ${libPath}/halChain.a ${libPath}/halLod.a ${libPath}/halLib.a ${basicLibs}
 
-test/blockVizBed: test/blockVizBed.c ${libPath}/halChain.a ${libPath}/halLod.a ${libPath}/halLib.a ${basicLib\
+test/blockVizMaf: test/blockVizMaf.c ${libPath}/halMaf.a ${libPath}/halChain.a ${libPath}/halLod.a ${libPath}/halLib.a ${basicLib\
 sDependencies}
-	${cpp} ${cppflags} -I inc -I impl -I ${libPath} -I impl -I tests -o test/blockVizBed test/blockVizBed.c ${libPath}/halChain.a ${libPath}/halLod.a ${libPath}/halLib.a ${basicLibs}
+	${cpp} ${cppflags} -I inc -I impl -I ${libPath} -I impl -I tests -o test/blockVizMaf test/blockVizMaf.c ${libPath}/halMaf.a ${libPath}/halChain.a ${libPath}/halLod.a ${libPath}/halLib.a ${basicLibs}
+
+test/blockVizBed: test/blockVizBed.c ${libPath}/halMaf.a ${libPath}/halChain.a ${libPath}/halLod.a ${libPath}/halLib.a ${basicLib\
+sDependencies}
+	${cpp} ${cppflags} -I inc -I impl -I ${libPath} -I impl -I tests -o test/blockVizBed test/blockVizBed.c ${libPath}/halMaf.a ${libPath}/halChain.a ${libPath}/halLod.a ${libPath}/halLib.a ${basicLibs}
 
-test/blockInterpolateTest: test/blockInterpolateTest.cpp ${libPath}/halChain.a ${libPath}/halLib.a ${basicLib\
+test/blockInterpolateTest: test/blockInterpolateTest.cpp ${libPath}/halMaf.a ${libPath}/halChain.a ${libPath}/halLib.a ${basicLib\
 sDependencies}
-	${cpp} ${cppflags} -I inc -I impl -I ${libPath} -I impl -I tests -o test/blockInterpolateTest test/blockInterpolateTest.cpp ${libPath}/halChain.a ${libPath}/halLib.a ${basicLibs}
+	${cpp} ${cppflags} -I inc -I impl -I ${libPath} -I impl -I tests -o test/blockInterpolateTest test/blockInterpolateTest.cpp ${libPath}/halMaf.a ${libPath}/halChain.a ${libPath}/halLib.a ${basicLibs}
 
 ${binPath}/halChainTests : ${libTests} ${libTestsHeaders} ${libTestsCommon} ${libTestsHeadersCommon} ${libSources} ${libHeaders} ${libInternalHeaders} ${libPath}/halLib.a ${basicLibsDependencies}
-	${cpp} ${cppflags} -I inc -I impl -I ${libPath} -I tests -I ../api/tests -o ${binPath}/halChainTests ${libHalTests} ${libTests} ${libPath}/halLib.a ${libPath}/halChain.a ${basicLibs}
+	${cpp} ${cppflags} -I inc -I impl -I ${libPath} -I tests -I ../api/tests -o ${binPath}/halChainTests ${libHalTests} ${libTests} ${libPath}/halLib.a ${libPath}/halMaf.a ${libPath}/halChain.a ${basicLibs}
 
 ${binPath}/hal2assemblyHub.py : hal2assemblyHub.py
 	cp hal2assemblyHub.py ${binPath}/hal2assemblyHub.py

diff --git a/chain/impl/halBlockViz.cpp b/chain/impl/halBlockViz.cpp
@@ -18,6 +18,7 @@
 #include "halBlockViz.h"
 #include "halBlockMapper.h"
 #include "halLodManager.h"
+#include "halMafExport.h"
 
 #ifdef ENABLE_UDC
 #include <pthread.h>
@@ -277,6 +278,86 @@ struct hal_block_results_t *halGetBlocksInTargetRange(int halHandle,
   return results;
 }
 
+extern "C" hal_int_t halGetMAF(FILE* outFile,
+                               int halHandle, 
+                               hal_species_t* qSpeciesNames,
+                               char* tSpecies,
+                               char* tChrom,
+                               hal_int_t tStart, 
+                               hal_int_t tEnd,
+                               int doDupes)
+{
+  HAL_LOCK
+  hal_int_t numBytes = 0;
+  try
+  {
+    hal_int_t rangeLength = tEnd - tStart;
+    if (rangeLength < 0)
+    {
+      stringstream ss;
+      ss << "Invalid query range [" << tStart << "," << tEnd << ").";
+      throw hal_exception(ss.str());
+    }
+    AlignmentConstPtr alignment = 
+       getExistingAlignment(halHandle, hal_size_t(0), true);
+
+    set<const Genome*> qGenomeSet;
+    for (hal_species_t* qSpecies = qSpeciesNames; qSpecies != NULL;
+         qSpecies = qSpecies->next)
+    {
+      checkGenomes(halHandle, alignment, qSpecies->name, tSpecies, tChrom);
+      const Genome* qGenome = alignment->openGenome(qSpecies->name);
+      qGenomeSet.insert(qGenome);
+    }
+
+    const Genome* tGenome = alignment->openGenome(tSpecies);
+    const Sequence* tSequence = tGenome->getSequence(tChrom);
+
+    hal_index_t myEnd = tEnd > 0 ? tEnd : tSequence->getSequenceLength();
+    hal_index_t absStart = tSequence->getStartPosition() + tStart;
+    hal_index_t absEnd = tSequence->getStartPosition() + myEnd - 1;
+    if (absStart > absEnd)
+    {
+      throw hal_exception("Invalid range");
+    }
+    if (absEnd > tSequence->getEndPosition())
+    {
+      throw hal_exception("Target end position outside of target sequence");
+    }
+
+    stringstream mafBuffer;
+    MafExport mafExport;
+    mafExport.setNoDupes(doDupes == 0);
+    mafExport.setUcscNames(true);
+    mafExport.convertSegmentedSequence(mafBuffer, alignment, tGenome, 
+                                       absStart, 1 + absEnd - absStart,
+                                       qGenomeSet);
+    // if these buffers are very big, it is my intuition that 
+    // chopping them up and, saw, only converting and writing a few kb
+    // at a time would be more efficient... but i think that's the least
+    // of our problems right now.  
+    string mafStringBuffer = mafBuffer.str();
+    if (mafStringBuffer.empty() == false)
+    {
+      numBytes = (hal_int_t)fwrite(mafStringBuffer.c_str(), 
+                                   mafStringBuffer.length(), 
+                                   sizeof(char), outFile);
+    }                 
+  }
+  catch(exception& e)
+  {
+    cerr << "Exception caught: " << e.what() << endl;
+    numBytes = -1;
+  }
+  catch(...)
+  {
+    cerr << "Error in hal maf query";
+    numBytes = -1;
+  }
+  return numBytes;
+  HAL_UNLOCK
+}
+
 extern "C" struct hal_species_t *halGetSpecies(int halHandle)
 {
   HAL_LOCK
@@ -858,7 +939,7 @@ void cleanTargetDupesList(vector<hal_target_dupe_list_t*>& dupeList)
   dupeList.resize(dupeList.size() - deadCount);
 }
 
-static void flipStrand(hal_block_t* firstBlock)
+void flipStrand(hal_block_t* firstBlock)
 {
   // below is stupid and naive and will not work for anything but 
   // dense mode (ie where no edges between blocks).  put off fixing

diff --git a/chain/inc/halBlockViz.h b/chain/inc/halBlockViz.h
@@ -152,7 +152,31 @@ struct hal_block_results_t *halGetBlocksInTargetRange(int halHandle,
                                                       hal_int_t tEnd,
                                                       int getSequenceString,
                                                       int doDupes);
-
+
+/** Read alignment into an output file in MAF format.  Interface very 
+ * similar to halGetBlocksInTargetRange except multiple query species 
+ * can be specified
+ *
+ * @param halHandle handle for the HAL alignment obtained from halOpen
+ * @param qSpeciesNames the names of the query species (no other information
+ * is read from the hal_species_t structure -- just the names).
+ * @param tSpecies the name of the reference species.
+ * @param tChrom name of the chromosome in reference.
+ * @param tStart start position in reference  
+ * @param tEnd last + 1 position in reference (if 0, then the size of the 
+ * chromosome is used). 
+ * @param doDupes create blocks for duplications if not 0.  When this 
+ * option is enabled, the same region can appear in more than one block.
+ * @return  number of bytes written
+ */
+hal_int_t halGetMAF(FILE* outFile,
+                    int halHandle, 
+                    hal_species_t* qSpeciesNames,
+                    char* tSpecies,
+                    char* tChrom,
+                    hal_int_t tStart, 
+                    hal_int_t tEnd,
+                    int doDupes); 
 
 /** Create a linked list of the species in the hal file.
  * @param halHandle handle for the HAL alignment obtained from halOpen