annotate average_fst.py @ 0:99a67ac88802 default tip

Uploaded
author rico
date Thu, 05 Apr 2012 14:22:50 -0400
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
99a67ac88802 Uploaded
rico
parents:
diff changeset
1 #!/usr/bin/env python
99a67ac88802 Uploaded
rico
parents:
diff changeset
2
99a67ac88802 Uploaded
rico
parents:
diff changeset
3 import sys
99a67ac88802 Uploaded
rico
parents:
diff changeset
4 import subprocess
99a67ac88802 Uploaded
rico
parents:
diff changeset
5 from Population import Population
99a67ac88802 Uploaded
rico
parents:
diff changeset
6
99a67ac88802 Uploaded
rico
parents:
diff changeset
7 ################################################################################
99a67ac88802 Uploaded
rico
parents:
diff changeset
8
99a67ac88802 Uploaded
rico
parents:
diff changeset
9 if len(sys.argv) < 12:
99a67ac88802 Uploaded
rico
parents:
diff changeset
10 print >> sys.stderr, "Usage"
99a67ac88802 Uploaded
rico
parents:
diff changeset
11 sys.exit(1)
99a67ac88802 Uploaded
rico
parents:
diff changeset
12
99a67ac88802 Uploaded
rico
parents:
diff changeset
13 input, p1_input, p2_input, data_source, min_total_count, discard_fixed, biased, output, shuffles, p0_input = sys.argv[1:11]
99a67ac88802 Uploaded
rico
parents:
diff changeset
14 individual_metadata = sys.argv[11:]
99a67ac88802 Uploaded
rico
parents:
diff changeset
15
99a67ac88802 Uploaded
rico
parents:
diff changeset
16 try:
99a67ac88802 Uploaded
rico
parents:
diff changeset
17 shuffle_count = int(shuffles)
99a67ac88802 Uploaded
rico
parents:
diff changeset
18 except:
99a67ac88802 Uploaded
rico
parents:
diff changeset
19 shuffle_count = 0
99a67ac88802 Uploaded
rico
parents:
diff changeset
20
99a67ac88802 Uploaded
rico
parents:
diff changeset
21 p_total = Population()
99a67ac88802 Uploaded
rico
parents:
diff changeset
22 p_total.from_tag_list(individual_metadata)
99a67ac88802 Uploaded
rico
parents:
diff changeset
23
99a67ac88802 Uploaded
rico
parents:
diff changeset
24 p1 = Population()
99a67ac88802 Uploaded
rico
parents:
diff changeset
25 p1.from_population_file(p1_input)
99a67ac88802 Uploaded
rico
parents:
diff changeset
26 if not p_total.is_superset(p1):
99a67ac88802 Uploaded
rico
parents:
diff changeset
27 print >> sys.stderr, 'There is an individual in population 1 that is not in the SNP table'
99a67ac88802 Uploaded
rico
parents:
diff changeset
28 sys.exit(1)
99a67ac88802 Uploaded
rico
parents:
diff changeset
29
99a67ac88802 Uploaded
rico
parents:
diff changeset
30 p2 = Population()
99a67ac88802 Uploaded
rico
parents:
diff changeset
31 p2.from_population_file(p2_input)
99a67ac88802 Uploaded
rico
parents:
diff changeset
32 if not p_total.is_superset(p2):
99a67ac88802 Uploaded
rico
parents:
diff changeset
33 print >> sys.stderr, 'There is an individual in population 2 that is not in the SNP table'
99a67ac88802 Uploaded
rico
parents:
diff changeset
34 sys.exit(1)
99a67ac88802 Uploaded
rico
parents:
diff changeset
35
99a67ac88802 Uploaded
rico
parents:
diff changeset
36 p0 = None
99a67ac88802 Uploaded
rico
parents:
diff changeset
37 if shuffle_count > 0:
99a67ac88802 Uploaded
rico
parents:
diff changeset
38 p0 = Population()
99a67ac88802 Uploaded
rico
parents:
diff changeset
39 p0.from_population_file(p0_input)
99a67ac88802 Uploaded
rico
parents:
diff changeset
40 if not p_total.is_superset(p0):
99a67ac88802 Uploaded
rico
parents:
diff changeset
41 print >> sys.stderr, 'There is an individual in population 0 that is not in the SNP table'
99a67ac88802 Uploaded
rico
parents:
diff changeset
42 sys.exit(1)
99a67ac88802 Uploaded
rico
parents:
diff changeset
43
99a67ac88802 Uploaded
rico
parents:
diff changeset
44 ################################################################################
99a67ac88802 Uploaded
rico
parents:
diff changeset
45
99a67ac88802 Uploaded
rico
parents:
diff changeset
46 prog = 'Fst_ave'
99a67ac88802 Uploaded
rico
parents:
diff changeset
47
99a67ac88802 Uploaded
rico
parents:
diff changeset
48 args = []
99a67ac88802 Uploaded
rico
parents:
diff changeset
49 args.append(prog)
99a67ac88802 Uploaded
rico
parents:
diff changeset
50 args.append(input)
99a67ac88802 Uploaded
rico
parents:
diff changeset
51 args.append(data_source)
99a67ac88802 Uploaded
rico
parents:
diff changeset
52 args.append(min_total_count)
99a67ac88802 Uploaded
rico
parents:
diff changeset
53 args.append(discard_fixed)
99a67ac88802 Uploaded
rico
parents:
diff changeset
54 args.append(biased)
99a67ac88802 Uploaded
rico
parents:
diff changeset
55 args.append(shuffles)
99a67ac88802 Uploaded
rico
parents:
diff changeset
56
99a67ac88802 Uploaded
rico
parents:
diff changeset
57 columns = p1.column_list()
99a67ac88802 Uploaded
rico
parents:
diff changeset
58 for column in columns:
99a67ac88802 Uploaded
rico
parents:
diff changeset
59 args.append('{0}:1'.format(column))
99a67ac88802 Uploaded
rico
parents:
diff changeset
60
99a67ac88802 Uploaded
rico
parents:
diff changeset
61 columns = p2.column_list()
99a67ac88802 Uploaded
rico
parents:
diff changeset
62 for column in columns:
99a67ac88802 Uploaded
rico
parents:
diff changeset
63 args.append('{0}:2'.format(column))
99a67ac88802 Uploaded
rico
parents:
diff changeset
64
99a67ac88802 Uploaded
rico
parents:
diff changeset
65 if p0 is not None:
99a67ac88802 Uploaded
rico
parents:
diff changeset
66 columns = p0.column_list()
99a67ac88802 Uploaded
rico
parents:
diff changeset
67 for column in columns:
99a67ac88802 Uploaded
rico
parents:
diff changeset
68 args.append('{0}:0'.format(column))
99a67ac88802 Uploaded
rico
parents:
diff changeset
69
99a67ac88802 Uploaded
rico
parents:
diff changeset
70 fh = open(output, 'w')
99a67ac88802 Uploaded
rico
parents:
diff changeset
71
99a67ac88802 Uploaded
rico
parents:
diff changeset
72 #print "args:", ' '.join(args)
99a67ac88802 Uploaded
rico
parents:
diff changeset
73 p = subprocess.Popen(args, bufsize=-1, stdin=None, stdout=fh, stderr=sys.stderr)
99a67ac88802 Uploaded
rico
parents:
diff changeset
74 rc = p.wait()
99a67ac88802 Uploaded
rico
parents:
diff changeset
75 fh.close()
99a67ac88802 Uploaded
rico
parents:
diff changeset
76
99a67ac88802 Uploaded
rico
parents:
diff changeset
77 sys.exit(0)
99a67ac88802 Uploaded
rico
parents:
diff changeset
78