/tools/filters/ucsc_gene_bed_to_intron_bed.py

https://bitbucket.org/cistrome/cistrome-harvard/ · Python · 98 lines · 49 code · 16 blank · 33 comment · 10 complexity · b77a60199c29e28eeb0da5b070d84ba9 MD5 · raw file

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