annotate pca.py @ 25:cba0d7a63b82

workaround for gd_genotype datatype admix shift int -> float
author Richard Burhans <burhans@bx.psu.edu>
date Wed, 29 May 2013 13:49:19 -0400
parents 248b06e86022
children 8997f2ca8c7a
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
1 #!/usr/bin/env python
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
2
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
3 import errno
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
4 import os
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
5 import shutil
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
6 import subprocess
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
7 import sys
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
8 from BeautifulSoup import BeautifulSoup
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
9 import gd_composite
24
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
10 import re
0
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
11
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
12 ################################################################################
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
13
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
14 def mkdir_p(path):
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
15 try:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
16 os.makedirs(path)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
17 except OSError, e:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
18 if e.errno <> errno.EEXIST:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
19 raise
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
20
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
21 ################################################################################
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
22
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
23 def run_program(prog, args, stdout_file=None):
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
24 #print "args: ", ' '.join(args)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
25 p = subprocess.Popen(args, bufsize=-1, executable=prog, stdin=None, stdout=subprocess.PIPE, stderr=subprocess.PIPE)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
26 (stdoutdata, stderrdata) = p.communicate()
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
27 rc = p.returncode
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
28
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
29 if stdout_file is not None:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
30 with open(stdout_file, 'w') as ofh:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
31 print >> ofh, stdoutdata
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
32
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
33 if rc != 0:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
34 print >> sys.stderr, "FAILED: rc={0}: {1}".format(rc, ' '.join(args))
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
35 print >> sys.stderr, stderrdata
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
36 sys.exit(1)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
37
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
38 ################################################################################
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
39
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
40 def do_ped2geno(input, output):
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
41 lines = []
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
42 with open(input) as fh:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
43 for line in fh:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
44 line = line.rstrip('\r\n')
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
45 lines.append(line.split())
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
46
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
47 pair_map = {
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
48 '0':{ '0':'9', '1':'9', '2':'9' },
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
49 '1':{ '0':'1', '1':'2', '2':'1' },
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
50 '2':{ '0':'1', '1':'1', '2':'0' }
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
51 }
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
52 with open(output, 'w') as ofh:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
53 for a_idx in xrange(6, len(lines[0]), 2):
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
54 b_idx = a_idx + 1
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
55 print >> ofh, ''.join(map(lambda line: pair_map[line[a_idx]][line[b_idx]], lines))
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
56
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
57 def do_map2snp(input, output):
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
58 with open(output, 'w') as ofh:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
59 with open(input) as fh:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
60 for line in fh:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
61 elems = line.split()
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
62 print >> ofh, ' {0} 11 0.002 2000 A T'.format(elems[1])
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
63
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
64 def make_ind_file(ind_file, input):
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
65 pops = []
24
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
66 name_map = []
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
67 name_idx = 0
0
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
68
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
69 ofh = open(ind_file, 'w')
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
70
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
71 with open(input) as fh:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
72 soup = BeautifulSoup(fh)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
73 misc = soup.find('div', {'id': 'gd_misc'})
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
74 populations = misc('ul')[0]
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
75
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
76 i = 0
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
77 for entry in populations:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
78 if i % 2 == 1:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
79 population_name = entry.contents[0].encode('utf8').strip().replace(' ', '_')
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
80 pops.append(population_name)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
81 individuals = entry.ol('li')
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
82 for individual in individuals:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
83 individual_name = individual.string.encode('utf8').strip()
24
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
84 name_map.append(individual_name)
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
85 print >> ofh, 'ind_%s' % name_idx, 'M', population_name
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
86 name_idx += 1
0
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
87 i += 1
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
88
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
89 ofh.close()
24
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
90 return pops, name_map
0
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
91
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
92 def make_par_file(par_file, geno_file, snp_file, ind_file, evec_file, eval_file):
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
93 with open(par_file, 'w') as fh:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
94 print >> fh, 'genotypename: {0}'.format(geno_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
95 print >> fh, 'snpname: {0}'.format(snp_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
96 print >> fh, 'indivname: {0}'.format(ind_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
97 print >> fh, 'evecoutname: {0}'.format(evec_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
98 print >> fh, 'evaloutname: {0}'.format(eval_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
99 print >> fh, 'altnormstyle: NO'
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
100 print >> fh, 'numoutevec: 2'
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
101
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
102 def do_smartpca(par_file):
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
103 prog = 'smartpca'
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
104
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
105 args = [ prog ]
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
106 args.append('-p')
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
107 args.append(par_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
108
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
109 #print "args: ", ' '.join(args)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
110 p = subprocess.Popen(args, bufsize=-1, stdin=None, stdout=subprocess.PIPE, stderr=sys.stderr)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
111 (stdoutdata, stderrdata) = p.communicate()
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
112 rc = p.returncode
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
113
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
114 if rc != 0:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
115 print >> sys.stderr, "FAILED: rc={0}: {1}".format(rc, ' '.join(args))
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
116 print >> sys.stderr, stderrdata
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
117 sys.exit(1)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
118
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
119 stats = []
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
120
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
121 save_line = False
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
122 for line in stdoutdata.split('\n'):
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
123 if line.startswith(('## Average divergence', '## Anova statistics', '## Statistical significance')):
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
124 stats.append('')
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
125 save_line = True
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
126 if line.strip() == '':
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
127 save_line = False
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
128 if save_line:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
129 stats.append(line)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
130
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
131 return '\n'.join(stats[1:])
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
132
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
133 def do_ploteig(evec_file, population_names):
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
134 prog = 'gd_ploteig'
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
135
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
136 args = [ prog ]
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
137 args.append('-i')
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
138 args.append(evec_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
139 args.append('-c')
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
140 args.append('1:2')
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
141 args.append('-p')
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
142 args.append(':'.join(population_names))
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
143 args.append('-x')
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
144
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
145 run_program(None, args)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
146
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
147 def do_eval2pct(eval_file, explained_file):
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
148 prog = 'eval2pct'
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
149
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
150 args = [ prog ]
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
151 args.append(eval_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
152
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
153 with open(explained_file, 'w') as ofh:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
154 #print "args:", ' '.join(args)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
155 p = subprocess.Popen(args, bufsize=-1, stdin=None, stdout=ofh, stderr=subprocess.PIPE)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
156 (stdoutdata, stderrdata) = p.communicate()
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
157 rc = p.returncode
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
158
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
159 if rc != 0:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
160 print >> sys.stderr, "FAILED: rc={0}: {1}".format(rc, ' '.join(args))
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
161 print >> sys.stderr, stderrdata
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
162 sys.exit(1)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
163
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
164 def do_coords2admix(coords_file):
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
165 prog = 'coords2admix'
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
166
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
167 args = [ prog ]
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
168 args.append(coords_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
169
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
170 with open('fake', 'w') as ofh:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
171 #print "args:", ' '.join(args)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
172 p = subprocess.Popen(args, bufsize=-1, stdin=None, stdout=ofh, stderr=subprocess.PIPE)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
173 (stdoutdata, stderrdata) = p.communicate()
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
174 rc = p.returncode
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
175
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
176 if rc != 0:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
177 print >> sys.stderr, "FAILED: rc={0}: {1}".format(rc, ' '.join(args))
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
178 print >> sys.stderr, stderrdata
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
179 sys.exit(1)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
180
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
181 shutil.copy2('fake', coords_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
182
24
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
183 ind_regex = re.compile('ind_([0-9]+)')
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
184
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
185 def fix_names(name_map, files):
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
186 for file in files:
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
187 tmp_filename = '%s.tmp' % file
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
188 with open(tmp_filename, 'w') as ofh:
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
189 with open(file) as fh:
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
190 for line in fh:
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
191 line = line.rstrip('\r\n')
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
192 match = ind_regex.search(line)
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
193 if match:
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
194 idx = int(match.group(1))
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
195 old = 'ind_%s' % idx
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
196 new = name_map[idx].replace(' ', '_')
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
197 line = line.replace(old, new)
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
198 print >> ofh, line
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
199
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
200 shutil.copy2(tmp_filename, file)
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
201 os.unlink(tmp_filename)
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
202
0
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
203 ################################################################################
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
204
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
205 if len(sys.argv) != 5:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
206 print "usage"
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
207 sys.exit(1)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
208
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
209 input, input_files_path, output, output_files_path = sys.argv[1:5]
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
210
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
211 mkdir_p(output_files_path)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
212
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
213 ped_file = os.path.join(input_files_path, 'admix.ped')
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
214 geno_file = os.path.join(output_files_path, 'admix.geno')
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
215 do_ped2geno(ped_file, geno_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
216
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
217 map_file = os.path.join(input_files_path, 'admix.map')
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
218 snp_file = os.path.join(output_files_path, 'admix.snp')
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
219 do_map2snp(map_file, snp_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
220
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
221 ind_file = os.path.join(output_files_path, 'admix.ind')
24
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
222 population_names, name_map = make_ind_file(ind_file, input)
0
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
223
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
224 par_file = os.path.join(output_files_path, 'par.admix')
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
225 evec_file = os.path.join(output_files_path, 'coordinates.txt')
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
226 eval_file = os.path.join(output_files_path, 'admix.eval')
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
227 make_par_file(par_file, geno_file, snp_file, ind_file, evec_file, eval_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
228
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
229 smartpca_stats = do_smartpca(par_file)
24
248b06e86022 Added gd_genotype datatype. Modified tools to support new datatype.
Richard Burhans <burhans@bx.psu.edu>
parents: 0
diff changeset
230 fix_names(name_map, [ind_file, evec_file])
0
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
231
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
232 do_ploteig(evec_file, population_names)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
233 plot_file = 'coordinates.txt.1:2.{0}.pdf'.format(':'.join(population_names))
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
234 output_plot_file = os.path.join(output_files_path, 'PCA.pdf')
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
235 shutil.copy2(plot_file, output_plot_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
236 os.unlink(plot_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
237
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
238 do_eval2pct(eval_file, os.path.join(output_files_path, 'explained.txt'))
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
239 os.unlink(eval_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
240
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
241 do_coords2admix(evec_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
242
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
243 ################################################################################
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
244
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
245 info_page = gd_composite.InfoPage()
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
246 info_page.set_title('PCA Galaxy Composite Dataset')
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
247
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
248 display_file = gd_composite.DisplayFile()
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
249 display_value = gd_composite.DisplayValue()
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
250
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
251 out_pdf = gd_composite.Parameter(name='PCA.pdf', value='PCA.pdf', display_type=display_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
252 out_evec = gd_composite.Parameter(name='coordinates.txt', value='coordinates.txt', display_type=display_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
253 out_explained = gd_composite.Parameter(name='explained.txt', value='explained.txt', display_type=display_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
254
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
255 evec_prefix = 'coordinates.txt.1:2.{0}'.format(':'.join(population_names))
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
256 ps_file = '{0}.ps'.format(evec_prefix)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
257 xtxt_file = '{0}.xtxt'.format(evec_prefix)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
258
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
259 os.unlink(os.path.join(output_files_path, ps_file))
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
260 os.unlink(os.path.join(output_files_path, xtxt_file))
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
261
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
262 info_page.add_output_parameter(out_pdf)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
263 info_page.add_output_parameter(out_evec)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
264 info_page.add_output_parameter(out_explained)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
265
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
266 in_admix = gd_composite.Parameter(name='par.admix', value='par.admix', display_type=display_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
267 in_geno = gd_composite.Parameter(name='admix.geno', value='admix.geno', display_type=display_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
268 in_snp = gd_composite.Parameter(name='admix.snp', value='admix.snp', display_type=display_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
269 in_ind = gd_composite.Parameter(name='admix.ind', value='admix.ind', display_type=display_file)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
270
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
271 info_page.add_input_parameter(in_admix)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
272 info_page.add_input_parameter(in_geno)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
273 info_page.add_input_parameter(in_snp)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
274 info_page.add_input_parameter(in_ind)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
275
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
276 misc_stats = gd_composite.Parameter(description='Stats<p/><pre>\n{0}\n</pre>'.format(smartpca_stats), display_type=display_value)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
277
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
278 info_page.add_misc(misc_stats)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
279
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
280 with open (output, 'w') as ofh:
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
281 print >> ofh, info_page.render()
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
282
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
283 sys.exit(0)
2c498d40ecde Uploaded
miller-lab
parents:
diff changeset
284