code: def prepare(args):
      
      p = OptionParser(prepare.__doc__)
      p.add_option("--accuracy", default=.85,
                   help="Sequencing per-base accuracy")
      opts, args = p.parse_args(args)
  
      if len(args) != 2:
          sys.exit(not p.print_help())
  
      vcffile, bamfile = args
      right = "{:.2f}".format(opts.accuracy)
      wrong = "{:.2f}".format(1 - opts.accuracy)
      vcf_reader = vcf.Reader(open(vcffile))
      variants = []
      for r in vcf_reader:
          v = CPRA(r)
          if not v.is_valid:
              continue
          variants.append(v)
  
      logging.debug("A total of {} bi-allelic SNVs imported from `{}`".\
                      format(len(variants), vcffile))
  
      bamfile = pysam.AlignmentFile(bamfile, "rb")
      for v in variants:
          pos = v.pos - 1
          for column in bamfile.pileup(v.chr, pos, pos + 1, truncate=True):
              for read in column.pileups:
                  query_position = read.query_position
                  if query_position is None:
                      continue
                  read_name = read.alignment.query_name
                  query_base = read.alignment.query_sequence[query_position]
                  a, b = v.alleles
                  if query_base == a:
                      other_base = b
                  elif query_base == b:
                      other_base = a
                  else:
                      continue
                  print(" ".join(str(x) for x in \
                             (v, read_name, query_base, right, other_base, wrong)))
docs: %prog prepare vcffile bamfile
  
      Convert vcf and bam to variant list. Inputs are:
      - vcffile: contains the positions of variants
      - bamfile: contains the reads that hold the variants
  
      Outputs:
      - reads_to_phase: phasing for each read
      - variants_to_phase: in format of phased vcf