/tools/filters/ucsc_gene_bed_to_exon_bed.py

https://bitbucket.org/cistrome/cistrome-harvard/ · Python · 152 lines · 114 code · 12 blank · 26 comment · 30 complexity · 33abf73608482a1de79d32bc4f4db63c 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, codon, intron, 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', 'intron', 'codon' ), "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. """
  48. # Read table and handle each gene
  49. for line in in_file:
  50. try:
  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. # Determine the subset of the transcribed region we are interested in
  63. if options.region == 'utr3':
  64. if strand == '-': region_start, region_end = tx_start, cds_start
  65. else: region_start, region_end = cds_end, tx_end
  66. elif options.region == 'utr5':
  67. if strand == '-': region_start, region_end = cds_end, tx_end
  68. else: region_start, region_end = tx_start, cds_start
  69. elif options.region == 'coding' or options.region == 'codon':
  70. region_start, region_end = cds_start, cds_end
  71. else:
  72. region_start, region_end = tx_start, tx_end
  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. # options.exons is always TRUE
  76. if options.exons:
  77. exon_starts = map( int, fields[11].rstrip( ',\n' ).split( ',' ) )
  78. exon_starts = map((lambda x: x + tx_start ), exon_starts)
  79. exon_ends = map( int, fields[10].rstrip( ',\n' ).split( ',' ) )
  80. exon_ends = map((lambda x, y: x + y ), exon_starts, exon_ends);
  81. #for Intron regions:
  82. if options.region == 'intron':
  83. i=0
  84. while i < len(exon_starts)-1:
  85. intron_starts = exon_ends[i]
  86. intron_ends = exon_starts[i+1]
  87. if strand: print_tab_sep(out_file, chrom, intron_starts, intron_ends, name, "0", strand )
  88. else: print_tab_sep(out_file, chrom, intron_starts, intron_ends )
  89. i+=1
  90. #for non-intron regions:
  91. else:
  92. for start, end in zip( exon_starts, exon_ends ):
  93. start = max( start, region_start )
  94. end = min( end, region_end )
  95. if start < end:
  96. if options.region == 'codon':
  97. start += (3 - ((start-region_start)%3))%3
  98. c_start = start
  99. while c_start+3 <= end:
  100. if strand:
  101. print_tab_sep(out_file, chrom, c_start, c_start+3, name, "0", strand )
  102. else:
  103. print_tab_sep(out_file, chrom, c_start, c_start+3)
  104. c_start += 3
  105. else:
  106. if strand:
  107. print_tab_sep(out_file, chrom, start, end, name, "0", strand )
  108. else:
  109. print_tab_sep(out_file, chrom, start, end )
  110. """
  111. else:
  112. if options.region == 'codon':
  113. c_start = start
  114. c_end = end
  115. if c_start > c_end:
  116. t = c_start
  117. c_start = c_end
  118. c_end = t
  119. while c_start+3 <= c_end:
  120. if strand:
  121. print_tab_sep(out_file, chrom, c_start, c_start+3, name, "0", strand )
  122. else:
  123. print_tab_sep(out_file, chrom, c_start, c_start+3)
  124. c_start += 3
  125. else:
  126. if strand:
  127. print_tab_sep(out_file, chrom, region_start, region_end, name, "0", strand )
  128. else:
  129. print_tab_sep(out_file, chrom, region_start, region_end )
  130. """
  131. except:
  132. continue
  133. def print_tab_sep(out_file, *args ):
  134. """Print items in `l` to stdout separated by tabs"""
  135. print >>out_file, string.join( [ str( f ) for f in args ], '\t' )
  136. if __name__ == "__main__": main()