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.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
if (len(sys.argv) == 1) or not os.path.exists(sys.argv):
raise Exception("The first and only argument is the reference fasta.")
allele0 = 
seqname = ""
with open(sys.argv) as f:
for line in f:
if line == ">":
seqname = line[1:].split(" ").strip()
for c in line.strip().upper():
allele1 = [x for x in allele0]
for line in sys.stdin:
if line[0:2] == "M1":
fields = line.strip().split()
if fields == seqname:
pos = int(fields) - 1
allele0[pos] = fields
allele1[pos] = fields