|
0
|
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()
|
