Mercurial > repos > xuebing > sharplabtool

comparison tools/rgenetics/rgGRR.py @ 0:9071e359b9a3

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author xuebing
date Fri, 09 Mar 2012 19:37:19 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:9071e359b9a3
1 """
2 # july 2009: Need to see outliers so need to draw them last?
3 # could use clustering on the zscores to guess real relationships for unrelateds
4 # but definitely need to draw last
5 # added MAX_SHOW_ROWS to limit the length of the main report page
6 # Changes for Galaxy integration
7 # added more robust knuth method for one pass mean and sd
8 # no difference really - let's use scipy.mean() and scipy.std() instead...
9 # fixed labels and changed to .xls for outlier reports so can open in excel
10 # interesting - with a few hundred subjects, 5k gives good resolution
11 # and 100k gives better but not by much
12 # TODO remove non autosomal markers
13 # TODO it would be best if label had the zmean and zsd as these are what matter for
14 # outliers rather than the group mean/sd
15 # mods to rgGRR.py from channing CVS which John Ziniti has rewritten to produce SVG plots
16 # to make a Galaxy tool - we need the table of mean and SD for interesting pairs, the SVG and the log
17 # so the result should be an HTML file
18
19 # rgIBS.py
20 # use a random subset of markers for a quick ibs
21 # to identify sample dups and closely related subjects
22 # try snpMatrix and plink and see which one works best for us?
23 # abecasis grr plots mean*sd for every subject to show clusters
24 # mods june 23 rml to avoid non-autosomal markers
25 # we seem to be distinguishing parent-child by gender - 2 clouds!
26
27
28 snpMatrix from David Clayton has:
29 ibs.stats function to calculate the identity-by-state stats of a group of samples
30 Description
31 Given a snp.matrix-class or a X.snp.matrix-class object with N samples, calculates some statistics
32 about the relatedness of every pair of samples within.
33
34 Usage
35 ibs.stats(x)
36 8 ibs.stats
37 Arguments
38 x a snp.matrix-class or a X.snp.matrix-class object containing N samples
39 Details
40 No-calls are excluded from consideration here.
41 Value
42 A data.frame containing N(N - 1)/2 rows, where the row names are the sample name pairs separated
43 by a comma, and the columns are:
44 Count count of identical calls, exclusing no-calls
45 Fraction fraction of identical calls comparied to actual calls being made in both samples
46 Warning
47 In some applications, it may be preferable to subset a (random) selection of SNPs first - the
48 calculation
49 time increases as N(N - 1)M/2 . Typically for N = 800 samples and M = 3000 SNPs, the
50 calculation time is about 1 minute. A full GWA scan could take hours, and quite unnecessary for
51 simple applications such as checking for duplicate or related samples.
52 Note
53 This is mostly written to find mislabelled and/or duplicate samples.
54 Illumina indexes their SNPs in alphabetical order so the mitochondria SNPs comes first - for most
55 purpose it is undesirable to use these SNPs for IBS purposes.
56 TODO: Worst-case S4 subsetting seems to make 2 copies of a large object, so one might want to
57 subset before rbind(), etc; a future version of this routine may contain a built-in subsetting facility
58 """
59 import sys,os,time,random,string,copy,optparse
60
61 try:
62 set
63 except NameError:
64 from Sets import Set as set
65
66 from rgutils import timenow,plinke
67
68 import plinkbinJZ
69
70
71 opts = None
72 verbose = False
73
74 showPolygons = False
75
76 class NullDevice:
77 def write(self, s):
78 pass
79
80 tempstderr = sys.stderr # save
81 #sys.stderr = NullDevice()
82 # need to avoid blather about deprecation and other strange stuff from scipy
83 # the current galaxy job runner assumes that
84 # the job is in error if anything appears on sys.stderr
85 # grrrrr. James wants to keep it that way instead of using the
86 # status flag for some strange reason. Presumably he doesn't use R or (in this case, scipy)
87 import numpy
88 import scipy
89 from scipy import weave
90
91
92 sys.stderr=tempstderr
93
94
95 PROGNAME = os.path.split(sys.argv[0])[-1]
96 X_AXIS_LABEL = 'Mean Alleles Shared'
97 Y_AXIS_LABEL = 'SD Alleles Shared'
98 LEGEND_ALIGN = 'topleft'
99 LEGEND_TITLE = 'Relationship'
100 DEFAULT_SYMBOL_SIZE = 1.0 # default symbol size
101 DEFAULT_SYMBOL_SIZE = 0.5 # default symbol size
102
103 ### Some colors for R/rpy
104 R_BLACK = 1
105 R_RED = 2
106 R_GREEN = 3
107 R_BLUE = 4
108 R_CYAN = 5
109 R_PURPLE = 6
110 R_YELLOW = 7
111 R_GRAY = 8
112
113 ### ... and some point-styles
114
115 ###
116 PLOT_HEIGHT = 600
117 PLOT_WIDTH = 1150
118
119
120 #SVG_COLORS = ('black', 'darkblue', 'blue', 'deepskyblue', 'firebrick','maroon','crimson')
121 #SVG_COLORS = ('cyan','dodgerblue','mediumpurple', 'fuchsia', 'red','gold','gray')
122 SVG_COLORS = ('cyan','dodgerblue','mediumpurple','forestgreen', 'lightgreen','gold','gray')
123 # dupe,parentchild,sibpair,halfsib,parents,unrel,unkn
124 #('orange', 'red', 'green', 'chartreuse', 'blue', 'purple', 'gray')
125
126 OUTLIERS_HEADER_list = ['Mean','Sdev','ZMean','ZSdev','FID1','IID1','FID2','IID2','RelMean_M','RelMean_SD','RelSD_M','RelSD_SD','PID1','MID1','PID2','MID2','Ped']
127 OUTLIERS_HEADER = '\t'.join(OUTLIERS_HEADER_list)
128 TABLE_HEADER='fid1_iid1\tfid2_iid2\tmean\tsdev\tzmean\tzsdev\tgeno\trelcode\tpid1\tmid1\tpid2\tmid2\n'
129
130
131 ### Relationship codes, text, and lookups/mappings
132 N_RELATIONSHIP_TYPES = 7
133 REL_DUPE, REL_PARENTCHILD, REL_SIBS, REL_HALFSIBS, REL_RELATED, REL_UNRELATED, REL_UNKNOWN = range(N_RELATIONSHIP_TYPES)
134 REL_LOOKUP = {
135 REL_DUPE: ('dupe', R_BLUE, 1),
136 REL_PARENTCHILD: ('parentchild', R_YELLOW, 1),
137 REL_SIBS: ('sibpairs', R_RED, 1),
138 REL_HALFSIBS: ('halfsibs', R_GREEN, 1),
139 REL_RELATED: ('parents', R_PURPLE, 1),
140 REL_UNRELATED: ('unrelated', R_CYAN, 1),
141 REL_UNKNOWN: ('unknown', R_GRAY, 1),
142 }
143 OUTLIER_STDEVS = {
144 REL_DUPE: 2,
145 REL_PARENTCHILD: 2,
146 REL_SIBS: 2,
147 REL_HALFSIBS: 2,
148 REL_RELATED: 2,
149 REL_UNRELATED: 3,
150 REL_UNKNOWN: 2,
151 }
152 # note now Z can be passed in
153
154 REL_STATES = [REL_LOOKUP[r][0] for r in range(N_RELATIONSHIP_TYPES)]
155 REL_COLORS = SVG_COLORS
156 REL_POINTS = [REL_LOOKUP[r][2] for r in range(N_RELATIONSHIP_TYPES)]
157
158 DEFAULT_MAX_SAMPLE_SIZE = 10000
159
160 REF_COUNT_HOM1 = 3
161 REF_COUNT_HET = 2
162 REF_COUNT_HOM2 = 1
163 MISSING = 0
164 MAX_SHOW_ROWS = 100 # framingham has millions - delays showing output page - so truncate and explain
165 MARKER_PAIRS_PER_SECOND_SLOW = 15000000.0
166 MARKER_PAIRS_PER_SECOND_FAST = 70000000.0
167
168
169 galhtmlprefix = """<?xml version="1.0" encoding="utf-8" ?>
170 <!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
171 <html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">
172 <head>
173 <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
174 <meta name="generator" content="Galaxy %s tool output - see http://g2.trac.bx.psu.edu/" />
175 <title></title>
176 <link rel="stylesheet" href="/static/style/base.css" type="text/css" />
177 </head>
178 <body>
179 <div class="document">
180 """
181
182
183 SVG_HEADER = '''<?xml version="1.0" standalone="no"?>
184 <!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.2//EN" "http://www.w3.org/Graphics/SVG/1.2/DTD/svg12.dtd">
185
186 <svg width="1280" height="800"
187 xmlns="http://www.w3.org/2000/svg" version="1.2"
188 xmlns:xlink="http://www.w3.org/1999/xlink" viewBox="0 0 1280 800" onload="init()">
189
190 <script type="text/ecmascript" xlink:href="/static/scripts/checkbox_and_radiobutton.js"/>
191 <script type="text/ecmascript" xlink:href="/static/scripts/helper_functions.js"/>
192 <script type="text/ecmascript" xlink:href="/static/scripts/timer.js"/>
193 <script type="text/ecmascript">
194 <![CDATA[
195 var checkBoxes = new Array();
196 var radioGroupBandwidth;
197 var colours = ['%s','%s','%s','%s','%s','%s','%s'];
198 function init() {
199 var style = {"font-family":"Arial,Helvetica", "fill":"black", "font-size":12};
200 var dist = 12;
201 var yOffset = 4;
202
203 //A checkBox for each relationship type dupe,parentchild,sibpair,halfsib,parents,unrel,unkn
204 checkBoxes["dupe"] = new checkBox("dupe","checkboxes",20,40,"cbRect","cbCross",true,"Duplicate",style,dist,yOffset,undefined,hideShowLayer);
205 checkBoxes["parentchild"] = new checkBox("parentchild","checkboxes",20,60,"cbRect","cbCross",true,"Parent-Child",style,dist,yOffset,undefined,hideShowLayer);
206 checkBoxes["sibpairs"] = new checkBox("sibpairs","checkboxes",20,80,"cbRect","cbCross",true,"Sib-pairs",style,dist,yOffset,undefined,hideShowLayer);
207 checkBoxes["halfsibs"] = new checkBox("halfsibs","checkboxes",20,100,"cbRect","cbCross",true,"Half-sibs",style,dist,yOffset,undefined,hideShowLayer);
208 checkBoxes["parents"] = new checkBox("parents","checkboxes",20,120,"cbRect","cbCross",true,"Parents",style,dist,yOffset,undefined,hideShowLayer);
209 checkBoxes["unrelated"] = new checkBox("unrelated","checkboxes",20,140,"cbRect","cbCross",true,"Unrelated",style,dist,yOffset,undefined,hideShowLayer);
210 checkBoxes["unknown"] = new checkBox("unknown","checkboxes",20,160,"cbRect","cbCross",true,"Unknown",style,dist,yOffset,undefined,hideShowLayer);
211
212 }
213
214 function hideShowLayer(id, status, label) {
215 var vis = "hidden";
216 if (status) {
217 vis = "visible";
218 }
219 document.getElementById(id).setAttributeNS(null, 'visibility', vis);
220 }
221
222 function showBTT(evt, rel, mm, dm, md, dd, n, mg, dg, lg, hg) {
223 var x = parseInt(evt.pageX)-250;
224 var y = parseInt(evt.pageY)-110;
225 switch(rel) {
226 case 0:
227 fill = colours[rel];
228 relt = "dupe";
229 break;
230 case 1:
231 fill = colours[rel];
232 relt = "parentchild";
233 break;
234 case 2:
235 fill = colours[rel];
236 relt = "sibpairs";
237 break;
238 case 3:
239 fill = colours[rel];
240 relt = "halfsibs";
241 break;
242 case 4:
243 fill = colours[rel];
244 relt = "parents";
245 break;
246 case 5:
247 fill = colours[rel];
248 relt = "unrelated";
249 break;
250 case 6:
251 fill = colours[rel];
252 relt = "unknown";
253 break;
254 default:
255 fill = "cyan";
256 relt = "ERROR_CODE: "+rel;
257 }
258
259 document.getElementById("btRel").textContent = "GROUP: "+relt;
260 document.getElementById("btMean").textContent = "mean="+mm+" +/- "+dm;
261 document.getElementById("btSdev").textContent = "sdev="+dm+" +/- "+dd;
262 document.getElementById("btPair").textContent = "npairs="+n;
263 document.getElementById("btGeno").textContent = "ngenos="+mg+" +/- "+dg+" (min="+lg+", max="+hg+")";
264 document.getElementById("btHead").setAttribute('fill', fill);
265
266 var tt = document.getElementById("btTip");
267 tt.setAttribute("transform", "translate("+x+","+y+")");
268 tt.setAttribute('visibility', 'visible');
269 }
270
271 function showOTT(evt, rel, s1, s2, mean, sdev, ngeno, rmean, rsdev) {
272 var x = parseInt(evt.pageX)-150;
273 var y = parseInt(evt.pageY)-180;
274
275 switch(rel) {
276 case 0:
277 fill = colours[rel];
278 relt = "dupe";
279 break;
280 case 1:
281 fill = colours[rel];
282 relt = "parentchild";
283 break;
284 case 2:
285 fill = colours[rel];
286 relt = "sibpairs";
287 break;
288 case 3:
289 fill = colours[rel];
290 relt = "halfsibs";
291 break;
292 case 4:
293 fill = colours[rel];
294 relt = "parents";
295 break;
296 case 5:
297 fill = colours[rel];
298 relt = "unrelated";
299 break;
300 case 6:
301 fill = colours[rel];
302 relt = "unknown";
303 break;
304 default:
305 fill = "cyan";
306 relt = "ERROR_CODE: "+rel;
307 }
308
309 document.getElementById("otRel").textContent = "PAIR: "+relt;
310 document.getElementById("otS1").textContent = "s1="+s1;
311 document.getElementById("otS2").textContent = "s2="+s2;
312 document.getElementById("otMean").textContent = "mean="+mean;
313 document.getElementById("otSdev").textContent = "sdev="+sdev;
314 document.getElementById("otGeno").textContent = "ngenos="+ngeno;
315 document.getElementById("otRmean").textContent = "relmean="+rmean;
316 document.getElementById("otRsdev").textContent = "relsdev="+rsdev;
317 document.getElementById("otHead").setAttribute('fill', fill);
318
319 var tt = document.getElementById("otTip");
320 tt.setAttribute("transform", "translate("+x+","+y+")");
321 tt.setAttribute('visibility', 'visible');
322 }
323
324 function hideBTT(evt) {
325 document.getElementById("btTip").setAttributeNS(null, 'visibility', 'hidden');
326 }
327
328 function hideOTT(evt) {
329 document.getElementById("otTip").setAttributeNS(null, 'visibility', 'hidden');
330 }
331
332 ]]>
333 </script>
334 <defs>
335 <!-- symbols for check boxes -->
336 <symbol id="cbRect" overflow="visible">
337 <rect x="-5" y="-5" width="10" height="10" fill="white" stroke="dimgray" stroke-width="1" cursor="pointer"/>
338 </symbol>
339 <symbol id="cbCross" overflow="visible">
340 <g pointer-events="none" stroke="black" stroke-width="1">
341 <line x1="-3" y1="-3" x2="3" y2="3"/>
342 <line x1="3" y1="-3" x2="-3" y2="3"/>
343 </g>
344 </symbol>
345 </defs>
346
347 <desc>Developer Works Dynamic Scatter Graph Scaling Example</desc>
348
349 <!-- Now Draw the main X and Y axis -->
350 <g style="stroke-width:1.0; stroke:black; shape-rendering:crispEdges">
351 <!-- X Axis top and bottom -->
352 <path d="M 100 100 L 1250 100 Z"/>
353 <path d="M 100 700 L 1250 700 Z"/>
354
355 <!-- Y Axis left and right -->
356 <path d="M 100 100 L 100 700 Z"/>
357 <path d="M 1250 100 L 1250 700 Z"/>
358 </g>
359
360 <g transform="translate(100,100)">
361
362 <!-- Grid Lines -->
363 <g style="fill:none; stroke:#dddddd; stroke-width:1; stroke-dasharray:2,2; text-anchor:end; shape-rendering:crispEdges">
364
365 <!-- Vertical grid lines -->
366 <line x1="125" y1="0" x2="115" y2="600" />
367 <line x1="230" y1="0" x2="230" y2="600" />
368 <line x1="345" y1="0" x2="345" y2="600" />
369 <line x1="460" y1="0" x2="460" y2="600" />
370 <line x1="575" y1="0" x2="575" y2="600" style="stroke-dasharray:none;" />
371 <line x1="690" y1="0" x2="690" y2="600" />
372 <line x1="805" y1="0" x2="805" y2="600" />
373 <line x1="920" y1="0" x2="920" y2="600" />
374 <line x1="1035" y1="0" x2="1035" y2="600" />
375
376 <!-- Horizontal grid lines -->
377 <line x1="0" y1="60" x2="1150" y2="60" />
378 <line x1="0" y1="120" x2="1150" y2="120" />
379 <line x1="0" y1="180" x2="1150" y2="180" />
380 <line x1="0" y1="240" x2="1150" y2="240" />
381 <line x1="0" y1="300" x2="1150" y2="300" style="stroke-dasharray:none;" />
382 <line x1="0" y1="360" x2="1150" y2="360" />
383 <line x1="0" y1="420" x2="1150" y2="420" />
384 <line x1="0" y1="480" x2="1150" y2="480" />
385 <line x1="0" y1="540" x2="1150" y2="540" />
386 </g>
387
388 <!-- Legend -->
389 <g style="fill:black; stroke:none" font-size="12" font-family="Arial" transform="translate(25,25)">
390 <rect width="160" height="270" style="fill:none; stroke:black; shape-rendering:crispEdges" />
391 <text x="5" y="20" style="fill:black; stroke:none;" font-size="13" font-weight="bold">Given Pair Relationship</text>
392 <rect x="120" y="35" width="10" height="10" fill="%s" stroke="%s" stroke-width="1" cursor="pointer"/>
393 <rect x="120" y="55" width="10" height="10" fill="%s" stroke="%s" stroke-width="1" cursor="pointer"/>
394 <rect x="120" y="75" width="10" height="10" fill="%s" stroke="%s" stroke-width="1" cursor="pointer"/>
395 <rect x="120" y="95" width="10" height="10" fill="%s" stroke="%s" stroke-width="1" cursor="pointer"/>
396 <rect x="120" y="115" width="10" height="10" fill="%s" stroke="%s" stroke-width="1" cursor="pointer"/>
397 <rect x="120" y="135" width="10" height="10" fill="%s" stroke="%s" stroke-width="1" cursor="pointer"/>
398 <rect x="120" y="155" width="10" height="10" fill="%s" stroke="%s" stroke-width="1" cursor="pointer"/>
399 <text x="15" y="195" style="fill:black; stroke:none" font-size="12" font-family="Arial" >Zscore gt 15</text>
400 <circle cx="125" cy="192" r="6" style="stroke:red; fill:gold; fill-opacity:1.0; stroke-width:1;"/>
401 <text x="15" y="215" style="fill:black; stroke:none" font-size="12" font-family="Arial" >Zscore 4 to 15</text>
402 <circle cx="125" cy="212" r="3" style="stroke:gold; fill:gold; fill-opacity:1.0; stroke-width:1;"/>
403 <text x="15" y="235" style="fill:black; stroke:none" font-size="12" font-family="Arial" >Zscore lt 4</text>
404 <circle cx="125" cy="232" r="2" style="stroke:gold; fill:gold; fill-opacity:1.0; stroke-width:1;"/>
405 <g id="checkboxes">
406 </g>
407 </g>
408
409
410 <g style='fill:black; stroke:none' font-size="17" font-family="Arial">
411 <!-- X Axis Labels -->
412 <text x="480" y="660">Mean Alleles Shared</text>
413 <text x="0" y="630" >1.0</text>
414 <text x="277" y="630" >1.25</text>
415 <text x="564" y="630" >1.5</text>
416 <text x="842" y="630" >1.75</text>
417 <text x="1140" y="630" >2.0</text>
418 </g>
419
420 <g transform="rotate(270)" style="fill:black; stroke:none" font-size="17" font-family="Arial">
421 <!-- Y Axis Labels -->
422 <text x="-350" y="-40">SD Alleles Shared</text>
423 <text x="-20" y="-10" >1.0</text>
424 <text x="-165" y="-10" >0.75</text>
425 <text x="-310" y="-10" >0.5</text>
426 <text x="-455" y="-10" >0.25</text>
427 <text x="-600" y="-10" >0.0</text>
428 </g>
429
430 <!-- Plot Title -->
431 <g style="fill:black; stroke:none" font-size="18" font-family="Arial">
432 <text x="425" y="-30">%s</text>
433 </g>
434
435 <!-- One group/layer of points for each relationship type -->
436 '''
437
438 SVG_FOOTER = '''
439 <!-- End of Data -->
440 </g>
441 <g id="btTip" visibility="hidden" style="stroke-width:1.0; fill:black; stroke:none;" font-size="10" font-family="Arial">
442 <rect width="250" height="110" style="fill:silver" rx="2" ry="2"/>
443 <rect id="btHead" width="250" height="20" rx="2" ry="2" />
444 <text id="btRel" y="14" x="85">unrelated</text>
445 <text id="btMean" y="40" x="4">mean=1.5 +/- 0.04</text>
446 <text id="btSdev" y="60" x="4">sdev=0.7 +/- 0.03</text>
447 <text id="btPair" y="80" x="4">npairs=1152</text>
448 <text id="btGeno" y="100" x="4">ngenos=4783 +/- 24 (min=1000, max=5000)</text>
449 </g>
450
451 <g id="otTip" visibility="hidden" style="stroke-width:1.0; fill:black; stroke:none;" font-size="10" font-family="Arial">
452 <rect width="150" height="180" style="fill:silver" rx="2" ry="2"/>
453 <rect id="otHead" width="150" height="20" rx="2" ry="2" />
454 <text id="otRel" y="14" x="40">sibpairs</text>
455 <text id="otS1" y="40" x="4">s1=fid1,iid1</text>
456 <text id="otS2" y="60" x="4">s2=fid2,iid2</text>
457 <text id="otMean" y="80" x="4">mean=1.82</text>
458 <text id="otSdev" y="100" x="4">sdev=0.7</text>
459 <text id="otGeno" y="120" x="4">ngeno=4487</text>
460 <text id="otRmean" y="140" x="4">relmean=1.85</text>
461 <text id="otRsdev" y="160" x="4">relsdev=0.65</text>
462 </g>
463 </svg>
464 '''
465
466
467 DEFAULT_MAX_SAMPLE_SIZE = 5000
468
469 REF_COUNT_HOM1 = 3
470 REF_COUNT_HET = 2
471 REF_COUNT_HOM2 = 1
472 MISSING = 0
473
474 MARKER_PAIRS_PER_SECOND_SLOW = 15000000
475 MARKER_PAIRS_PER_SECOND_FAST = 70000000
476
477 POLYGONS = {
478 REL_UNRELATED: ((1.360, 0.655), (1.385, 0.730), (1.620, 0.575), (1.610, 0.505)),
479 REL_HALFSIBS: ((1.630, 0.500), (1.630, 0.550), (1.648, 0.540), (1.648, 0.490)),
480 REL_SIBS: ((1.660, 0.510), (1.665, 0.560), (1.820, 0.410), (1.820, 0.390)),
481 REL_PARENTCHILD: ((1.650, 0.470), (1.650, 0.490), (1.750, 0.440), (1.750, 0.420)),
482 REL_DUPE: ((1.970, 0.000), (1.970, 0.150), (2.000, 0.150), (2.000, 0.000)),
483 }
484
485 def distance(point1, point2):
486 """ Calculate the distance between two points
487 """
488 (x1,y1) = [float(d) for d in point1]
489 (x2,y2) = [float(d) for d in point2]
490 dx = abs(x1 - x2)
491 dy = abs(y1 - y2)
492 return math.sqrt(dx**2 + dy**2)
493
494 def point_inside_polygon(x, y, poly):
495 """ Determine if a point (x,y) is inside a given polygon or not
496 poly is a list of (x,y) pairs.
497
498 Taken from: http://www.ariel.com.au/a/python-point-int-poly.html
499 """
500
501 n = len(poly)
502 inside = False
503
504 p1x,p1y = poly[0]
505 for i in range(n+1):
506 p2x,p2y = poly[i % n]
507 if y > min(p1y,p2y):
508 if y <= max(p1y,p2y):
509 if x <= max(p1x,p2x):
510 if p1y != p2y:
511 xinters = (y-p1y)*(p2x-p1x)/(p2y-p1y)+p1x
512 if p1x == p2x or x <= xinters:
513 inside = not inside
514 p1x,p1y = p2x,p2y
515 return inside
516
517 def readMap(pedfile):
518 """
519 """
520 mapfile = pedfile.replace('.ped', '.map')
521 marker_list = []
522 if os.path.exists(mapfile):
523 print 'readMap: %s' % (mapfile)
524 fh = file(mapfile, 'r')
525 for line in fh:
526 marker_list.append(line.strip().split())
527 fh.close()
528 print 'readMap: %s markers' % (len(marker_list))
529 return marker_list
530
531 def calcMeanSD(useme):
532 """
533 A numerically stable algorithm is given below. It also computes the mean.
534 This algorithm is due to Knuth,[1] who cites Welford.[2]
535 n = 0
536 mean = 0
537 M2 = 0
538
539 foreach x in data:
540 n = n + 1
541 delta = x - mean
542 mean = mean + delta/n
543 M2 = M2 + delta*(x - mean) // This expression uses the new value of mean
544 end for
545
546 variance_n = M2/n
547 variance = M2/(n - 1)
548 """
549 mean = 0.0
550 M2 = 0.0
551 sd = 0.0
552 n = len(useme)
553 if n > 1:
554 for i,x in enumerate(useme):
555 delta = x - mean
556 mean = mean + delta/(i+1) # knuth uses n+=1 at start
557 M2 = M2 + delta*(x - mean) # This expression uses the new value of mean
558 variance = M2/(n-1) # assume is sample so lose 1 DOF
559 sd = pow(variance,0.5)
560 return mean,sd
561
562
563 def doIBSpy(ped=None,basename='',outdir=None,logf=None,
564 nrsSamples=10000,title='title',pdftoo=0,Zcutoff=2.0):
565 #def doIBS(pedName, title, nrsSamples=None, pdftoo=False):
566 """ started with snpmatrix but GRR uses actual IBS counts and sd's
567 """
568 repOut = [] # text strings to add to the html display
569 refallele = {}
570 tblf = '%s_table.xls' % (title)
571 tbl = file(os.path.join(outdir,tblf), 'w')
572 tbl.write(TABLE_HEADER)
573 svgf = '%s.svg' % (title)
574 svg = file(os.path.join(outdir,svgf), 'w')
575
576 nMarkers = len(ped._markers)
577 if nMarkers < 5:
578 print sys.stderr, '### ERROR - %d is too few markers for reliable estimation in %s - terminating' % (nMarkers,PROGNAME)
579 sys.exit(1)
580 nSubjects = len(ped._subjects)
581 nrsSamples = min(nMarkers, nrsSamples)
582 if opts and opts.use_mito:
583 markers = range(nMarkers)
584 nrsSamples = min(len(markers), nrsSamples)
585 sampleIndexes = sorted(random.sample(markers, nrsSamples))
586 else:
587 autosomals = ped.autosomal_indices()
588 nrsSamples = min(len(autosomals), nrsSamples)
589 sampleIndexes = sorted(random.sample(autosomals, nrsSamples))
590
591 print ''
592 print 'Getting random.sample of %s from %s total' % (nrsSamples, nMarkers)
593 npairs = (nSubjects*(nSubjects-1))/2 # total rows in table
594 newfiles=[svgf,tblf]
595 explanations = ['rgGRR Plot (requires SVG)','Mean by SD alleles shared - %d rows' % npairs]
596 # these go with the output file links in the html file
597 s = 'Reading genotypes for %s subjects and %s markers\n' % (nSubjects, nrsSamples)
598 logf.write(s)
599 minUsegenos = nrsSamples/2 # must have half?
600 nGenotypes = nSubjects*nrsSamples
601 stime = time.time()
602 emptyRows = set()
603 genos = numpy.zeros((nSubjects, nrsSamples), dtype=int)
604 for s in xrange(nSubjects):
605 nValid = 0
606 #getGenotypesByIndices(self, s, mlist, format)
607 genos[s] = ped.getGenotypesByIndices(s, sampleIndexes, format='ref')
608 nValid = sum([1 for g in genos[s] if g])
609 if not nValid:
610 emptyRows.add(s)
611 sub = ped.getSubject(s)
612 print 'All missing for row %d (%s)' % (s, sub)
613 logf.write('All missing for row %d (%s)\n' % (s, sub))
614 rtime = time.time() - stime
615 if verbose:
616 print '@@Read %s genotypes in %s seconds' % (nGenotypes, rtime)
617
618
619 ### Now the expensive part. For each pair of subjects, we get the mean number
620 ### and standard deviation of shared alleles over all of the markers where both
621 ### subjects have a known genotype. Identical subjects should have mean shared
622 ### alleles very close to 2.0 with a standard deviation very close to 0.0.
623 tot = nSubjects*(nSubjects-1)/2
624 nprog = tot/10
625 nMarkerpairs = tot * nrsSamples
626 estimatedTimeSlow = nMarkerpairs/MARKER_PAIRS_PER_SECOND_SLOW
627 estimatedTimeFast = nMarkerpairs/MARKER_PAIRS_PER_SECOND_FAST
628
629 pairs = []
630 pair_data = {}
631 means = [] ## Mean IBS for each pair
632 ngenoL = [] ## Count of comparable genotypes for each pair
633 sdevs = [] ## Standard dev for each pair
634 rels = [] ## A relationship code for each pair
635 zmeans = [0.0 for x in xrange(tot)] ## zmean score for each pair for the relgroup
636 zstds = [0.0 for x in xrange(tot)] ## zstd score for each pair for the relgrp
637 skip = set()
638 ndone = 0 ## How many have been done so far
639
640 logf.write('Calculating %d pairs...\n' % (tot))
641 logf.write('Estimated time is %2.2f to %2.2f seconds ...\n' % (estimatedTimeFast, estimatedTimeSlow))
642
643 t1sum = 0
644 t2sum = 0
645 t3sum = 0
646 now = time.time()
647 scache = {}
648 _founder_cache = {}
649 C_CODE = """
650 #include "math.h"
651 int i;
652 int sumibs = 0;
653 int ssqibs = 0;
654 int ngeno = 0;
655 float mean = 0;
656 float M2 = 0;
657 float delta = 0;
658 float sdev=0;
659 float variance=0;
660 for (i=0; i<nrsSamples; i++) {
661 int a1 = g1[i];
662 int a2 = g2[i];
663 if (a1 != 0 && a2 != 0) {
664 ngeno += 1;
665 int shared = 2-abs(a1-a2);
666 delta = shared - mean;
667 mean = mean + delta/ngeno;
668 M2 += delta*(shared-mean);
669 // yes that second time, the updated mean is used see calcmeansd above;
670 //printf("%d %d %d %d %d %d\\n", i, a1, a2, ngeno, shared, squared);
671 }
672 }
673 if (ngeno > 1) {
674 variance = M2/(ngeno-1);
675 sdev = sqrt(variance);
676 //printf("OK: %d %3.2f %3.2f\\n", ngeno, mean, sdev);
677 }
678 //printf("%d %d %d %1.2f %1.2f\\n", ngeno, sumibs, ssqibs, mean, sdev);
679 result[0] = ngeno;
680 result[1] = mean;
681 result[2] = sdev;
682 return_val = ngeno;
683 """
684 started = time.time()
685 for s1 in xrange(nSubjects):
686 if s1 in emptyRows:
687 continue
688 (fid1,iid1,did1,mid1,sex1,phe1,iid1,d_sid1,m_sid1) = scache.setdefault(s1, ped.getSubject(s1))
689
690 isFounder1 = _founder_cache.setdefault(s1, (did1==mid1))
691 g1 = genos[s1]
692
693 for s2 in xrange(s1+1, nSubjects):
694 if s2 in emptyRows:
695 continue
696 t1s = time.time()
697
698 (fid2,iid2,did2,mid2,sex2,phe2,iid2,d_sid2,m_sid2) = scache.setdefault(s2, ped.getSubject(s2))
699
700 g2 = genos[s2]
701 isFounder2 = _founder_cache.setdefault(s2, (did2==mid2))
702
703 # Determine the relationship for this pair
704 relcode = REL_UNKNOWN
705 if (fid2 == fid1):
706 if iid1 == iid2:
707 relcode = REL_DUPE
708 elif (did2 == did1) and (mid2 == mid1) and did1 != mid1:
709 relcode = REL_SIBS
710 elif (iid1 == mid2) or (iid1 == did2) or (iid2 == mid1) or (iid2 == did1):
711 relcode = REL_PARENTCHILD
712 elif (str(did1) != '0' and (did2 == did1)) or (str(mid1) != '0' and (mid2 == mid1)):
713 relcode = REL_HALFSIBS
714 else:
715 # People in the same family should be marked as some other
716 # form of related. In general, these people will have a
717 # pretty random spread of similarity. This distinction is
718 # probably not very useful most of the time
719 relcode = REL_RELATED
720 else:
721 ### Different families
722 relcode = REL_UNRELATED
723
724 t1e = time.time()
725 t1sum += t1e-t1s
726
727
728 ### Calculate sum(2-abs(a1-a2)) and sum((2-abs(a1-a2))**2) and count
729 ### the number of contributing genotypes. These values are not actually
730 ### calculated here, but instead are looked up in a table for speed.
731 ### FIXME: This is still too slow ...
732 result = [0.0, 0.0, 0.0]
733 ngeno = weave.inline(C_CODE, ['g1', 'g2', 'nrsSamples', 'result'])
734 if ngeno >= minUsegenos:
735 _, mean, sdev = result
736 means.append(mean)
737 sdevs.append(sdev)
738 ngenoL.append(ngeno)
739 pairs.append((s1, s2))
740 rels.append(relcode)
741 else:
742 skip.add(ndone) # signal no comparable genotypes for this pair
743 ndone += 1
744 t2e = time.time()
745 t2sum += t2e-t1e
746 t3e = time.time()
747 t3sum += t3e-t2e
748
749 logme = [ 'T1: %s' % (t1sum), 'T2: %s' % (t2sum), 'T3: %s' % (t3sum),'TOT: %s' % (t3e-now),
750 '%s pairs with no (or not enough) comparable genotypes (%3.1f%%)' % (len(skip),
751 float(len(skip))/float(tot)*100)]
752 logf.write('%s\n' % '\t'.join(logme))
753 ### Calculate mean and standard deviation of scores on a per relationship
754 ### type basis, allowing us to flag outliers for each particular relationship
755 ### type
756 relstats = {}
757 relCounts = {}
758 outlierFiles = {}
759 for relCode, relInfo in REL_LOOKUP.items():
760 relName, relColor, relStyle = relInfo
761 useme = [means[x] for x in xrange(len(means)) if rels[x] == relCode]
762 relCounts[relCode] = len(useme)
763 mm = scipy.mean(useme)
764 ms = scipy.std(useme)
765 useme = [sdevs[x] for x in xrange(len(sdevs)) if rels[x] == relCode]
766 sm = scipy.mean(useme)
767 ss = scipy.std(useme)
768 relstats[relCode] = {'sd':(sm,ss), 'mean':(mm,ms)}
769 s = 'Relstate %s (n=%d): mean(mean)=%3.2f sdev(mean)=%3.2f, mean(sdev)=%3.2f sdev(sdev)=%3.2f\n' % \
770 (relName,relCounts[relCode], mm, ms, sm, ss)
771 logf.write(s)
772
773 ### now fake z scores for each subject like abecasis recommends max(|zmu|,|zsd|)
774 ### within each group, for each pair, z=(groupmean-pairmean)/groupsd
775 available = len(means)
776 logf.write('%d pairs are available of %d\n' % (available, tot))
777 ### s = '\nOutliers:\nrelationship\tzmean\tzsd\tped1\tped2\tmean\tsd\trmeanmean\trmeansd\trsdmean\trsdsd\n'
778 ### logf.write(s)
779 pairnum = 0
780 offset = 0
781 nOutliers = 0
782 cexs = []
783 outlierRecords = dict([(r, []) for r in range(N_RELATIONSHIP_TYPES)])
784 zsdmax = 0
785 for s1 in range(nSubjects):
786 if s1 in emptyRows:
787 continue
788 (fid1,iid1,did1,mid1,sex1,aff1,ok1,d_sid1,m_sid1) = scache[s1]
789 for s2 in range(s1+1, nSubjects):
790 if s2 in emptyRows:
791 continue
792 if pairnum not in skip:
793 ### Get group stats for this relationship
794 (fid2,iid2,did2,mid2,sex2,aff2,ok2,d_sid2,m_sid2) = scache[s2]
795 try:
796 r = rels[offset]
797 except IndexError:
798 logf.write('###OOPS offset %d available %d pairnum %d len(rels) %d', offset, available, pairnum, len(rels))
799 notfound = ('?',('?','0','0'))
800 relInfo = REL_LOOKUP.get(r,notfound)
801 relName, relColor, relStyle = relInfo
802 rmm,rmd = relstats[r]['mean'] # group mean, group meansd alleles shared
803 rdm,rdd = relstats[r]['sd'] # group sdmean, group sdsd alleles shared
804
805 try:
806 zsd = (sdevs[offset] - rdm)/rdd # distance from group mean in group sd units
807 except:
808 zsd = 1
809 if abs(zsd) > zsdmax:
810 zsdmax = zsd # keep for sort scaling
811 try:
812 zmean = (means[offset] - rmm)/rmd # distance from group mean
813 except:
814 zmean = 1
815 zmeans[offset] = zmean
816 zstds[offset] = zsd
817 pid=(s1,s2)
818 zrad = max(zsd,zmean)
819 if zrad < 4:
820 zrad = 2
821 elif 4 < zrad < 15:
822 zrad = 3 # to 9
823 else: # > 15 6=24+
824 zrad=zrad/4
825 zrad = min(zrad,6) # scale limit
826 zrad = max(2,max(zsd,zmean)) # as > 2, z grows
827 pair_data[pid] = (zmean,zsd,r,zrad)
828 if max(zsd,zmean) > Zcutoff: # is potentially interesting
829 mean = means[offset]
830 sdev = sdevs[offset]
831 outlierRecords[r].append((mean, sdev, zmean, zsd, fid1, iid1, fid2, iid2, rmm, rmd, rdm, rdd,did1,mid1,did2,mid2))
832 nOutliers += 1
833 tbl.write('%s_%s\t%s_%s\t%f\t%f\t%f\t%f\t%d\t%s\t%s\t%s\t%s\t%s\n' % \
834 (fid1, iid1, fid2, iid2, mean, sdev, zmean,zsd, ngeno, relName, did1,mid1,did2,mid2))
835 offset += 1
836 pairnum += 1
837 logf.write( 'Outliers: %s\n' % (nOutliers))
838
839 ### Write outlier files for each relationship type
840 repOut.append('<h2>Outliers in tab delimited files linked above are also listed below</h2>')
841 lzsd = round(numpy.log10(zsdmax)) + 1
842 scalefactor = 10**lzsd
843 for relCode, relInfo in REL_LOOKUP.items():
844 relName, _, _ = relInfo
845 outliers = outlierRecords[relCode]
846 if not outliers:
847 continue
848 outliers = [(scalefactor*int(abs(x[3]))+ int(abs(x[2])),x) for x in outliers] # decorate
849 outliers.sort()
850 outliers.reverse() # largest deviation first
851 outliers = [x[1] for x in outliers] # undecorate
852 nrows = len(outliers)
853 truncated = 0
854 if nrows > MAX_SHOW_ROWS:
855 s = '<h3>%s outlying pairs (top %d of %d) from %s</h3><table border="0" cellpadding="3">' % \
856 (relName,MAX_SHOW_ROWS,nrows,title)
857 truncated = nrows - MAX_SHOW_ROWS
858 else:
859 s = '<h3>%s outlying pairs (n=%d) from %s</h3><table border="0" cellpadding="3">' % (relName,nrows,title)
860 repOut.append(s)
861 fhname = '%s_rgGRR_%s_outliers.xls' % (title, relName)
862 fhpath = os.path.join(outdir,fhname)
863 fh = open(fhpath, 'w')
864 newfiles.append(fhname)
865 explanations.append('%s Outlier Pairs %s, N=%d, Cutoff SD=%f' % (relName,title,len(outliers),Zcutoff))
866 fh.write(OUTLIERS_HEADER)
867 s = ''.join(['<th>%s</th>' % x for x in OUTLIERS_HEADER_list])
868 repOut.append('<tr align="center">%s</tr>' % s)
869 for n,rec in enumerate(outliers):
870 #(mean, sdev, zmean, zsd, fid1, iid1, fid2, iid2, rmm, rmd, rdm, rdd) = rec
871 s = '%f\t%f\t%f\t%f\t%s\t%s\t%s\t%s\t%f\t%f\t%f\t%f\t%s\t%s\t%s\t%s\t' % tuple(rec)
872 fh.write('%s%s\n' % (s,relName))
873 # (mean, sdev, zmean, zsd, fid1, iid1, fid2, iid2, rmm, rmd, rdm, rdd, did1,mid1,did2,mid2))
874 s = '''<td>%f</td><td>%f</td><td>%f</td><td>%f</td><td>%s</td><td>%s</td>
875 <td>%s</td><td>%s</td><td>%f</td><td>%f</td><td>%f</td><td>%f</td><td>%s</td><td>%s</td><td>%s</td><td>%s</td>''' % tuple(rec)
876 s = '%s<td>%s</td>' % (s,relName)
877 if n < MAX_SHOW_ROWS:
878 repOut.append('<tr align="center">%s</tr>' % s)
879 if truncated > 0:
880 repOut.append('<H2>WARNING: %d rows truncated - see outlier file for all %d rows</H2>' % (truncated,
881 nrows))
882 fh.close()
883 repOut.append('</table><p>')
884
885 ### Now, draw the plot in jpeg and svg formats, and optionally in the PDF format
886 ### if requested
887 logf.write('Plotting ...')
888 pointColors = [REL_COLORS[rel] for rel in rels]
889 pointStyles = [REL_POINTS[rel] for rel in rels]
890
891 mainTitle = '%s (%s subjects, %d snp)' % (title, nSubjects, nrsSamples)
892 svg.write(SVG_HEADER % (SVG_COLORS[0],SVG_COLORS[1],SVG_COLORS[2],SVG_COLORS[3],SVG_COLORS[4],
893 SVG_COLORS[5],SVG_COLORS[6],SVG_COLORS[0],SVG_COLORS[0],SVG_COLORS[1],SVG_COLORS[1],
894 SVG_COLORS[2],SVG_COLORS[2],SVG_COLORS[3],SVG_COLORS[3],SVG_COLORS[4],SVG_COLORS[4],
895 SVG_COLORS[5],SVG_COLORS[5],SVG_COLORS[6],SVG_COLORS[6],mainTitle))
896 #rpy.r.jpeg(filename='%s.jpg' % (title), width=1600, height=1200, pointsize=12, quality=100, bg='white')
897 #rpy.r.par(mai=(1,1,1,0.5))
898 #rpy.r('par(xaxs="i",yaxs="i")')
899 #rpy.r.plot(means, sdevs, main=mainTitle, ylab=Y_AXIS_LABEL, xlab=X_AXIS_LABEL, cex=cexs, col=pointColors, pch=pointStyles, xlim=(0,2), ylim=(0,2))
900 #rpy.r.legend(LEGEND_ALIGN, legend=REL_STATES, pch=REL_POINTS, col=REL_COLORS, title=LEGEND_TITLE)
901 #rpy.r.grid(nx=10, ny=10, col='lightgray', lty='dotted')
902 #rpy.r.dev_off()
903
904 ### We will now go through each relationship type to partition plot points
905 ### into "bulk" and "outlier" groups. Bulk points will represent common
906 ### mean/sdev pairs and will cover the majority of the points in the plot --
907 ### they will use generic tooltip informtion about all of the pairs
908 ### represented by that point. "Outlier" points will be uncommon pairs,
909 ### with very specific information in their tooltips. It would be nice to
910 ### keep hte total number of plotted points in the SVG representation to
911 ### ~10000 (certainly less than 100000?)
912 pointMap = {}
913 orderedRels = [y[1] for y in reversed(sorted([(relCounts.get(x, 0),x) for x in REL_LOOKUP.keys()]))]
914 # do we really want this? I want out of zone points last and big
915 for relCode in orderedRels:
916 svgColor = SVG_COLORS[relCode]
917 relName, relColor, relStyle = REL_LOOKUP[relCode]
918 svg.write('<g id="%s" style="stroke:%s; fill:%s; fill-opacity:1.0; stroke-width:1;" cursor="pointer">\n' % (relName, svgColor, svgColor))
919 pMap = pointMap.setdefault(relCode, {})
920 nPoints = 0
921 rpairs=[]
922 rgenos=[]
923 rmeans=[]
924 rsdevs=[]
925 rz = []
926 for x,rel in enumerate(rels): # all pairs
927 if rel == relCode:
928 s1,s2 = pairs[x]
929 pid=(s1,s2)
930 zmean,zsd,r,zrad = pair_data[pid][:4]
931 rpairs.append(pairs[x])
932 rgenos.append(ngenoL[x])
933 rmeans.append(means[x])
934 rsdevs.append(sdevs[x])
935 rz.append(zrad)
936 ### Now add the svg point group for this relationship to the svg file
937 for x in range(len(rmeans)):
938 svgX = '%d' % ((rmeans[x] - 1.0) * PLOT_WIDTH) # changed so mean scale is 1-2
939 svgY = '%d' % (PLOT_HEIGHT - (rsdevs[x] * PLOT_HEIGHT)) # changed so sd scale is 0-1
940 s1, s2 = rpairs[x]
941 (fid1,uid1,did1,mid1,sex1,phe1,iid1,d_sid1,m_sid1) = scache[s1]
942 (fid2,uid2,did2,mid2,sex2,phe2,iid2,d_sid2,m_sid2) = scache[s2]
943 ngenos = rgenos[x]
944 nPoints += 1
945 point = pMap.setdefault((svgX, svgY), [])
946 point.append((rmeans[x], rsdevs[x], fid1, iid1, did1, mid1, fid2, iid2, did2, mid2, ngenos,rz[x]))
947 for (svgX, svgY) in pMap:
948 points = pMap[(svgX, svgY)]
949 svgX = int(svgX)
950 svgY = int(svgY)
951 if len(points) > 1:
952 mmean,dmean = calcMeanSD([p[0] for p in points])
953 msdev,dsdev = calcMeanSD([p[1] for p in points])
954 mgeno,dgeno = calcMeanSD([p[-1] for p in points])
955 mingeno = min([p[-1] for p in points])
956 maxgeno = max([p[-1] for p in points])
957 svg.write("""<circle cx="%d" cy="%d" r="2"
958 onmouseover="showBTT(evt, %d, %1.2f, %1.2f, %1.2f, %1.2f, %d, %d, %d, %d, %d)"
959 onmouseout="hideBTT(evt)" />\n""" % (svgX, svgY, relCode, mmean, dmean, msdev, dsdev, len(points), mgeno, dgeno, mingeno, maxgeno))
960 else:
961 mean, sdev, fid1, iid1, did1, mid1, fid2, iid2, did2, mid2, ngenos, zrad = points[0][:12]
962 rmean = float(relstats[relCode]['mean'][0])
963 rsdev = float(relstats[relCode]['sd'][0])
964 if zrad < 4:
965 zrad = 2
966 elif 4 < zrad < 9:
967 zrad = 3 # to 9
968 else: # > 9 5=15+
969 zrad=zrad/3
970 zrad = min(zrad,5) # scale limit
971 if zrad <= 3:
972 svg.write('<circle cx="%d" cy="%d" r="%s" onmouseover="showOTT(evt, %d, \'%s,%s,%s,%s\', \'%s,%s,%s,%s\', %1.2f, %1.2f, %s, %1.2f, %1.2f)" onmouseout="hideOTT(evt)" />\n' % (svgX, svgY, zrad, relCode, fid1, iid1, did1, mid1, fid2, iid2, did2, mid2, mean, sdev, ngenos, rmean, rsdev))
973 else: # highlight pairs a long way from expectation by outlining circle in red
974 svg.write("""<circle cx="%d" cy="%d" r="%s" style="stroke:red; fill:%s; fill-opacity:1.0; stroke-width:1;"
975 onmouseover="showOTT(evt, %d, \'%s,%s,%s,%s\', \'%s,%s,%s,%s\', %1.2f, %1.2f, %s, %1.2f, %1.2f)"
976 onmouseout="hideOTT(evt)" />\n""" % \
977 (svgX, svgY, zrad, svgColor, relCode, fid1, iid1, did1, mid1, fid2, iid2, did2, mid2, mean, sdev, ngenos, rmean, rsdev))
978 svg.write('</g>\n')
979
980 ### Create a pdf as well if indicated on the command line
981 ### WARNING! for framingham share, with about 50M pairs, this is a 5.5GB pdf!
982 ## if pdftoo:
983 ## pdfname = '%s.pdf' % (title)
984 ## rpy.r.pdf(pdfname, 6, 6)
985 ## rpy.r.par(mai=(1,1,1,0.5))
986 ## rpy.r('par(xaxs="i",yaxs="i")')
987 ## rpy.r.plot(means, sdevs, main='%s, %d snp' % (title, nSamples), ylab=Y_AXIS_LABEL, xlab=X_AXIS_LABEL, cex=cexs, col=pointColors, pch=pointStyles, xlim=(0,2), ylim=(0,2))
988 ## rpy.r.legend(LEGEND_ALIGN, legend=REL_STATES, pch=REL_POINTS, col=REL_COLORS, title=LEGEND_TITLE)
989 ## rpy.r.grid(nx=10, ny=10, col='lightgray', lty='dotted')
990 ## rpy.r.dev_off()
991
992 ### Draw polygons
993 if showPolygons:
994 svg.write('<g id="polygons" cursor="pointer">\n')
995 for rel, poly in POLYGONS.items():
996 points = ' '.join(['%s,%s' % ((p[0]-1.0)*float(PLOT_WIDTH), (PLOT_HEIGHT - p[1]*PLOT_HEIGHT)) for p in poly])
997 svg.write('<polygon points="%s" fill="transparent" style="stroke:%s; stroke-width:1"/>\n' % (points, SVG_COLORS[rel]))
998 svg.write('</g>\n')
999
1000
1001 svg.write(SVG_FOOTER)
1002 svg.close()
1003 return newfiles,explanations,repOut
1004
1005 def doIBS(n=100):
1006 """parse parameters from galaxy
1007 expect 'input pbed path' 'basename' 'outpath' 'title' 'logpath' 'n'
1008 <command interpreter="python">
1009 rgGRR.py $i.extra_files_path/$i.metadata.base_name "$i.metadata.base_name"
1010 '$out_file1' '$out_file1.files_path' "$title1" '$n' '$Z'
1011 </command>
1012
1013 """
1014 u="""<command interpreter="python">
1015 rgGRR.py $i.extra_files_path/$i.metadata.base_name "$i.metadata.base_name"
1016 '$out_file1' '$out_file1.files_path' "$title1" '$n' '$Z'
1017 </command>
1018 """
1019
1020
1021 if len(sys.argv) < 7:
1022 print >> sys.stdout, 'Need pbed inpath, basename, out_htmlname, outpath, title, logpath, nSNP, Zcutoff on command line please'
1023 print >> sys.stdout, u
1024 sys.exit(1)
1025 ts = '%s%s' % (string.punctuation,string.whitespace)
1026 ptran = string.maketrans(ts,'_'*len(ts))
1027 inpath = sys.argv[1]
1028 ldinpath = os.path.split(inpath)[0]
1029 basename = sys.argv[2]
1030 outhtml = sys.argv[3]
1031 newfilepath = sys.argv[4]
1032 title = sys.argv[5].translate(ptran)
1033 logfname = 'Log_%s.txt' % title
1034 logpath = os.path.join(newfilepath,logfname) # log was a child - make part of html extra_files_path zoo
1035 n = int(sys.argv[6])
1036 try:
1037 Zcutoff = float(sys.argv[7])
1038 except:
1039 Zcutoff = 2.0
1040 try:
1041 os.makedirs(newfilepath)
1042 except:
1043 pass
1044 logf = file(logpath,'w')
1045 efp,ibase_name = os.path.split(inpath) # need to use these for outputs in files_path
1046 ped = plinkbinJZ.BPed(inpath)
1047 ped.parse(quick=True)
1048 if ped == None:
1049 print >> sys.stderr, '## doIBSpy problem - cannot open %s or %s - cannot run' % (ldreduced,basename)
1050 sys.exit(1)
1051 newfiles,explanations,repOut = doIBSpy(ped=ped,basename=basename,outdir=newfilepath,
1052 logf=logf,nrsSamples=n,title=title,pdftoo=0,Zcutoff=Zcutoff)
1053 logf.close()
1054 logfs = file(logpath,'r').readlines()
1055 lf = file(outhtml,'w')
1056 lf.write(galhtmlprefix % PROGNAME)
1057 # this is a mess. todo clean up - should each datatype have it's own directory? Yes
1058 # probably. Then titles are universal - but userId libraries are separate.
1059 s = '<div>Output from %s run at %s<br>\n' % (PROGNAME,timenow())
1060 lf.write('<h4>%s</h4>\n' % s)
1061 fixed = ["'%s'" % x for x in sys.argv] # add quotes just in case
1062 s = 'If you need to rerun this analysis, the command line was\n<pre>%s</pre>\n</div>' % (' '.join(fixed))
1063 lf.write(s)
1064 # various ways of displaying svg - experiments related to missing svg mimetype on test (!)
1065 #s = """<object data="%s" type="image/svg+xml" width="%d" height="%d">
1066 # <embed src="%s" type="image/svg+xml" width="%d" height="%d" />
1067 # </object>""" % (newfiles[0],PLOT_WIDTH,PLOT_HEIGHT,newfiles[0],PLOT_WIDTH,PLOT_HEIGHT)
1068 s = """ <embed src="%s" type="image/svg+xml" width="%d" height="%d" />""" % (newfiles[0],PLOT_WIDTH,PLOT_HEIGHT)
1069 #s = """ <iframe src="%s" type="image/svg+xml" width="%d" height="%d" />""" % (newfiles[0],PLOT_WIDTH,PLOT_HEIGHT)
1070 lf.write(s)
1071 lf.write('<div><h4>Click the links below to save output files and plots</h4><br><ol>\n')
1072 for i in range(len(newfiles)):
1073 if i == 0:
1074 lf.write('<li><a href="%s" type="image/svg+xml" >%s</a></li>\n' % (newfiles[i],explanations[i]))
1075 else:
1076 lf.write('<li><a href="%s">%s</a></li>\n' % (newfiles[i],explanations[i]))
1077 flist = os.listdir(newfilepath)
1078 for fname in flist:
1079 if not fname in newfiles:
1080 lf.write('<li><a href="%s">%s</a></li>\n' % (fname,fname))
1081 lf.write('</ol></div>')
1082 lf.write('<div>%s</div>' % ('\n'.join(repOut))) # repOut is a list of tables
1083 lf.write('<div><hr><h3>Log from this job (also stored in %s)</h3><pre>%s</pre><hr></div>' % (logfname,''.join(logfs)))
1084 lf.write('</body></html>\n')
1085 lf.close()
1086 logf.close()
1087
1088 if __name__ == '__main__':
1089 doIBS()