1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
|
# Copyright (C) University of Tennessee Health Science Center, Memphis, TN.
#
# This program is free software: you can redistribute it and/or modify it
# under the terms of the GNU Affero General Public License
# as published by the Free Software Foundation, either version 3 of the
# License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
# See the GNU Affero General Public License for more details.
#
# This program is available from Source Forge: at GeneNetwork Project
# (sourceforge.net/projects/genenetwork/).
#
# Contact Drs. Robert W. Williams and Xiaodong Zhou (2010)
# at rwilliams@uthsc.edu and xzhou15@uthsc.edu
#
#
#
# This module is used by GeneNetwork project (www.genenetwork.org)
#
# Created by GeneNetwork Core Team 2010/08/10
#
# Last updated by GeneNetwork Core Team 2010/10/20
import string
import piddle as pid
import os
from htmlgen import HTMLgen2 as HT
import reaper
from utility import Plot
from base.templatePage import templatePage
from base import webqtlConfig
from utility import webqtlUtil
class CompositeMarkerRegressionPage(templatePage):
def __init__(self, fd):
templatePage.__init__(self, fd)
if not fd.genotype:
fd.readData()
fd.parentsf14regression = fd.formdata.getvalue('parentsf14regression')
weightedRegression = fd.formdata.getvalue('applyVarianceSE')
if fd.parentsf14regression and fd.genotype_2:
_genotype = fd.genotype_2
else:
_genotype = fd.genotype_1
_strains, _vals, _vars, N = fd.informativeStrains(_genotype.prgy, weightedRegression)
self.data = fd
if self.data.identification:
heading2 = HT.Paragraph('Trait ID: %s' % self.data.identification)
heading2.__setattr__("class","subtitle")
self.dict['title'] = '%s: Composite Regression' % self.data.identification
else:
heading2 = ""
self.dict['title'] = 'Composite Regression'
if self.data.traitInfo:
symbol,chromosome,MB = string.split(fd.traitInfo,'\t')
heading3 = HT.Paragraph('[ ',HT.Strong(HT.Italic('%s' % symbol,id="green")),' on Chr %s @ %s Mb ]' % (chromosome,MB))
else:
heading3 = ""
if N < webqtlConfig.KMININFORMATIVE:
heading = "Composite Regression"
detail = ['Fewer than %d strain data were entered for %s data set. No mapping attempted.' % (webqtlConfig.KMININFORMATIVE, self.data.RISet)]
self.error(heading=heading,detail=detail)
return
else:
heading = HT.Paragraph('Trait Data Entered for %s Set' % self.data.RISet)
heading.__setattr__("class","title")
tt = HT.TableLite()
for ii in range(N/2):
tt.append(HT.TR(HT.TD(_strains[2*ii],nowrap="yes"),HT.TD(width=10), HT.TD(_vals[2*ii], nowrap="yes"), \
HT.TD(width=20), HT.TD(_strains[2*ii+1],nowrap="yes"),HT.TD(width=10), HT.TD(_vals[2*ii+1],nowrap="yes")))
if N % 2:
tt.append(HT.TR(HT.TD(_strains[N-1],nowrap="yes"),HT.TD(width=10), HT.TD(_vals[N-1],nowrap="yes"), \
HT.TD(width=20), HT.TD("",nowrap="yes"),HT.TD(width=10), HT.TD("",nowrap="yes")))
indata = tt
mean, median, var, stdev, sem, N = reaper.anova(_vals)
stats = HT.Paragraph('Number of entered values = %d ' % N,HT.BR(),\
'Mean value = %8.3f ' % mean, HT.BR(), \
'Median value = %8.3f ' % median, HT.BR(), \
'Variance = %8.3f ' % var, HT.BR(), \
'Standard Deviation = %8.3f ' % stdev, HT.BR(), \
'Standard Error = %8.3f ' % sem)
self.controlLocus = fd.formdata.getvalue('controlLocus')
heading4 = HT.Blockquote('Control Background Selected for %s Data Set:' % self.data.RISet)
heading4.__setattr__("class","subtitle")
datadiv = HT.TD(heading, HT.Center(heading2,heading3,indata, stats, heading4,HT.Center(self.controlLocus)), width='45%',valign='top', bgColor='#eeeeee')
resultstable = self.GenReport(fd, _genotype, _strains, _vals, _vars)
self.dict['body'] = str(datadiv)+str(resultstable)
def GenReport(self, fd, _genotype, _strains, _vals, _vars= []):
'Create an HTML division which reports any loci which are significantly associated with the submitted trait data.'
if webqtlUtil.ListNotNull(_vars):
qtlresults = _genotype.regression(strains = _strains, trait = _vals, variance = _vars, control = self.controlLocus)
LRSArray = _genotype.permutation(strains = _strains, trait = _vals, variance = _vars, nperm=fd.nperm)
else:
qtlresults = _genotype.regression(strains = _strains, trait = _vals, control = self.controlLocus)
LRSArray = _genotype.permutation(strains = _strains, trait = _vals,nperm=fd.nperm)
myCanvas = pid.PILCanvas(size=(400,300))
#plotBar(myCanvas,10,10,390,290,LRSArray,XLabel='LRS',YLabel='Frequency',title=' Histogram of Permutation Test',identification=fd.identification)
Plot.plotBar(myCanvas, LRSArray,XLabel='LRS',YLabel='Frequency',title=' Histogram of Permutation Test')
filename= webqtlUtil.genRandStr("Reg_")
myCanvas.save(webqtlConfig.IMGDIR+filename, format='gif')
img=HT.Image('/image/'+filename+'.gif',border=0,alt='Histogram of Permutation Test')
if fd.suggestive == None:
fd.suggestive = LRSArray[int(fd.nperm*0.37-1)]
else:
fd.suggestive = float(fd.suggestive)
if fd.significance == None:
fd.significance = LRSArray[int(fd.nperm*0.95-1)]
else:
fd.significance = float(fd.significance)
#########################################
# Permutation Graph
#########################################
permutationHeading = HT.Paragraph('Histogram of Permutation Test')
permutationHeading.__setattr__("class","title")
lrs = HT.Blockquote('Total of %d permutations' % fd.nperm,HT.P(),'Suggestive LRS = %2.2f' % LRSArray[int(fd.nperm*0.37-1)],\
HT.BR(),'Significant LRS = %2.2f' % LRSArray[int(fd.nperm*0.95-1)],HT.BR(),'Highly Significant LRS =%2.2f' % LRSArray[int(fd.nperm*0.99-1)])
permutation = HT.TableLite()
permutation.append(HT.TR(HT.TD(img)),HT.TR(HT.TD(lrs)))
_dispAllLRS = 0
if fd.formdata.getvalue('displayAllLRS'):
_dispAllLRS = 1
qtlresults2 = []
if _dispAllLRS:
filtered = qtlresults[:]
else:
filtered = filter(lambda x, y=fd.suggestive: x.lrs > y, qtlresults)
if len(filtered) == 0:
qtlresults2 = qtlresults[:]
qtlresults2.sort()
filtered = qtlresults2[-10:]
#########################################
# Marker regression report
#########################################
locusFormName = webqtlUtil.genRandStr("fm_")
locusForm = HT.Form(cgi = os.path.join(webqtlConfig.CGIDIR, webqtlConfig.SCRIPTFILE), \
enctype='multipart/form-data', name=locusFormName, submit=HT.Input(type='hidden'))
hddn = {'FormID':'showDatabase','ProbeSetID':'_','database':fd.RISet+"Geno",'CellID':'_', \
'RISet':fd.RISet, 'incparentsf1':'on'}
for key in hddn.keys():
locusForm.append(HT.Input(name=key, value=hddn[key], type='hidden'))
regressionHeading = HT.Paragraph('Marker Regression Report')
regressionHeading.__setattr__("class","title")
if qtlresults2 != []:
report = HT.Blockquote(HT.Font('No association ',color="#FF0000"),HT.Font('with a likelihood ratio statistic greater than %3.1f was found. Here are the top 10 LRSs.' % fd.suggestive,color="#000000"))
else:
report = HT.Blockquote('The following loci in the %s data set have associations with the above trait data.\n' % fd.RISet, HT.P())
report.__setattr__("class","normalsize")
fpText = open('%s.txt' % (webqtlConfig.TMPDIR+filename), 'wb')
textUrl = HT.Href(text = 'Download', url= '/tmp/'+filename+'.txt', target = "_blank", Class='fs12 fwn')
bottomInfo = HT.Paragraph(textUrl, ' result in tab-delimited text format.', HT.BR(), HT.BR(),'LRS values marked with',HT.Font(' * ',color="red"), 'are greater than the significance threshold (specified by you or by permutation test). ' , HT.BR(), HT.BR(), HT.Strong('Additive Effect'), ' is half the difference in the mean phenotype of all cases that are homozygous for one parental allel at this marker minus the mean of all cases that are homozygous for the other parental allele at this marker. ','In the case of %s strains, for example,' % fd.RISet,' A positive additive effect indicates that %s alleles increase trait values. Negative additive effect indicates that %s alleles increase trait values.'% (fd.ppolar,fd.mpolar),Class="fs12 fwn")
c1 = HT.TD('LRS',Class="fs14 fwb ffl b1 cw cbrb")
c2 = HT.TD('Chr',Class="fs14 fwb ffl b1 cw cbrb")
c3 = HT.TD('Mb',Class="fs14 fwb ffl b1 cw cbrb")
c4 = HT.TD('Locus',Class="fs14 fwb ffl b1 cw cbrb")
c5 = HT.TD('Additive Effect',Class="fs14 fwb ffl b1 cw cbrb")
fpText.write('LRS\tChr\tMb\tLocus\tAdditive Effect\n')
hr = HT.TR(c1, c2, c3, c4, c5)
tbl = HT.TableLite(border=0, width="90%", cellpadding=0, cellspacing=0, Class="collap")
tbl.append(hr)
for ii in filtered:
#add by NL 06-22-2011: set LRS to 460 when LRS is infinite,
if ii.lrs==float('inf') or ii.lrs>webqtlConfig.MAXLRS:
LRS=webqtlConfig.MAXLRS #maximum LRS value
else:
LRS=ii.lrs
fpText.write('%2.3f\t%s\t%s\t%s\t%2.3f\n' % (LRS, ii.locus.chr, ii.locus.Mb, ii.locus.name, ii.additive))
if LRS > fd.significance:
c1 = HT.TD('%3.3f*' % LRS, Class="fs13 b1 cbw cr")
else:
c1 = HT.TD('%3.3f' % LRS,Class="fs13 b1 cbw c222")
tbl.append(HT.TR(c1, HT.TD(ii.locus.chr,Class="fs13 b1 cbw c222"), HT.TD(ii.locus.Mb,Class="fs13 b1 cbw c222"), HT.TD(HT.Href(text=ii.locus.name, url = "javascript:showTrait('%s','%s');" % (locusFormName, ii.locus.name), Class='normalsize'), Class="fs13 b1 cbw c222"), HT.TD('%3.3f' % ii.additive,Class="fs13 b1 cbw c222"),bgColor='#eeeeee'))
locusForm.append(tbl)
tbl2 = HT.TableLite(border=0, cellspacing=0, cellpadding=0,width="90%")
tbl2.append(HT.TR(HT.TD(bottomInfo)))
rv=HT.TD(permutationHeading,HT.Center(permutation),regressionHeading,report, HT.Center(locusForm,HT.P(),tbl2,HT.P()),width='55%',valign='top', bgColor='#eeeeee')
return rv
|
