How to generate two fasta files from samtools phase output? 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?
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 &gt; 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] == "&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))