PedPhase.py

#!/usr/bin/python

"""
#
# Phase individuals genotype based only on the pedigree information. 
#
# Usage: PedPhase.py -g <genotype_table> -p <pedigree_info> -o <output_prefix>
#
# Where:
# genotype_table = table generated by the VariantsToTable tool from GATK, containing at least the CHROM,POS,REF,ALT and .GT fields and only one alternative haplotype (https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_variantutils_VariantsToTable.php)
# pedigree_info = a file containing the pedigree information for all individuals of interest (format file as described in http://zzz.bwh.harvard.edu/plink/data.shtml#ped)
# output_prefix = prefix name to save the output
#
# Options: -h for usage help
#          -clean for a clean run, where all intermediate files are erased. Default = keep all files 
#
# Outputs: One table named <output_prefix>_final.txt.
#			This tab limited file will contain the phase information in the following format:
#			CHROM	POS 	IND1   IND2 	...
# 			Where: CHROM and POS, as in genotype_table; IND(N), the phased or unphased genotype in 0:1 mode (0= REF; 1= ALT) of each analyzed individual.
#
"""

import re, sys, getopt

import os

clean = False

# Check for the arguments, open the inputs and print useful help messages

try:
    opts, args = getopt.getopt(sys.argv[1:],"hg:p:o:",["clean"])
except getopt.GetoptError:
    print '\n', '####     Invalid use     ####', '\n'
    print 'Usage: PedPhase.py -g <genotype_table> -p <pedigree_info> -o <output_prefix>'
    print 'For help use PedPhase.py -h'
    sys.exit(99)
    

for opt, arg in opts:
    if opt == '-h':
        print '\n', 'Phase individuals genotype based only on the pedigree information. ', '\n'
        print 'Usage: PedPhase.py -g <genotype_table> -p <pedigree_info> -o <output_prefix>', '\n'
        print 'Where: genotype_table = table generated by the VariantsToTable tool from GATK, containing at least the CHROM,POS,REF,ALT and .GT fields and only one alternative haplotype (https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_variantutils_VariantsToTable.php)'
        print 'pedigree_info = a file containing the pedigree information for all individuals of interest (format file as described in http://zzz.bwh.harvard.edu/plink/data.shtml#ped)'
        print 'output_prefix = prefix name to save the output'
        print 'Options: -h for help'
        print '--clean for a clean run, where all intermediate files are erased. Default = keep all files'
        
        sys.exit()

    elif opt in ("-g"):
        gen = arg
    elif opt in ("-p"):
        ped = arg
    elif opt in ("-o"):
        out_prefix = arg
    elif opt in ("--clean"):
        clean = True

    else:
        assert False, "unhandled option"
        
        
        
### Create useful variables

id = ""
father = ""
mother = ""
off_phase = ""
father_phase = "" 
mother_phase = ""	


h = 0
id_pos = 0
f_pos = 0
m_pos = 0
par_pos = 0
chr_pos = 0
pos_pos = 0
ref_pos = 0
alt_pos = 0


mama = False
papa = False

dict_gametes= dict()


cwd = os.getcwd()

################################

### defining functions

def update_gen(inp, reference, alternative): # change the genotype mode of input
	inp = inp.replace(reference, "0")
	inp = inp.replace(alternative, "1")
	return inp
	
	

########################################################



## Start the first phasing step - for all individuals in the pedigree that have both parents phase info based on the parental genotype
## Intermediate output (*_phase1.txt) will generated for each individual


print "Starting the first phasing step - Phase individuals based on both parents info.\n"

pedigree = open(ped)

