Mercurial > repos > miller-lab > genome_diversity

comparison average_fst.py @ 12:4b6590dd7250

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author miller-lab
date Wed, 12 Sep 2012 17:10:26 -0400
parents 2c498d40ecde
children f04f40a36cc8
comparison
equal deleted inserted replaced
11:d4ec09e8079f 12:4b6590dd7250
1 #!/usr/bin/env python
2
3 import sys
4 import subprocess
5 from Population import Population
6
7 ################################################################################
8
9 if len(sys.argv) < 12:
10 print >> sys.stderr, "Usage"
11 sys.exit(1)
12
13 input, p1_input, p2_input, data_source, min_total_count, discard_fixed, biased, output, shuffles, p0_input = sys.argv[1:11]
14 individual_metadata = sys.argv[11:]
15
16 try:
17 shuffle_count = int(shuffles)
18 except:
19 shuffle_count = 0
20
21 p_total = Population()
22 p_total.from_tag_list(individual_metadata)
23
24 p1 = Population()
25 p1.from_population_file(p1_input)
26 if not p_total.is_superset(p1):
27 print >> sys.stderr, 'There is an individual in population 1 that is not in the SNP table'
28 sys.exit(1)
29
30 p2 = Population()
31 p2.from_population_file(p2_input)
32 if not p_total.is_superset(p2):
33 print >> sys.stderr, 'There is an individual in population 2 that is not in the SNP table'
34 sys.exit(1)
35
36 p0 = None
37 if shuffle_count > 0:
38 p0 = Population()
39 p0.from_population_file(p0_input)
40 if not p_total.is_superset(p0):
41 print >> sys.stderr, 'There is an individual in population 0 that is not in the SNP table'
42 sys.exit(1)
43
44 ################################################################################
45
46 prog = 'Fst_ave'
47
48 args = []
49 args.append(prog)
50 args.append(input)
51 args.append(data_source)
52 args.append(min_total_count)
53 args.append(discard_fixed)
54 args.append(biased)
55 args.append(shuffles)
56
57 columns = p1.column_list()
58 for column in columns:
59 args.append('{0}:1'.format(column))
60
61 columns = p2.column_list()
62 for column in columns:
63 args.append('{0}:2'.format(column))
64
65 if p0 is not None:
66 columns = p0.column_list()
67 for column in columns:
68 args.append('{0}:0'.format(column))
69
70 fh = open(output, 'w')
71
72 #print "args:", ' '.join(args)
73 p = subprocess.Popen(args, bufsize=-1, stdin=None, stdout=fh, stderr=sys.stderr)
74 rc = p.wait()
75 fh.close()
76
77 sys.exit(0)
78