|
0
|
1 """
|
|
|
2 fakephe.py
|
|
|
3 ross lazarus sept 30 2007
|
|
|
4 This is available under the LGPL as defined then.
|
|
|
5
|
|
|
6 use the pedigree data for ids
|
|
|
7
|
|
|
8 use pythons generators to literally generate a bunch of random phenotype measurements
|
|
|
9
|
|
|
10 Plink format at http://pngu.mgh.harvard.edu/~purcell/plink/data.shtml#pheno
|
|
|
11 is
|
|
|
12
|
|
|
13 To specify an alternate phenotype for analysis, i.e. other than the one in the *.ped file (or, if using a binary fileset, the
|
|
|
14 *.fam file), use the --pheno option:
|
|
|
15 plink --file mydata --pheno pheno.txt
|
|
|
16
|
|
|
17 where pheno.txt is a file that contains 3 columns (one row per individual):
|
|
|
18
|
|
|
19 Family ID
|
|
|
20 Individual ID
|
|
|
21 Phenotype
|
|
|
22
|
|
|
23 NOTE The original PED file must still contain a phenotype in column 6 (even if this is a dummy phenotype, e.g. all missing),
|
|
|
24 unless the --no-pheno flag is given.
|
|
|
25
|
|
|
26 If an individual is in the original file but not listed in the alternate phenotype file, that person's phenotype will be set to
|
|
|
27 missing. If a person is in the alternate phenotype file but not in the original file, that entry will be ignored. The order of
|
|
|
28 the alternate phenotype file need not be the same as for the original file.
|
|
|
29
|
|
|
30 If the phenotype file contains more than one phenotype, then use the --mpheno N option to specify the Nth phenotype is the one
|
|
|
31 to be used:
|
|
|
32 plink --file mydata --pheno pheno2.txt --mpheno 4
|
|
|
33
|
|
|
34 where pheno2.txt contains 5 different phenotypes (i.e. 7 columns in total), this command will use the 4th for analysis
|
|
|
35 (phenotype D):
|
|
|
36
|
|
|
37 Family ID
|
|
|
38 Individual ID
|
|
|
39 Phenotype A
|
|
|
40 Phenotype B
|
|
|
41 Phenotype C
|
|
|
42 Phenotype D
|
|
|
43 Phenotype E
|
|
|
44
|
|
|
45 Alternatively, your alternate phenotype file can have a header row, in which case you can use variable names to specify which
|
|
|
46 phenotype to use. If you have a header row, the first two variables must be labelled FID and IID. All subsequent variable names
|
|
|
47 cannot have any whitespace in them. For example,
|
|
|
48
|
|
|
49 FID IID qt1 bmi site
|
|
|
50 F1 1110 2.3 22.22 2
|
|
|
51 F2 2202 34.12 18.23 1
|
|
|
52 ...
|
|
|
53
|
|
|
54 then
|
|
|
55 plink --file mydata --pheno pheno2.txt --pheno-name bmi --assoc
|
|
|
56
|
|
|
57 will select the second phenotype labelled "bmi", for analysis
|
|
|
58
|
|
|
59 Finally, if there is more than one phenotype, then for basic association tests, it is possible to specify that all phenotypes
|
|
|
60 be tested, sequentially, with the output sent to different files: e.g. if bigpheno.raw contains 10,000 phenotypes, then
|
|
|
61 plink --bfile mydata --assoc --pheno bigpheno.raw --all-pheno
|
|
|
62
|
|
|
63 will loop over all of these, one at a time testing for association with SNP, generating a lot of output. You might want to use
|
|
|
64 the --pfilter command in this case, to only report results with a p-value less than a certain value, e.g. --pfilter 1e-3.
|
|
|
65
|
|
|
66 WARNING Currently, all phenotypes must be numerically coded, including missing values, in the alternate phenotype file. The
|
|
|
67 default missing value is -9, change this with --missing-phenotype, but it must be a numeric value still (in contrast to the
|
|
|
68 main phenotype in the PED/FAM file. This issue will be fixed in future releases.
|
|
|
69 Covariate files
|
|
|
70
|
|
|
71 ===========================
|
|
|
72 rgfakePhe.xml
|
|
|
73 <tool id="fakePhe1" name="Fake phenos">
|
|
|
74 <description>for multiple null fake phenotype</description>
|
|
|
75 <command interpreter="python2.4">rgfakePhe.py $input1 '$title1' $out_file1 $log_file1 $script_file</command>
|
|
|
76 <inputs>
|
|
|
77 <page>
|
|
|
78 <param name="input1"
|
|
|
79 type="library" format="lped"
|
|
|
80 label="Pedigree from Dataset"/>
|
|
|
81 <param name="title1" type="text"
|
|
|
82 value="My fake phenos" size="60"
|
|
|
83 label="Title for outputs"/>
|
|
|
84 </page>
|
|
|
85 <page>
|
|
|
86 <repeat name="fakePhe" title="Phenotypes to Fake">
|
|
|
87 <param name="pName" type="text" label="Phenotype Name">
|
|
|
88 </param>
|
|
|
89 <conditional name="series">
|
|
|
90 <param name="phetype" type="select" label="Phenotype Type">
|
|
|
91 <option value="rnorm" selected="true">Random normal variate</option>
|
|
|
92 <option value="cat">Random categorical</option>
|
|
|
93 </param>
|
|
|
94 <when value="rnorm">
|
|
|
95 <param name="Mean" type="float" value="0.0" label="Mean">
|
|
|
96 </param>
|
|
|
97 <param name="SD" type="float" label="SD" value="1.0">
|
|
|
98 </param>
|
|
|
99 </when>
|
|
|
100 <when value="cat">
|
|
|
101 <param name="values" type="text" value="1,2,3,fiddle" label="comma separated values to choose from">
|
|
|
102 </param>
|
|
|
103 </when>
|
|
|
104 </conditional>
|
|
|
105 </repeat>
|
|
|
106 </page>
|
|
|
107 </inputs>
|
|
|
108 <configfiles>
|
|
|
109 <configfile name="script_file">
|
|
|
110 #for $n, $f in enumerate($fakePhe)
|
|
|
111 #if $f.series.phetype=='rnorm'
|
|
|
112 {'pN':'$f.pName','pT':'$f.series.phetype','pP':"{'Mean':'$f.series.Mean', 'SD':'$f.series.SD'}"}
|
|
|
113 #elif $f.series.phetype=='cat'
|
|
|
114 {'pN':'$f.pName','pT':'$f.series.phetype','pP':"{'values':'$f.series.values'}"}
|
|
|
115 #end if
|
|
|
116 #end for
|
|
|
117 </configfile>
|
|
|
118 </configfiles>
|
|
|
119
|
|
|
120 <outputs>
|
|
|
121 <data format="pphe" name="out_file1" />
|
|
|
122 <data format="text" name="log_file1" parent="out_file1"/>
|
|
|
123 </outputs>
|
|
|
124
|
|
|
125 <help>
|
|
|
126 .. class:: infomark
|
|
|
127
|
|
|
128 This tool allows you to generate an arbitrary (sort of)
|
|
|
129 synthetic phenotype file with measurements drawn from normal,
|
|
|
130 gamma, or categorical distributions. These are for testing under
|
|
|
131 the null hypothesis of no association - the values are random but
|
|
|
132 from user specified distributions.
|
|
|
133
|
|
|
134 -----
|
|
|
135
|
|
|
136 .. class:: warningmark
|
|
|
137
|
|
|
138 This tool is very experimental
|
|
|
139
|
|
|
140 -----
|
|
|
141
|
|
|
142 - **Pedigree** is a library pedigree file - the id's will be used in the synthetic null phenotypes
|
|
|
143 - **Title** is a name to give to the output phenotype file
|
|
|
144
|
|
|
145 On the next page, you can add an unlimited number of various kinds of phenotypes including choices for
|
|
|
146 categorical ones or distributions with specific parameters
|
|
|
147
|
|
|
148 Just keep adding new ones until you're done then use the Execute button to run the generation
|
|
|
149
|
|
|
150
|
|
|
151
|
|
|
152
|
|
|
153 **Example from another tool to keep you busy and in awe**
|
|
|
154
|
|
|
155 Input file::
|
|
|
156
|
|
|
157 1 68 4.1
|
|
|
158 2 71 4.6
|
|
|
159 3 62 3.8
|
|
|
160 4 75 4.4
|
|
|
161 5 58 3.2
|
|
|
162 6 60 3.1
|
|
|
163 7 67 3.8
|
|
|
164 8 68 4.1
|
|
|
165 9 71 4.3
|
|
|
166 10 69 3.7
|
|
|
167
|
|
|
168 Create a two series XY plot on the above data:
|
|
|
169
|
|
|
170 - Series 1: Red Dashed-Line plot between columns 1 and 2
|
|
|
171 - Series 2: Blue Circular-Point plot between columns 3 and 2
|
|
|
172
|
|
|
173 .. image:: ./static/images/xy_example.jpg
|
|
|
174 </help>
|
|
|
175 </tool>
|
|
|
176
|
|
|
177
|
|
|
178
|
|
|
179 """
|
|
|
180
|
|
|
181 import random,copy,sys,os,time,string,math
|
|
|
182
|
|
|
183 doFbatphe = False
|
|
|
184
|
|
|
185 galhtmlprefix = """<?xml version="1.0" encoding="utf-8" ?>
|
|
|
186 <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
|
|
|
187 <html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">
|
|
|
188 <head>
|
|
|
189 <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
|
|
|
190 <meta name="generator" content="Galaxy %s tool output - see http://g2.trac.bx.psu.edu/" />
|
|
|
191 <title></title>
|
|
|
192 <link rel="stylesheet" href="/static/style/base.css" type="text/css" />
|
|
|
193 </head>
|
|
|
194 <body>
|
|
|
195 <div class="document">
|
|
|
196 """
|
|
|
197
|
|
|
198
|
|
|
199 def timenow():
|
|
|
200 """return current time as a string
|
|
|
201 """
|
|
|
202 return time.strftime('%d/%m/%Y %H:%M:%S', time.localtime(time.time()))
|
|
|
203
|
|
|
204
|
|
|
205 def poisson(lamb=2):
|
|
|
206 """
|
|
|
207 algorithm poisson random number (Knuth):
|
|
|
208 init:
|
|
|
209 Let L = e^-lamb, k = 0 and p = 1.
|
|
|
210 do:
|
|
|
211 k = k + 1.
|
|
|
212 Generate uniform random number u in [0,1] and let p = p u.
|
|
|
213 while p >= L.
|
|
|
214 return k - 1.
|
|
|
215 """
|
|
|
216 lamb = float(lamb)
|
|
|
217 l = math.e**(-lamb)
|
|
|
218 k=0
|
|
|
219 p=1
|
|
|
220 while 1:
|
|
|
221 while p >= l:
|
|
|
222 k += 1
|
|
|
223 u = random.uniform(0.0,1.0)
|
|
|
224 p *= u
|
|
|
225 yield '%e' % (k - 1)
|
|
|
226
|
|
|
227 def gammaPhe(alpha=1,beta=1):
|
|
|
228 """generator for random values from a gamma
|
|
|
229 """
|
|
|
230 while 1: # an iterator for a random phenotype
|
|
|
231 dat = random.gammavariate(float(alpha),float(beta))
|
|
|
232 yield '%e' % dat
|
|
|
233
|
|
|
234 def weibullPhe(alpha=1,beta=1):
|
|
|
235 """generator for random values from a weibull distribution
|
|
|
236 """
|
|
|
237 while 1: # an iterator for a random phenotype
|
|
|
238 dat = random.weibullvariate(float(alpha),float(beta))
|
|
|
239 yield '%e' % dat
|
|
|
240
|
|
|
241 def normPhe(mean=0,sd=1):
|
|
|
242 """generator for random values from a normal distribution
|
|
|
243 """
|
|
|
244 while 1:# an iterator for a random phenotype
|
|
|
245 dat = random.normalvariate(float(mean),float(sd))
|
|
|
246 yield '%e' % dat
|
|
|
247
|
|
|
248 def expoPhe(mean=1):
|
|
|
249 """generator for random values from an exponential distribution
|
|
|
250 """
|
|
|
251 lamb = 1.0/float(mean)
|
|
|
252 while 1:# an iterator for a random phenotype
|
|
|
253 dat = random.expovariate(lamb)
|
|
|
254 yield '%e' % dat
|
|
|
255
|
|
|
256
|
|
|
257 def catPhe(values='1,2,3'):
|
|
|
258 """ Schrodinger's of course.
|
|
|
259 """
|
|
|
260 v = values.split(',')
|
|
|
261 while 1:# an iterator for a random phenotype
|
|
|
262 dat = random.choice(v)
|
|
|
263 yield dat
|
|
|
264
|
|
|
265 def uniPhe(low=0.0,hi=1.0):
|
|
|
266 """generator for a uniform distribution
|
|
|
267 what if low=-5 and hi=-2
|
|
|
268 """
|
|
|
269 low = float(low)
|
|
|
270 hi = float(hi)
|
|
|
271 while 1: # unif phenotype
|
|
|
272 v = random.uniform(low,hi) # 0-1
|
|
|
273 yield '%e' % v
|
|
|
274
|
|
|
275 def getIds(indir='foo'):
|
|
|
276 """read identifiers - first 2 cols from a pedigree file or fam file
|
|
|
277 """
|
|
|
278 inpath = os.path.split(indir)[0] # get root
|
|
|
279 infam = '%s.fam' % indir
|
|
|
280 inped = '%s.ped' % indir # if there's a ped
|
|
|
281 flist = os.listdir(inpath)
|
|
|
282 if len(flist) == 0:
|
|
|
283 print >> sys.stderr, '### Error - input path %s is empty' % indir
|
|
|
284 sys.exit(1)
|
|
|
285 if os.path.exists(infam):
|
|
|
286 pfname = infam
|
|
|
287 elif os.path.exists(inped):
|
|
|
288 pfname = inped
|
|
|
289 else:
|
|
|
290 print >> sys.stderr, '### Error - input path %s contains no ped or fam' % indir
|
|
|
291 sys.exit(1)
|
|
|
292 f = file(pfname,'r')
|
|
|
293 ids = []
|
|
|
294 for l in f:
|
|
|
295 ll = l.split()
|
|
|
296 if len(ll) > 5: # ok line?
|
|
|
297 ids.append(ll[:2])
|
|
|
298 return ids
|
|
|
299
|
|
|
300 def makePhe(phes = [],ids=[]):
|
|
|
301 """Create the null phenotype values and append them to the case id
|
|
|
302 phes is the (generator, headername) for each column
|
|
|
303 for a phe file, ids are the case identifiers for the phenotype file
|
|
|
304 res contains the final strings for the file
|
|
|
305 each column is derived by iterating
|
|
|
306 over the generator actions set up by makePhes
|
|
|
307 """
|
|
|
308 header = ['FID','IID'] # for plink
|
|
|
309 res = copy.copy(ids)
|
|
|
310 for (f,fname) in phes:
|
|
|
311 header.append(fname)
|
|
|
312 for n,subject in enumerate(ids):
|
|
|
313 res[n].append(f.next()) # generate the right value
|
|
|
314 res.insert(0,header)
|
|
|
315 res = [' '.join(x) for x in res] # must be space delim for fbat
|
|
|
316 return res
|
|
|
317
|
|
|
318 def makePhes(pheTypes=[], pheNames=[], pheParams=[]):
|
|
|
319 """set up phes for makePhe
|
|
|
320 each has to have an iterator (action) and a name
|
|
|
321 """
|
|
|
322 action = None
|
|
|
323 phes = []
|
|
|
324 for n,pt in enumerate(pheTypes):
|
|
|
325 name = pheNames[n]
|
|
|
326 if pt == 'rnorm':
|
|
|
327 m = pheParams[n].get('Mean',0.0)
|
|
|
328 s = pheParams[n].get('SD',1.0)
|
|
|
329 action = normPhe(mean=m,sd=s) # so we can just iterate over action
|
|
|
330 elif pt == 'rgamma':
|
|
|
331 m = pheParams[n].get('Alpha',0.0)
|
|
|
332 s = pheParams[n].get('Beta',1.0)
|
|
|
333 action = gammaPhe(alpha=m,beta=s)
|
|
|
334 if pt == 'exponential':
|
|
|
335 m = pheParams[n].get('Mean',1.0)
|
|
|
336 action = expoPhe(mean=m) # so we can just iterate over action
|
|
|
337 elif pt == 'weibull':
|
|
|
338 m = pheParams[n].get('Alpha',0.0)
|
|
|
339 s = pheParams[n].get('Beta',1.0)
|
|
|
340 action = weibullPhe(alpha=m,beta=s)
|
|
|
341 elif pt == 'cat':
|
|
|
342 v = pheParams[n].get('values',['?',])
|
|
|
343 action = catPhe(values=v)
|
|
|
344 elif pt == 'unif':
|
|
|
345 low = pheParams[n].get('low',0.0)
|
|
|
346 hi = pheParams[n].get('hi',1.0)
|
|
|
347 action = uniPhe(low=low,hi=hi)
|
|
|
348 elif pt == 'poisson':
|
|
|
349 lamb = pheParams[n].get('lamb',1.0)
|
|
|
350 action = poisson(lamb=lamb)
|
|
|
351 phes.append((action,name))
|
|
|
352 return phes
|
|
|
353
|
|
|
354 def doImport(outfile='test',flist=[],expl='',mylog=[]):
|
|
|
355 """ import into one of the new html composite data types for Rgenetics
|
|
|
356 Dan Blankenberg with mods by Ross Lazarus
|
|
|
357 October 2007
|
|
|
358 """
|
|
|
359 outf = open(outfile,'w')
|
|
|
360 progname = os.path.basename(sys.argv[0])
|
|
|
361 outf.write(galhtmlprefix % progname)
|
|
|
362 if len(flist) > 0:
|
|
|
363 outf.write('<ol>\n')
|
|
|
364 for i, data in enumerate( flist ):
|
|
|
365 outf.write('<li><a href="%s">%s %s</a></li>\n' % (os.path.split(data)[-1],os.path.split(data)[-1],expl))
|
|
|
366 outf.write('</ol>\n')
|
|
|
367 else:
|
|
|
368 outf.write('No files found')
|
|
|
369 outf.write('<hr/><h4>Log of run follows:</h4><br/><pre>%s\n</pre><br/><hr/>' % ('\n'.join(mylog)))
|
|
|
370 outf.write("</div></body></html>\n")
|
|
|
371 outf.close()
|
|
|
372
|
|
|
373
|
|
|
374 def test():
|
|
|
375 """test case
|
|
|
376 need to get these out of a galaxy form - series of pages - get types
|
|
|
377 on first screen, names on second, params on third?
|
|
|
378 holy shit. this actually works I think
|
|
|
379 """
|
|
|
380 pT = ['rnorm','rnorm','rnorm','rnorm','cat','unif']
|
|
|
381 pN = ['SysBP','DiaBP','HtCM','WtKG','Race','age']
|
|
|
382 pP = [{'Mean':'120','SD':'10'},{'Mean':'90','SD':'15'},{'Mean':'160','SD':'20'},{'Mean':'60','SD':'20'}, \
|
|
|
383 {'values':'Blink,What,Yours,green'},{'low':16,'hi':99}]
|
|
|
384 phes = makePhes(pheTypes=pT, pheNames=pN, pheParams=pP)
|
|
|
385 ids = []
|
|
|
386 for i in range(10):
|
|
|
387 ids.append(['fid%d' % i,'iid%d' % i])
|
|
|
388 pheres = makePhe(phes=phes,ids=ids)
|
|
|
389 res = [''.join(x) for x in pheres]
|
|
|
390 print '\n'.join(res)
|
|
|
391
|
|
|
392
|
|
|
393
|
|
|
394 if __name__ == "__main__":
|
|
|
395 """
|
|
|
396 <command interpreter="python">rgfakePhe.py '$infile1.extra_files_path/$infile1.metadata.base_name'
|
|
|
397 "$title1" '$ppheout' '$ppheout.files_path' '$script_file '
|
|
|
398 </command>
|
|
|
399 The xml file for this tool is complex, and allows any arbitrary number of
|
|
|
400 phenotype columns to be specified from a couple of optional types - rnorm, cat
|
|
|
401 are working now.
|
|
|
402
|
|
|
403 Note that we create new files in their respective library directories and link to them in the output file
|
|
|
404 so they can be displayed and downloaded separately
|
|
|
405
|
|
|
406 """
|
|
|
407 killme = string.punctuation + string.whitespace
|
|
|
408 trantab = string.maketrans(killme,'_'*len(killme))
|
|
|
409 progname = os.path.basename(sys.argv[0])
|
|
|
410 cl = '## at %s, %s got cl= %s' % (timenow(),progname,' '.join(sys.argv))
|
|
|
411 print >> sys.stdout,cl
|
|
|
412 if len(sys.argv) < 5:
|
|
|
413 test()
|
|
|
414 else:
|
|
|
415 inped = sys.argv[1]
|
|
|
416 title = sys.argv[2].translate(trantab)
|
|
|
417 ppheout = sys.argv[3]
|
|
|
418 pphe_path = sys.argv[4]
|
|
|
419 scriptfile = sys.argv[5]
|
|
|
420 ind = file(scriptfile,'r').readlines()
|
|
|
421 mylog = []
|
|
|
422 s = '## %s starting at %s<br/>\n' % (progname,timenow())
|
|
|
423 mylog.append(s)
|
|
|
424 mylog.append(cl)
|
|
|
425 s = '## params = %s<br/>\n' % (' '.join(sys.argv[1:]))
|
|
|
426 mylog.append(s)
|
|
|
427 s = '\n'.join(ind)
|
|
|
428 mylog.append('Script file %s contained %s<br/>\n' % (scriptfile,s))
|
|
|
429 pT = []
|
|
|
430 pN = []
|
|
|
431 pP = []
|
|
|
432 for l in ind:
|
|
|
433 l = l.strip()
|
|
|
434 if len(l) > 1:
|
|
|
435 adict = eval(l)
|
|
|
436 pT.append(adict.get('pT',None))
|
|
|
437 pN.append(adict.get('pN',None))
|
|
|
438 pP.append(eval(adict.get('pP',None)))
|
|
|
439 s = '## pt,pn,pp=%s,%s,%s<br/>\n' % (str(pT),str(pN),str(pP))
|
|
|
440 mylog.append(s)
|
|
|
441 phes = makePhes(pheTypes=pT, pheNames=pN, pheParams=pP)
|
|
|
442 ids = getIds(indir=inped) # for labelling rows
|
|
|
443 pheres = makePhe(phes=phes,ids=ids) # random values from given distributions
|
|
|
444 try:
|
|
|
445 os.makedirs(pphe_path)
|
|
|
446 except:
|
|
|
447 pass
|
|
|
448 outname = os.path.join(pphe_path,title)
|
|
|
449 pphefname = '%s.pphe' % outname
|
|
|
450 f = file(pphefname, 'w')
|
|
|
451 f.write('\n'.join(pheres))
|
|
|
452 f.write('\n')
|
|
|
453 f.close()
|
|
|
454 if doFbatphe:
|
|
|
455 try:
|
|
|
456 os.makedirs(fphe_path)
|
|
|
457 except:
|
|
|
458 pass
|
|
|
459 outname = os.path.join(fphe_path,title)
|
|
|
460 fphefname = '%s.phe' % outname
|
|
|
461 f = file(fphefname, 'w')
|
|
|
462 header = pheres[0].split()
|
|
|
463 pheres[0] = ' '.join(header[2:])# remove iid fid from header for fbat
|
|
|
464 f.write('\n'.join(pheres))
|
|
|
465 f.close()
|
|
|
466 doImport(outfile=fpheout,flist=[fphefname,],expl='(FBAT phenotype format)',mylog=mylog)
|
|
|
467 #doImport(outfile='test',flist=[],expl='',mylog=[]):
|
|
|
468 doImport(outfile=ppheout,flist=[pphefname,],expl='(Plink phenotype format)',mylog=mylog)
|
|
|
469
|
