code: def count(args):
      
      p = OptionParser(count.__doc__)
      p.add_option("--type", default="exon",
                   help="Only count feature type")
      p.set_cpus(cpus=8)
      opts, args = p.parse_args(args)
  
      if len(args) != 2:
          sys.exit(not p.print_help())
  
      bamfile, gtf = args
      cpus = opts.cpus
      pf = bamfile.split(".")[0]
      countfile = pf + ".count"
      if not need_update(bamfile, countfile):
          return
  
      nsorted = pf + "_nsorted"
      nsortedbam, nsortedsam = nsorted + ".bam", nsorted + ".sam"
      if need_update(bamfile, nsortedsam):
          cmd = "samtools sort -@ {0} -n {1} {2}".format(cpus, bamfile, nsorted)
          sh(cmd)
          cmd = "samtools view -@ {0} -h {1}".format(cpus, nsortedbam)
          sh(cmd, outfile=nsortedsam)
  
      if need_update(nsortedsam, countfile):
          cmd = "htseq-count --stranded=no --minaqual=10"
          cmd += " -t {0}".format(opts.type)
          cmd += " {0} {1}".format(nsortedsam, gtf)
          sh(cmd, outfile=countfile)
docs: %prog count bamfile gtf
  
      Count the number of reads mapped using `htseq-count`.