diff gops_complement.py @ 0:d958d5a0d1e8

Imported from capsule None

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/gops_complement.py	Tue Apr 01 10:51:48 2014 -0400
@@ -0,0 +1,96 @@
+#!/usr/bin/env python
+"""
+Complement regions.
+
+usage: %prog in_file out_file
+    -1, --cols1=N,N,N,N: Columns for chrom, start, end, strand in file
+    -l, --lengths=N: Filename of .len file for species (chromosome lengths)
+    -a, --all: Complement all chromosomes (Genome-wide complement)
+"""
+
+import sys, traceback, fileinput
+from warnings import warn
+from bx.intervals import *
+from bx.intervals.io import *
+from bx.intervals.operations.complement import complement
+from bx.intervals.operations.subtract import subtract
+from bx.cookbook import doc_optparse
+from galaxy.tools.util.galaxyops import *
+
+assert sys.version_info[:2] >= ( 2, 4 )
+
+def main():
+    allchroms = False
+    upstream_pad = 0
+    downstream_pad = 0
+
+    options, args = doc_optparse.parse( __doc__ )
+    try:
+        chr_col_1, start_col_1, end_col_1, strand_col_1 = parse_cols_arg( options.cols1 )
+        lengths = options.lengths
+        if options.all: allchroms = True
+        in_fname, out_fname = args
+    except:
+        doc_optparse.exception()
+
+    g1 = NiceReaderWrapper( fileinput.FileInput( in_fname ),
+                            chrom_col=chr_col_1,
+                            start_col=start_col_1,
+                            end_col=end_col_1,
+                            strand_col=strand_col_1,
+                            fix_strand=True )
+
+    lens = dict()
+    chroms = list()
+    # dbfile is used to determine the length of each chromosome.  The lengths
+    # are added to the lens dict and passed copmlement operation code in bx.
+    dbfile = fileinput.FileInput( lengths )
+    
+    if dbfile:
+        if not allchroms:
+            try:
+                for line in dbfile:
+                    fields = line.split("\t")
+                    lens[fields[0]] = int(fields[1])
+            except:
+                # assume LEN doesn't exist or is corrupt somehow
+                pass
+        elif allchroms:
+            try:
+                for line in dbfile:
+                    fields = line.split("\t")
+                    end = int(fields[1])
+                    chroms.append("\t".join([fields[0],"0",str(end)]))
+            except:
+                pass
+
+    # Safety...if the dbfile didn't exist and we're on allchroms, then
+    # default to generic complement
+    if allchroms and len(chroms) == 0:
+        allchroms = False
+
+    if allchroms:
+        chromReader = GenomicIntervalReader(chroms)
+        generator = subtract([chromReader, g1])
+    else:
+        generator = complement(g1, lens)
+
+    out_file = open( out_fname, "w" )
+
+    try:
+        for interval in generator:
+            if type( interval ) is GenomicInterval:
+                out_file.write( "%s\n" % "\t".join( interval ) )
+            else:
+                out_file.write( "%s\n" % interval )
+    except ParseError, exc:
+        out_file.close()
+        fail( "Invalid file format: %s" % str( exc ) )
+
+    out_file.close()
+
+    if g1.skipped > 0:
+        print skipped( g1, filedesc="" )
+
+if __name__ == "__main__":
+    main()