Mercurial > repos > xuebing > sharplabtool
diff tools/new_operations/gops_complement.py @ 0:9071e359b9a3
Uploaded
author | xuebing |
---|---|
date | Fri, 09 Mar 2012 19:37:19 -0500 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tools/new_operations/gops_complement.py Fri Mar 09 19:37:19 2012 -0500 @@ -0,0 +1,98 @@ +#!/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) +""" +from galaxy import eggs +import pkg_resources +pkg_resources.require( "bx-python" ) +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()