text: user2085: samtools phase is an easy-to-use tool to phase variants from a bam file and a reference.
  For a sequence that has two haplotypes, like this:
     position: 123456789
  haplotype 0: GACCATATT
         SNPs:  **    *
  haplotype 1: GTACATACT
  
  The output looks like this:
  CC  M?   chr    PS     pos  allele0 allele1 hetIndex #supports0 #errors0 #supp1 #err1
  CC
  CC
  PS      frame-1   2    9
  M1      frame-1  60    2      A       T       1       3       0       23      0
  M1      frame-1  60    3      C       A       2       35      1       24      0
  M1      frame-1  60    8      T       C       3       46      0       30      0
  
  How can I take the reference fasta file and the output of samtools phase to generate one fasta file for each of the haplotypes?
  
  user2085: Answer January 20th, 2020. If you run samtools phase with the option -b my_prefix it outputs a bam file for both haplotypes, named my_prefix.0.bam and my_prefix.1.bam. If the sequence has been correctly phased, use these bam files to correct the reference using pilon. This is preferable to the method described in the January 14th, 2020 answer because pilon will correct SNPs and indels, whereas the python program below will only correct the specific SNPs that samtools phase found. samtools phase does not find indels.
  Answer January 14th, 2020. This python program below takes a fasta file as an argument and the output of samtools phase piped in, then prints out a fasta file of the two haplotypes.
  Run the program like this: 
  samtools phase reads.bam | phase2fasta.py ref.fa > haplotypes.fasta.
  Here is the code of phase2fasta.py:
  #!/usr/bin/env python
  
  import os
  import sys
  
  if (len(sys.argv) == 1) or not os.path.exists(sys.argv[1]):
      raise Exception("The first and only argument is the reference fasta.")
  allele0 = []
  seqname = ""
  with open(sys.argv[1]) as f:
      for line in f:
          if line[0] == ">":
              seqname = line[1:].split(" ")[0].strip()
          else:
              for c in line.strip().upper():
                  allele0.append(c)
  allele1 = [x for x in allele0]
  for line in sys.stdin:
      if line[0:2] == "M1":
          fields = line.strip().split()
          if fields[1] == seqname:
              pos = int(fields[3]) - 1
              allele0[pos] = fields[4]
              allele1[pos] = fields[5]
  print(">{}_allele0".format(seqname))
  print("".join(allele0))
  print(">{}_allele1".format(seqname))
  print("".join(allele1))
  
  
original_text: 2085: <p><a href="http://www.htslib.org/doc/samtools-phase.1.html" rel="nofollow noreferrer">samtools phase</a> is an easy-to-use tool to phase variants from a bam file and a reference.</p>
  
  <p>For a sequence that has two haplotypes, like this:</p>
  
  <pre><code>   position: 123456789
  haplotype 0: GACCATATT
         SNPs:  **    *
  haplotype 1: GTACATACT
  </code></pre>
  
  <p>The output looks like this:</p>
  
  <pre><code>CC  M?   chr    PS     pos  allele0 allele1 hetIndex #supports0 #errors0 #supp1 #err1
  CC
  CC
  PS      frame-1   2    9
  M1      frame-1  60    2      A       T       1       3       0       23      0
  M1      frame-1  60    3      C       A       2       35      1       24      0
  M1      frame-1  60    8      T       C       3       46      0       30      0
  </code></pre>
  
  <p>How can I take the reference fasta file and the output of <code>samtools phase</code> to generate one <code>fasta</code> file for each of the haplotypes?</p>
   2085: <p><strong>Answer January 20th, 2020.</strong> If you run <code>samtools phase</code> with the option <code>-b my_prefix</code> it outputs a bam file for both haplotypes, named <code>my_prefix.0.bam</code> and <code>my_prefix.1.bam</code>. If the sequence has been correctly phased, use these bam files to correct the reference using <a href="https://github.com/broadinstitute/pilon" rel="nofollow noreferrer">pilon</a>. This is preferable to the method described in the <strong>January 14th, 2020</strong> answer because pilon will correct SNPs and indels, whereas the python program below will only correct the specific SNPs that <code>samtools phase</code> found. <code>samtools phase</code> does not find indels.</p>
  
  <p><strong>Answer January 14th, 2020.</strong> This python program below takes a fasta file as an argument and the output of samtools phase piped in, then prints out a fasta file of the two haplotypes.</p>
  
  <p>Run the program like this: </p>
  
  <p><code>samtools phase reads.bam | phase2fasta.py ref.fa &gt; haplotypes.fasta</code>.</p>
  
  <p>Here is the code of <code>phase2fasta.py</code>:</p>
  
  <pre class="lang-py prettyprint-override"><code>#!/usr/bin/env python
  
  import os
  import sys
  
  if (len(sys.argv) == 1) or not os.path.exists(sys.argv[1]):
      raise Exception("The first and only argument is the reference fasta.")
  allele0 = []
  seqname = ""
  with open(sys.argv[1]) as f:
      for line in f:
          if line[0] == "&gt;":
              seqname = line[1:].split(" ")[0].strip()
          else:
              for c in line.strip().upper():
                  allele0.append(c)
  allele1 = [x for x in allele0]
  for line in sys.stdin:
      if line[0:2] == "M1":
          fields = line.strip().split()
          if fields[1] == seqname:
              pos = int(fields[3]) - 1
              allele0[pos] = fields[4]
              allele1[pos] = fields[5]
  print("&gt;{}_allele0".format(seqname))
  print("".join(allele0))
  print("&gt;{}_allele1".format(seqname))
  print("".join(allele1))
  </code></pre>
  