/tools/filters/ucsc_gene_bed_to_exon_bed.py

  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, codon, intron, 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', 'intron', 'codon' ), "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    
 59    # Read table and handle each gene
 60    for line in in_file:
 61        try:
 62            if line[0:1] == "#":
 63                continue
 64            # Parse fields from gene tabls
 65            fields = line.split( '\t' )
 66            chrom     = fields[0]
 67            tx_start  = int( fields[1] )
 68            tx_end    = int( fields[2] )
 69            name      = fields[3]
 70            strand    = fields[5].replace(" ","_")
 71            cds_start = int( fields[6] )
 72            cds_end   = int( fields[7] )
 73
 74            # Determine the subset of the transcribed region we are interested in
 75            if options.region == 'utr3':
 76                if strand == '-': region_start, region_end = tx_start, cds_start
 77                else: region_start, region_end = cds_end, tx_end 
 78            elif options.region == 'utr5':
 79                if strand == '-': region_start, region_end = cds_end, tx_end
 80                else: region_start, region_end = tx_start, cds_start
 81            elif options.region == 'coding' or options.region == 'codon':
 82                region_start, region_end = cds_start, cds_end
 83            else:
 84                region_start, region_end = tx_start, tx_end
 85
 86            # If only interested in exons, print the portion of each exon overlapping
 87            # the region of interest, otherwise print the span of the region
 88        # options.exons is always TRUE
 89            if options.exons:
 90                exon_starts = map( int, fields[11].rstrip( ',\n' ).split( ',' ) )
 91                exon_starts = map((lambda x: x + tx_start ), exon_starts)
 92                exon_ends = map( int, fields[10].rstrip( ',\n' ).split( ',' ) )
 93                exon_ends = map((lambda x, y: x + y ), exon_starts, exon_ends);
 94
 95        #for Intron regions:
 96            if options.region == 'intron':
 97                i=0
 98                while i < len(exon_starts)-1:
 99                    intron_starts = exon_ends[i]
100                    intron_ends = exon_starts[i+1]
101                    if strand: print_tab_sep(out_file, chrom, intron_starts, intron_ends, name, "0", strand )
102                    else: print_tab_sep(out_file, chrom, intron_starts, intron_ends )
103                    i+=1
104        #for non-intron regions:
105            else:
106                for start, end in zip( exon_starts, exon_ends ):
107                    start = max( start, region_start )
108                    end = min( end, region_end )
109                    if start < end:
110                        if options.region == 'codon':
111                            start += (3 - ((start-region_start)%3))%3
112                            c_start = start 
113                            while c_start+3 <= end:
114                                if strand:
115                                    print_tab_sep(out_file, chrom, c_start, c_start+3, name, "0", strand )
116                                else:
117                                    print_tab_sep(out_file, chrom, c_start, c_start+3)
118                                c_start += 3
119                        else:
120                            if strand:
121                                print_tab_sep(out_file, chrom, start, end, name, "0", strand )
122                            else: 
123                                print_tab_sep(out_file, chrom, start, end )
124                    """
125                    else:
126                        if options.region == 'codon':
127                            c_start = start
128                            c_end = end
129                            if c_start > c_end:
130                                t = c_start
131                                c_start = c_end
132                                c_end = t
133                            while c_start+3 <= c_end:
134                                if strand:
135                                    print_tab_sep(out_file, chrom, c_start, c_start+3, name, "0", strand )
136                                else:
137                                    print_tab_sep(out_file, chrom, c_start, c_start+3)
138                                c_start += 3
139                        else:
140                            if strand:
141                                print_tab_sep(out_file, chrom, region_start, region_end, name, "0", strand )
142                            else: 
143                                print_tab_sep(out_file, chrom, region_start, region_end )
144                    """
145        except:
146            continue
147
148def print_tab_sep(out_file, *args ):
149    """Print items in `l` to stdout separated by tabs"""
150    print >>out_file, string.join( [ str( f ) for f in args ], '\t' )
151
152if __name__ == "__main__": main()