Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-01-15T12:40:48.690Z
Answers.Text: Oh hey, thanks for this! We need this sort of self-answered post!
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-01-14T23:17:51.827Z
Answers.LastActivityDate: 2020-01-20T12:51:38.150Z
Answers.LastEditDate: 2020-01-20T12:51:38.150Z
Answers.UserDisplayName: conchoecia
Body: <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>
  
Tags: <sam><fasta><samtools><phase>
Title: How to generate two fasta files from samtools phase output?
ThreadText: # How to generate two fasta files from samtools phase output?
  
  (Asked by: username_0 on 2020-01-14)
  
  [samtools phase](http://www.htslib.org/doc/samtools-phase.1.html) 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?
  
  **Tags:** sam, fasta, samtools, phase
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_0 on 2020-01-14. Score: 3 | Accepted answer
  
  **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](https://github.com/broadinstitute/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`:
  
  ```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))
  
  ```
  
  > **Comments:**
  > **username_1**: Oh hey, thanks for this! We need this sort of self-answered post!  
  
  ---
  
  