Mercurial > repos > miller-lab > genome_diversity

comparison aggregate_gd_indivs.py @ 27:8997f2ca8c7a

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Update to Miller Lab devshed revision bae0d3306d3b
author Richard Burhans <burhans@bx.psu.edu>
date Mon, 15 Jul 2013 10:47:35 -0400
parents 91e835060ad2
children
comparison
equal deleted inserted replaced
26:91e835060ad2 27:8997f2ca8c7a
1 #!/usr/bin/env python 1 #!/usr/bin/env python
2 2
3 import gd_util
3 import sys 4 import sys
4 import subprocess
5 from Population import Population 5 from Population import Population
6 6
7 ################################################################################ 7 ################################################################################
8 8
9 if len(sys.argv) < 6: 9 if len(sys.argv) != 6:
10 print >> sys.stderr, "Usage" 10 gd_util.dir('Usage')
11 sys.exit(1)
12 11
13 input, p1_input, output, input_type = sys.argv[1:5] 12 input, p1_input, output, input_type, ind_arg = sys.argv[1:]
14 individual_metadata = sys.argv[5:]
15 13
16 p_total = Population() 14 p_total = Population()
17 p_total.from_tag_list(individual_metadata) 15 p_total.from_wrapped_dict(ind_arg)
18 16
19 p1 = Population() 17 p1 = Population()
20 p1.from_population_file(p1_input) 18 p1.from_population_file(p1_input)
21 19
22 if not p_total.is_superset(p1): 20 if not p_total.is_superset(p1):
23 print >> sys.stderr, 'There is an individual in the population that is not in the SNP table' 21 gd_util.die('There is an individual in the population that is not in the SNP table')
24 sys.exit(1)
25 22
26 ################################################################################ 23 ################################################################################
27 24
28 prog = 'aggregate' 25 prog = 'aggregate'
29 26
30 args = [] 27 args = [ prog ]
31 args.append(prog)
32 args.append(input) 28 args.append(input)
33 29
34 if input_type == 'gd_snp': 30 if input_type == 'gd_snp':
35 args.append('1') 31 args.append(1)
36 elif input_type == 'gd_genotype': 32 elif input_type == 'gd_genotype':
37 args.append('0') 33 args.append(0)
38 else: 34 else:
39 print >> sys.stderr, "unknown input type:", input_type 35 die('unknown input type: {0}'.format(input_type))
40 sys.exit(1)
41 36
42 columns = p1.column_list() 37 columns = p1.column_list()
43 38
44 for column in sorted(columns): 39 for column in sorted(columns):
45 if input_type == 'gd_genotype': 40 if input_type == 'gd_genotype':
46 column = str(int(column) - 2) 41 column = str(int(column) - 2)
47 args.append(column) 42 args.append(column)
48 43
49 fh = open(output, 'w') 44 with open(output, 'w') as fh:
50 45 gd_util.run_program(prog, args, stdout=fh)
51 #print "args:", ' '.join(args)
52 p = subprocess.Popen(args, bufsize=-1, stdin=None, stdout=fh, stderr=sys.stderr)
53 rc = p.wait()
54 fh.close()
55 46
56 sys.exit(0) 47 sys.exit(0)
57 48