Mercurial > repos > miller-lab > genome_diversity

comparison add_fst_column.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 248b06e86022
children
comparison
equal deleted inserted replaced
26:91e835060ad2 27:8997f2ca8c7a
1 #!/usr/bin/env python 1 #!/usr/bin/env python
2 2
3 # <command interpreter="python"> 3 import gd_util
4 # add_fst_column.py "$input" "$p1_input" "$p2_input" "$data_source.choice" "$data_source.min_value" "$retain" "$discard_fixed" "$biased" "$output"
5 # #for $individual, $individual_col in zip($input.dataset.metadata.individual_names, $input.dataset.metadata.individual_columns)
6 # #set $arg = '%s:%s' % ($individual_col, $individual)
7 # "$arg"
8 # #end for
9 # </command>
10
11 import sys 4 import sys
12 import subprocess
13 from Population import Population 5 from Population import Population
14 6
15 ################################################################################ 7 ################################################################################
16 8
17 if len(sys.argv) < 12: 9 if len(sys.argv) != 13:
18 print >> sys.stderr, "Usage" 10 gd_util.die('Usage')
19 sys.exit(1)
20 11
21 input, p1_input, p2_input, input_type, genotypes, min_reads, min_qual, retain, discard_fixed, biased, output = sys.argv[1:12] 12 input, p1_input, p2_input, input_type, data_source, min_reads, min_qual, retain, discard_fixed, biased, output, ind_arg = sys.argv[1:]
22 individual_metadata = sys.argv[12:]
23 13
24 p_total = Population() 14 p_total = Population()
25 p_total.from_tag_list(individual_metadata) 15 p_total.from_wrapped_dict(ind_arg)
26 16
27 p1 = Population() 17 p1 = Population()
28 p1.from_population_file(p1_input) 18 p1.from_population_file(p1_input)
29 if not p_total.is_superset(p1): 19 if not p_total.is_superset(p1):
30 print >> sys.stderr, 'There is an individual in population 1 that is not in the SNP table' 20 gd_util.die('There is an individual in population 1 that is not in the SNP table')
31 sys.exit(1)
32 21
33 p2 = Population() 22 p2 = Population()
34 p2.from_population_file(p2_input) 23 p2.from_population_file(p2_input)
35 if not p_total.is_superset(p2): 24 if not p_total.is_superset(p2):
36 print >> sys.stderr, 'There is an individual in population 2 that is not in the SNP table' 25 gd_util.die('There is an individual in population 2 that is not in the SNP table')
37 sys.exit(1)
38 26
39 ################################################################################ 27 ################################################################################
40 28
41 prog = 'Fst_column' 29 prog = 'Fst_column'
42 30
43 args = [] 31 args = [ prog ]
44 args.append(prog)
45 args.append(input) 32 args.append(input)
46 args.append(genotypes) 33 args.append(data_source)
47 args.append(min_reads) 34 args.append(min_reads)
48 args.append(min_qual) 35 args.append(min_qual)
49 args.append(retain) 36 args.append(retain)
50 args.append(discard_fixed) 37 args.append(discard_fixed)
51 args.append(biased) 38 args.append(biased)
60 for column in columns: 47 for column in columns:
61 if input_type == 'gd_genotype': 48 if input_type == 'gd_genotype':
62 column = int(column) - 2 49 column = int(column) - 2
63 args.append('{0}:2'.format(column)) 50 args.append('{0}:2'.format(column))
64 51
65 fh = open(output, 'w') 52 with open(output, 'w') as fh:
66 53 gd_util.run_program(prog, args, stdout=fh)
67 #print "args:", ' '.join(args)
68 p = subprocess.Popen(args, bufsize=-1, stdin=None, stdout=fh, stderr=sys.stderr)
69 rc = p.wait()
70 fh.close()
71 54
72 sys.exit(0) 55 sys.exit(0)
73 56