code: def read_alignment(out_sam, loci, seqs, out_file):
      
      hits = defaultdict(list)
      with open(out_file, "w") as out_handle:
          samfile = pysam.Samfile(out_sam, "r")
          for a in samfile.fetch():
              if not a.is_unmapped:
                  nm = int([t[1] for t in a.tags if t[0] == "NM"][0])
                  a = makeBED(a)
                  if not a:
                      continue
                  ref, locus = get_loci(samfile.getrname(int(a.chr)), loci)
                  hits[a.name].append((nm, "%s %s %s %s %s %s" % (a.name, a.name.split("-")[0], locus, ref, a.start, a.end)))
          for hit in hits.values():
              nm = hit[0][0]
              for l in hit:
                  if nm == l[0]:
                      print(l[1], file=out_handle)
      return out_file
docs: read which seqs map to which loci and
      return a tab separated file