annotate tools/rgenetics/rgutils.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 # utilities for rgenetics
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
2 #
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
3 # copyright 2009 ross lazarus
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
4 # released under the LGPL
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
5 #
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
6
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
7 import subprocess, os, sys, time, tempfile,string,plinkbinJZ
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
8 import datetime
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
9
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
10 galhtmlprefix = """<?xml version="1.0" encoding="utf-8" ?>
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
11 <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
12 <html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
13 <head>
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
14 <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
15 <meta name="generator" content="Galaxy %s tool output - see http://g2.trac.bx.psu.edu/" />
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
16 <title></title>
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
17 <link rel="stylesheet" href="/static/style/base.css" type="text/css" />
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
18 </head>
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
19 <body>
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
20 <div class="document">
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
21 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
22 galhtmlattr = """<h3><a href="http://rgenetics.org">Rgenetics</a> tool %s run at %s</h3>"""
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
23 galhtmlpostfix = """</div></body></html>\n"""
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
24
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
25 plinke = 'plink' # changed jan 2010 - all exes must be on path
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
26 rexe = 'R' # to avoid cluster/platform dependencies
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
27 smartpca = 'smartpca.perl'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
28
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
29 def timenow():
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
30 """return current time as a string
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
31 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
32 return time.strftime('%d/%m/%Y %H:%M:%S', time.localtime(time.time()))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
33
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
34 def timestamp():
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
35 return datetime.datetime.now().strftime('%Y%m%d%H%M%S')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
36
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
37 def fail( message ):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
38 print >> sys.stderr, message
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
39 return -1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
40
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
41 def whereis(program):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
42 for path in os.environ.get('PATH', '').split(':'):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
43 if os.path.exists(os.path.join(path, program)) and \
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
44 not os.path.isdir(os.path.join(path, program)):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
45 return os.path.join(path, program)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
46 return None
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
47
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
48
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
49 def bedToPicInterval(infile=None):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
50 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
51 Picard tools requiring targets want
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
52 a sam style header which incidentally, MUST be sorted in natural order - not lexicographic order:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
53
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
54 @SQ SN:chrM LN:16571
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
55 @SQ SN:chr1 LN:247249719
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
56 @SQ SN:chr2 LN:242951149
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
57 @SQ SN:chr3 LN:199501827
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
58 @SQ SN:chr4 LN:191273063
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
59 added to the start of what looks like a bed style file
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
60 chr1 67052400 67052451 - CCDS635.1_cds_0_0_chr1_67052401_r
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
61 chr1 67060631 67060788 - CCDS635.1_cds_1_0_chr1_67060632_r
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
62 chr1 67065090 67065317 - CCDS635.1_cds_2_0_chr1_67065091_r
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
63 chr1 67066082 67066181 - CCDS635.1_cds_3_0_chr1_67066083_r
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
64
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
65
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
66 see http://genome.ucsc.edu/FAQ/FAQtracks.html#tracks1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
67 we need to add 1 to start coordinates on the way through - but length calculations are easier
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
68 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
69 # bedToPicard.py
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
70 # ross lazarus October 2010
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
71 # LGPL
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
72 # for Rgenetics
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
73
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
74 def getFlen(bedfname=None):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
75 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
76 find all features in a BED file and sum their lengths
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
77 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
78 features = {}
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
79 try:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
80 infile = open(bedfname,'r')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
81 except:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
82 print '###ERROR: getFlen unable to open bedfile %s' % bedfname
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
83 sys.exit(1)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
84 for i,row in enumerate(infile):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
85 if row[0] == '@': # shouldn't happen given a bed file!
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
86 print 'row %d=%s - should NOT start with @!' % (i,row)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
87 sys.exit(1)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
88 row = row.strip()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
89 if len(row) > 0:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
90 srow = row.split('\t')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
91 f = srow[0]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
92 spos = srow[1] # zero based from UCSC so no need to add 1 - eg 0-100 is 100 bases numbered 0-99 (!)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
93 epos = srow[2] # see http://genome.ucsc.edu/FAQ/FAQtracks.html#tracks1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
94 flen = int(epos) - int(spos)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
95 features.setdefault(f,0)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
96 features[f] += flen
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
97 infile.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
98 return features
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
99
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
100 def keynat(string):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
101 '''
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
102 borrowed from http://code.activestate.com/recipes/285264-natural-string-sorting/
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
103 A natural sort helper function for sort() and sorted()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
104 without using regular expressions or exceptions.
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
105
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
106 >>> items = ('Z', 'a', '10th', '1st', '9')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
107 >>> sorted(items)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
108 ['10th', '1st', '9', 'Z', 'a']
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
109 >>> sorted(items, key=keynat)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
110 ['1st', '9', '10th', 'a', 'Z']
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
111 '''
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
112 it = type(1)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
113 r = []
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
114 for c in string:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
115 if c.isdigit():
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
116 d = int(c)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
117 if r and type( r[-1] ) == it:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
118 r[-1] = r[-1] * 10 + d
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
119 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
120 r.append(d)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
121 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
122 r.append(c.lower())
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
123 return r
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
124
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
125 def writePic(outfname=None,bedfname=None):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
126 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
127 collect header info and rewrite bed with header for picard
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
128 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
129 featlen = getFlen(bedfname=bedfname)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
130 try:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
131 outf = open(outfname,'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
132 except:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
133 print '###ERROR: writePic unable to open output picard file %s' % outfname
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
134 sys.exit(1)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
135 infile = open(bedfname,'r') # already tested in getFlen
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
136 k = featlen.keys()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
137 fk = sorted(k, key=keynat)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
138 header = ['@SQ\tSN:%s\tLN:%d' % (x,featlen[x]) for x in fk]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
139 outf.write('\n'.join(header))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
140 outf.write('\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
141 for row in infile:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
142 row = row.strip()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
143 if len(row) > 0: # convert zero based start coordinate to 1 based
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
144 srow = row.split('\t')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
145 srow[1] = '%d' % (int(srow[1])+1) # see http://genome.ucsc.edu/FAQ/FAQtracks.html#tracks1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
146 outf.write('\t'.join(srow))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
147 outf.write('\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
148 outf.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
149 infile.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
150
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
151
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
152
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
153 # bedToPicInterval starts here
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
154 fd,outf = tempfile.mkstemp(prefix='rgPicardHsMetrics')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
155 writePic(outfname=outf,bedfname=infile)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
156 return outf
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
157
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
158
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
159 def getFileString(fpath, outpath):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
160 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
161 format a nice file size string
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
162 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
163 size = ''
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
164 fp = os.path.join(outpath, fpath)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
165 s = '? ?'
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
166 if os.path.isfile(fp):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
167 n = float(os.path.getsize(fp))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
168 if n > 2**20:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
169 size = ' (%1.1f MB)' % (n/2**20)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
170 elif n > 2**10:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
171 size = ' (%1.1f KB)' % (n/2**10)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
172 elif n > 0:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
173 size = ' (%d B)' % (int(n))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
174 s = '%s %s' % (fpath, size)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
175 return s
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
176
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
177
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
178 def fixPicardOutputs(tempout=None,output_dir=None,log_file=None,html_output=None,progname=None,cl=[],transpose=True):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
179 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
180 picard produces long hard to read tab header files
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
181 make them available but present them transposed for readability
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
182 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
183 rstyle="""<style type="text/css">
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
184 tr.d0 td {background-color: oldlace; color: black;}
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
185 tr.d1 td {background-color: aliceblue; color: black;}
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
186 </style>"""
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
187 cruft = []
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
188 dat = []
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
189 try:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
190 r = open(tempout,'r').readlines()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
191 except:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
192 r = []
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
193 for row in r:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
194 if row.strip() > '':
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
195 srow = row.split('\t')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
196 if row[0] == '#':
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
197 cruft.append(row.strip()) # want strings
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
198 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
199 dat.append(srow) # want lists
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
200
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
201 res = [rstyle,]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
202 res.append(galhtmlprefix % progname)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
203 res.append(galhtmlattr % (progname,timenow()))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
204 flist = os.listdir(output_dir)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
205 pdflist = [x for x in flist if os.path.splitext(x)[-1].lower() == '.pdf']
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
206 if len(pdflist) > 0: # assumes all pdfs come with thumbnail .jpgs
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
207 for p in pdflist:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
208 imghref = '%s.jpg' % os.path.splitext(p)[0] # removes .pdf
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
209 res.append('<table cellpadding="10"><tr><td>\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
210 res.append('<a href="%s"><img src="%s" alt="Click thumbnail to download %s" hspace="10" align="middle"></a>\n' % (p,imghref,p))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
211 res.append('</tr></td></table>\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
212 res.append('<b>Your job produced the following output files.</b><hr/>\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
213 res.append('<table>\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
214 for i,f in enumerate(flist):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
215 fn = os.path.split(f)[-1]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
216 res.append('<tr><td><a href="%s">%s</a></td></tr>\n' % (fn,fn))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
217 res.append('</table><p/>\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
218 if len(cruft) + len(dat) > 0:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
219 res.append('<b>Picard on line resources</b><ul>\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
220 res.append('<li><a href="http://picard.sourceforge.net/index.shtml">Click here for Picard Documentation</a></li>\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
221 res.append('<li><a href="http://picard.sourceforge.net/picard-metric-definitions.shtml">Click here for Picard Metrics definitions</a></li></ul><hr/>\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
222 if transpose:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
223 res.append('<b>Picard output (transposed for readability)</b><hr/>\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
224 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
225 res.append('<b>Picard output</b><hr/>\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
226 res.append('<table cellpadding="3" >\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
227 if len(cruft) > 0:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
228 cres = ['<tr class="d%d"><td>%s</td></tr>' % (i % 2,x) for i,x in enumerate(cruft)]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
229 res += cres
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
230 if len(dat) > 0:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
231 maxrows = 100
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
232 if transpose:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
233 tdat = map(None,*dat) # transpose an arbitrary list of lists
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
234 missing = len(tdat) - maxrows
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
235 tdat = ['<tr class="d%d"><td>%s</td><td>%s</td></tr>\n' % ((i+len(cruft)) % 2,x[0],x[1]) for i,x in enumerate(tdat) if i < maxrows]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
236 if len(tdat) > maxrows:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
237 tdat.append('<tr><td colspan="2">...WARNING: %d rows deleted for sanity...see raw files for all rows</td></tr>' % missing)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
238 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
239 tdat = ['<tr class="d%d"><td>%s</td></tr>\n' % ((i+len(cruft)) % 2,x) for i,x in enumerate(dat) if i < maxrows]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
240 if len(dat) > maxrows:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
241 missing = len(dat) - maxrows
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
242 tdat.append('<tr><td>...WARNING: %d rows deleted for sanity...see raw files for all rows</td></tr>' % missing)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
243 res += tdat
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
244 res.append('</table>\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
245 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
246 res.append('<b>No Picard output found - please consult the Picard log above for an explanation</b>')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
247 l = open(log_file,'r').readlines()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
248 if len(l) > 0:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
249 res.append('<b>Picard log</b><hr/>\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
250 rlog = ['<pre>',]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
251 rlog += l
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
252 rlog.append('</pre>')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
253 res += rlog
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
254 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
255 res.append("Odd, Picard left no log file %s - must have really barfed badly?" % log_file)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
256 res.append('<hr/>The freely available <a href="http://picard.sourceforge.net/command-line-overview.shtml">Picard software</a> \n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
257 res.append( 'generated all outputs reported here, using this command line:<br/>\n<pre>%s</pre>\n' % ''.join(cl))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
258 res.append(galhtmlpostfix)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
259 outf = open(html_output,'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
260 outf.write(''.join(res))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
261 outf.write('\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
262 outf.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
263
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
264 def keynat(string):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
265 '''
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
266 borrowed from http://code.activestate.com/recipes/285264-natural-string-sorting/
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
267 A natural sort helper function for sort() and sorted()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
268 without using regular expressions or exceptions.
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
269
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
270 >>> items = ('Z', 'a', '10th', '1st', '9')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
271 >>> sorted(items)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
272 ['10th', '1st', '9', 'Z', 'a']
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
273 >>> sorted(items, key=keynat)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
274 ['1st', '9', '10th', 'a', 'Z']
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
275 '''
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
276 it = type(1)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
277 r = []
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
278 for c in string:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
279 if c.isdigit():
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
280 d = int(c)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
281 if r and type( r[-1] ) == it:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
282 r[-1] = r[-1] * 10 + d
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
283 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
284 r.append(d)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
285 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
286 r.append(c.lower())
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
287 return r
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
288
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
289 def getFlen(bedfname=None):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
290 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
291 find all features in a BED file and sum their lengths
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
292 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
293 features = {}
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
294 otherHeaders = []
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
295 try:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
296 infile = open(bedfname,'r')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
297 except:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
298 print '###ERROR: getFlen unable to open bedfile %s' % bedfname
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
299 sys.exit(1)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
300 for i,row in enumerate(infile):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
301 if row.startswith('@'): # add to headers if not @SQ
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
302 if not row.startswith('@SQ'):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
303 otherHeaders.append(row)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
304 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
305 row = row.strip()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
306 if row.startswith('#') or row.lower().startswith('browser') or row.lower().startswith('track'):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
307 continue # ignore headers
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
308 srow = row.split('\t')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
309 if len(srow) > 3:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
310 srow = row.split('\t')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
311 f = srow[0]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
312 spos = srow[1] # zero based from UCSC so no need to add 1 - eg 0-100 is 100 bases numbered 0-99 (!)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
313 epos = srow[2] # see http://genome.ucsc.edu/FAQ/FAQtracks.html#tracks1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
314 flen = int(epos) - int(spos)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
315 features.setdefault(f,0)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
316 features[f] += flen
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
317 infile.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
318 fk = features.keys()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
319 fk = sorted(fk, key=keynat)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
320 return features,fk,otherHeaders
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
321
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
322 def bedToPicInterval(infile=None,outfile=None):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
323 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
324 Picard tools requiring targets want
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
325 a sam style header which incidentally, MUST be sorted in natural order - not lexicographic order:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
326
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
327 @SQ SN:chrM LN:16571
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
328 @SQ SN:chr1 LN:247249719
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
329 @SQ SN:chr2 LN:242951149
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
330 @SQ SN:chr3 LN:199501827
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
331 @SQ SN:chr4 LN:191273063
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
332 added to the start of what looks like a bed style file
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
333 chr1 67052400 67052451 - CCDS635.1_cds_0_0_chr1_67052401_r
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
334 chr1 67060631 67060788 - CCDS635.1_cds_1_0_chr1_67060632_r
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
335 chr1 67065090 67065317 - CCDS635.1_cds_2_0_chr1_67065091_r
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
336 chr1 67066082 67066181 - CCDS635.1_cds_3_0_chr1_67066083_r
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
337
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
338 see http://genome.ucsc.edu/FAQ/FAQtracks.html#tracks1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
339 we need to add 1 to start coordinates on the way through - but length calculations are easier
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
340 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
341 # bedToPicard.py
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
342 # ross lazarus October 2010
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
343 # LGPL
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
344 # for Rgenetics
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
345 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
346 collect header info and rewrite bed with header for picard
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
347 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
348 featlen,fk,otherHeaders = getFlen(bedfname=infile)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
349 try:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
350 outf = open(outfile,'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
351 except:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
352 print '###ERROR: writePic unable to open output picard file %s' % outfile
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
353 sys.exit(1)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
354 inf = open(infile,'r') # already tested in getFlen
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
355 header = ['@SQ\tSN:%s\tLN:%d' % (x,featlen[x]) for x in fk]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
356 if len(otherHeaders) > 0:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
357 header += otherHeaders
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
358 outf.write('\n'.join(header))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
359 outf.write('\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
360 for row in inf:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
361 row = row.strip()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
362 if len(row) > 0: # convert zero based start coordinate to 1 based
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
363 if row.startswith('@'):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
364 continue
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
365 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
366 srow = row.split('\t')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
367 srow[1] = '%d' % (int(srow[1])+1) # see http://genome.ucsc.edu/FAQ/FAQtracks.html#tracks1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
368 outf.write('\t'.join(srow))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
369 outf.write('\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
370 outf.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
371 inf.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
372
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
373
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
374 def oRRun(rcmd=[],outdir=None,title='myR',rexe='R'):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
375 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
376 run an r script, lines in rcmd,
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
377 in a temporary directory
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
378 move everything, r script and all back to outdir which will be an html file
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
379
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
380
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
381 # test
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
382 RRun(rcmd=['print("hello cruel world")','q()'],title='test')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
383 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
384 rlog = []
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
385 print '### rexe = %s' % rexe
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
386 assert os.path.isfile(rexe)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
387 rname = '%s.R' % title
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
388 stoname = '%s.R.log' % title
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
389 rfname = rname
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
390 stofname = stoname
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
391 if outdir: # want a specific path
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
392 rfname = os.path.join(outdir,rname)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
393 stofname = os.path.join(outdir,stoname)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
394 try:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
395 os.makedirs(outdir) # might not be there yet...
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
396 except:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
397 pass
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
398 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
399 outdir = tempfile.mkdtemp(prefix=title)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
400 rfname = os.path.join(outdir,rname)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
401 stofname = os.path.join(outdir,stoname)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
402 rmoutdir = True
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
403 f = open(rfname,'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
404 if type(rcmd) == type([]):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
405 f.write('\n'.join(rcmd))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
406 else: # string
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
407 f.write(rcmd)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
408 f.write('\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
409 f.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
410 sto = file(stofname,'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
411 vcl = [rexe,"--vanilla --slave", '<', rfname ]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
412 x = subprocess.Popen(' '.join(vcl),shell=True,stderr=sto,stdout=sto,cwd=outdir)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
413 retval = x.wait()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
414 sto.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
415 rlog = file(stofname,'r').readlines()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
416 rlog.insert(0,'## found R at %s' % rexe)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
417 if outdir <> None:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
418 flist = os.listdir(outdir)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
419 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
420 flist = os.listdir('.')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
421 flist.sort
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
422 flist = [(x,x) for x in flist]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
423 for i,x in enumerate(flist):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
424 if x == rname:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
425 flist[i] = (x,'R script for %s' % title)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
426 elif x == stoname:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
427 flist[i] = (x,'R log for %s' % title)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
428 if False and rmoutdir:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
429 os.removedirs(outdir)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
430 return rlog,flist # for html layout
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
431
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
432
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
433
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
434
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
435 def RRun(rcmd=[],outdir=None,title='myR',tidy=True):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
436 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
437 run an r script, lines in rcmd,
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
438 in a temporary directory
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
439 move everything, r script and all back to outdir which will be an html file
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
440
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
441
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
442 # test
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
443 RRun(rcmd=['print("hello cruel world")','q()'],title='test')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
444 echo "a <- c(5, 5); b <- c(0.5, 0.5)" | cat - RScript.R | R --slave \ --vanilla
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
445 suggested by http://tolstoy.newcastle.edu.au/R/devel/05/09/2448.html
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
446 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
447 killme = string.punctuation + string.whitespace
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
448 trantab = string.maketrans(killme,'_'*len(killme))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
449 title = title.translate(trantab)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
450 rlog = []
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
451 tempout=False
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
452 rname = '%s.R' % title
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
453 stoname = '%s.R.log' % title
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
454 cwd = os.getcwd()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
455 if outdir: # want a specific path
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
456 try:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
457 os.makedirs(outdir) # might not be there yet...
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
458 except:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
459 pass
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
460 os.chdir(outdir)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
461 if type(rcmd) == type([]):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
462 script = '\n'.join(rcmd)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
463 else: # string
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
464 script = rcmd
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
465 sto = file(stoname,'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
466 rscript = file(rname,'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
467 rscript.write(script)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
468 rscript.write('\n#R script autogenerated by rgenetics/rgutils.py on %s\n' % timenow())
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
469 rscript.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
470 vcl = '%s --slave --vanilla < %s' % (rexe,rname)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
471 if outdir:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
472 x = subprocess.Popen(vcl,shell=True,stderr=sto,stdout=sto,cwd=outdir)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
473 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
474 x = subprocess.Popen(vcl,shell=True,stderr=sto,stdout=sto)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
475 retval = x.wait()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
476 sto.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
477 rlog = file(stoname,'r').readlines()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
478 if retval <> 0:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
479 rlog.insert(0,'Nonzero exit code = %d' % retval) # indicate failure
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
480 if outdir:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
481 flist = os.listdir(outdir)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
482 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
483 flist = os.listdir(os.getcwd())
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
484 flist.sort
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
485 flist = [(x,x) for x in flist]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
486 for i,x in enumerate(flist):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
487 if x == rname:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
488 flist[i] = (x,'R script for %s' % title)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
489 elif x == stoname:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
490 flist[i] = (x,'R log for %s' % title)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
491 if outdir:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
492 os.chdir(cwd)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
493 return rlog,flist # for html layout
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
494
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
495 def runPlink(bfn='bar',ofn='foo',logf=None,plinktasks=[],cd='./',vclbase = []):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
496 """run a series of plink tasks and append log results to stdout
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
497 vcl has a list of parameters for the spawnv
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
498 common settings can all go in the vclbase list and are added to each plinktask
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
499 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
500 # root for all
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
501 fplog,plog = tempfile.mkstemp()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
502 if type(logf) == type(' '): # open otherwise assume is file - ugh I'm in a hurry
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
503 mylog = file(logf,'a+')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
504 else:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
505 mylog = logf
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
506 mylog.write('## Rgenetics: http://rgenetics.org Galaxy Tools rgQC.py Plink runner\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
507 for task in plinktasks: # each is a list
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
508 vcl = vclbase + task
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
509 sto = file(plog,'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
510 x = subprocess.Popen(' '.join(vcl),shell=True,stdout=sto,stderr=sto,cwd=cd)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
511 retval = x.wait()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
512 sto.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
513 try:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
514 lplog = file(plog,'r').read()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
515 mylog.write(lplog)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
516 os.unlink(plog) # no longer needed
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
517 except:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
518 mylog.write('### %s Strange - no std out from plink when running command line\n%s' % (timenow(),' '.join(vcl)))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
519
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
520 def pruneLD(plinktasks=[],cd='./',vclbase = []):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
521 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
522 plink blathers when doing pruning - ignore
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
523 Linkage disequilibrium based SNP pruning
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
524 if a million snps in 3 billion base pairs, have mean 3k spacing
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
525 assume 40-60k of ld in ceu, a window of 120k width is about 40 snps
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
526 so lots more is perhaps less efficient - each window computational cost is
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
527 ON^2 unless the code is smart enough to avoid unecessary computation where
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
528 allele frequencies make it impossible to see ld > the r^2 cutoff threshold
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
529 So, do a window and move forward 20?
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
530 The fine Plink docs at http://pngu.mgh.harvard.edu/~purcell/plink/summary.shtml#prune
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
531 reproduced below
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
532
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
533 Sometimes it is useful to generate a pruned subset of SNPs that are in approximate linkage equilibrium with each other. This can be achieved via two commands:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
534 --indep which prunes based on the variance inflation factor (VIF), which recursively removes SNPs within a sliding window; second, --indep-pairwise which is
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
535 similar, except it is based only on pairwise genotypic correlation.
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
536
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
537 Hint The output of either of these commands is two lists of SNPs: those that are pruned out and those that are not. A separate command using the --extract or
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
538 --exclude option is necessary to actually perform the pruning.
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
539
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
540 The VIF pruning routine is performed:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
541 plink --file data --indep 50 5 2
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
542
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
543 will create files
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
544
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
545 plink.prune.in
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
546 plink.prune.out
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
547
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
548 Each is a simlpe list of SNP IDs; both these files can subsequently be specified as the argument for
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
549 a --extract or --exclude command.
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
550
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
551 The parameters for --indep are: window size in SNPs (e.g. 50), the number of SNPs to shift the
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
552 window at each step (e.g. 5), the VIF threshold. The VIF is 1/(1-R^2) where R^2 is the multiple correlation coefficient for a SNP being regressed on all other
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
553 SNPs simultaneously. That is, this considers the correlations between SNPs but also between linear combinations of SNPs. A VIF of 10 is often taken to represent
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
554 near collinearity problems in standard multiple regression analyses (i.e. implies R^2 of 0.9). A VIF of 1 would imply that the SNP is completely independent of
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
555 all other SNPs. Practically, values between 1.5 and 2 should probably be used; particularly in small samples, if this threshold is too low and/or the window
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
556 size is too large, too many SNPs may be removed.
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
557
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
558 The second procedure is performed:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
559 plink --file data --indep-pairwise 50 5 0.5
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
560
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
561 This generates the same output files as the first version; the only difference is that a
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
562 simple pairwise threshold is used. The first two parameters (50 and 5) are the same as above (window size and step); the third parameter represents the r^2
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
563 threshold. Note: this represents the pairwise SNP-SNP metric now, not the multiple correlation coefficient; also note, this is based on the genotypic
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
564 correlation, i.e. it does not involve phasing.
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
565
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
566 To give a concrete example: the command above that specifies 50 5 0.5 would a) consider a
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
567 window of 50 SNPs, b) calculate LD between each pair of SNPs in the window, b) remove one of a pair of SNPs if the LD is greater than 0.5, c) shift the window 5
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
568 SNPs forward and repeat the procedure.
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
569
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
570 To make a new, pruned file, then use something like (in this example, we also convert the
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
571 standard PED fileset to a binary one):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
572 plink --file data --extract plink.prune.in --make-bed --out pruneddata
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
573 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
574 fplog,plog = tempfile.mkstemp()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
575 alog = []
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
576 alog.append('## Rgenetics: http://rgenetics.org Galaxy Tools rgQC.py Plink pruneLD runner\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
577 for task in plinktasks: # each is a list
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
578 vcl = vclbase + task
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
579 sto = file(plog,'w')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
580 x = subprocess.Popen(' '.join(vcl),shell=True,stdout=sto,stderr=sto,cwd=cd)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
581 retval = x.wait()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
582 sto.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
583 try:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
584 lplog = file(plog,'r').readlines()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
585 lplog = [x for x in lplog if x.find('Pruning SNP') == -1]
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
586 alog += lplog
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
587 alog.append('\n')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
588 os.unlink(plog) # no longer needed
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
589 except:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
590 alog.append('### %s Strange - no std out from plink when running command line\n%s\n' % (timenow(),' '.join(vcl)))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
591 return alog
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
592
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
593 def readMap(mapfile=None,allmarkers=False,rsdict={},c=None,spos=None,epos=None):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
594 """abstract out - keeps reappearing
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
595 """
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
596 mfile = open(mapfile, 'r')
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
597 markers = []
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
598 snpcols = {}
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
599 snpIndex = 0 # in case empty or comment lines
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
600 for rownum,row in enumerate(mfile):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
601 line = row.strip()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
602 if not line or line[0]=='#': continue
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
603 chrom, snp, genpos, abspos = line.split()[:4] # just in case more cols
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
604 try:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
605 abspos = int(abspos)
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
606 except:
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
607 abspos = 0 # stupid framingham data grumble grumble
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
608 if allmarkers or rsdict.get(snp,None) or (chrom == c and (spos <= abspos <= epos)):
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
609 markers.append((chrom,abspos,snp)) # decorate for sort into genomic
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
610 snpcols[snp] = snpIndex # so we know which col to find genos for this marker
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
611 snpIndex += 1
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
612 markers.sort()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
613 rslist = [x[2] for x in markers] # drop decoration
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
614 rsdict = dict(zip(rslist,rslist))
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
615 mfile.close()
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
616 return markers,snpcols,rslist,rsdict
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
617
9071e359b9a3 Uploaded
xuebing
parents:
diff changeset
618