/tools/new_operations/gops_coverage.py

https://bitbucket.org/cistrome/cistrome-harvard/ · Python · 68 lines · 50 code · 9 blank · 9 comment · 9 complexity · b599448b791491cf8f279659003959eb MD5 · raw file

  1. #!/usr/bin/env python
  2. """
  3. Calculate coverage of one query on another, and append the coverage to
  4. the last two columns as bases covered and percent coverage.
  5. usage: %prog bed_file_1 bed_file_2 out_file
  6. -1, --cols1=N,N,N,N: Columns for start, end, strand in first file
  7. -2, --cols2=N,N,N,N: Columns for start, end, strand in second file
  8. """
  9. from galaxy import eggs
  10. import pkg_resources
  11. pkg_resources.require( "bx-python" )
  12. import sys, traceback, fileinput
  13. from warnings import warn
  14. from bx.intervals import *
  15. from bx.intervals.io import *
  16. from bx.intervals.operations.coverage import *
  17. from bx.cookbook import doc_optparse
  18. from galaxy.tools.util.galaxyops import *
  19. assert sys.version_info[:2] >= ( 2, 4 )
  20. def main():
  21. upstream_pad = 0
  22. downstream_pad = 0
  23. options, args = doc_optparse.parse( __doc__ )
  24. try:
  25. chr_col_1, start_col_1, end_col_1, strand_col_1 = parse_cols_arg( options.cols1 )
  26. chr_col_2, start_col_2, end_col_2, strand_col_2 = parse_cols_arg( options.cols2 )
  27. in_fname, in2_fname, out_fname = args
  28. except:
  29. doc_optparse.exception()
  30. g1 = NiceReaderWrapper( fileinput.FileInput( in_fname ),
  31. chrom_col=chr_col_1,
  32. start_col=start_col_1,
  33. end_col=end_col_1,
  34. strand_col=strand_col_1,
  35. fix_strand=True )
  36. g2 = NiceReaderWrapper( fileinput.FileInput( in2_fname ),
  37. chrom_col=chr_col_2,
  38. start_col=start_col_2,
  39. end_col=end_col_2,
  40. strand_col=strand_col_2,
  41. fix_strand=True )
  42. out_file = open( out_fname, "w" )
  43. try:
  44. for line in coverage( [g1,g2] ):
  45. if type( line ) is GenomicInterval:
  46. out_file.write( "%s\n" % "\t".join( line.fields ) )
  47. else:
  48. out_file.write( "%s\n" % line )
  49. except ParseError, exc:
  50. out_file.close()
  51. fail( "Invalid file format: %s" % str( exc ) )
  52. out_file.close()
  53. if g1.skipped > 0:
  54. print skipped( g1, filedesc=" of 1st dataset" )
  55. if g2.skipped > 0:
  56. print skipped( g2, filedesc=" of 2nd dataset" )
  57. if __name__ == "__main__":
  58. main()