Convert PEDMAP to beagle markers user Kantale

From PyPedia
Jump to: navigation, search

Contents

[edit] Documentation

Convert from PED and MAP (http://pngu.mgh.harvard.edu/~purcell/plink/dataman.shtml) to beagle's "markers" format.


[edit] Parameters

<inputs>
</inputs>


[edit] Return

[edit] See also

Convert_list_PEDMAP_to_beagle_markers_user_Kantale

[edit] Code

def Convert_PEDMAP_to_beagle_markers_user_Kantale(
	input_PED_filename = None,
	input_MAP_filename = None,
	output_markers_filename = None,
):

	output_markers_file = open(output_markers_filename, "w")
	
	input_PED = [x.replace("\n", "").split() for x in open(input_PED_filename).readlines()]
	samples = len(input_PED)
	snps = (len(input_PED[0]) - 6) / 2

	input_MAP_file = open(input_MAP_filename)
	problems_single_allele = 0
	problems_zero_alleles = 0
	for snp in xrange(snps):
		line = input_MAP_file.readline()
		if not line:
			print "MAP contains less lines than PED SNPs"
			break

		line_s = line.replace("\n", "").split()

		line_to_print = [line_s[1], line_s[3]]

		allele_1 = None
		allele_2 = None
		for sample in xrange(samples):

			current_allele_1 = input_PED[sample][(snp * 2) + 6]
			current_allele_2 = input_PED[sample][((snp * 2) + 1) + 6]

			if current_allele_1 != "0" and not allele_1:
				allele_1 = current_allele_1

			if current_allele_1 != allele_1 and current_allele_1 != "0" and allele_1:
				allele_2 = current_allele_1
				break

			if current_allele_2 != allele_1 and current_allele_1 != "0" and allele_1:
				allele_2 = current_allele_2
				break

		if allele_2:
			line_to_print += [allele_1, allele_2]
		elif allele_1:
			line_to_print += [allele_1, allele_1]
			print "Warning, reference:", line_s[1], "   pos:", line_s[3], "  Has a single allele:", allele_1
			problems_single_allele += 1
		else:
			line_to_print += ["A", "A"] #Put a random allele
			print "Warning, reference:", line_s[1], "   pos:", line_s[3], "  Has zero alleles:"
			problems_zero_alleles += 1
		output_markers_file.write(str.join("\t", line_to_print) + "\n")

	while True:
		line = input_MAP_file.readline()
		if line:
			print "Extra line in MAP file:", line
		else:
			break

	output_markers_file.close()
	input_MAP_file.close()

	print "SNPs with single alleles:", problems_single_allele 
	print "SNPs with zero alleles:", problems_zero_alleles

	print "Output file:", output_markers_filename

[edit] Unit Tests

def uni1():
	return True

[edit] Development Code

def Convert_PEDMAP_to_beagle_markers_user_Kantale():
	pass

[edit] Permissions

[edit] Documentation Permissions

Kantale

[edit] Code Permissions

Kantale

[edit] Unit Tests Permissions

Kantale

[edit] Permissions Permissions

Kantale

Personal tools
Namespaces

Variants
Actions
Navigation
Toolbox