From 987e2c7c4e5797858a8ee65e16bec009b488ec9c Mon Sep 17 00:00:00 2001
From: Guruprasad Anada <gua110@bx.psu.edu>
Date: Thu, 17 Jan 2008 19:24:34 +0000
Subject: [PATCH] Adding liftOver tool. Test section is commented as the
 functinal test module doesn't handle multiple outputs yet. I've downloaded
 liftOver mapping files from UCSC for all the genomes for which we have
 sequences. I've added a function in get_ucsc_scores.pl file (which contains
 functions to fetch sequences and alignments) to fetch liftOver mappings too.

---
 tool_conf.xml.sample                   |  3 ++
 tools/extract/liftOver_wrapper.py      | 44 ++++++++++++++++++++++
 tools/extract/liftOver_wrapper.xml     | 52 ++++++++++++++++++++++++++
 tools/extract/liftOver_wrapper_code.py |  6 +++
 4 files changed, 105 insertions(+)
 create mode 100644 tools/extract/liftOver_wrapper.py
 create mode 100644 tools/extract/liftOver_wrapper.xml
 create mode 100644 tools/extract/liftOver_wrapper_code.py

diff --git a/tool_conf.xml.sample b/tool_conf.xml.sample
index 912c755c015..75163d8c41d 100644
--- a/tool_conf.xml.sample
+++ b/tool_conf.xml.sample
@@ -27,6 +27,9 @@
     <tool file="encode/gencode_partition.xml" />
     <tool file="encode/random_intervals.xml" />
   </section>
+  <section name="LiftOver" id="liftOver">
+    <tool file="extract/liftOver_wrapper.xml" />
+  </section>
   <section name="Text Manipulation" id="textutil">
     <tool file="filters/fixedValueColumn.xml" />
     <tool file="stats/column_maker.xml" />
diff --git a/tools/extract/liftOver_wrapper.py b/tools/extract/liftOver_wrapper.py
new file mode 100644
index 00000000000..d072101bd21
--- /dev/null
+++ b/tools/extract/liftOver_wrapper.py
@@ -0,0 +1,44 @@
+#!/usr/bin/env python2.4
+#Guru
+"""
+Converts coordinates from one build/assembly to another using liftOver binary and mapping files downloaded from UCSC.
+"""
+
+import sys, os, string
+
+def stop_err(msg):
+    sys.stderr.write(msg)
+    sys.exit()
+    
+if len(sys.argv) != 6:
+    stop_error("USAGE: prog input out_file1 out_file2 input_dbkey output_dbkey")
+
+infile = sys.argv[1]
+outfile1 = sys.argv[2]
+outfile2 = sys.argv[3]
+in_dbkey = sys.argv[4]
+out_dbkey = sys.argv[5]
+
+#ensure dbkeys are set
+if in_dbkey == "?": 
+    stop_err("You must specify a build to the input dataset in order to covert genome coordinates.")
+if out_dbkey == "?": 
+    stop_err("Please specify a build for the output dataset using the 'To' dropdown menu.")
+
+#Check if the apping file exists    
+#example file path: hg18ToHg17.over.chain
+mapfilename = in_dbkey + "To" + out_dbkey[0].capitalize() + out_dbkey[1:] + ".over.chain"
+mapfilepath = "/depot/data2/galaxy/" + in_dbkey + "/liftOver/" + mapfilename
+try:
+    open(mapfilepath, 'r')
+except Exception, ex:
+    stop_err("Mapping information from %s to %s is unavailable." %(in_dbkey, out_dbkey))
+    
+print "Mapping from %s to %s" %(in_dbkey,out_dbkey)
+try:
+    cmd_line = "liftOver " + infile + " " + mapfilepath + " " + outfile1 + " " + outfile2 + "  > /dev/null 2>&1"
+    os.system(cmd_line)
+except Exception, exc:
+    print >>sys.stderr, exc
+  
+    
diff --git a/tools/extract/liftOver_wrapper.xml b/tools/extract/liftOver_wrapper.xml
new file mode 100644
index 00000000000..bc7d4b49470
--- /dev/null
+++ b/tools/extract/liftOver_wrapper.xml
@@ -0,0 +1,52 @@
+<tool id="liftOver1" name="Convert genome coordinates">
+  <description> between assemblies and genomes</description>
+  <command interpreter="python">liftOver_wrapper.py $input "$out_file1" "$out_file2" $dbkey $to_dbkey</command>
+  <inputs>
+    <param format="bed" name="input" type="data" label="Convert coordinates of"/>
+    <param name="to_dbkey" type="genomebuild" label="To" />
+  </inputs>
+  <outputs>
+    <data format="bed" name="out_file1" />
+    <data format="bed" name="out_file2" />
+  </outputs>
+  <!--
+  <tests>
+    <test>
+      <param name="input" value="5.bed" dbkey="hg18" ftype="bed" />
+      <param name="to_dbkey" value="panTro2" />
+      <output name="out_file1" file="5_liftover_mapped.bed"/>
+      <output name="out_file2" file="5_liftover_unmapped.bed"/>
+    </test>
+  </tests>
+  -->
+  <help>
+
+.. class:: warningmark
+
+Make sure that the genome build is specified for the input dataset (click the pencil icon if it is not specified).
+
+-----
+
+**What it does**
+
+This tool converts genome coordinates and annotation files between assemblies and genomes. It outputs 2 files, one containing all the mapped coordinates and one containing the unmapped coordinates (if any). 
+
+-----
+
+**Example**
+
+Running this tool to convert the following hg16 intervals into hg18 intervals::
+
+    chrX  85170   112199  AK002185  0  +
+    chrX  110458  112199  AK097346  0  +
+    chrX  112203  121212  AK074528  0  -
+
+will return the following hg18 intervals::
+
+    chrX  132991  160020  AK002185  0  +
+    chrX  158279  160020  AK097346  0  +
+    chrX  160024  169033  AK074528  0  -
+
+</help>
+<code file="liftOver_wrapper_code.py"/>
+</tool>
diff --git a/tools/extract/liftOver_wrapper_code.py b/tools/extract/liftOver_wrapper_code.py
new file mode 100644
index 00000000000..0d7ec90629a
--- /dev/null
+++ b/tools/extract/liftOver_wrapper_code.py
@@ -0,0 +1,6 @@
+def exec_before_job(app, inp_data, out_data, param_dict, tool):
+    out_data['out_file1'].set_dbkey(param_dict['to_dbkey'])
+    out_data['out_file2'].set_dbkey(param_dict['to_dbkey'])
+    out_data['out_file1'].name = out_data['out_file1'].name + " [ MAPPED COORDINATES ]"
+    out_data['out_file2'].name = out_data['out_file2'].name + " [ UNMAPPED COORDINATES ]"
+