PageRenderTime 33ms CodeModel.GetById 15ms app.highlight 14ms RepoModel.GetById 1ms app.codeStats 0ms

/tools/filters/ucsc_gene_table_to_intervals.py

https://bitbucket.org/cistrome/cistrome-harvard/
Python | 106 lines | 71 code | 13 blank | 22 comment | 17 complexity | 6026a813e2982105c780e308163b6f69 MD5 | raw file
  1#!/usr/bin/env python
  2
  3"""
  4Read a table dump in the UCSC gene table format and print a tab separated
  5list of intervals corresponding to requested features of each gene.
  6
  7usage: ucsc_gene_table_to_intervals.py [options]
  8
  9options:
 10  -h, --help                  show this help message and exit
 11  -rREGION, --region=REGION
 12                              Limit to region: one of coding, utr3, utr5, transcribed [default]
 13  -e, --exons                 Only print intervals overlapping an exon
 14  -i, --input=inputfile       input file
 15  -o, --output=outputfile     output file
 16"""
 17
 18import optparse, string, sys
 19
 20assert sys.version_info[:2] >= ( 2, 4 )
 21
 22def main():
 23
 24    # Parse command line    
 25    parser = optparse.OptionParser( usage="%prog [options] " )
 26    parser.add_option( "-r", "--region", dest="region", default="transcribed",
 27                       help="Limit to region: one of coding, utr3, utr5, transcribed [default]" )
 28    parser.add_option( "-e", "--exons",  action="store_true", dest="exons",
 29                       help="Only print intervals overlapping an exon" )
 30    parser.add_option( "-s", "--strand",  action="store_true", dest="strand",
 31                       help="Print strand after interval" )
 32    parser.add_option( "-i", "--input",  dest="input",  default=None,
 33                       help="Input file" )
 34    parser.add_option( "-o", "--output", dest="output", default=None,
 35                       help="Output file" )
 36    options, args = parser.parse_args()
 37    assert options.region in ( 'coding', 'utr3', 'utr5', 'transcribed' ), "Invalid region argument"
 38    
 39    try:
 40        out_file = open (options.output,"w")
 41    except:
 42        print >> sys.stderr, "Bad output file."
 43        sys.exit(0)
 44    
 45    try:
 46        in_file = open (options.input)
 47    except:
 48        print >> sys.stderr, "Bad input file."
 49        sys.exit(0)
 50    
 51    print "Region:", options.region+";"
 52    print "Only overlap with Exons:",
 53    if options.exons:
 54        print "Yes"
 55    else:
 56        print "No"
 57    
 58    # Read table and handle each gene
 59    for line in in_file:
 60        try:
 61            if line[0:1] == "#":
 62                continue
 63            # Parse fields from gene tabls
 64            fields = line.split( '\t' )
 65            name = fields[0]
 66            chrom = fields[1]
 67            strand = fields[2].replace(" ","_")
 68            tx_start = int( fields[3] )
 69            tx_end = int( fields[4] )
 70            cds_start = int( fields[5] )
 71            cds_end = int( fields[6] )
 72
 73            # Determine the subset of the transcribed region we are interested in
 74            if options.region == 'utr3':
 75                if strand == '-': region_start, region_end = tx_start, cds_start
 76                else: region_start, region_end = cds_end, tx_end 
 77            elif options.region == 'utr5':
 78                if strand == '-': region_start, region_end = cds_end, tx_end
 79                else: region_start, region_end = tx_start, cds_start
 80            elif options.region == 'coding':
 81                region_start, region_end = cds_start, cds_end
 82            else:
 83                region_start, region_end = tx_start, tx_end
 84
 85            # If only interested in exons, print the portion of each exon overlapping
 86            # the region of interest, otherwise print the span of the region
 87            if options.exons:
 88                exon_starts = map( int, fields[8].rstrip( ',\n' ).split( ',' ) )
 89                exon_ends = map( int, fields[9].rstrip( ',\n' ).split( ',' ) )
 90                for start, end in zip( exon_starts, exon_ends ):
 91                    start = max( start, region_start )
 92                    end = min( end, region_end )
 93                    if start < end:
 94                        if strand: print_tab_sep(out_file, chrom, start, end, name, "0", strand )
 95                        else: print_tab_sep(out_file, chrom, start, end )
 96            else:
 97                if strand: print_tab_sep(out_file, chrom, region_start, region_end, name, "0", strand )
 98                else: print_tab_sep(out_file, chrom, region_start, region_end )
 99        except:
100            continue
101
102def print_tab_sep(out_file, *args ):
103    """Print items in `l` to stdout separated by tabs"""
104    print >>out_file, string.join( [ str( f ) for f in args ], '\t' )
105
106if __name__ == "__main__": main()