<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>

CC BY-SA 4.0
2020-01-15T12:40:48.690Z
Oh hey, thanks for this! We need this sort of self-answered post!
terdon
CC BY-SA 4.0
2020-01-14T23:17:51.827Z
2020-01-20T12:51:38.150Z
2020-01-20T12:51:38.150Z
conchoecia
<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>

<sam><fasta><samtools><phase>
How to generate two fasta files from samtools phase output?
# How to generate two fasta files from samtools phase output?

(Asked by: username_0 on 2020-01-14)

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?

**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. 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))

> **Comments:**
> **username_1**: Oh hey, thanks for this! We need this sort of self-answered post!  

---