for entry in pedigree:
	if entry.startswith("#"): # skip the header
		continue
	entry = entry.split("\t")
	id = str(entry[1])
	father = str(entry[2])
	mother = str(entry[3])
	if father != "0"  and mother != "0": # check if individual have mother and father info to phase
		
		# start saving the first output
		outname = out_prefix + id + "_phase1.txt"
		out1 = open(outname,"w")
		out1.write("CHROM\tPOS\tGT\tPHASE\tMP_PHASE\tFP_PHASE\n")
		
		# convert the ids for the format in the table
		id2 = id + ".GT" 
		father2 = father + ".GT" 
		mother2 = mother + ".GT" 
		h = 0
		genotype = open(gen)
		
		for line in genotype: 
			header = line.split("\t")
			if h == 0: # get the position containing info for the desired information
				id_pos = header.index(id2)
				f_pos = header.index(father2)
				m_pos = header.index(mother2)
				chr_pos = header.index("CHROM")
				pos_pos = header.index("POS")
				ref_pos = header.index("REF")
				alt_pos = header.index("ALT")
				h = 1
				continue
			else:
				# change genotype mode
				ref = header[ref_pos]
				alt = header[alt_pos]				
				off_gt = header[id_pos]
				off_gt = update_gen(off_gt, ref, alt)
				father_gt = header[f_pos]
				father_gt = update_gen(father_gt, ref, alt)
				mother_gt = header[m_pos]
				mother_gt = update_gen(mother_gt, ref, alt)
				
				# check the phase and phase if posible, if deletion (*) or more than one ALT genotype in the position the phasing wont work
							
				## Phase homozygous offspring
				if off_gt == "0/0": 
					off_phase = "0:0" 
					father_phase = "0" 
					mother_phase = "0"
				if off_gt == "1/1": 
					off_phase = "1:1" 
					father_phase = "1" 
					mother_phase = "1"
			
				## Phase homozygous parents
				elif father_gt == "0/0" and mother_gt == "1/1":
					if off_gt == "0/1" or off_gt == "1/0":   # check if there is a genotype error before phasing
						off_phase = "0:1"
						father_phase = "0" 
						mother_phase = "1"
					else:
						off_phase = ""
						father_phase = "" 
						mother_phase = ""				
				elif father_gt == "1/1" and mother_gt == "0/0":
					if off_gt == "0/1" or off_gt == "1/0":   # check if there is a genotype error before phasing
						off_phase = "1:0" 
						father_phase = "1" 
						mother_phase = "0"
					else:
						off_phase = ""
						father_phase = "" 
						mother_phase = ""				
				
				## Phase one homozygous and one heterozygous parent
				elif father_gt == "0/0" and (mother_gt == "0/1" or mother_gt == "1/0"):
					if off_gt == "0/1" or off_gt == "1/0":  
						off_phase = "0:1" 
						father_phase = "0" 
						mother_phase = "1"
					elif off_gt == "0/0":  # homozygous offspring
						off_phase = "0:0" 
						father_phase = "0" 
						mother_phase = "0"
					else: # it is an error
						off_phase = ""
						father_phase = "" 
						mother_phase = ""		
				elif (father_gt == "0/1" or father_gt == "1/0") and mother_gt == "0/0":
					if off_gt == "0/1" or off_gt == "1/0":  
						off_phase = "1:0" 
						father_phase = "1" 
						mother_phase = "0"
					elif off_gt == "0/0":  # homozygous offspring
						off_phase = "0:0" 
						father_phase = "0" 
						mother_phase = "0"
					else: # it is an error
						off_phase = ""
						father_phase = "" 
						mother_phase = ""
				
				elif father_gt == "1/1" and (mother_gt == "0/1" or mother_gt == "1/0"):
					if off_gt == "0/1" or off_gt == "1/0":  
						off_phase = "1:0" 
						father_phase = "1" 
						mother_phase = "0"
					elif off_gt == "1/1":  # homozygous offspring
						off_phase = "1:1" 
						father_phase = "1" 
						mother_phase = "1"
					else: # it is an error
						off_phase = ""
						father_phase = "" 
						mother_phase = ""		
				elif (father_gt == "0/1" or father_gt == "1/0") and mother_gt == "1/1":
					if off_gt == "0/1" or off_gt == "1/0":  
						off_phase = "0:1" 
						father_phase = "0" 
						mother_phase = "1"
					elif off_gt == "1/1":  # homozygous offspring
						off_phase = "1:1" 
						father_phase = "1" 
						mother_phase = "1"
					else: # it is an error
						off_phase = ""
						father_phase = "" 
						mother_phase = ""
	
				
				else: # it is not possible to phase / possible error
					off_phase = ""
					father_phase = "" 
					mother_phase = ""
				
				# save the information of this position in the output
				chr = header[chr_pos]
				pos = header[pos_pos]		
				out1.write(str(chr) + "\t" + str(pos) + "\t" + str(off_gt) + "\t" + str(off_phase) + "\t" + str(father_phase) + "\t" + str(mother_phase) + "\n")
		out1.close()





## Start the second phasing step - for all individuals with only one parental
## Intermediate output (*_phase2.txt) will generated for each new individual phased

print "Starting the second phasing step - Phase individuals that have only one parent sequenced.\n"

pedigree = open(ped)

