annotate tools/rgenetics/rgTDT.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 # hack to run and process a plink tdt
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
3 # expects args as
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
4 # bfilepath outname jobname outformat (wig,xls)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
5 # ross lazarus
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
6 # for wig files, we need annotation so look for map file or complain
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
7
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
8 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
9 Parameters for wiggle track definition lines
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
10 All options are placed in a single line separated by spaces:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
11
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
12 track type=wiggle_0 name=track_label description=center_label \
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
13 visibility=display_mode color=r,g,b altColor=r,g,b \
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
14 priority=priority autoScale=on|off \
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
15 gridDefault=on|off maxHeightPixels=max:default:min \
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
16 graphType=bar|points viewLimits=lower:upper \
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
17 yLineMark=real-value yLineOnOff=on|off \
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
18 windowingFunction=maximum|mean|minimum smoothingWindow=off|2-16
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
19 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
20
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
21 import sys,math,shutil,subprocess,os,time,tempfile,shutil,string
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
22 from os.path import abspath
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
23 from optparse import OptionParser
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
24 from rgutils import timenow, plinke
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
25 myversion = 'v0.003 January 2010'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
26 verbose = False
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
27
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
28
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
29
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
30 def makeGFF(resf='',outfname='',logf=None,twd='.',name='track name',description='track description',topn=1000):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
31 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
32 score must be scaled to 0-1000
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
33
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
34 Want to make some wig tracks from each analysis
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
35 Best n -log10(p). Make top hit the window.
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
36 we use our tab output which has
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
37 rs chrom offset ADD_stat ADD_p ADD_log10p
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
38 rs3094315 1 792429 1.151 0.2528 0.597223
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
39
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
40 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
41
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
42 def is_number(s):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
43 try:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
44 float(s)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
45 return True
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
46 except ValueError:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
47 return False
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
48 header = 'track name=%s description="%s" visibility=2 useScore=1 color=0,60,120\n' % (name,description)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
49 column_names = [ 'Seqname', 'Source', 'Feature', 'Start', 'End', 'Score', 'Strand', 'Frame', 'Group' ]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
50 halfwidth=100
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
51 resfpath = os.path.join(twd,resf)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
52 resf = open(resfpath,'r')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
53 resfl = resf.readlines() # dumb but convenient for millions of rows
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
54 resfl = [x.split() for x in resfl]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
55 headl = resfl[0]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
56 resfl = resfl[1:]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
57 headl = [x.strip().upper() for x in headl]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
58 headIndex = dict(zip(headl,range(0,len(headl))))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
59 # s = 'rs\tchrom\toffset\ta1\ta2\ttransmitted\tuntransmitted\tTDTChiSq\tTDTp\t-log10TDTp\tAbsTDTOR\tTDTOR'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
60 chrpos = headIndex.get('CHROM',None)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
61 rspos = headIndex.get('RS',None)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
62 offspos = headIndex.get('OFFSET',None)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
63 ppos = headIndex.get('-LOG10TDTP',None)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
64 wewant = [chrpos,rspos,offspos,ppos]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
65 if None in wewant: # missing something
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
66 logf.write('### Error missing a required header in makeGFF - headIndex=%s\n' % headIndex)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
67 return
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
68 resfl = [x for x in resfl if x[ppos] > '']
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
69 resfl = [(float(x[ppos]),x) for x in resfl] # decorate
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
70 resfl.sort()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
71 resfl.reverse() # using -log10 so larger is better
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
72 resfl = resfl[:topn] # truncate
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
73 pvals = [x[0] for x in resfl] # need to scale
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
74 resfl = [x[1] for x in resfl] # drop decoration
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
75 maxp = max(pvals) # need to scale
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
76 minp = min(pvals)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
77 prange = abs(maxp-minp) + 0.5 # fudge
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
78 scalefact = 1000.0/prange
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
79 logf.write('###maxp=%f,minp=%f,prange=%f,scalefact=%f\n' % (maxp,minp,prange,scalefact))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
80 for i,row in enumerate(resfl):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
81 row[ppos] = '%d' % (int(scalefact*pvals[i]))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
82 resfl[i] = row # replace
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
83 outf = file(outfname,'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
84 outf.write(header)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
85 outres = [] # need to resort into chrom offset order
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
86 for i,lrow in enumerate(resfl):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
87 chrom,snp,offset,p, = [lrow[x] for x in wewant]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
88 gff = ('chr%s' % chrom,'rgTDT','variation','%d' % (int(offset)-halfwidth),
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
89 '%d' % (int(offset)+halfwidth),p,'.','.','%s logp=%1.2f' % (snp,pvals[i]))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
90 outres.append(gff)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
91 outres = [(x[0],int(x[3]),x) for x in outres] # decorate
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
92 outres.sort() # into chrom offset
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
93 outres=[x[2] for x in outres] # undecorate
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
94 outres = ['\t'.join(x) for x in outres]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
95 outf.write('\n'.join(outres))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
96 outf.write('\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
97 outf.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
98
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
99
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
100
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
101 def xformTDT(infname='',resf='',outfname='',name='foo',mapf='/usr/local/galaxy/data/rg/lped/x.bim'):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
102 """munge a plink .tdt file into either a ucsc track or an xls file
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
103 CHR SNP A1:A2 T:U_TDT OR_TDT CHISQ_TDT P_TDT
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
104 0 MitoT217C 2:3 0:0 NA NA NA
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
105 0 MitoG228A 1:4 0:0 NA NA NA
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
106 0 MitoT250C 2:3 0:0 NA NA NA
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
107 map file has
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
108 1 rs4378174 0 003980745
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
109 1 rs10796404 0 005465256
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
110 1 rs2697965 0 014023092
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
111
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
112 grrr!
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
113 Changed in 1.01 to
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
114 [rerla@hg fresh]$ head database/job_working_directory/445/rgTDT.tdt
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
115 CHR SNP BP A1 A2 T U OR CHISQ P
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
116 1 rs12562034 758311 1 3 71 79 0.8987 0.4267 0.5136
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
117 1 rs3934834 995669 4 2 98 129 0.7597 4.233 0.03963
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
118
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
119
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
120 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
121 if verbose:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
122 print 'Rgenetics Galaxy Tools, rgTDT.py.xformTDT got resf=%s, outtype=%s, outfname=%s' % (resf,outtype,outfname)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
123 wewantcols = ['SNP','CHR','BP','A1','A2','T','U','OR','CHISQ','P']
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
124 res = []
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
125 s = 'rs\tchrom\toffset\ta1\ta2\ttransmitted\tuntransmitted\tTDTChiSq\tTDTp\t-log10TDTp\tAbsTDTOR\tTDTOR' # header
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
126 res.append(s)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
127 rsdict = {}
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
128 if not mapf:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
129 sys.stderr.write('rgTDT called but no map file provided - cannot determine locations')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
130 sys.exit(1)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
131 map = file(mapf,'r')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
132 for l in map: # plink map
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
133 ll = l.strip().split()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
134 if len(ll) >= 3:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
135 rs=ll[1].strip()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
136 chrom = ll[0]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
137 if chrom.lower() == 'x':
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
138 chrom = '23'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
139 if chrom.lower() == 'y':
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
140 chrom = '24'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
141 if chrom.lower() == 'mito':
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
142 chrom = '25'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
143 offset = ll[3]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
144 rsdict[rs] = (chrom,offset)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
145 f = open(resf,'r')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
146 headl = f.next().strip()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
147 headl = headl.split()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
148 wewant = [headl.index(x) for x in wewantcols]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
149 llen = len(headl)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
150 lnum = anum = 0
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
151 for l in f:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
152 lnum += 1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
153 ll = l.strip().split()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
154 if len(ll) >= llen: # valid line
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
155 snp,chrom,offset,a1,a2,t,u,orat,chisq,p = [ll[x] for x in wewant]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
156 if chisq == 'NA' or p == 'NA' or orat == 'NA':
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
157 continue # can't use these lines - gg gets unhappy
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
158 snp = snp.strip()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
159 lp = '0.0'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
160 fp = '1.0'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
161 fakeorat = '1.0'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
162 if p.upper().strip() <> 'NA':
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
163 try:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
164 fp = float(p)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
165 if fp <> 0:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
166 lp = '%6f' % -math.log10(fp)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
167 fp = '%6f' % fp
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
168 except:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
169 pass
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
170 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
171 p = '1.0'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
172 if orat.upper().strip() <> 'NA':
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
173 try:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
174 fakeorat = orat
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
175 if float(orat) < 1.0:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
176 fakeorat = '%6f' % (1.0/float(orat)) # invert so large values big
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
177 except:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
178 pass
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
179 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
180 orat = '1.0'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
181 outl = '\t'.join([snp,chrom,offset,a1,a2,t,u,chisq,p,lp,fakeorat,orat])
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
182 res.append(outl)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
183 f = file(outfname,'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
184 res.append('')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
185 f.write('\n'.join(res))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
186 f.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
187
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
188
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
189 if __name__ == "__main__":
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
190 """ called as
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
191 <command interpreter="python">
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
192 rgTDT.py -i '$infile.extra_files_path/$infile.metadata.base_name' -o '$title' -f '$outformat' -r '$out_file1' -l '$logf' -x '${GALAXY_DATA_INDEX_DIR}/rg/bin/pl$
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
193
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
194 </command>
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
195
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
196 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
197 u = """ called in xml as
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
198 <command interpreter="python2.4">
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
199 rgTDT.py -i $i -o $out_prefix -f $outformat -r $out_file1 -l $logf
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
200 </command>
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
201 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
202 if len(sys.argv) < 6:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
203 s = '## Error rgTDT.py needs 5 command line params - got %s \n' % (sys.argv)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
204 if verbose:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
205 print >> sys.stdout, s
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
206 sys.exit(0)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
207 parser = OptionParser(usage=u, version="%prog 0.01")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
208 a = parser.add_option
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
209 a("-i","--infile",dest="bfname")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
210 a("-o","--oprefix",dest="oprefix")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
211 a("-f","--formatOut",dest="outformat")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
212 a("-r","--results",dest="outfname")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
213 a("-l","--logfile",dest="logf")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
214 a("-d","--du",dest="uId")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
215 a("-e","--email",dest="uEmail")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
216 a("-g","--gff",dest="gffout",default="")
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
217 (options,args) = parser.parse_args()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
218 killme = string.punctuation + string.whitespace
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
219 trantab = string.maketrans(killme,'_'*len(killme))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
220 title = options.oprefix
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
221 title = title.translate(trantab)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
222 map_file = '%s.bim' % (options.bfname) #
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
223 me = sys.argv[0]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
224 alogf = options.logf # absolute paths
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
225 od = os.path.split(alogf)[0]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
226 try:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
227 os.path.makedirs(od)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
228 except:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
229 pass
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
230 aoutf = options.outfname # absolute paths
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
231 od = os.path.split(aoutf)[0]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
232 try:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
233 os.path.makedirs(od)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
234 except:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
235 pass
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
236 vcl = [plinke,'--noweb', '--bfile',options.bfname,'--out',title,'--mind','0.5','--tdt']
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
237 logme = []
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
238 if verbose:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
239 s = 'Rgenetics %s http://rgenetics.org Galaxy Tools rgTDT.py started %s\n' % (myversion,timenow())
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
240 print >> sys.stdout,s
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
241 logme.append(s)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
242 s ='rgTDT.py: bfname=%s, logf=%s, argv = %s\n' % (options.bfname,alogf, sys.argv)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
243 print >> sys.stdout,s
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
244 logme.append(s)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
245 s = 'rgTDT.py: vcl=%s\n' % (' '.join(vcl))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
246 print >> sys.stdout,s
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
247 logme.append(s)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
248 twd = tempfile.mkdtemp(suffix='rgTDT') # make sure plink doesn't spew log file into the root!
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
249 tname = os.path.join(twd,title)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
250 p=subprocess.Popen(' '.join(vcl),shell=True,cwd=twd)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
251 retval = p.wait()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
252 shutil.copy('%s.log' % tname,alogf)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
253 sto = file(alogf,'a')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
254 sto.write('\n'.join(logme))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
255 resf = '%s.tdt' % tname # plink output is here we hope
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
256 xformTDT(options.bfname,resf,aoutf,title,map_file) # leaves the desired summary file
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
257 gffout = options.gffout
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
258 if gffout > '':
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
259 makeGFF(resf=aoutf,outfname=gffout,logf=sto,twd='.',name='rgTDT_Top_Table',description=title,topn=1000)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
260 shutil.rmtree(twd)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
261 sto.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
262
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
263
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
264