code: def coverage(args):
      
      p = OptionParser(coverage.__doc__)
      p.add_option("--format", default="bigwig",
                   choices=("bedgraph", "bigwig", "coverage"),
                   help="Output format")
      p.add_option("--nosort", default=False, action="store_true",
                   help="Do not sort BAM")
      p.set_outfile()
      opts, args = p.parse_args(args)
  
      if len(args) != 2:
          sys.exit(not p.print_help())
  
      fastafile, bamfile = args
      format = opts.format
      if opts.nosort:
          logging.debug("BAM sorting skipped")
      else:
          bamfile = index([bamfile, "--fasta={0}".format(fastafile)])
  
      pf = bamfile.rsplit(".", 2)[0]
      sizesfile = Sizes(fastafile).filename
      cmd = "genomeCoverageBed -ibam {0} -g {1}".format(bamfile, sizesfile)
      if format in ("bedgraph", "bigwig"):
          cmd += " -bg"
          bedgraphfile = pf + ".bedgraph"
          sh(cmd, outfile=bedgraphfile)
  
          if format == "bedgraph":
              return bedgraphfile
  
          bigwigfile = pf + ".bigwig"
          cmd = "bedGraphToBigWig {0} {1} {2}".\
                      format(bedgraphfile, sizesfile, bigwigfile)
          sh(cmd)
          return bigwigfile
  
      coveragefile = pf + ".coverage"
      if need_update(fastafile, coveragefile):
          sh(cmd, outfile=coveragefile)
  
      gcf = GenomeCoverageFile(coveragefile)
      fw = must_open(opts.outfile, "w")
      for seqid, cov in gcf.iter_coverage_seqid():
          print("\t".join((seqid, "{0:.1f}".format(cov))), file=fw)
      fw.close()
docs: %prog coverage fastafile bamfile
  
      Calculate coverage for BAM file. BAM file will be sorted unless with
      --nosort.