/tools/maf/maf_stats.py

https://bitbucket.org/cistrome/cistrome-harvard/ · Python · 110 lines · 90 code · 9 blank · 11 comment · 29 complexity · fd1069da7e6fd3dc809ae4343f843f6b MD5 · raw file

  1. #!/usr/bin/env python
  2. #Dan Blankenberg
  3. """
  4. Reads a list of intervals and a maf. Outputs a new set of intervals with statistics appended.
  5. """
  6. import sys
  7. from galaxy import eggs
  8. import pkg_resources; pkg_resources.require( "bx-python" )
  9. import bx.intervals.io
  10. from bx.bitset import BitSet
  11. from galaxy.tools.util import maf_utilities
  12. assert sys.version_info[:2] >= ( 2, 4 )
  13. def __main__():
  14. maf_source_type = sys.argv.pop( 1 )
  15. input_maf_filename = sys.argv[1].strip()
  16. input_interval_filename = sys.argv[2].strip()
  17. output_filename = sys.argv[3].strip()
  18. dbkey = sys.argv[4].strip()
  19. try:
  20. chr_col = int( sys.argv[5].strip() ) - 1
  21. start_col = int( sys.argv[6].strip() ) - 1
  22. end_col = int( sys.argv[7].strip() ) - 1
  23. except:
  24. print >>sys.stderr, "You appear to be missing metadata. You can specify your metadata by clicking on the pencil icon associated with your interval file."
  25. sys.exit()
  26. summary = sys.argv[8].strip()
  27. if summary.lower() == "true": summary = True
  28. else: summary = False
  29. mafIndexFile = "%s/maf_index.loc" % sys.argv[9]
  30. try:
  31. maf_index_filename = sys.argv[10].strip()
  32. except:
  33. maf_index_filename = None
  34. index = index_filename = None
  35. if maf_source_type == "user":
  36. #index maf for use here
  37. index, index_filename = maf_utilities.open_or_build_maf_index( input_maf_filename, maf_index_filename, species = [dbkey] )
  38. if index is None:
  39. print >>sys.stderr, "Your MAF file appears to be malformed."
  40. sys.exit()
  41. elif maf_source_type == "cached":
  42. #access existing indexes
  43. index = maf_utilities.maf_index_by_uid( input_maf_filename, mafIndexFile )
  44. if index is None:
  45. print >> sys.stderr, "The MAF source specified (%s) appears to be invalid." % ( input_maf_filename )
  46. sys.exit()
  47. else:
  48. print >>sys.stdout, 'Invalid source type specified: %s' % maf_source_type
  49. sys.exit()
  50. out = open(output_filename, 'w')
  51. num_region = None
  52. num_bad_region = 0
  53. species_summary = {}
  54. total_length = 0
  55. #loop through interval file
  56. for num_region, region in enumerate( bx.intervals.io.NiceReaderWrapper( open( input_interval_filename, 'r' ), chrom_col = chr_col, start_col = start_col, end_col = end_col, fix_strand = True, return_header = False, return_comments = False ) ):
  57. src = "%s.%s" % ( dbkey, region.chrom )
  58. region_length = region.end - region.start
  59. if region_length < 1:
  60. num_bad_region += 1
  61. continue
  62. total_length += region_length
  63. coverage = { dbkey: BitSet( region_length ) }
  64. for block in index.get_as_iterator( src, region.start, region.end ):
  65. for spec in maf_utilities.get_species_in_block( block ):
  66. if spec not in coverage: coverage[spec] = BitSet( region_length )
  67. for block in maf_utilities.iter_blocks_split_by_species( block ):
  68. if maf_utilities.component_overlaps_region( block.get_component_by_src( src ), region ):
  69. #need to chop and orient the block
  70. block = maf_utilities.orient_block_by_region( maf_utilities.chop_block_by_region( block, src, region ), src, region, force_strand = '+' )
  71. start_offset, alignment = maf_utilities.reduce_block_by_primary_genome( block, dbkey, region.chrom, region.start )
  72. for i in range( len( alignment[dbkey] ) ):
  73. for spec, text in alignment.items():
  74. if text[i] != '-':
  75. coverage[spec].set( start_offset + i )
  76. if summary:
  77. #record summary
  78. for key in coverage.keys():
  79. if key not in species_summary: species_summary[key] = 0
  80. species_summary[key] = species_summary[key] + coverage[key].count_range()
  81. else:
  82. #print coverage for interval
  83. coverage_sum = coverage[dbkey].count_range()
  84. out.write( "%s\t%s\t%s\t%s\n" % ( "\t".join( region.fields ), dbkey, coverage_sum, region_length - coverage_sum ) )
  85. keys = coverage.keys()
  86. keys.remove( dbkey )
  87. keys.sort()
  88. for key in keys:
  89. coverage_sum = coverage[key].count_range()
  90. out.write( "%s\t%s\t%s\t%s\n" % ( "\t".join( region.fields ), key, coverage_sum, region_length - coverage_sum ) )
  91. if summary:
  92. out.write( "#species\tnucleotides\tcoverage\n" )
  93. for spec in species_summary:
  94. out.write( "%s\t%s\t%.4f\n" % ( spec, species_summary[spec], float( species_summary[spec] ) / total_length ) )
  95. out.close()
  96. if num_region is not None:
  97. print "%i regions were processed with a total length of %i." % ( num_region + 1, total_length )
  98. if num_bad_region:
  99. print "%i regions were invalid." % ( num_bad_region )
  100. maf_utilities.remove_temp_index_file( index_filename )
  101. if __name__ == "__main__": __main__()