diff options
Diffstat (limited to 'web/webqtl/markerRegression/CompositeMarkerRegressionPage.py')
| -rwxr-xr-x | web/webqtl/markerRegression/CompositeMarkerRegressionPage.py | 211 |
1 files changed, 211 insertions, 0 deletions
diff --git a/web/webqtl/markerRegression/CompositeMarkerRegressionPage.py b/web/webqtl/markerRegression/CompositeMarkerRegressionPage.py new file mode 100755 index 00000000..6cd8c53a --- /dev/null +++ b/web/webqtl/markerRegression/CompositeMarkerRegressionPage.py @@ -0,0 +1,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 + |