annotate tools/rgenetics/rgGLM.py @ 0:9071e359b9a3

Uploaded
author xuebing
date Fri, 09 Mar 2012 19:37:19 -0500
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
1 #!/usr/local/bin/python
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
2 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
3 # added most of the available options for linear models
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
4 # june 2009 rml
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
5 # hack to run and process a plink quantitative trait
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
6 #
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
7
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
8 This is a wrapper for Shaun Purcell's Plink linear/logistic models for
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
9 traits, covariates and genotypes.
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
10
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
11 It requires some judgement to interpret the findings
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
12 We need some better visualizations - manhattan plots are good.
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
13 svg with rs numbers for top 1%?
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
14
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
15 toptable tools - truncate a gg file down to some low percentile
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
16
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
17 intersect with other tables - eg gene expression regressions on snps
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
18
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
19
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
20
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
21 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
22
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
23 import sys,math,shutil,subprocess,os,string,tempfile,shutil,commands
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
24 from rgutils import plinke
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
25
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
26 def makeGFF(resf='',outfname='',logf=None,twd='.',name='track name',description='track description',topn=1000):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
27 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
28 score must be scaled to 0-1000
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
29
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
30 Want to make some wig tracks from each analysis
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
31 Best n -log10(p). Make top hit the window.
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
32 we use our tab output which has
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
33 rs chrom offset ADD_stat ADD_p ADD_log10p
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
34 rs3094315 1 792429 1.151 0.2528 0.597223
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
35
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
36 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
37
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
38 def is_number(s):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
39 try:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
40 float(s)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
41 return True
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
42 except ValueError:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
43 return False
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
44 header = 'track name=%s description="%s" visibility=2 useScore=1 color=0,60,120\n' % (name,description)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
45 column_names = [ 'Seqname', 'Source', 'Feature', 'Start', 'End', 'Score', 'Strand', 'Frame', 'Group' ]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
46 halfwidth=100
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
47 resfpath = os.path.join(twd,resf)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
48 resf = open(resfpath,'r')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
49 resfl = resf.readlines() # dumb but convenient for millions of rows
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
50 resfl = [x.split() for x in resfl]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
51 headl = resfl[0]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
52 resfl = resfl[1:]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
53 headl = [x.strip().upper() for x in headl]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
54 headIndex = dict(zip(headl,range(0,len(headl))))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
55 chrpos = headIndex.get('CHROM',None)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
56 rspos = headIndex.get('RS',None)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
57 offspos = headIndex.get('OFFSET',None)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
58 ppos = headIndex.get('ADD_LOG10P',None)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
59 wewant = [chrpos,rspos,offspos,ppos]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
60 if None in wewant: # missing something
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
61 logf.write('### Error missing a required header in makeGFF - headIndex=%s\n' % headIndex)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
62 return
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
63 resfl = [x for x in resfl if x[ppos] > '']
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
64 resfl = [(float(x[ppos]),x) for x in resfl] # decorate
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
65 resfl.sort()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
66 resfl.reverse() # using -log10 so larger is better
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
67 resfl = resfl[:topn] # truncate
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
68 pvals = [x[0] for x in resfl] # need to scale
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
69 resfl = [x[1] for x in resfl] # drop decoration
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
70 if len(pvals) == 0:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
71 logf.write('### no pvalues found in resfl - %s' % (resfl[:3]))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
72 sys.exit(1)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
73 maxp = max(pvals) # need to scale
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
74 minp = min(pvals)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
75 prange = abs(maxp-minp) + 0.5 # fudge
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
76 scalefact = 1000.0/prange
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
77 logf.write('###maxp=%f,minp=%f,prange=%f,scalefact=%f\n' % (maxp,minp,prange,scalefact))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
78 for i,row in enumerate(resfl):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
79 row[ppos] = '%d' % (int(scalefact*pvals[i]))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
80 resfl[i] = row # replace
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
81 outf = file(outfname,'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
82 outf.write(header)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
83 outres = [] # need to resort into chrom offset order
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
84 for i,lrow in enumerate(resfl):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
85 chrom,snp,offset,p, = [lrow[x] for x in wewant]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
86 gff = ('chr%s' % chrom,'rgGLM','variation','%d' % (int(offset)-halfwidth),
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
87 '%d' % (int(offset)+halfwidth),p,'.','.','%s logp=%1.2f' % (snp,pvals[i]))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
88 outres.append(gff)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
89 outres = [(x[0],int(x[3]),x) for x in outres] # decorate
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
90 outres.sort() # into chrom offset
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
91 outres=[x[2] for x in outres] # undecorate
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
92 outres = ['\t'.join(x) for x in outres]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
93 outf.write('\n'.join(outres))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
94 outf.write('\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
95 outf.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
96
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
97
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
98
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
99 def xformQassoc(resf='',outfname='',logf=None,twd='.'):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
100 """ plink.assoc.linear to gg file
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
101 from the docs
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
102 The output per each SNP might look something like:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
103
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
104 CHR SNP BP A1 TEST NMISS OR STAT P
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
105 5 rs000001 10001 A ADD 664 0.7806 -1.942 0.05216
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
106 5 rs000001 10001 A DOMDEV 664 0.9395 -0.3562 0.7217
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
107 5 rs000001 10001 A COV1 664 0.9723 -0.7894 0.4299
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
108 5 rs000001 10001 A COV2 664 1.159 0.5132 0.6078
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
109 5 rs000001 10001 A GENO_2DF 664 NA 5.059 0.0797
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
110 need to transform into gg columns for each distinct test
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
111 or bed for tracks?
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
112
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
113 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
114 logf.write('xformQassoc got resf=%s, outfname=%s\n' % (resf,outfname))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
115 resdict = {}
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
116 rsdict = {}
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
117 markerlist = []
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
118 # plink is "clever" - will run logistic if only 2 categories such as gender
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
119 resfs = resf.split('.')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
120 if resfs[-1] == 'logistic':
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
121 resfs[-1] = 'linear'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
122 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
123 resfs[-1] = 'logistic'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
124 altresf = '.'.join(resfs)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
125
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
126 altresfpath = os.path.join(twd,altresf)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
127 resfpath = os.path.join(twd,resf)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
128 try:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
129 resf = open(resfpath,'r')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
130 except:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
131 try:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
132 resf = open(altresfpath,'r')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
133 except:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
134 print >> sys.stderr, '## error - no file plink output %s or %s found - cannot continue' % (resfpath, altresfpath)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
135 sys.exit(1)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
136 for lnum,row in enumerate(resf):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
137 if lnum == 0:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
138 headl = row.split()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
139 headl = [x.strip().upper() for x in headl]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
140 headIndex = dict(zip(headl,range(0,len(headl))))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
141 chrpos = headIndex.get('CHR',None)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
142 rspos = headIndex.get('SNP',None)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
143 offspos = headIndex.get('BP',None)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
144 nmisspos = headIndex.get('NMISS',None)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
145 testpos = headIndex.get('TEST',None)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
146 ppos = headIndex.get('P',None)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
147 coeffpos = headIndex.get('OR',None)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
148 if not coeffpos:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
149 coeffpos = headIndex.get('BETA',None)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
150 apos = headIndex.get('A1',None)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
151 statpos = headIndex.get('STAT',None)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
152 wewant = [chrpos,rspos,offspos,testpos,statpos,ppos,coeffpos,apos]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
153 if None in wewant: # missing something
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
154 logf.write('missing a required header in xformQassoc - headIndex=%s\n' % headIndex)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
155 return
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
156 llen = len(headl)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
157 else: # no Nones!
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
158 ll = row.split()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
159 if len(ll) >= llen: # valid line
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
160 chrom,snp,offset,test,stat,p,coeff,allele = [ll[x] for x in wewant]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
161 snp = snp.strip()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
162 if p <> 'NA' :
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
163 try:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
164 ffp = float(p)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
165 if ffp <> 0:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
166 lp = -math.log10(ffp)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
167 except:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
168 lp = 0.0
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
169 resdict.setdefault(test,{})
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
170 resdict[test][snp] = (stat,p,'%f' % lp)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
171 if rsdict.get(snp,None) == None:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
172 rsdict[snp] = (chrom,offset)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
173 markerlist.append(snp)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
174 # now have various tests indexed by rs
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
175 tk = resdict.keys()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
176 tk.sort() # tests
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
177 ohead = ['rs','chrom','offset']
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
178 for t in tk: # add headers
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
179 ohead.append('%s_stat' % t)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
180 ohead.append('%s_p' % t)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
181 ohead.append('%s_log10p' % t)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
182 oheads = '\t'.join(ohead)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
183 res = [oheads,]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
184 for snp in markerlist: # retain original order
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
185 chrom,offset = rsdict[snp]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
186 outl = [snp,chrom,offset]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
187 for t in tk:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
188 outl += resdict[t][snp] # add stat,p for this test
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
189 outs = '\t'.join(outl)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
190 res.append(outs)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
191 f = file(outfname,'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
192 res.append('')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
193 f.write('\n'.join(res))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
194 f.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
195
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
196
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
197 if __name__ == "__main__":
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
198 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
199
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
200 <command interpreter="python">
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
201 rgGLM.py '$i.extra_files_path/$i.metadata.base_name' '$phef.extra_files_path/$phef.metadata.base_name'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
202 "$title1" '$predvar' '$covar' '$out_file1' '$logf' '$i.metadata.base_name'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
203 '$inter' '$cond' '$gender' '$mind' '$geno' '$maf' '$logistic' '$wigout'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
204 </command>
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
205 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
206 topn = 1000
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
207 killme = string.punctuation+string.whitespace
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
208 trantab = string.maketrans(killme,'_'*len(killme))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
209 if len(sys.argv) < 17:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
210 s = 'rgGLM.py needs 17 params - got %s \n' % (sys.argv)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
211 sys.stderr.write(s) # print >>,s would probably also work?
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
212 sys.exit(0)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
213 blurb = 'rgGLM.py called with %s' % sys.argv
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
214 print >> sys.stdout,blurb
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
215 bfname = sys.argv[1]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
216 phename = sys.argv[2]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
217 title = sys.argv[3]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
218 title.translate(trantab)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
219 predvar = sys.argv[4]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
220 covar = sys.argv[5].strip()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
221 outfname = sys.argv[6]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
222 logfname = sys.argv[7]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
223 op = os.path.split(logfname)[0]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
224 try: # for test - needs this done
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
225 os.makedirs(op)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
226 except:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
227 pass
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
228 basename = sys.argv[8].translate(trantab)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
229 inter = sys.argv[9] == '1'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
230 cond = sys.argv[10].strip()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
231 if cond == 'None':
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
232 cond = ''
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
233 gender = sys.argv[11] == '1'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
234 mind = sys.argv[12]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
235 geno = sys.argv[13]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
236 maf = sys.argv[14]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
237 logistic = sys.argv[15].strip()=='1'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
238 gffout = sys.argv[16]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
239 me = sys.argv[0]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
240 phepath = '%s.pphe' % phename
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
241 twd = tempfile.mkdtemp(suffix='rgGLM') # make sure plink doesn't spew log file into the root!
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
242 tplog = os.path.join(twd,'%s.log' % basename) # should be path to plink log
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
243 vcl = [plinke,'--noweb','--bfile',bfname,'--pheno-name','"%s"' % predvar,'--pheno',
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
244 phepath,'--out',basename,'--mind %s' % mind, '--geno %s' % geno,
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
245 '--maf %s' % maf]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
246 if logistic:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
247 vcl.append('--logistic')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
248 resf = '%s.assoc.logistic' % basename # plink output is here we hope
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
249 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
250 vcl.append('--linear')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
251 resf = '%s.assoc.linear' % basename # plink output is here we hope
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
252 resf = os.path.join(twd,resf)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
253 if gender:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
254 vcl.append('--sex')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
255 if inter:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
256 vcl.append('--interaction')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
257 if covar > 'None':
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
258 vcl += ['--covar',phepath,'--covar-name',covar] # comma sep list of covariates
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
259 tcfile = None
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
260 if len(cond) > 0: # plink wants these in a file..
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
261 dummy,tcfile = tempfile.mkstemp(suffix='condlist') #
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
262 f = open(tcfile,'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
263 cl = cond.split()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
264 f.write('\n'.join(cl))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
265 f.write('\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
266 f.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
267 vcl.append('--condition-list %s' % tcfile)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
268 p=subprocess.Popen(' '.join(vcl),shell=True,cwd=twd)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
269 retval = p.wait()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
270 if tcfile:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
271 os.unlink(tcfile)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
272 plinklog = file(tplog,'r').read()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
273 logf = file(logfname,'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
274 logf.write(blurb)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
275 logf.write('\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
276 logf.write('vcl=%s\n' % vcl)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
277 xformQassoc(resf=resf,outfname=outfname,logf=logf,twd=twd) # leaves the desired summary file
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
278 makeGFF(resf=outfname,outfname=gffout,logf=logf,twd=twd,name='rgGLM_TopTable',description=title,topn=topn)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
279 logf.write('\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
280 logf.write(plinklog)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
281 logf.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
282 #shutil.rmtree(twd) # clean up
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
283
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
284
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
285
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
286
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
287