docs: def quantify(**kwargs):
      """
      Quantify expected read counts
  
      :param alnfile: alignment incidence file (h5)
      :param grpfile: gene ID to isoform ID mapping info (tsv)
      :param lenfile: transcript lengths (tsv)
      :param multiread_model: emase model (default: 4)
      :param read_length: read length (default: 100)
      :param pseudocount: prior read count (default: 0.0)
      :param tolerance: tolerance for EM termination (default: 0.0001 in TPM)
      :param max_iters: maximum iterations for EM iteration
      :param report_alignment_counts: whether to report alignment counts (default: False)
      :param report_posterior:
      :return: Expected read counts (tsv)
      """
      alnfile = kwargs.get('alnfile')
      grpfile = kwargs.get('grpfile')
      if grpfile is None:
          grpfile2chk = os.path.join(DATA_DIR, 'ref.gene2transcripts.tsv')
          if os.path.exists(grpfile2chk):
              grpfile = grpfile2chk
          else:
              print >> sys.stderr, '[gbrs::quantify] A group file is not given. Group-level results will not be reported.'
  
      outbase = kwargs.get('outbase')
      gtypefile = kwargs.get('gtypefile')
      pseudocount = kwargs.get('pseudocount')
      lenfile = kwargs.get('lenfile')
      if lenfile is None:
          lenfile2chk = os.path.join(DATA_DIR, 'gbrs.hybridized.targets.info')
          if os.path.exists(lenfile2chk):
              lenfile = lenfile2chk
          else:
              print >> sys.stderr, '[gbrs::quantify] A length file is not given. Transcript length adjustment will *not* be performed.'
  
      read_length = kwargs.get('read_length')
      multiread_model = kwargs.get('multiread_model')
      tolerance = kwargs.get('tolerance')
      max_iters = kwargs.get('max_iters')
      report_group_counts = grpfile is not None  # If grpfile exist, always report groupwise results too
      report_alignment_counts = kwargs.get('report_alignment_counts')
      report_posterior = kwargs.get('report_posterior')
  
      # Load alignment incidence matrix ('alnfile' is assumed to be in multiway transcriptome)
      alnmat = emase.AlignmentPropertyMatrix(h5file=alnfile, grpfile=grpfile)
  
      # Load genotype calls
      if gtypefile is not None:  # Genotype calls are at the gene level
          outbase = outbase + '.diploid'
          hid = dict(zip(alnmat.hname, np.arange(alnmat.num_haplotypes)))
          gid = dict(zip(alnmat.gname, np.arange(len(alnmat.gname))))
          gtmask = np.zeros((alnmat.num_haplotypes, alnmat.num_loci))
          gtcall_g = dict.fromkeys(alnmat.gname)
          gtcall_t = dict.fromkeys(alnmat.lname)
          with open(gtypefile) as fh:
              for curline in dropwhile(is_comment, fh):
                  item = curline.rstrip().split("\t")
                  g, gt = item[:2]
                  gtcall_g[g] = gt
                  hid2set = np.array([hid[c] for c in gt])
                  tid2set = np.array(alnmat.groups[gid[g]])
                  gtmask[np.meshgrid(hid2set, tid2set)] = 1.0
                  for t in tid2set:
                      gtcall_t[alnmat.lname[t]] = gt
          alnmat.multiply(gtmask, axis=2)
          for h in xrange(alnmat.num_haplotypes):
              alnmat.data[h].eliminate_zeros()
      else:
          outbase = outbase + ".multiway"
          gtcall_g = None
          gtcall_t = None
  
      # Run emase
      em_factory = emase.EMfactory(alnmat)
      em_factory.prepare(pseudocount=pseudocount, lenfile=lenfile, read_length=read_length)
      em_factory.run(model=multiread_model, tol=tolerance, max_iters=max_iters, verbose=True)
      em_factory.report_depths(filename="%s.isoforms.tpm" % outbase, tpm=True, notes=gtcall_t)
      em_factory.report_read_counts(filename="%s.isoforms.expected_read_counts" % outbase, notes=gtcall_t)
      if report_posterior:
          em_factory.export_posterior_probability(filename="%s.posterior.h5" % outbase)
      if report_group_counts:
          em_factory.report_depths(filename="%s.genes.tpm" % outbase, tpm=True, grp_wise=True, notes=gtcall_g)
          em_factory.report_read_counts(filename="%s.genes.expected_read_counts" % outbase, grp_wise=True, notes=gtcall_g)
      if report_alignment_counts:
          alnmat = emase.AlignmentPropertyMatrix(h5file=alnfile, grpfile=grpfile)
          alnmat.report_alignment_counts(filename="%s.isoforms.alignment_counts" % outbase)
          if report_group_counts:
              alnmat._bundle_inline(reset=True)
              alnmat.report_alignment_counts(filename="%s.genes.alignment_counts" % outbase)