docs: def tagcount(sam, out, genemap, output_evidence_table, positional, minevidence,
               cb_histogram, cb_cutoff, no_scale_evidence, subsample, sparse,
               parse_tags, gene_tags):
      ''' Count up evidence for tagged molecules
      '''
      from pysam import AlignmentFile
  
      from io import StringIO
      import pandas as pd
  
      from utils import weigh_evidence
  
      logger.info('Reading optional files')
  
      gene_map = None
      if genemap:
          with open(genemap) as fh:
              try:
                  gene_map = dict(p.strip().split() for p in fh)
              except ValueError:
                  logger.error('Incorrectly formatted gene_map, need to be tsv.')
                  sys.exit()
  
      if positional:
          tuple_template = '{0},{1},{2},{3}'
      else:
          tuple_template = '{0},{1},{3}'
  
      if not cb_cutoff:
          cb_cutoff = 0
  
      if cb_histogram and cb_cutoff == "auto":
          cb_cutoff = guess_depth_cutoff(cb_histogram)
  
      cb_cutoff = int(cb_cutoff)
  
      cb_hist = None
      filter_cb = False
      if cb_histogram:
          cb_hist = pd.read_csv(cb_histogram, index_col=0, header=-1, squeeze=True, sep="\t")
          total_num_cbs = cb_hist.shape[0]
          cb_hist = cb_hist[cb_hist > cb_cutoff]
          logger.info('Keeping {} out of {} cellular barcodes.'.format(cb_hist.shape[0], total_num_cbs))
          filter_cb = True
  
      parser_re = re.compile('.*:CELL_(?P<CB>.*):UMI_(?P<MB>.*)')
  
      if subsample:
          logger.info('Creating reservoir of subsampled reads ({} per cell)'.format(subsample))
          start_sampling  = time.time()
  
          reservoir = collections.defaultdict(list)
          cb_hist_sampled = 0 * cb_hist
          cb_obs = 0 * cb_hist
  
          track = stream_bamfile(sam)
          current_read = 'none_observed_yet'
          for i, aln in enumerate(track):
              if aln.qname == current_read:
                  continue
  
              current_read = aln.qname
  
              if parse_tags:
                  CB = aln.get_tag('CR')
              else:
                  match = parser_re.match(aln.qname)
                  CB = match.group('CB')
  
              if CB not in cb_hist.index:
                  continue
  
              cb_obs[CB] += 1
              if len(reservoir[CB]) < subsample:
                  reservoir[CB].append(i)
                  cb_hist_sampled[CB] += 1
              else:
                  s = pd.np.random.randint(0, cb_obs[CB])
                  if s < subsample:
                      reservoir[CB][s] = i
  
          index_filter = set(itertools.chain.from_iterable(reservoir.values()))
          sam_file.close()
          sampling_time = time.time() - start_sampling
          logger.info('Sampling done - {:.3}s'.format(sampling_time))
  
      evidence = collections.defaultdict(int)
  
      logger.info('Tallying evidence')
      start_tally = time.time()
  
      sam_mode = 'r' if sam.endswith(".sam") else 'rb'
      sam_file = AlignmentFile(sam, mode=sam_mode)
      targets = [x["SN"] for x in sam_file.header["SQ"]]
      track = sam_file.fetch(until_eof=True)
      count = 0
      unmapped = 0
      kept = 0
      nomatchcb = 0
      current_read = 'none_observed_yet'
      count_this_read = True
      missing_transcripts = set()
      for i, aln in enumerate(track):
          if count and not count % 1000000:
              logger.info("Processed %d alignments, kept %d." % (count, kept))
              logger.info("%d were filtered for being unmapped." % unmapped)
              if filter_cb:
                  logger.info("%d were filtered for not matching known barcodes."
                              % nomatchcb)
          count += 1
  
          if aln.is_unmapped:
              unmapped += 1
              continue
  
          if gene_tags and not aln.has_tag('GX'):
              unmapped += 1
              continue
  
          if aln.qname != current_read:
              current_read = aln.qname
              if subsample and i not in index_filter:
                  count_this_read = False
                  continue
              else:
                  count_this_read = True
          else:
              if not count_this_read:
                  continue
  
          if parse_tags:
              CB = aln.get_tag('CR')
          else:
              match = parser_re.match(aln.qname)
              CB = match.group('CB')
  
          if filter_cb:
              if CB not in cb_hist.index:
                  nomatchcb += 1
                  continue
  
          if parse_tags:
              MB = aln.get_tag('UM')
          else:
              MB = match.group('MB')
  
          if gene_tags:
              target_name = aln.get_tag('GX').split(',')[0]
          else:
              txid = sam_file.getrname(aln.reference_id)
              if gene_map:
                  if txid in gene_map:
                      target_name = gene_map[txid]
                  else:
                      missing_transcripts.add(txid)
                      target_name = txid
              else:
                  target_name = txid
  
          e_tuple = tuple_template.format(CB, target_name, aln.pos, MB)
  
          # Scale evidence by number of hits
          if no_scale_evidence:
              evidence[e_tuple] += 1.0
          else:
              evidence[e_tuple] += weigh_evidence(aln.tags)
  
          kept += 1
  
      tally_time = time.time() - start_tally
      if missing_transcripts:
          logger.warn('The following transcripts were missing gene_ids, so we added them as the transcript ids: %s' % str(missing_transcripts))
      logger.info('Tally done - {:.3}s, {:,} alns/min'.format(tally_time, int(60. * count / tally_time)))
      logger.info('Collapsing evidence')
  
      logger.info('Writing evidence')
      with tempfile.NamedTemporaryFile('w+t') as out_handle:
          for key in evidence:
              line = '{},{}\n'.format(key, evidence[key])
              out_handle.write(line)
  
          out_handle.flush()
          out_handle.seek(0)
  
          evidence_table = pd.read_csv(out_handle, header=None)
  
      del evidence
  
      evidence_query = 'evidence >= %f' % minevidence
      if positional:
          evidence_table.columns=['cell', 'gene', 'umi', 'pos', 'evidence']
          collapsed = evidence_table.query(evidence_query).groupby(['cell', 'gene'])['umi', 'pos'].size()
  
      else:
          evidence_table.columns=['cell', 'gene', 'umi', 'evidence']
          collapsed = evidence_table.query(evidence_query).groupby(['cell', 'gene'])['umi'].size()
  
      expanded = collapsed.unstack().T
  
      if gene_map:
          # This Series is just for sorting the index
          genes = pd.Series(index=set(gene_map.values()))
          genes = genes.sort_index()
          # Now genes is assigned to a DataFrame
          genes = expanded.ix[genes.index]
  
      elif gene_tags:
          expanded.sort_index()
          genes = expanded
  
      else:
          # make data frame have a complete accounting of transcripts
          targets = pd.Series(index=set(targets))
          targets = targets.sort_index()
          expanded = expanded.reindex(targets.index.values, fill_value=0)
          genes = expanded
  
      genes.fillna(0, inplace=True)
      genes = genes.astype(int)
      genes.index.name = "gene"
  
      logger.info('Output results')
  
      if subsample:
          cb_hist_sampled.to_csv('ss_{}_'.format(subsample) + os.path.basename(cb_histogram), sep='\t')
  
      if output_evidence_table:
          import shutil
          buf.seek(0)
          with open(output_evidence_table, 'w') as etab_fh:
              shutil.copyfileobj(buf, etab_fh)
  
      if sparse:
          pd.Series(genes.index).to_csv(out + ".rownames", index=False, header=False)
          pd.Series(genes.columns.values).to_csv(out + ".colnames", index=False, header=False)
          with open(out, "w+b") as out_handle:
              scipy.io.mmwrite(out_handle, scipy.sparse.csr_matrix(genes))
      else:
          genes.to_csv(out)