for entry in pedigree:
	if entry.startswith("#"): # skip the header
		continue
	entry = entry.split("\t")
	id = str(entry[1])
	father = str(entry[2])
	mother = str(entry[3])
	mama = False
	papa = False
	if (father == "0" and mother != "0") or (father != "0" and mother == "0"): # check if individual have only one of the parents, i.e. was not phased before
		
		# start saving the second output
		outname = out_prefix + id + "_phase2.txt"
		out2 = open(outname,"w")
		out2.write("CHROM\tPOS\tGT\tPHASE\tMP_PHASE\tFP_PHASE\n")
		
		
		# convert the ids for the format in the table
		id2 = id + ".GT" 
		father2 = father + ".GT" 
		mother2 = mother + ".GT" 
		h = 0
		genotype = open(gen)
		
		for line in genotype: 
			header = line.split("\t")
			if h == 0: # get the position containing info for the desired information
				id_pos = header.index(id2)
				chr_pos = header.index("CHROM")
				pos_pos = header.index("POS")
				ref_pos = header.index("REF")
				alt_pos = header.index("ALT")
				h = 1
				
				# take only the info from parent that is sequenced
				if (father == "0"):
					par_pos = header.index(mother2)
					mama = True
				if (mother == "0"):
					par_pos = header.index(father2)	
					papa = True		 
				continue
			else:
				# change genotype mode
				ref = header[ref_pos]
				alt = header[alt_pos]				
				off_gt = header[id_pos]
				off_gt = update_gen(off_gt, ref, alt)
				parental_gt = header[par_pos]
				parental_gt = update_gen(parental_gt, ref, alt)
				
				
				#Phase, if possible: loci is homozygous in the parent and heterozygous in offspring
				if off_gt == "0/1" or off_gt == "1/0":
		
					# if mother is homozygous
					if (parental_gt == "0/0" and mama == True):			
						off_phase = "1:0"
						father_phase = "1" 
						mother_phase = "0" 
					elif (parental_gt == "1/1" and mama == True):			
						off_phase = "0:1"
						father_phase = "0" 
						mother_phase = "1"
					
					# if father is homozygous
					elif (parental_gt == "0/0" and papa == True):			
						off_phase = "0:1"
						father_phase = "0" 
						mother_phase = "1"
					elif (parental_gt == "1/1" and papa == True):			
						off_phase = "1:0"
						father_phase = "1" 
						mother_phase = "0"
					
					elif (parental_gt != "1/1" and parental_gt != "0/0"): # it is not possible to phase 
						off_phase = ""
						father_phase = "" 
						mother_phase = ""
				
				
				# Phase homozygous
				elif off_gt == "0/0": 
					off_phase = "0:0"
					father_phase = "0" 
					mother_phase = "0"
					
				elif off_gt == "1/1": 
					off_phase = "1:1" 
					father_phase = "1" 
					mother_phase = "1"
				
				# save the information of this position in the output
				chr = header[chr_pos]
				pos = header[pos_pos]			
				out2.write(str(chr) + "\t" + str(pos) + "\t" + str(off_gt) + "\t" + str(off_phase) + "\t" + str(father_phase) + "\t" + str(mother_phase) + "\n")
		out2.close()




############################################

print "Saving the final output...\n"

## save the output final file with the gamete info (*_final.txt)
## CHROM	POS 	IND1  IND2 	...



## Final output (*_final.txt) will be generated for each individual phased in previous steps

phased = list()

g = 0

father_gam = ""
mother_gam = ""


## get the list of phased individuals in previous steps


for file in os.listdir(cwd):
    if file.endswith("_phase1.txt") or file.endswith("_phase2.txt"):
        phased_ind = file.split(out_prefix )
        phased_ind = str(phased_ind[1])
        phased_ind = phased_ind.split("_")
        phased_ind = str(phased_ind[0])
        phased.append(phased_ind)


## Creating the header and the final output file

outname = out_prefix + "_final.txt"
out3 = open(outname,"w")
out3.write("CHROM\tPOS")
header = list()

for ind in phased:
	header.append(ind)
	out3.write("\t" + ind)

out3.write("\n")



## Save all genotypes (phased when available) for all individuals

all_variants = {}
dict_alleles = dict()

for id in header: # the order of columns in saved in header
	ind_file_name = out_prefix + id
			
	for file in os.listdir(cwd): # open the file with phase info
		if (file.startswith(ind_file_name)) and (file.endswith("_phase1.txt") or file.endswith("_phase2.txt")):
			ind_phase_file = open(file)
			
			for entry in ind_phase_file:
				entry = entry.split("\t")
				position = str(entry[1])
				phase = str(entry[3])
				if len(phase) > 1:
					dict_alleles[position] = phase
				else:
					unphased = str(entry[2])
					dict_alleles[position] = unphased
					
	all_variants[id] = dict_alleles
	dict_alleles = dict()
		


## Save the genotype information for each individual concatenating phased and unphased data when not phased


position1 = ""


par_phase_file = ""
ind_phase_file = ""


chr_pos = 0
pos_pos = 0


genotype = open(gen)
		
for line in genotype: 
	if line.startswith("CHROM"): # skip the header and get info
		head = line.split("\t")
		chr_pos = head.index("CHROM")
		pos_pos = head.index("POS")
		continue
		
	else: # write the position info
		line = line.split("\t")
		chrom = str(line[chr_pos])
		position1 = str(line[pos_pos])
		out3.write(chrom + "\t" + position1)
		
		
		for id in header: # the order in which the individuals appear in columns is saved in header
			pos_dict = all_variants[id]
			phase = str(pos_dict[position1])
			out3.write("\t" + phase)
			
		out3.write("\n")
	
out3.close()		



#########################
## Clean intermediate outputs if the -clean option is used


if clean is True:
	print "Phasing completed. Removing intermediated outputs as requested \n"
	
	for file in os.listdir(cwd): 
		if file.endswith("_phase1.txt") or file.endswith("_phase2.txt"):
			os.remove(file)
					
#######################


print "Phasing program complete! =D See final output in", outname, "\n"