diff options
author | Zachary Sloan | 2013-01-03 18:15:32 -0600 |
---|---|---|
committer | Zachary Sloan | 2013-01-03 18:15:32 -0600 |
commit | 91ed29ef68e8ad29b728f7f574ccc83730d9f7ab (patch) | |
tree | 14ab5e3014f5c0a1d8197d2a5587b627e97133ea | |
parent | fb1c44d09f47fa2e72a90df286ff8bb2a3107f2c (diff) | |
download | genenetwork2-91ed29ef68e8ad29b728f7f574ccc83730d9f7ab.tar.gz |
Began working on marker_regression.py and created Chromosomes class
in species.py
-rwxr-xr-x | wqflask/base/data_set.py | 6 | ||||
-rw-r--r-- | wqflask/base/species.py | 52 | ||||
-rwxr-xr-x | wqflask/wqflask/interval_analyst/GeneUtil.py | 124 | ||||
-rwxr-xr-x | wqflask/wqflask/interval_analyst/IntervalAnalystPage.py | 405 | ||||
-rw-r--r-- | wqflask/wqflask/interval_analyst/__init__.py | 0 | ||||
-rw-r--r-- | wqflask/wqflask/marker_regression/MarkerRegressionPage.py | 1648 | ||||
-rw-r--r-- | wqflask/wqflask/marker_regression/__init__.py | 0 | ||||
-rwxr-xr-x | wqflask/wqflask/marker_regression/marker_regression.py | 1648 |
8 files changed, 3870 insertions, 13 deletions
diff --git a/wqflask/base/data_set.py b/wqflask/base/data_set.py index 36d4acaf..50ef8f57 100755 --- a/wqflask/base/data_set.py +++ b/wqflask/base/data_set.py @@ -85,8 +85,8 @@ class DatasetGroup(object): self.f1list = None self.parlist = None self.allsamples = None - - + + #def read_genotype(self): # self.read_genotype_file() # @@ -158,8 +158,8 @@ class DataSet(object): self.retrieve_other_names() - self.species = species.TheSpecies(self) self.group = DatasetGroup(self) # sets self.group and self.group_id and gets genotype + self.species = species.TheSpecies(self) diff --git a/wqflask/base/species.py b/wqflask/base/species.py index 98941ce5..1fd76772 100644 --- a/wqflask/base/species.py +++ b/wqflask/base/species.py @@ -1,16 +1,48 @@ -from __future__ import print_function, division +from __future__ import absolute_import, print_function, division +import collections + +from flask import Flask, g + +#from MySQLdb import escape_string as escape + +from pprint import pformat as pf class TheSpecies(object): def __init__(self, dataset): self.dataset = dataset + print("self.dataset is:", pf(self.dataset.__dict__)) + self.chromosomes = Chromosomes(self.dataset.group.name) + + #@property + #def chromosomes(self): + # chromosomes = [("All", -1)] + # + # for counter, genotype in enumerate(self.dataset.group.genotype): + # if len(genotype) > 1: + # chromosomes.append((genotype.name, counter)) + # + # print("chromosomes is: ", pf(chromosomes)) + # + # return chromosomes + + + +class Chromosomes(object): + def __init__(self, group_name): + self.chromosomes = collections.OrderedDict() + + results = g.db.execute(""" + Select + Chr_Length.Name, Length from Chr_Length, InbredSet + where + Chr_Length.SpeciesId = InbredSet.SpeciesId AND + InbredSet.Name = %s + Order by OrderId + """, group_name).fetchall() + print("bike:", results) + + for item in results: + self.chromosomes[item.Name] = item.Length - @property - def chromosomes(self): - chromosomes = [("All", -1)] - - for counter, genotype in enumerate(self.dataset.group.genotype): - if len(genotype) > 1: - chromosomes.append((genotype.name, counter)) - - return chromosomes + print("self.chromosomes:", self.chromosomes) diff --git a/wqflask/wqflask/interval_analyst/GeneUtil.py b/wqflask/wqflask/interval_analyst/GeneUtil.py new file mode 100755 index 00000000..43008ecf --- /dev/null +++ b/wqflask/wqflask/interval_analyst/GeneUtil.py @@ -0,0 +1,124 @@ +# 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 + +#Just return a list of dictionaries +#each dictionary contains sub-dictionary +def loadGenes(cursor, chrName, diffCol, startMb, endMb, webqtlDb =None, species='mouse'): + #cursor.execute("desc GeneList") + #results = cursor.fetchall() + #fetchFields = map(lambda X:X[0], results) + fetchFields = ['SpeciesId', 'Id', 'GeneSymbol', 'GeneDescription', 'Chromosome', 'TxStart', 'TxEnd', + 'Strand', 'GeneID', 'NM_ID', 'kgID', 'GenBankID', 'UnigenID', 'ProteinID', 'AlignID', + 'exonCount', 'exonStarts', 'exonEnds', 'cdsStart', 'cdsEnd'] + + ##List All Species in the Gene Table + speciesDict = {} + cursor.execute("select Species.Name, GeneList.SpeciesId from Species, GeneList where \ + GeneList.SpeciesId = Species.Id group by GeneList.SpeciesId") + results = cursor.fetchall() + for item in results: + speciesDict[item[0]] = item[1] + + ##List current Species and other Species + speciesId = speciesDict[species] + otherSpecies = map(lambda X: [X, speciesDict[X]], speciesDict.keys()) + otherSpecies.remove([species, speciesId]) + + cursor.execute("""SELECT %s from GeneList + where + SpeciesId = %d AND Chromosome = '%s' AND + ((TxStart > %f and TxStart <= %f) OR (TxEnd > %f and TxEnd <= %f)) + order by txStart + """ + % (string.join(fetchFields, ", "), speciesId, chrName, startMb, endMb, startMb, endMb)) + results = cursor.fetchall() + GeneList = [] + + if results: + for result in results: + newdict = {} + for j, item in enumerate(fetchFields): + newdict[item] = result[j] + #count SNPs if possible + if diffCol and species=='mouse': + cursor.execute(""" + select + count(*) from BXDSnpPosition + where + Chr = '%s' AND Mb >= %2.6f AND Mb < %2.6f AND + StrainId1 = %d AND StrainId2 = %d + """ % (chrName, newdict["TxStart"], newdict["TxEnd"], diffCol[0], diffCol[1])) + newdict["snpCount"] = cursor.fetchone()[0] + newdict["snpDensity"] = newdict["snpCount"]/(newdict["TxEnd"]-newdict["TxStart"])/1000.0 + else: + newdict["snpDensity"] = newdict["snpCount"] = 0 + + try: + newdict['GeneLength'] = 1000.0*(newdict['TxEnd'] - newdict['TxStart']) + except: + pass + + #load gene from other Species by the same name + for item in otherSpecies: + othSpec, othSpecId = item + newdict2 = {} + + cursor.execute("SELECT %s from GeneList where SpeciesId = %d and geneSymbol= '%s' limit 1" % + (string.join(fetchFields, ", "), othSpecId, newdict["GeneSymbol"])) + resultsOther = cursor.fetchone() + if resultsOther: + for j, item in enumerate(fetchFields): + newdict2[item] = resultsOther[j] + + #count SNPs if possible, could be a separate function + if diffCol and othSpec == 'mouse': + cursor.execute(""" + select + count(*) from BXDSnpPosition + where + Chr = '%s' AND Mb >= %2.6f AND Mb < %2.6f AND + StrainId1 = %d AND StrainId2 = %d + """ % (chrName, newdict["TxStart"], newdict["TxEnd"], diffCol[0], diffCol[1])) + + newdict2["snpCount"] = cursor.fetchone()[0] + newdict2["snpDensity"] = newdict2["snpCount"]/(newdict2["TxEnd"]-newdict2["TxStart"])/1000.0 + else: + newdict2["snpDensity"] = newdict2["snpCount"] = 0 + + try: + newdict2['GeneLength'] = 1000.0*(newdict2['TxEnd'] - newdict2['TxStart']) + except: + pass + + newdict['%sGene' % othSpec] = newdict2 + + GeneList.append(newdict) + + return GeneList + + diff --git a/wqflask/wqflask/interval_analyst/IntervalAnalystPage.py b/wqflask/wqflask/interval_analyst/IntervalAnalystPage.py new file mode 100755 index 00000000..ec9aa29c --- /dev/null +++ b/wqflask/wqflask/interval_analyst/IntervalAnalystPage.py @@ -0,0 +1,405 @@ +# 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 + +from mod_python import apache, util, Cookie +import os +import time +import pyXLWriter as xl + +from htmlgen import HTMLgen2 as HT + +import GeneUtil +from base.templatePage import templatePage +from utility import webqtlUtil +from base import webqtlConfig + + +class IntervalAnalystPage(templatePage): + filename = webqtlUtil.genRandStr("Itan_") + + _scriptfile = "main.py?FormID=intervalAnalyst" + + #A dictionary that lets us map the html form names "txStart_mm6" -> "Mb Start (mm8)" + #the first item is the short name (column headers) and the second item is the long name (dropdown list) + # [short name, long name, category] + columnNames = {"GeneSymbol" : ["Gene", "Gene Name", 'gene'], + "GeneDescription" : ["Description", "Gene Description", 'species'], + 'GeneNeighborsCount' : ["Neighbors", "Gene Neighbors", 'gene'], + 'GeneNeighborsRange' : ["Neighborhood", "Gene Neighborhood (Mb)", 'gene'], + 'GeneNeighborsDensity' : ["Gene Density", "Gene Density (Neighbors/Mb)", 'gene'], + "ProteinID" : ["Prot ID", "Protein ID", 'protein'], + "Chromosome" : ["Chr", "Chromosome", 'species'], + "TxStart" : ["Start", "Mb Start", 'species'], + "TxEnd" : ["End", "Mb End", 'species'], + "GeneLength" : ["Length", "Kb Length", 'species'], + "cdsStart" : ["CDS Start", "Mb CDS Start", 'species'], + "cdsEnd" : ["CDS End", "Mb CDS End", 'species'], + "exonCount" : ["Num Exons", "Exon Count", 'species'], + "exonStarts" : ["Exon Starts", "Exon Starts", 'species'], + "exonEnds" : ["Exon Ends", "Exon Ends", 'species'], + "Strand" : ["Strand", "Strand", 'species'], + "GeneID" : ["Gene ID", "Gene ID", 'species'], + "GenBankID" : ["GenBank", "GenBank ID", 'species'], + "UnigenID" : ["Unigen", "Unigen ID", 'species'], + "NM_ID" : ["NM ID", "NM ID", 'species'], + "kgID" : ["kg ID", "kg ID", 'species'], + "snpCount" : ["SNPs", "SNP Count", 'species'], + "snpDensity" : ["SNP Density", "SNP Density", 'species'], + "lrs" : ["LRS", "Likelihood Ratio Statistic", 'misc'], + "lod" : ["LOD", "Likelihood Odds Ratio", 'misc'], + "pearson" : ["Pearson", "Pearson Product Moment", 'misc'], + "literature" : ["Lit Corr", "Literature Correlation", 'misc'], + } + + ###Species Freeze + speciesFreeze = {'mouse':'mm9', 'rat':'rn3', 'human':'hg19'} + for key in speciesFreeze.keys(): + speciesFreeze[speciesFreeze[key]] = key + + def __init__(self, fd): + + templatePage.__init__(self, fd) + + fd.formdata['remote_ip'] = fd.remote_ip + if not self.openMysql(): + return + + self.species = fd.formdata.getvalue("species", "mouse") + try: + self.startMb = float(fd.formdata.getvalue("startMb")) + except: + self.startMb = 10 + try: + self.endMb = float(fd.formdata.getvalue("endMb")) + except: + self.endMb = self.startMb + 10 + + self.Chr = fd.formdata.getvalue("chromosome", "1") + self.xls = fd.formdata.getvalue("xls", "1") + try: + s1 = int(fd.formdata.getvalue("s1")) + s2 = int(fd.formdata.getvalue("s2")) + self.diffColDefault = self.diffCol = [s1, s2] + except: + self.diffColDefault = self.diffCol = [] + if self.species != 'mouse': + self.diffColDefault = [2, 3]#default is B6 and D2 for other species + + controlFrm, dispFields = self.genControlForm(fd) + geneTable, filename = self.genGeneTable(fd, dispFields) + + infoTD = HT.TD(width=400, valign= "top") + infoTD.append(HT.Paragraph("Interval Analyst : Chr %s" % self.Chr, Class="title"), + HT.Strong("Species : "), self.species.title(), HT.BR(), + HT.Strong("Database : "), "UCSC %s" % self.speciesFreeze[self.species], HT.BR(), + HT.Strong("Range : "), "%2.6f Mb - %2.6f Mb" % (self.startMb, self.endMb), HT.BR(), + ) + if filename: + infoTD.append(HT.BR(), HT.BR(), HT.Href(text="Download", url = "/tmp/" + filename, Class="normalsize") + , " output in MS excel format.") + + mainTable = HT.TableLite(HT.TR(infoTD, HT.TD(controlFrm, Class="doubleBorder", width=400), HT.TD(" ", width="")), cellpadding=10) + mainTable.append(HT.TR(HT.TD(geneTable, colspan=3))) + self.dict['body'] = HT.TD(mainTable) + self.dict['title'] = "Interval Analyst" + + def genGeneTable(self, fd, dispFields): + filename = "" + if self.xls: + #import pyXLWriter as xl + filename = "IntAn_Chr%s_%2.6f-%2.6f" % (self.Chr, self.startMb, self.endMb) + filename += ".xls" + + # Create a new Excel workbook + workbook = xl.Writer(os.path.join(webqtlConfig.TMPDIR, filename)) + worksheet = workbook.add_worksheet() + titleStyle = workbook.add_format(align = 'left', bold = 0, size=18, border = 1, border_color="gray") + headingStyle = workbook.add_format(align = 'center', bold = 1, size=13, fg_color = 0x1E, color="white", border = 1, border_color="gray") + + ##Write title Info + worksheet.write([0, 0], "GeneNetwork Interval Analyst Table", titleStyle) + worksheet.write([1, 0], "%s%s" % (webqtlConfig.PORTADDR, os.path.join(webqtlConfig.CGIDIR, self._scriptfile))) + # + worksheet.write([2, 0], "Date : %s" % time.strftime("%B %d, %Y", time.gmtime())) + worksheet.write([3, 0], "Time : %s GMT" % time.strftime("%H:%M ", time.gmtime())) + worksheet.write([4, 0], "Search by : %s" % fd.formdata['remote_ip']) + worksheet.write([5, 0], "view region : Chr %s %2.6f - %2.6f Mb" % (self.Chr, self.startMb, self.endMb)) + nTitleRow = 7 + + geneTable = HT.TableLite(Class="collap", cellpadding=5) + headerRow = HT.TR(HT.TD(" ", Class="fs13 fwb ffl b1 cw cbrb", width="1")) + if self.xls: + worksheet.write([nTitleRow, 0], "Index", headingStyle) + + for ncol, column in enumerate(dispFields): + if len(column) == 1: + headerRow.append(HT.TD(self.columnNames[column[0]][0], Class="fs13 fwb ffl b1 cw cbrb", NOWRAP=1,align="Center")) + if self.xls: + colTitle = self.columnNames[column[0]][0] + worksheet.write([nTitleRow, ncol+1], colTitle, headingStyle) + worksheet.set_column([ncol+1, ncol+1], 2*len(colTitle)) + else: + headerRow.append(HT.TD(self.columnNames[column[0]][0], HT.BR(), " (%s)" % self.speciesFreeze[column[1]], + Class="fs13 fwb ffl b1 cw cbrb", NOWRAP=1, align="Center")) + if self.xls: + colTitle = self.columnNames[column[0]][0] + " (%s)" % self.speciesFreeze[column[1]] + worksheet.write([nTitleRow, ncol+1], colTitle, headingStyle) + worksheet.set_column([ncol+1, ncol+1], 2*len(colTitle)) + #headerRow.append(HT.TD(self.columnNames[column[0]][0], HT.BR(), + # "(%s %s)" % (column[1].title(), self.speciesFreeze[column[1]]), + # Class="colorBlue", NOWRAP=1, align="Center")) + geneTable.append(headerRow) + + geneCol = GeneUtil.loadGenes(self.cursor, self.Chr, self.diffColDefault, self.startMb, self.endMb, species=self.species) + for gIndex, theGO in enumerate(geneCol): + geneRow = HT.TR(HT.TD(gIndex+1, Class="fs12 fwn b1", align="right")) + if self.xls: + nTitleRow += 1 + worksheet.write([nTitleRow, 0], gIndex + 1) + + for ncol, column in enumerate(dispFields): + if len(column) == 1 or column[1]== self.species: + keyValue = "" + fieldName = column[0] + curSpecies = self.species + curGO = theGO + if theGO.has_key(fieldName): + keyValue = theGO[fieldName] + else: + fieldName , othSpec = column + curSpecies = othSpec + subGO = '%sGene' % othSpec + keyValue = "" + curGO = theGO[subGO] + if theGO[subGO].has_key(fieldName): + keyValue = theGO[subGO][fieldName] + + if self.xls: + worksheet.write([nTitleRow, ncol+1], keyValue) + geneRow.append(self.formatTD(keyValue, fieldName, curSpecies, curGO)) + + geneTable.append(geneRow) + + if self.xls: + workbook.close() + return geneTable, filename + + def formatTD(self, keyValue, fieldName, Species, theGO): + if keyValue is None: + keyValue = "" + if keyValue != "": + if fieldName in ("exonStarts", "exonEnds"): + keyValue = string.replace(keyValue, ',', ' ') + return HT.TD(HT.Span(keyValue, Class="code", Id="green"), width=350, Class="fs12 fwn b1") + elif fieldName in ("GeneDescription"): + if keyValue == "---": + keyValue = "" + return HT.TD(keyValue, Class="fs12 fwn b1", width=300) + elif fieldName in ("GeneSymbol"): + webqtlLink = HT.Href("./%s?cmd=sch&gene=%s&alias=1&species=%s" % (webqtlConfig.SCRIPTFILE, keyValue, Species), + HT.Image("/images/webqtl_search.gif", border=0, valign="top"), target="_blank") + if theGO['GeneID']: + geneSymbolLink = HT.Href(webqtlConfig.NCBI_LOCUSID % theGO['GeneID'], keyValue, Class="normalsize", target="_blank") + else: + geneSymbolLink = keyValue + return HT.TD(webqtlLink, geneSymbolLink, Class="fs12 fwn b1",NOWRAP=1) + elif fieldName == 'UnigenID': + try: + gurl = HT.Href(webqtlConfig.UNIGEN_ID % tuple(string.split(keyValue,'.')[:2]), keyValue, Class="normalsize", target="_blank") + except: + gurl = keyValue + return HT.TD(gurl, Class="fs12 fwn b1",NOWRAP=1) + elif fieldName in ("exonCount", "Chromosome"): + return HT.TD(keyValue, Class="fs12 fwn b1",align="right") + elif fieldName in ("snpCount"): + if keyValue: + snpString = HT.Href(url="%s&chr=%s&start=%s&end=%s&geneName=%s&s1=%d&s2=%d" % (os.path.join(webqtlConfig.CGIDIR, 'main.py?FormID=snpBrowser'), + theGO["Chromosome"], theGO["TxStart"], theGO["TxEnd"], theGO["GeneSymbol"], self.diffColDefault[0], self.diffColDefault[1]), + text=theGO["snpCount"], target="_blank", Class="normalsize") + else: + snpString = keyValue + return HT.TD(snpString, Class="fs12 fwn b1",align="right") + elif fieldName in ("snpDensity", "GeneLength"): + if keyValue: keyValue = "%2.3f" % keyValue + else: keyValue = "" + return HT.TD(keyValue, Class="fs12 fwn b1",align="right") + elif fieldName in ("TxStart", "TxEnd"): + return HT.TD("%2.6f" % keyValue, Class="fs12 fwn b1",align="right") + else: + return HT.TD(keyValue, Class="fs12 fwn b1",NOWRAP=1) + else: + return HT.TD(keyValue, Class="fs12 fwn b1",NOWRAP=1,align="right") + + def genControlForm(self, fd): + ##desc GeneList + self.cursor.execute("Desc GeneList") + GeneListFields = self.cursor.fetchall() + GeneListFields = map(lambda X: X[0], GeneListFields) + + #group columns by category--used for creating the dropdown list of possible columns + categories = {} + for item in self.columnNames.keys(): + category = self.columnNames[item] + if category[-1] not in categories.keys(): + categories[category[-1]] = [item ] + else: + categories[category[-1]] = categories[category[-1]]+[item] + + ##List All Species in the Gene Table + speciesDict = {} + self.cursor.execute("select Species.Name, GeneList.SpeciesId from Species, GeneList where \ + GeneList.SpeciesId = Species.Id group by GeneList.SpeciesId order by Species.Id") + results = self.cursor.fetchall() + speciesField = categories.pop('species', []) + categoriesOrder = ['gene', 'protein'] + for item in results: + specName, specId = item + categoriesOrder.append(specName) + speciesDict[specName] = specId + AppliedField = [] + for item2 in speciesField: + if item2 in GeneListFields: + self.cursor.execute("select %s from GeneList where SpeciesId = %d and %s is not NULL limit 1 " % (item2, specId, item2)) + columnApply = self.cursor.fetchone() + if not columnApply: + continue + elif specName != 'mouse' and item2 in ('snpCount', 'snpDensity'): + continue + else: + pass + AppliedField.append(item2) + categories[specName] = AppliedField + + categoriesOrder += ['misc'] + + ############################################################ + ## Create the list of possible columns for the dropdown list + ############################################################ + allColumnsList = HT.Select(name="allColumns", Class="snpBrowserDropBox") + + for category in categoriesOrder: + allFields = categories[category] + if allFields: + geneOpt = HT.Optgroup(label=category.title()) + for item in allFields: + if category in self.speciesFreeze.keys(): + geneOpt.append(("%s (%s %s)" % (self.columnNames[item][1], category.title(), self.speciesFreeze[category]), + "%s__%s" % (item, self.speciesFreeze[category]))) + else: + geneOpt.append((self.columnNames[item][1], item)) + geneOpt.sort() + allColumnsList.append(geneOpt) + + ###################################### + ## Create the list of selected columns + ###################################### + + #cols contains the value of all the selected columns + submitCols = cols = fd.formdata.getvalue("columns", "default") + + if cols == "default": + if self.species=="mouse": #these are the same columns that are shown on intervalPage.py + cols = ['GeneSymbol', 'GeneDescription', 'Chromosome', 'TxStart', 'Strand', 'GeneLength', 'GeneID', 'NM_ID', 'snpCount', 'snpDensity'] + elif self.species=="rat": + cols = ['GeneSymbol', 'GeneDescription', 'Chromosome', 'TxStart', 'GeneLength', 'Strand', 'GeneID', 'UnigenID'] + else: + #should not happen + cols = [] + else: + if type(cols)==type(""): + cols = [cols] + + colsLst = [] + dispFields = [] + for column in cols: + if submitCols == "default" and column not in ('GeneSymbol') and (column in GeneListFields or column in speciesField): + colsLst.append(("%s (%s %s)" % (self.columnNames[column][1], self.species.title(), self.speciesFreeze[self.species]), + "%s__%s" % (column, self.speciesFreeze[self.species]))) + dispFields.append([column, self.species]) + else: + column2 = column.split("__") + if len(column2) == 1: + colsLst.append((self.columnNames[column2[0]][1], column)) + dispFields.append([column]) + else: + thisSpecies = self.speciesFreeze[column2[1]] + colsLst.append(("%s (%s %s)" % (self.columnNames[column2[0]][1], thisSpecies.title(), column2[1]), + column)) + dispFields.append((column2[0], thisSpecies)) + selectedColumnsList = HT.Select(name="columns", Class="snpBrowserSelectBox", multiple="true", data=colsLst, size=6) + + ########################## + ## Create the columns form + ########################## + columnsForm = HT.Form(name="columnsForm", submit=HT.Input(type='hidden'), cgi=os.path.join(webqtlConfig.CGIDIR, self._scriptfile), enctype="multipart/form-data") + columnsForm.append(HT.Input(type="hidden", name="fromdatabase", value= fd.formdata.getvalue("fromdatabase", "unknown"))) + columnsForm.append(HT.Input(type="hidden", name="species", value=self.species)) + if self.diffCol: + columnsForm.append(HT.Input(type="hidden", name="s1", value=self.diffCol[0])) + columnsForm.append(HT.Input(type="hidden", name="s2", value=self.diffCol[1])) + startBox = HT.Input(type="text", name="startMb", value=self.startMb, size=10) + endBox = HT.Input(type="text", name="endMb", value=self.endMb, size=10) + addButton = HT.Input(type="button", name="add", value="Add", Class="button", onClick="addToList(this.form.allColumns.options[this.form.allColumns.selectedIndex].text, this.form.allColumns.options[this.form.allColumns.selectedIndex].value, this.form.columns)") + removeButton = HT.Input(type="button", name="remove", value="Remove", Class="button", onClick="removeFromList(this.form.columns.selectedIndex, this.form.columns)") + upButton = HT.Input(type="button", name="up", value="Up", Class="button", onClick="swapOptions(this.form.columns.selectedIndex, this.form.columns.selectedIndex-1, this.form.columns)") + downButton = HT.Input(type="button", name="down", value="Down", Class="button", onClick="swapOptions(this.form.columns.selectedIndex, this.form.columns.selectedIndex+1, this.form.columns)") + clearButton = HT.Input(type="button", name="clear", value="Clear", Class="button", onClick="deleteAllElements(this.form.columns)") + submitButton = HT.Input(type="submit", value="Refresh", Class="button", onClick="selectAllElements(this.form.columns)") + + selectChrBox = HT.Select(name="chromosome") + self.cursor.execute(""" + Select + Chr_Length.Name, Length from Chr_Length, Species + where + Chr_Length.SpeciesId = Species.Id AND + Species.Name = '%s' + Order by + Chr_Length.OrderId + """ % self.species) + + results = self.cursor.fetchall() + for chrInfo in results: + selectChrBox.append((chrInfo[0], chrInfo[0])) + selectChrBox.selected.append(self.Chr) + + innerColumnsTable = HT.TableLite(border=0, Class="collap", cellpadding = 2) + innerColumnsTable.append(HT.TR(HT.TD(selectedColumnsList)), + HT.TR(HT.TD(clearButton, removeButton, upButton, downButton))) + columnsTable = HT.TableLite(border=0, cellpadding=2, cellspacing=0) + columnsTable.append(HT.TR(HT.TD(HT.Font("Chr: ", size=-1)), + HT.TD(selectChrBox, submitButton)), + HT.TR(HT.TD(HT.Font("View: ", size=-1)), + HT.TD(startBox, HT.Font("Mb to ", size=-1), endBox, HT.Font("Mb", size=-1))), + HT.TR(HT.TD(HT.Font("Show: ", size=-1)), + HT.TD(allColumnsList, addButton)), + HT.TR(HT.TD(""), + HT.TD(innerColumnsTable))) + columnsForm.append(columnsTable) + #columnsForm.append(HT.Input(type="hidden", name="sort", value=diffCol), + # HT.Input(type="hidden", name="identification", value=identification), + # HT.Input(type="hidden", name="traitInfo", value=traitInfo)) + + return columnsForm, dispFields diff --git a/wqflask/wqflask/interval_analyst/__init__.py b/wqflask/wqflask/interval_analyst/__init__.py new file mode 100644 index 00000000..e69de29b --- /dev/null +++ b/wqflask/wqflask/interval_analyst/__init__.py diff --git a/wqflask/wqflask/marker_regression/MarkerRegressionPage.py b/wqflask/wqflask/marker_regression/MarkerRegressionPage.py new file mode 100644 index 00000000..d02d80b3 --- /dev/null +++ b/wqflask/wqflask/marker_regression/MarkerRegressionPage.py @@ -0,0 +1,1648 @@ +# 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 time +import string +import math +from math import * +import piddle as pid +import sys,os +import httplib, urllib + +from htmlgen import HTMLgen2 as HT +from utility import Plot +from intervalAnalyst import GeneUtil +from base.webqtlTrait import webqtlTrait +from base.templatePage import templatePage +from utility import webqtlUtil +from base import webqtlConfig +from dbFunction import webqtlDatabaseFunction +from base.GeneralObject import GeneralObject + +import reaper +import cPickle +from utility.THCell import THCell +from utility.TDCell import TDCell + +class MarkerRegressionPage(templatePage): + + def __init__(self, fd): + + templatePage.__init__(self, fd) + + if not self.openMysql(): + return + + self.initializeParameters(fd) + + filename= webqtlUtil.genRandStr("Itvl_") + ChrList,ChrNameOrderIdDict,ChrOrderIdNameDict,ChrLengthMbList= self.getChrNameOrderIdLength(RISet=fd.RISet) + + if self.mappingMethodId == '4': # For PLINK + + traitInfoList = string.split(string.strip(fd.identification),':') + probesetName = string.strip(traitInfoList[-1]) + plinkOutputFileName= webqtlUtil.genRandStr("%s_%s_"%(fd.RISet,probesetName)) + + # get related values from fd.allTraitData; the format of 'allTraitValueDict'is {strainName1: value=-0.2...} + fd.readData() + allTraitValueDict = fd.allTraitData + + #automatically generate pheno txt file for PLINK + self.genPhenoTxtFileForPlink(phenoFileName=plinkOutputFileName,RISetName=fd.RISet,probesetName=probesetName, valueDict=allTraitValueDict) + # os.system full path is required for input and output files; specify missing value is -9999 + plink_command = '%splink/plink --noweb --ped %splink/%s.ped --no-fid --no-parents --no-sex --no-pheno --map %splink/%s.map --pheno %s/%s.txt --pheno-name %s --missing-phenotype -9999 --out %s%s --assoc ' % (webqtlConfig.HTMLPATH, webqtlConfig.HTMLPATH, fd.RISet, webqtlConfig.HTMLPATH, fd.RISet, webqtlConfig.TMPDIR, plinkOutputFileName, probesetName, webqtlConfig.TMPDIR, plinkOutputFileName) + + os.system(plink_command) + + if fd.identification: + heading2 = HT.Paragraph('Trait ID: %s' % fd.identification) + heading2.__setattr__("class","subtitle") + self.dict['title'] = '%s: Genome Association' % fd.identification + else: + heading2 = "" + self.dict['title'] = 'Genome Association' + + if fd.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 = "" + + heading = HT.Paragraph('Trait Data Entered for %s Set' % fd.RISet) + heading.__setattr__("class","title") + + # header info part:Trait Data Entered for HLC Set & Trait ID: + headerdiv = HT.TR(HT.TD(heading, heading2,heading3, width='45%',valign='top', align='left', bgColor='#eeeeee')) + + self.ChrList=ChrList # get chr name from '1' to 'X' + self.ChrLengthMbList = ChrLengthMbList + + # build plink result dict based on chr, key is chr name, value is in list type including Snpname, bp and pvalue info + plinkResultDict={} + count,minPvalue,plinkResultDict =self.getPlinkResultDict(outputFileName=plinkOutputFileName,thresholdPvalue=self.pValue,ChrOrderIdNameDict=ChrOrderIdNameDict) + + # if can not find results which are matched with assigned p-value, system info will show up + if count >0: + + #for genome association report table + reportTable="" + # sortable table object + resultstable,tblobj,bottomInfo = self.GenReportForPLINK(ChrNameOrderIdDict=ChrNameOrderIdDict, RISet=fd.RISet,plinkResultDict=plinkResultDict,thresholdPvalue=self.pValue,chrList=self.ChrList) + + # creat object for result table for sort function + objfile = open('%s.obj' % (webqtlConfig.TMPDIR+filename), 'wb') + cPickle.dump(tblobj, objfile) + objfile.close() + + sortby = ("Index", "up") + reportTable =HT.Div(webqtlUtil.genTableObj(tblobj=tblobj, file=filename, sortby=sortby, tableID = "sortable", addIndex = "0"), Id="sortable") + + descriptionTable = HT.TableLite(border=0, cellpadding=0, cellspacing=0) + descriptionTable.append(HT.TR(HT.TD(reportTable, colspan=3))) + descriptionTable.append(HT.TR(HT.TD(HT.BR(),HT.BR()))) + descriptionTable.append(bottomInfo) + + # get each chr's length + self.ChrLengthMbList = map(lambda x: x/1000000.0, self.ChrLengthMbList) # change unit from bp to mb + self.ChrLengthMbSum = reduce(lambda x, y:x+y, self.ChrLengthMbList, 0.0)# get total length of all chrs + if self.ChrLengthMbList: + self.GraphInterval = self.ChrLengthMbSum/(len(self.ChrLengthMbList)*12) #Empirical Mb interval + else: + self.GraphInterval = 1 + + # for human data, there's no CM value + self.ChrLengthCMList = [] + self.ChrLengthCMSum = 0 + + # begin: common part with human data + intCanvas = pid.PILCanvas(size=(self.graphWidth,self.graphHeight)) + gifmap = self.plotIntMappingForPLINK(fd, intCanvas, startMb = self.startMb, endMb = self.endMb, plinkResultDict=plinkResultDict) + + intCanvas.save(os.path.join(webqtlConfig.IMGDIR, filename), format='png') + intImg=HT.Image('/image/'+filename+'.png', border=0, usemap='#WebQTLImageMap') + + TD_LR = HT.TR(HT.TD(HT.Blockquote(gifmap,intImg, HT.P()), bgColor='#eeeeee', height = 200)) + self.dict['body'] = str(headerdiv)+str(TD_LR)+str(resultstable)+str(HT.TR(HT.TD(descriptionTable))) + + else: + heading = "Genome Association" + detail = ['There is no association with marker that meets this criteria. Please provide a less stringend threshold. The minimun p-value is %s.'%minPvalue] + self.error(heading=heading,detail=detail) + return + + elif self.mappingMethodId == '1': # QTLreaper result + 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) + + if fd.identification: + heading2 = HT.Paragraph('Trait ID: %s' % fd.identification) + heading2.__setattr__("class","subtitle") + self.dict['title'] = '%s: Genome Association' % fd.identification + else: + heading2 = "" + self.dict['title'] = 'Genome Association' + + if fd.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 = "Genome Association" + detail = ['Fewer than %d strain data were entered for %s data set. No mapping attempted.' % (webqtlConfig.KMININFORMATIVE, fd.RISet)] + self.error(heading=heading,detail=detail) + return + else: + heading = HT.Paragraph('Trait Data Entered for %s Set' % fd.RISet) + heading.__setattr__("class","title") + + datadiv = HT.TD(heading, heading2,heading3, width='45%',valign='top', align='left', bgColor='#eeeeee') + resultstable,tblobj,bottomInfo = self.GenReport(ChrNameOrderIdDict,fd, _genotype, _strains, _vals, _vars) + #resultstable = self.GenReport(fd, _genotype, _strains, _vals, _vars) + + # creat object for result table for sort function + objfile = open('%s.obj' % (webqtlConfig.TMPDIR+filename), 'wb') + cPickle.dump(tblobj, objfile) + objfile.close() + + sortby = ("Index", "up") + reportTable =HT.Div(webqtlUtil.genTableObj(tblobj=tblobj, file=filename, sortby=sortby, tableID = "sortable", addIndex = "0"), Id="sortable") + + descriptionTable = HT.TableLite(border=0, cellpadding=0, cellspacing=0) + descriptionTable.append(HT.TR(HT.TD(reportTable, colspan=3))) + descriptionTable.append(HT.TR(HT.TD(HT.BR(),HT.BR()))) + descriptionTable.append(bottomInfo) + + self.traitList=_vals + + ##########################plot####################### + + ################################################################ + # Generate Chr list and Retrieve Length Information + ################################################################ + self.genotype= _genotype + self.ChrList = [("All", -1)] + + for i, indChr in enumerate(self.genotype): + self.ChrList.append((indChr.name, i)) + + self.cursor.execute(""" + Select + Length from Chr_Length, InbredSet + where + Chr_Length.SpeciesId = InbredSet.SpeciesId AND + InbredSet.Name = '%s' AND + Chr_Length.Name in (%s) + Order by + OrderId + """ % (fd.RISet, string.join(map(lambda X: "'%s'" % X[0], self.ChrList[1:]), ", "))) + + self.ChrLengthMbList = self.cursor.fetchall() + self.ChrLengthMbList = map(lambda x: x[0]/1000000.0, self.ChrLengthMbList) + self.ChrLengthMbSum = reduce(lambda x, y:x+y, self.ChrLengthMbList, 0.0) + if self.ChrLengthMbList: + self.MbGraphInterval = self.ChrLengthMbSum/(len(self.ChrLengthMbList)*12) #Empirical Mb interval + else: + self.MbGraphInterval = 1 + + self.ChrLengthCMList = [] + for i, _chr in enumerate(self.genotype): + self.ChrLengthCMList.append(_chr[-1].cM - _chr[0].cM) + self.ChrLengthCMSum = reduce(lambda x, y:x+y, self.ChrLengthCMList, 0.0)# used for calculate plot scale + + self.GraphInterval = self.MbGraphInterval #Mb + + # begin: common part with human data + intCanvas = pid.PILCanvas(size=(self.graphWidth,self.graphHeight)) + gifmap = self.plotIntMapping(fd, intCanvas, startMb = self.startMb, endMb = self.endMb, showLocusForm= "") + filename= webqtlUtil.genRandStr("Itvl_") + intCanvas.save(os.path.join(webqtlConfig.IMGDIR, filename), format='png') + intImg=HT.Image('/image/'+filename+'.png', border=0, usemap='#WebQTLImageMap') + + ################################################################ + # footnote goes here + ################################################################ + btminfo = HT.Paragraph(Id="smallsize") #Small('More information about this graph is available here.') + + if (self.additiveChecked): + btminfo.append(HT.BR(), 'A positive additive coefficient (', HT.Font('green', color='green'), ' line) indicates that %s alleles increase trait values. In contrast, a negative additive coefficient (' % fd.ppolar, HT.Font('red', color='red'), ' line) indicates that %s alleles increase trait values.' % fd.mpolar) + + + TD_LR = HT.TR(HT.TD(HT.Blockquote(gifmap,intImg, HT.P()), bgColor='#eeeeee', height = 200)) + + self.dict['body'] = str(datadiv)+str(TD_LR)+str(resultstable)+str(HT.TR(HT.TD(descriptionTable))) + + # end: common part with human data + + else: + pass + + + # add by NL 10-2-2011 + def initializeParameters(self, fd): + """ + Initializes all of the MarkerRegressionPage class parameters, + acquiring most values from the formdata (fd) + """ + ################################### + # manhattam plot parameters + ################################### + + self.graphHeight = 600 + self.graphWidth = 1280 + self.plotScale = 'physic' + self.selectedChr = -1 + self.GRAPH_BACK_DARK_COLOR = pid.HexColor(0xF1F1F9) + self.GRAPH_BACK_LIGHT_COLOR = pid.HexColor(0xFBFBFF) + self.LRS_COLOR = pid.HexColor(0x0000FF) + self.LRS_LOD ='LRS' + self.lrsMax = float(fd.formdata.getvalue('lrsMax', 0)) + self.startMb = fd.formdata.getvalue('startMb', "-1") + self.endMb = fd.formdata.getvalue('endMb', "-1") + self.mappingMethodId = fd.formdata.getvalue('mappingMethodId', "0") + self.permChecked=True + self.multipleInterval=False + self.SIGNIFICANT_WIDTH = 5 + self.SUGGESTIVE_WIDTH = 5 + self.SIGNIFICANT_COLOR = pid.HexColor(0xEBC7C7) + self.SUGGESTIVE_COLOR = pid.gainsboro + self.colorCollection = [self.LRS_COLOR] + self.additiveChecked= True + self.ADDITIVE_COLOR_POSITIVE = pid.green + self.legendChecked =False + self.pValue=float(fd.formdata.getvalue('pValue',-1)) + + # allow user to input p-value greater than 1, + # in this case, the value will be treated as -lgP value. so the input value needs to be transferred to power of 10 format + if self.pValue >1: + self.pValue =10**-(self.pValue) + + try: + self.startMb = float(self.startMb) + self.endMb = float(self.endMb) + if self.startMb > self.endMb: + temp = self.startMb + self.startMb = self.endMb + self.endMb = temp + #minimal distance 10bp + if self.endMb - self.startMb < 0.00001: + self.endMb = self.startMb + 0.00001 + except: + self.startMb = self.endMb = -1 + + def GenReportForPLINK(self, ChrNameOrderIdDict={},RISet='',plinkResultDict= {},thresholdPvalue=-1,chrList=[]): + + 'Create an HTML division which reports any loci which are significantly associated with the submitted trait data.' + ######################################### + # Genome Association 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':RISet+"Geno",'CellID':'_', \ + 'RISet':RISet, 'incparentsf1':'on'} + for key in hddn.keys(): + locusForm.append(HT.Input(name=key, value=hddn[key], type='hidden')) + + regressionHeading = HT.Paragraph('Genome Association Report') + regressionHeading.__setattr__("class","title") + + filename= webqtlUtil.genRandStr("GenomeAsscociation_") + fpText = open('%s.txt' % (webqtlConfig.TMPDIR+filename), 'wb') + fpText.write('The loci meet the criteria of P-Value <= %3.6f.\n'%thresholdPvalue) + pValueInfo =HT.Paragraph('The loci meet the criteria of P-Value <= %3.6f.\n'%thresholdPvalue) + + textUrl = HT.Href(text = 'Download', url= '/tmp/'+filename+'.txt', target = "_blank", Class='fs12 fwn') + bottomInfo = HT.TR(HT.TD(HT.Paragraph(textUrl, ' result in tab-delimited text format.', HT.BR(), HT.BR(),Class="fs12 fwn"), colspan=3)) + + tblobj={} # build dict for genTableObj function; keys include header and body + tblobj_header = [] # value of key 'header' + tblobj_body=[] # value of key 'body' + reportHeaderRow=[] # header row list for tblobj_header (html part) + headerList=['Index','SNP Name','Chr','Mb','-log(P)'] + headerStyle="fs14 fwb ffl b1 cw cbrb" # style of the header + cellColorStyle = "fs13 b1 fwn c222" # style of the cells + + if headerList: + for ncol, item in enumerate(headerList): + reportHeaderRow.append(THCell(HT.TD(item, Class=headerStyle, valign='bottom',nowrap='ON'),text=item, idx=ncol)) + #download file for table headers' names + fpText.write('SNP_Name\tChromosome\tMb\t-log(P)\n') + + tblobj_header.append(reportHeaderRow) + tblobj['header']=tblobj_header + + index=1 + for chr in chrList: + + if plinkResultDict.has_key(chr): + if chr in ChrNameOrderIdDict.keys(): + chrOrderId =ChrNameOrderIdDict[chr] + else: + chrOrderId=chr + + valueList=plinkResultDict[chr] + + for value in valueList: + reportBodyRow=[] # row list for tblobj_body (html part) + snpName=value[0] + bp=value[1] + mb=int(bp)/1000000.0 + + try: + pValue =float(value[2]) + except: + pValue =1 + formattedPvalue = -math.log10(pValue) + + formattedPvalue = webqtlUtil.SciFloat(formattedPvalue) + dbSnprs=snpName.replace('rs','') + SnpHref = HT.Href(text=snpName, url="http://www.ncbi.nlm.nih.gov/projects/SNP/snp_ref.cgi?rs=%s"%dbSnprs, target="_blank") + + selectCheck=HT.Input(type="checkbox", Class="checkbox", name="index",value=index, onClick="highlight(this)") + reportBodyRow.append(TDCell(HT.TD(str(index),selectCheck, align='right',Class=cellColorStyle,nowrap='ON'),str(index),index)) + reportBodyRow.append(TDCell(HT.TD(SnpHref, Class=cellColorStyle,nowrap='ON'),snpName, snpName)) + reportBodyRow.append(TDCell(HT.TD(chr, Class=cellColorStyle, align="center",nowrap='ON'),chr, chrOrderId)) + reportBodyRow.append(TDCell(HT.TD('%3.6f'%mb, Class=cellColorStyle, align="center",nowrap='ON'),mb, mb)) + reportBodyRow.append(TDCell(HT.TD(formattedPvalue, Class=cellColorStyle, align="center",nowrap='ON'),formattedPvalue, float(formattedPvalue))) + + fpText.write('%s\t%s\t%3.6f\t%s\n' % (snpName, str(chr), mb, formattedPvalue)) + index+=1 + + tblobj_body.append(reportBodyRow) + + tblobj['body']=tblobj_body + rv=HT.TR(HT.TD(regressionHeading,pValueInfo, locusForm, HT.P(), width='55%',valign='top', align='left',bgColor='#eeeeee')) + + return rv, tblobj,bottomInfo + + + def GenReport(self, ChrNameOrderIdDict,fd, _genotype, _strains, _vals, _vars= []): + 'Create an HTML division which reports any loci which are significantly associated with the submitted trait data.' + #calculate QTL for each trait + self.qtlresults = [] + if webqtlUtil.ListNotNull(_vars): + qtlresults = _genotype.regression(strains = _strains, trait = _vals, variance = _vars) + LRSArray = _genotype.permutation(strains = _strains, trait = _vals, variance = _vars, nperm=fd.nperm) + else: + qtlresults = _genotype.regression(strains = _strains, trait = _vals) + LRSArray = _genotype.permutation(strains = _strains, trait = _vals,nperm=fd.nperm) + + self.qtlresults.append(qtlresults) + + filename= webqtlUtil.genRandStr("GenomeAsscociation_") + + # set suggestive, significant and highly significant LRS + 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) + + self.significance =fd.significance + self.suggestive = fd.suggestive + self.highlysignificant = LRSArray[int(fd.nperm*0.99-1)] + _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:] + + ######################################### + # Permutation Graph + ######################################### + 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) + + permutationHeading = HT.Paragraph('Histogram of Permutation Test') + permutationHeading.__setattr__("class","title") + + permutation = HT.TableLite() + permutation.append(HT.TR(HT.TD(img))) + + + ######################################### + # Genome Association 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('Genome Association Report') + regressionHeading.__setattr__("class","title") + # report is the info part above report table + 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') + fpText.write('Suggestive LRS =%3.2f\n'%self.suggestive) + fpText.write('Significant LRS =%3.2f\n'%self.significance) + fpText.write('Highly Significant LRS =%3.2f\n'%self.highlysignificant) + LRSInfo =HT.Paragraph(' Suggestive LRS =%3.2f\n'%fd.suggestive, HT.BR(), ' Significant LRS =%3.2f\n'%fd.significance,HT.BR(),' Highly Significant LRS =%3.2f\n' % self.highlysignificant) + + textUrl = HT.Href(text = 'Download', url= '/tmp/'+filename+'.txt', target = "_blank", Class='fs12 fwn') + + bottomInfo = HT.TR(HT.TD(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"))) + + tblobj={} # build dict for genTableObj function; keys include header and body + tblobj_header = [] # value of key 'header' + tblobj_body=[] # value of key 'body' + reportHeaderRow=[] # header row list for tblobj_header (html part) + headerStyle="fs14 fwb ffl b1 cw cbrb" # style of the header + cellColorStyle = "fs13 b1 fwn c222" # style of the cells + + headerList=['Index','LRS','Chr','Mb','Locus','Additive Effect'] + for ncol, item in enumerate(headerList): + reportHeaderRow.append(THCell(HT.TD(item, Class=headerStyle, valign='bottom',nowrap='ON'),text=item, idx=ncol)) + + if fd.genotype.type == 'intercross': + ncol =len(headerList) + reportHeaderRow.append(THCell(HT.TD('Dominance Effect', Class=headerStyle, valign='bottom',nowrap='ON'),text='Dominance Effect', idx=ncol)) + + #download file for table headers' names + fpText.write('LRS\tChromosome\tMb\tLocus\tAdditive Effect\tDominance Effect\n') + + index=1 + for ii in filtered: + #add by NL 06-20-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 + + if LRS > fd.significance: + lrs = HT.TD(HT.Font('%3.3f*' % LRS, color='#FF0000'),Class=cellColorStyle) + else: + lrs = HT.TD('%3.3f' % LRS,Class=cellColorStyle) + + if ii.locus.chr in ChrNameOrderIdDict.keys(): + chrOrderId =ChrNameOrderIdDict[ii.locus.chr] + else: + chrOrderId=ii.locus.chr + + reportBodyRow=[] # row list for tblobj_body (html part) + selectCheck=HT.Input(type="checkbox", Class="checkbox", name="index",value=index, onClick="highlight(this)") + reportBodyRow.append(TDCell(HT.TD(str(index),selectCheck, align='right',Class=cellColorStyle,nowrap='ON'),str(index),index)) + reportBodyRow.append(TDCell(lrs,LRS, LRS)) + reportBodyRow.append(TDCell(HT.TD(ii.locus.chr, Class=cellColorStyle, align="center",nowrap='ON'),ii.locus.chr, chrOrderId)) + reportBodyRow.append(TDCell(HT.TD('%3.6f'%ii.locus.Mb, Class=cellColorStyle, align="center",nowrap='ON'),ii.locus.Mb, ii.locus.Mb)) + reportBodyRow.append(TDCell(HT.TD(HT.Href(text=ii.locus.name, url = "javascript:showTrait('%s','%s');" % (locusFormName, ii.locus.name), Class='normalsize'), Class=cellColorStyle, align="center",nowrap='ON'),ii.locus.name, ii.locus.name)) + reportBodyRow.append(TDCell(HT.TD('%3.3f' % ii.additive, Class=cellColorStyle, align="center",nowrap='ON'),ii.additive, ii.additive)) + reportBodyRow.append(TDCell(HT.TD('%3.3f' % ii.dominance, Class=cellColorStyle, align="center",nowrap='ON'),ii.dominance, ii.dominance)) + + fpText.write('%2.3f\t%s\t%3.6f\t%s\t%2.3f\t%2.3f\n' % (LRS, ii.locus.chr, ii.locus.Mb, ii.locus.name, ii.additive, ii.dominance)) + index+=1 + tblobj_body.append(reportBodyRow) + else: + #download file for table headers' names + fpText.write('LRS\tChromosome\tMb\tLocus\tAdditive Effect\n') + + index=1 + for ii in filtered: + #add by NL 06-20-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 + + if LRS > fd.significance: + lrs = HT.TD(HT.Font('%3.3f*' % LRS, color='#FF0000'),Class=cellColorStyle) + else: + lrs = HT.TD('%3.3f' % LRS,Class=cellColorStyle) + + if ii.locus.chr in ChrNameOrderIdDict.keys(): + chrOrderId =ChrNameOrderIdDict[ii.locus.chr] + else: + chrOrderId=ii.locus.chr + + reportBodyRow=[] # row list for tblobj_body (html part) + selectCheck=HT.Input(type="checkbox", Class="checkbox", name="index",value=index, onClick="highlight(this)") + reportBodyRow.append(TDCell(HT.TD(str(index),selectCheck, align='right',Class=cellColorStyle,nowrap='ON'),str(index),index)) + reportBodyRow.append(TDCell(lrs,LRS, LRS)) + reportBodyRow.append(TDCell(HT.TD(ii.locus.chr, Class=cellColorStyle, align="center",nowrap='ON'),ii.locus.chr, chrOrderId)) + reportBodyRow.append(TDCell(HT.TD('%3.6f'%ii.locus.Mb, Class=cellColorStyle, align="center",nowrap='ON'),ii.locus.Mb, ii.locus.Mb)) + reportBodyRow.append(TDCell(HT.TD(HT.Href(text=ii.locus.name, url = "javascript:showTrait('%s','%s');" % (locusFormName, ii.locus.name), Class='normalsize'), Class=cellColorStyle, align="center",nowrap='ON'),ii.locus.name, ii.locus.name)) + reportBodyRow.append(TDCell(HT.TD('%3.3f' % ii.additive, Class=cellColorStyle, align="center",nowrap='ON'),ii.additive, ii.additive)) + + fpText.write('%2.3f\t%s\t%3.6f\t%s\t%2.3f\n' % (LRS, ii.locus.chr, ii.locus.Mb, ii.locus.name, ii.additive)) + index+=1 + tblobj_body.append(reportBodyRow) + + tblobj_header.append(reportHeaderRow) + tblobj['header']=tblobj_header + tblobj['body']=tblobj_body + + rv=HT.TD(regressionHeading,LRSInfo,report, locusForm, HT.P(),width='55%',valign='top', align='left', bgColor='#eeeeee') + if fd.genotype.type == 'intercross': + bottomInfo.append(HT.BR(), HT.BR(), HT.Strong('Dominance Effect'),' is the difference between the mean trait value of cases heterozygous at a marker and the average mean for the two groups homozygous at this marker: e.g., BD - (BB+DD)/2]. A positive dominance effect indicates that the average phenotype of BD heterozygotes exceeds the mean of BB and DD homozygotes. No dominance deviation can be computed for a set of recombinant inbred strains or for a backcross.') + return rv,tblobj,bottomInfo + + return rv,tblobj,bottomInfo + + def plotIntMappingForPLINK(self, fd, canvas, offset= (80, 120, 20, 80), zoom = 1, startMb = None, endMb = None, showLocusForm = "",plinkResultDict={}): + #calculating margins + xLeftOffset, xRightOffset, yTopOffset, yBottomOffset = offset + + fontZoom = zoom + if zoom == 2: + fontZoom = 1.5 + + xLeftOffset = int(xLeftOffset*fontZoom) + xRightOffset = int(xRightOffset*fontZoom) + yBottomOffset = int(yBottomOffset*fontZoom) + + cWidth = canvas.size[0] + cHeight = canvas.size[1] + plotWidth = cWidth - xLeftOffset - xRightOffset + plotHeight = cHeight - yTopOffset - yBottomOffset + startPixelX = xLeftOffset + endPixelX = (xLeftOffset + plotWidth) + + #Drawing Area Height + drawAreaHeight = plotHeight + if self.plotScale == 'physic' and self.selectedChr > -1: # for single chr + drawAreaHeight -= self.ENSEMBL_BAND_HEIGHT + self.UCSC_BAND_HEIGHT+ self.WEBQTL_BAND_HEIGHT + 3*self.BAND_SPACING+ 10*zoom + if self.geneChecked: + drawAreaHeight -= self.NUM_GENE_ROWS*self.EACH_GENE_HEIGHT + 3*self.BAND_SPACING + 10*zoom + else: + if self.selectedChr > -1: + drawAreaHeight -= 20 + else:# for all chrs + drawAreaHeight -= 30 + + #Image map + gifmap = HT.Map(name='WebQTLImageMap') + + newoffset = (xLeftOffset, xRightOffset, yTopOffset, yBottomOffset) + # Draw the alternating-color background first and get plotXScale + plotXScale = self.drawGraphBackgroundForPLINK(canvas, gifmap, offset=newoffset, zoom= zoom, startMb=startMb, endMb = endMb,plinkResultDict=plinkResultDict) + + # Draw X axis + self.drawXAxisForPLINK(fd, canvas, drawAreaHeight, gifmap, plotXScale, showLocusForm, offset=newoffset, zoom= zoom, startMb=startMb, endMb = endMb) + # Draw manhattam plot + self.drawManhattanPlotForPLINK(canvas, drawAreaHeight, gifmap, plotXScale, offset=newoffset, zoom= zoom, startMb=startMb, endMb = endMb,plinkResultDict=plinkResultDict,thresholdPvalue=self.pValue) + + return gifmap + + + def plotIntMapping(self, fd, canvas, offset= (80, 120, 20, 80), zoom = 1, startMb = None, endMb = None, showLocusForm = ""): + #calculating margins + xLeftOffset, xRightOffset, yTopOffset, yBottomOffset = offset + + fontZoom = zoom + if zoom == 2: + fontZoom = 1.5 + + xLeftOffset = int(xLeftOffset*fontZoom) + xRightOffset = int(xRightOffset*fontZoom) + yBottomOffset = int(yBottomOffset*fontZoom) + + cWidth = canvas.size[0] + cHeight = canvas.size[1] + plotWidth = cWidth - xLeftOffset - xRightOffset + plotHeight = cHeight - yTopOffset - yBottomOffset + startPixelX = xLeftOffset + endPixelX = (xLeftOffset + plotWidth) + + #Drawing Area Height + drawAreaHeight = plotHeight + if self.plotScale == 'physic' and self.selectedChr > -1: # for single chr + drawAreaHeight -= self.ENSEMBL_BAND_HEIGHT + self.UCSC_BAND_HEIGHT+ self.WEBQTL_BAND_HEIGHT + 3*self.BAND_SPACING+ 10*zoom + if self.geneChecked: + drawAreaHeight -= self.NUM_GENE_ROWS*self.EACH_GENE_HEIGHT + 3*self.BAND_SPACING + 10*zoom + else:# for all chrs + if self.selectedChr > -1: + drawAreaHeight -= 20 + else: + drawAreaHeight -= 30 + + #Image map + gifmap = HT.Map(name='WebQTLImageMap') + + newoffset = (xLeftOffset, xRightOffset, yTopOffset, yBottomOffset) + # Draw the alternating-color background first and get plotXScale + plotXScale = self.drawGraphBackground(canvas, gifmap, offset=newoffset, zoom= zoom, startMb=startMb, endMb = endMb) + + # Draw X axis + self.drawXAxis(fd, canvas, drawAreaHeight, gifmap, plotXScale, showLocusForm, offset=newoffset, zoom= zoom, startMb=startMb, endMb = endMb) + # Draw QTL curve + self.drawQTL(canvas, drawAreaHeight, gifmap, plotXScale, offset=newoffset, zoom= zoom, startMb=startMb, endMb = endMb) + + #draw legend + if self.multipleInterval: + self.drawMultiTraitName(fd, canvas, gifmap, showLocusForm, offset=newoffset) + elif self.legendChecked: + self.drawLegendPanel(fd, canvas, offset=newoffset) + else: + pass + + #draw position, no need to use a separate function + if fd.genotype.Mbmap: + self.drawProbeSetPosition(canvas, plotXScale, offset=newoffset) + + return gifmap + + + # functions for manhattam plot of markers + def drawManhattanPlotForPLINK(self, canvas, drawAreaHeight, gifmap, plotXScale, offset= (40, 120, 80, 10), zoom = 1, startMb = None, endMb = None,plinkResultDict={},thresholdPvalue=-1): + + xLeftOffset, xRightOffset, yTopOffset, yBottomOffset = offset + plotWidth = canvas.size[0] - xLeftOffset - xRightOffset + plotHeight = canvas.size[1] - yTopOffset - yBottomOffset + fontZoom = zoom + if zoom == 2: + fontZoom = 1.5 + + # INTERCROSS = (self.genotype.type=="intercross") + INTERCROSS ='' #?????? + + ChrLengthDistList = self.ChrLengthMbList + drawRegionDistance = self.ChrLengthMbSum + GraphInterval=self.GraphInterval + pvalueHeightThresh = drawAreaHeight - 80 #ZS: Otherwise the plot gets very close to the chromosome labels + + #draw the pvalue scale + #We first determine whether or not we are using a sliding scale. + #If so, we need to compute the maximum pvalue value to determine where the max y-value should be, and call this pvalueMax. + #pvalueTop is then defined to be above the pvalueMax by enough to add one additional pvalueScale increment. + #if we are using a set-scale, then we set pvalueTop to be the user's value, and pvalueMax doesn't matter. + + # for human data we use p value instead of lrs + pValueList=[] + for key in plinkResultDict: + valueList = plinkResultDict[key] + for item in valueList: + pValue = item[-1] + pValueList.append(pValue) + + formattedPValueList=[] + for pValue in pValueList: + try: + pValue=float(pValue) + except: + pValue =1 + formattedpValue = -math.log10(pValue) + formattedPValueList.append(formattedpValue) + + #sliding scale + pvalueMax = max(formattedPValueList) + #pvalueMax =pvalueMax +1 + # no permutation result for plink func: GenReport() + pvalueMin = int(-math.log10(thresholdPvalue)) + + if pvalueMax> 100: + pvalueScale = 20.0 + elif pvalueMax > 20: + pvalueScale = 5.0 + elif pvalueMax > 7.5: + pvalueScale = 2.5 + else: + pvalueScale = 1.0 + + # the base line for x-axis is -log(thresholdPvalue) + pvalueAxisList = Plot.frange(pvalueMin, pvalueMax, pvalueScale) + #make sure the user's value appears on the y-axis + #ZS: There is no way to do this without making the position of the points not directly proportional to a given distance on the y-axis + #tempPvalueMax=round(pvalueMax) + tempPvalueMax = pvalueAxisList[len(pvalueAxisList)-1] + pvalueScale + pvalueAxisList.append(tempPvalueMax) + + #ZS: I don't understand this; the if statement will be true for any number that isn't exactly X.5. + #if abs(tempPvalueMax-pvalueMax) <0.5: + # tempPvalueMax=tempPvalueMax+1 + # pvalueAxisList.append(tempPvalueMax) + + #draw the "pvalue" string to the left of the axis + pvalueScaleFont=pid.Font(ttf="verdana", size=14*fontZoom, bold=0) + pvalueLODFont=pid.Font(ttf="verdana", size=14*zoom*1.5, bold=0) + yZero = yTopOffset + plotHeight + + #yAxis label display area + yAxis_label ='-log(P)' + canvas.drawString(yAxis_label, xLeftOffset - canvas.stringWidth("999.99", font=pvalueScaleFont) - 10*zoom, \ + yZero - 150, font=pvalueLODFont, color=pid.black, angle=90) + + for i,item in enumerate(pvalueAxisList): + ypvalue = yZero - (float(i)/float(len(pvalueAxisList) - 1)) * pvalueHeightThresh + canvas.drawLine(xLeftOffset, ypvalue, xLeftOffset - 4, ypvalue, color=self.LRS_COLOR, width=1*zoom) + scaleStr = "%2.1f" % item + #added by NL 6-24-2011:Y-axis scale display + canvas.drawString(scaleStr, xLeftOffset-4-canvas.stringWidth(scaleStr, font=pvalueScaleFont)-5, ypvalue+3, font=pvalueScaleFont, color=self.LRS_COLOR) + + ChrList=self.ChrList + startPosX = xLeftOffset + + for i, chr in enumerate(ChrList): + + if plinkResultDict.has_key(chr): + plinkresultList = plinkResultDict[chr] + + m = 0 + #add by NL 06-24-2011: for mahanttam plot + symbolFont = pid.Font(ttf="fnt_bs", size=5,bold=0) + # color for point in each chr + chrCount=len(ChrList) + chrColorDict =self.getColorForMarker(chrCount=chrCount,flag=1) + for j, item in enumerate(plinkresultList): + try : + mb=float(item[1])/1000000.0 + except: + mb=0 + + try : + pvalue =float(item[-1]) + except: + pvalue =1 + + try: + snpName = item[0] + except: + snpName='' + + formattedPvalue = -math.log10(pvalue) + + Xc = startPosX + (mb-startMb)*plotXScale + Yc = yZero - (formattedPvalue-pvalueMin)*pvalueHeightThresh/(tempPvalueMax - pvalueMin) + canvas.drawString("5", Xc-canvas.stringWidth("5",font=symbolFont)/2+1,Yc+2,color=chrColorDict[i], font=symbolFont) + m += 1 + + startPosX += (ChrLengthDistList[i]+GraphInterval)*plotXScale + + canvas.drawLine(xLeftOffset, yZero, xLeftOffset, yTopOffset, color=self.LRS_COLOR, width=1*zoom) #the blue line running up the y axis + + def drawQTL(self, canvas, drawAreaHeight, gifmap, plotXScale, offset= (40, 120, 80, 10), zoom = 1, startMb = None, endMb = None): + + xLeftOffset, xRightOffset, yTopOffset, yBottomOffset = offset + plotWidth = canvas.size[0] - xLeftOffset - xRightOffset + plotHeight = canvas.size[1] - yTopOffset - yBottomOffset + fontZoom = zoom + if zoom == 2: + fontZoom = 1.5 + + INTERCROSS = (self.genotype.type=="intercross") + + ChrLengthDistList = self.ChrLengthMbList + GraphInterval=self.GraphInterval + LRSHeightThresh = drawAreaHeight + AdditiveHeightThresh = drawAreaHeight/2 + DominanceHeightThresh = drawAreaHeight/2 + + #draw the LRS scale + #We first determine whether or not we are using a sliding scale. + #If so, we need to compute the maximum LRS value to determine where the max y-value should be, and call this LRSMax. + #LRSTop is then defined to be above the LRSMax by enough to add one additional LRSScale increment. + #if we are using a set-scale, then we set LRSTop to be the user's value, and LRSMax doesn't matter. + + if self.LRS_LOD == 'LOD': + lodm = self.LODFACTOR + else: + lodm = 1.0 + + if self.lrsMax <= 0: #sliding scale + LRSMax = max(map(max, self.qtlresults)).lrs + #genotype trait will give infinite LRS + LRSMax = min(LRSMax, webqtlConfig.MAXLRS) + LRSMax = max(self.significance, LRSMax) + else: + LRSMax = self.lrsMax*lodm + + if LRSMax/lodm > 100: + LRSScale = 20.0 + elif LRSMax/lodm > 20: + LRSScale = 5.0 + elif LRSMax/lodm > 7.5: + LRSScale = 2.5 + else: + LRSScale = 1.0 + + LRSAxisList = Plot.frange(LRSScale, LRSMax/lodm, LRSScale) + #make sure the user's value appears on the y-axis + #update by NL 6-21-2011: round the LOD value to 100 when LRSMax is equal to 460 + LRSAxisList.append(round(LRSMax/lodm)) + + #draw the "LRS" or "LOD" string to the left of the axis + LRSScaleFont=pid.Font(ttf="verdana", size=14*fontZoom, bold=0) + LRSLODFont=pid.Font(ttf="verdana", size=14*zoom*1.5, bold=0) + yZero = yTopOffset + plotHeight + + #yAxis label display area + canvas.drawString(self.LRS_LOD, xLeftOffset - canvas.stringWidth("999.99", font=LRSScaleFont) - 10*zoom, \ + yZero - 150, font=LRSLODFont, color=pid.black, angle=90) + + for item in LRSAxisList: + yLRS = yZero - (item*lodm/LRSMax) * LRSHeightThresh + canvas.drawLine(xLeftOffset, yLRS, xLeftOffset - 4, yLRS, color=self.LRS_COLOR, width=1*zoom) + scaleStr = "%2.1f" % item + #added by NL 6-24-2011:Y-axis scale display + canvas.drawString(scaleStr, xLeftOffset-4-canvas.stringWidth(scaleStr, font=LRSScaleFont)-5, yLRS+3, font=LRSScaleFont, color=self.LRS_COLOR) + + + #"Significant" and "Suggestive" Drawing Routine + # ======= Draw the thick lines for "Significant" and "Suggestive" ===== (crowell: I tried to make the SNPs draw over these lines, but piddle wouldn't have it...) + if self.permChecked and not self.multipleInterval: + significantY = yZero - self.significance*LRSHeightThresh/LRSMax + suggestiveY = yZero - self.suggestive*LRSHeightThresh/LRSMax + + + startPosX = xLeftOffset + for i, _chr in enumerate(self.genotype): + rightEdge = int(startPosX + self.ChrLengthDistList[i]*plotXScale - self.SUGGESTIVE_WIDTH/1.5) + #added by NL 6-24-2011:draw suggestive line (grey one) + canvas.drawLine(startPosX+self.SUGGESTIVE_WIDTH/1.5, suggestiveY, rightEdge, suggestiveY, color=self.SUGGESTIVE_COLOR, + width=self.SUGGESTIVE_WIDTH*zoom, clipX=(xLeftOffset, xLeftOffset + plotWidth-2)) + #added by NL 6-24-2011:draw significant line (pink one) + canvas.drawLine(startPosX+self.SUGGESTIVE_WIDTH/1.5, significantY, rightEdge, significantY, color=self.SIGNIFICANT_COLOR, + width=self.SIGNIFICANT_WIDTH*zoom, clipX=(xLeftOffset, xLeftOffset + plotWidth-2)) + sugg_coords = "%d, %d, %d, %d" % (startPosX, suggestiveY-2, rightEdge + 2*zoom, suggestiveY+2) + sig_coords = "%d, %d, %d, %d" % (startPosX, significantY-2, rightEdge + 2*zoom, significantY+2) + if self.LRS_LOD == 'LRS': + sugg_title = "Suggestive LRS = %0.2f" % self.suggestive + sig_title = "Significant LRS = %0.2f" % self.significance + else: + sugg_title = "Suggestive LOD = %0.2f" % (self.suggestive/4.61) + sig_title = "Significant LOD = %0.2f" % (self.significance/4.61) + Areas1 = HT.Area(shape='rect',coords=sugg_coords,title=sugg_title) + Areas2 = HT.Area(shape='rect',coords=sig_coords,title=sig_title) + gifmap.areas.append(Areas1) + gifmap.areas.append(Areas2) + + startPosX += (self.ChrLengthDistList[i]+self.GraphInterval)*plotXScale + + + if self.multipleInterval: + lrsEdgeWidth = 1 + else: + additiveMax = max(map(lambda X : abs(X.additive), self.qtlresults[0])) + if INTERCROSS: + dominanceMax = max(map(lambda X : abs(X.dominance), self.qtlresults[0])) + else: + dominanceMax = -1 + lrsEdgeWidth = 2 + for i, qtlresult in enumerate(self.qtlresults): + m = 0 + startPosX = xLeftOffset + thisLRSColor = self.colorCollection[i] + + #add by NL 06-24-2011: for mahanttam plot + symbolFont = pid.Font(ttf="fnt_bs", size=5,bold=0) + + for j, _chr in enumerate(self.genotype): + chrCount=len(self.genotype) + chrColorDict =self.getColorForMarker(chrCount=chrCount,flag=1) + LRSCoordXY = [] + AdditiveCoordXY = [] + DominanceCoordXY = [] + for k, _locus in enumerate(_chr): + if self.plotScale == 'physic': + Xc = startPosX + (_locus.Mb-startMb)*plotXScale + else: + Xc = startPosX + (_locus.cM-_chr[0].cM)*plotXScale + # updated by NL 06-18-2011: + # fix the over limit LRS graph issue since genotype trait may give infinite LRS; + # for any lrs is over than 460(LRS max in this system), it will be reset to 460 + if qtlresult[m].lrs> 460 or qtlresult[m].lrs=='inf': + Yc = yZero - webqtlConfig.MAXLRS*LRSHeightThresh/LRSMax + else: + Yc = yZero - qtlresult[m].lrs*LRSHeightThresh/LRSMax + + LRSCoordXY.append((Xc, Yc)) + #add by NL 06-24-2011: for mahanttam plot + #self.significance/4.61 consider chr and LOD + # significantY = yZero - self.significance*LRSHeightThresh/LRSMax + # if Yc >significantY: + # canvas.drawString(":", Xc-canvas.stringWidth(":",font=symbolFont)/2+1,Yc+2,color=pid.black, font=symbolFont) + # else: + # canvas.drawString(":", Xc-canvas.stringWidth(":",font=symbolFont)/2+1,Yc+2,color=pid.black, font=symbolFont) + + # add by NL 06-27-2011: eliminate imputed value when locus name is equal to '-' + if (qtlresult[m].locus.name) and (qtlresult[m].locus.name!=' - '): + canvas.drawString("5", Xc-canvas.stringWidth("5",font=symbolFont)/2+1,Yc+2,color=chrColorDict[j], font=symbolFont) + + if not self.multipleInterval and self.additiveChecked: + Yc = yZero - qtlresult[m].additive*AdditiveHeightThresh/additiveMax + AdditiveCoordXY.append((Xc, Yc)) + if not self.multipleInterval and INTERCROSS and self.additiveChecked: + Yc = yZero - qtlresult[m].dominance*DominanceHeightThresh/dominanceMax + DominanceCoordXY.append((Xc, Yc)) + m += 1 + + startPosX += (ChrLengthDistList[j]+GraphInterval)*plotXScale + + + ###draw additive scale + if not self.multipleInterval and self.additiveChecked: + additiveScaleFont=pid.Font(ttf="verdana",size=12*fontZoom,bold=0) + additiveScale = Plot.detScaleOld(0,additiveMax) + additiveStep = (additiveScale[1]-additiveScale[0])/additiveScale[2] + additiveAxisList = Plot.frange(0, additiveScale[1], additiveStep) + maxAdd = additiveScale[1] + addPlotScale = AdditiveHeightThresh/additiveMax + + additiveAxisList.append(additiveScale[1]) + for item in additiveAxisList: + additiveY = yZero - item*addPlotScale + canvas.drawLine(xLeftOffset + plotWidth,additiveY,xLeftOffset+4+ plotWidth,additiveY,color=self.ADDITIVE_COLOR_POSITIVE, width=1*zoom) + scaleStr = "%2.3f" % item + canvas.drawString(scaleStr,xLeftOffset + plotWidth +6,additiveY+5,font=additiveScaleFont,color=self.ADDITIVE_COLOR_POSITIVE) + + canvas.drawLine(xLeftOffset+plotWidth,additiveY,xLeftOffset+plotWidth,yZero,color=self.ADDITIVE_COLOR_POSITIVE, width=1*zoom) + + canvas.drawLine(xLeftOffset, yZero, xLeftOffset, yTopOffset, color=self.LRS_COLOR, width=1*zoom) #the blue line running up the y axis + + def drawGraphBackgroundForPLINK(self, canvas, gifmap, offset= (80, 120, 80, 50), zoom = 1, startMb = None, endMb = None,plinkResultDict={} ): + + xLeftOffset, xRightOffset, yTopOffset, yBottomOffset = offset + plotWidth = canvas.size[0] - xLeftOffset - xRightOffset + plotHeight = canvas.size[1] - yTopOffset - yBottomOffset + fontZoom = zoom + if zoom == 2: + fontZoom = 1.5 + + #calculate plot scale + #XZ: all of these global variables should be passed from function signiture + ChrLengthDistList = self.ChrLengthMbList + drawRegionDistance = self.ChrLengthMbSum + GraphInterval=self.GraphInterval + ChrList =self.ChrList + + #multiple chromosome view + plotXScale = plotWidth / ((len(ChrList)-1)*GraphInterval + drawRegionDistance) + + startPosX = xLeftOffset + chrLabelFont=pid.Font(ttf="verdana",size=24*fontZoom,bold=0) + + for i, _chr in enumerate(ChrList): + + if (i % 2 == 0): + theBackColor = self.GRAPH_BACK_DARK_COLOR + else: + theBackColor = self.GRAPH_BACK_LIGHT_COLOR + # NL:resize chr width for drawing + if float(ChrLengthDistList[i])<90: + ChrLengthDistList[i]=90 + #draw the shaded boxes and the sig/sug thick lines + canvas.drawRect(startPosX, yTopOffset, startPosX + ChrLengthDistList[i]*plotXScale, \ + yTopOffset+plotHeight, edgeColor=pid.gainsboro,fillColor=theBackColor) + + chrNameWidth = canvas.stringWidth(_chr, font=chrLabelFont) + chrStartPix = startPosX + (ChrLengthDistList[i]*plotXScale -chrNameWidth)/2 + chrEndPix = startPosX + (ChrLengthDistList[i]*plotXScale +chrNameWidth)/2 + + canvas.drawString(_chr, chrStartPix, yTopOffset +20,font = chrLabelFont,color=pid.dimgray) + COORDS = "%d,%d,%d,%d" %(chrStartPix, yTopOffset, chrEndPix,yTopOffset +20) + + #add by NL 09-03-2010 + HREF = "javascript:changeView(%d,%s);" % (i,ChrLengthDistList) + Areas = HT.Area(shape='rect',coords=COORDS,href=HREF) + gifmap.areas.append(Areas) + startPosX += (ChrLengthDistList[i]+GraphInterval)*plotXScale + + return plotXScale + + + def drawGraphBackground(self, canvas, gifmap, offset= (80, 120, 80, 50), zoom = 1, startMb = None, endMb = None): + ##conditions + ##multiple Chromosome view + ##single Chromosome Physical + ##single Chromosome Genetic + xLeftOffset, xRightOffset, yTopOffset, yBottomOffset = offset + plotWidth = canvas.size[0] - xLeftOffset - xRightOffset + plotHeight = canvas.size[1] - yTopOffset - yBottomOffset + fontZoom = zoom + if zoom == 2: + fontZoom = 1.5 + + #calculate plot scale + if self.plotScale != 'physic': + self.ChrLengthDistList = self.ChrLengthCMList + drawRegionDistance = self.ChrLengthCMSum + else: + self.ChrLengthDistList = self.ChrLengthMbList + drawRegionDistance = self.ChrLengthMbSum + + if self.selectedChr > -1: #single chromosome view + spacingAmt = plotWidth/13.5 + i = 0 + for startPix in Plot.frange(xLeftOffset, xLeftOffset+plotWidth, spacingAmt): + if (i % 2 == 0): + theBackColor = self.GRAPH_BACK_DARK_COLOR + else: + theBackColor = self.GRAPH_BACK_LIGHT_COLOR + i += 1 + canvas.drawRect(startPix, yTopOffset, min(startPix+spacingAmt, xLeftOffset+plotWidth), \ + yTopOffset+plotHeight, edgeColor=theBackColor, fillColor=theBackColor) + + drawRegionDistance = self.ChrLengthDistList[self.selectedChr] + self.ChrLengthDistList = [drawRegionDistance] + if self.plotScale == 'physic': + plotXScale = plotWidth / (endMb-startMb) + else: + plotXScale = plotWidth / drawRegionDistance + + else: #multiple chromosome view + plotXScale = plotWidth / ((len(self.genotype)-1)*self.GraphInterval + drawRegionDistance) + + startPosX = xLeftOffset + chrLabelFont=pid.Font(ttf="verdana",size=24*fontZoom,bold=0) + + for i, _chr in enumerate(self.genotype): + + if (i % 2 == 0): + theBackColor = self.GRAPH_BACK_DARK_COLOR + else: + theBackColor = self.GRAPH_BACK_LIGHT_COLOR + + #draw the shaded boxes and the sig/sug thick lines + canvas.drawRect(startPosX, yTopOffset, startPosX + self.ChrLengthDistList[i]*plotXScale, \ + yTopOffset+plotHeight, edgeColor=pid.gainsboro,fillColor=theBackColor) + + chrNameWidth = canvas.stringWidth(_chr.name, font=chrLabelFont) + chrStartPix = startPosX + (self.ChrLengthDistList[i]*plotXScale -chrNameWidth)/2 + chrEndPix = startPosX + (self.ChrLengthDistList[i]*plotXScale +chrNameWidth)/2 + + canvas.drawString(_chr.name, chrStartPix, yTopOffset +20,font = chrLabelFont,color=pid.dimgray) + COORDS = "%d,%d,%d,%d" %(chrStartPix, yTopOffset, chrEndPix,yTopOffset +20) + + #add by NL 09-03-2010 + HREF = "javascript:changeView(%d,%s);" % (i,self.ChrLengthMbList) + Areas = HT.Area(shape='rect',coords=COORDS,href=HREF) + gifmap.areas.append(Areas) + startPosX += (self.ChrLengthDistList[i]+self.GraphInterval)*plotXScale + + return plotXScale + + # XZ: The only difference of function drawXAxisForPLINK and function drawXAxis are the function name and the self.plotScale condition. + def drawXAxisForPLINK(self, fd, canvas, drawAreaHeight, gifmap, plotXScale, showLocusForm, offset= (40, 120, 80, 10), zoom = 1, startMb = None, endMb = None): + xLeftOffset, xRightOffset, yTopOffset, yBottomOffset = offset + plotWidth = canvas.size[0] - xLeftOffset - xRightOffset + plotHeight = canvas.size[1] - yTopOffset - yBottomOffset + yZero = canvas.size[1] - yBottomOffset + fontZoom = zoom + if zoom == 2: + fontZoom = 1.5 + + #Parameters + ChrLengthDistList = self.ChrLengthMbList + GraphInterval=self.GraphInterval + + NUM_MINOR_TICKS = 5 # Number of minor ticks between major ticks + X_MAJOR_TICK_THICKNESS = 2 + X_MINOR_TICK_THICKNESS = 1 + X_AXIS_THICKNESS = 1*zoom + + # ======= Alex: Draw the X-axis labels (megabase location) + MBLabelFont = pid.Font(ttf="verdana", size=12*fontZoom, bold=0) + xMajorTickHeight = 15 # How high the tick extends below the axis + xMinorTickHeight = 5*zoom + xAxisTickMarkColor = pid.black + xAxisLabelColor = pid.black + fontHeight = 12*fontZoom # How tall the font that we're using is + spacingFromLabelToAxis = 20 + spacingFromLineToLabel = 3 + + if self.plotScale == 'physic': + strYLoc = yZero + spacingFromLabelToAxis + canvas.fontHeight(MBLabelFont) + ###Physical single chromosome view + if self.selectedChr > -1: + graphMbWidth = endMb - startMb + XScale = Plot.detScale(startMb, endMb) + XStart, XEnd, XStep = XScale + if XStep < 8: + XStep *= 2 + spacingAmtX = spacingAmt = (XEnd-XStart)/XStep + + j = 0 + while abs(spacingAmtX -int(spacingAmtX)) >= spacingAmtX/100.0 and j < 6: + j += 1 + spacingAmtX *= 10 + + formatStr = '%%2.%df' % j + + for counter, _Mb in enumerate(Plot.frange(XStart, XEnd, spacingAmt / NUM_MINOR_TICKS)): + if _Mb < startMb or _Mb > endMb: + continue + Xc = xLeftOffset + plotXScale*(_Mb - startMb) + if counter % NUM_MINOR_TICKS == 0: # Draw a MAJOR mark, not just a minor tick mark + canvas.drawLine(Xc, yZero, Xc, yZero+xMajorTickHeight, color=xAxisTickMarkColor, width=X_MAJOR_TICK_THICKNESS) # Draw the MAJOR tick mark + labelStr = str(formatStr % _Mb) # What Mbase location to put on the label + strWidth = canvas.stringWidth(labelStr, font=MBLabelFont) + drawStringXc = (Xc - (strWidth / 2.0)) + canvas.drawString(labelStr, drawStringXc, strYLoc, font=MBLabelFont, color=xAxisLabelColor, angle=0) + else: + canvas.drawLine(Xc, yZero, Xc, yZero+xMinorTickHeight, color=xAxisTickMarkColor, width=X_MINOR_TICK_THICKNESS) # Draw the MINOR tick mark + # end else + + ###Physical genome wide view + else: + distScale = 0 + startPosX = xLeftOffset + for i, distLen in enumerate(ChrLengthDistList): + if distScale == 0: #universal scale in whole genome mapping + if distLen > 75: + distScale = 25 + elif distLen > 30: + distScale = 10 + else: + distScale = 5 + for tickdists in range(distScale, ceil(distLen), distScale): + canvas.drawLine(startPosX + tickdists*plotXScale, yZero, startPosX + tickdists*plotXScale, yZero + 7, color=pid.black, width=1*zoom) + canvas.drawString(str(tickdists), startPosX+tickdists*plotXScale, yZero + 10*zoom, color=pid.black, font=MBLabelFont, angle=270) + startPosX += (ChrLengthDistList[i]+GraphInterval)*plotXScale + + megabaseLabelFont = pid.Font(ttf="verdana", size=14*zoom*1.5, bold=0) + canvas.drawString("Megabases", xLeftOffset + (plotWidth -canvas.stringWidth("Megabases", font=megabaseLabelFont))/2, + strYLoc + canvas.fontHeight(MBLabelFont) + 5*zoom, font=megabaseLabelFont, color=pid.black) + pass + + canvas.drawLine(xLeftOffset, yZero, xLeftOffset+plotWidth, yZero, color=pid.black, width=X_AXIS_THICKNESS) # Draw the X axis itself + + def drawXAxis(self, fd, canvas, drawAreaHeight, gifmap, plotXScale, showLocusForm, offset= (40, 120, 80, 10), zoom = 1, startMb = None, endMb = None): + xLeftOffset, xRightOffset, yTopOffset, yBottomOffset = offset + plotWidth = canvas.size[0] - xLeftOffset - xRightOffset + plotHeight = canvas.size[1] - yTopOffset - yBottomOffset + yZero = canvas.size[1] - yBottomOffset + fontZoom = zoom + if zoom == 2: + fontZoom = 1.5 + + #Parameters + NUM_MINOR_TICKS = 5 # Number of minor ticks between major ticks + X_MAJOR_TICK_THICKNESS = 2 + X_MINOR_TICK_THICKNESS = 1 + X_AXIS_THICKNESS = 1*zoom + + # ======= Alex: Draw the X-axis labels (megabase location) + MBLabelFont = pid.Font(ttf="verdana", size=12*fontZoom, bold=0) + xMajorTickHeight = 15 # How high the tick extends below the axis + xMinorTickHeight = 5*zoom + xAxisTickMarkColor = pid.black + xAxisLabelColor = pid.black + fontHeight = 12*fontZoom # How tall the font that we're using is + spacingFromLabelToAxis = 20 + spacingFromLineToLabel = 3 + + if self.plotScale == 'physic': + strYLoc = yZero + spacingFromLabelToAxis + canvas.fontHeight(MBLabelFont) + ###Physical single chromosome view + if self.selectedChr > -1: + graphMbWidth = endMb - startMb + XScale = Plot.detScale(startMb, endMb) + XStart, XEnd, XStep = XScale + if XStep < 8: + XStep *= 2 + spacingAmtX = spacingAmt = (XEnd-XStart)/XStep + + j = 0 + while abs(spacingAmtX -int(spacingAmtX)) >= spacingAmtX/100.0 and j < 6: + j += 1 + spacingAmtX *= 10 + + formatStr = '%%2.%df' % j + + for counter, _Mb in enumerate(Plot.frange(XStart, XEnd, spacingAmt / NUM_MINOR_TICKS)): + if _Mb < startMb or _Mb > endMb: + continue + Xc = xLeftOffset + plotXScale*(_Mb - startMb) + if counter % NUM_MINOR_TICKS == 0: # Draw a MAJOR mark, not just a minor tick mark + canvas.drawLine(Xc, yZero, Xc, yZero+xMajorTickHeight, color=xAxisTickMarkColor, width=X_MAJOR_TICK_THICKNESS) # Draw the MAJOR tick mark + labelStr = str(formatStr % _Mb) # What Mbase location to put on the label + strWidth = canvas.stringWidth(labelStr, font=MBLabelFont) + drawStringXc = (Xc - (strWidth / 2.0)) + canvas.drawString(labelStr, drawStringXc, strYLoc, font=MBLabelFont, color=xAxisLabelColor, angle=0) + else: + canvas.drawLine(Xc, yZero, Xc, yZero+xMinorTickHeight, color=xAxisTickMarkColor, width=X_MINOR_TICK_THICKNESS) # Draw the MINOR tick mark + # end else + + ###Physical genome wide view + else: + distScale = 0 + startPosX = xLeftOffset + for i, distLen in enumerate(self.ChrLengthDistList): + if distScale == 0: #universal scale in whole genome mapping + if distLen > 75: + distScale = 25 + elif distLen > 30: + distScale = 10 + else: + distScale = 5 + for tickdists in range(distScale, ceil(distLen), distScale): + canvas.drawLine(startPosX + tickdists*plotXScale, yZero, startPosX + tickdists*plotXScale, yZero + 7, color=pid.black, width=1*zoom) + canvas.drawString(str(tickdists), startPosX+tickdists*plotXScale, yZero + 10*zoom, color=pid.black, font=MBLabelFont, angle=270) + startPosX += (self.ChrLengthDistList[i]+self.GraphInterval)*plotXScale + + megabaseLabelFont = pid.Font(ttf="verdana", size=14*zoom*1.5, bold=0) + canvas.drawString("Megabases", xLeftOffset + (plotWidth -canvas.stringWidth("Megabases", font=megabaseLabelFont))/2, + strYLoc + canvas.fontHeight(MBLabelFont) + 5*zoom, font=megabaseLabelFont, color=pid.black) + pass + else: + ChrAInfo = [] + preLpos = -1 + distinctCount = 0.0 + if len(self.genotype) > 1: + for i, _chr in enumerate(self.genotype): + thisChr = [] + Locus0CM = _chr[0].cM + nLoci = len(_chr) + if nLoci <= 8: + for _locus in _chr: + if _locus.name != ' - ': + if _locus.cM != preLpos: + distinctCount += 1 + preLpos = _locus.cM + thisChr.append([_locus.name, _locus.cM-Locus0CM]) + else: + for j in (0, nLoci/4, nLoci/2, nLoci*3/4, -1): + while _chr[j].name == ' - ': + j += 1 + if _chr[j].cM != preLpos: + distinctCount += 1 + preLpos = _chr[j].cM + thisChr.append([_chr[j].name, _chr[j].cM-Locus0CM]) + ChrAInfo.append(thisChr) + else: + for i, _chr in enumerate(self.genotype): + thisChr = [] + Locus0CM = _chr[0].cM + for _locus in _chr: + if _locus.name != ' - ': + if _locus.cM != preLpos: + distinctCount += 1 + preLpos = _locus.cM + thisChr.append([_locus.name, _locus.cM-Locus0CM]) + ChrAInfo.append(thisChr) + + stepA = (plotWidth+0.0)/distinctCount + + LRectWidth = 10 + LRectHeight = 3 + offsetA = -stepA + lineColor = pid.lightblue + startPosX = xLeftOffset + for j, ChrInfo in enumerate(ChrAInfo): + preLpos = -1 + for i, item in enumerate(ChrInfo): + Lname,Lpos = item + if Lpos != preLpos: + offsetA += stepA + differ = 1 + else: + differ = 0 + preLpos = Lpos + Lpos *= plotXScale + if self.selectedChr > -1: + Zorder = i % 5 + else: + Zorder = 0 + if differ: + canvas.drawLine(startPosX+Lpos,yZero,xLeftOffset+offsetA,\ + yZero+25, color=lineColor) + canvas.drawLine(xLeftOffset+offsetA,yZero+25,xLeftOffset+offsetA,\ + yZero+40+Zorder*(LRectWidth+3),color=lineColor) + rectColor = pid.orange + else: + canvas.drawLine(xLeftOffset+offsetA, yZero+40+Zorder*(LRectWidth+3)-3,\ + xLeftOffset+offsetA, yZero+40+Zorder*(LRectWidth+3),color=lineColor) + rectColor = pid.deeppink + canvas.drawRect(xLeftOffset+offsetA, yZero+40+Zorder*(LRectWidth+3),\ + xLeftOffset+offsetA-LRectHeight,yZero+40+Zorder*(LRectWidth+3)+LRectWidth,\ + edgeColor=rectColor,fillColor=rectColor,edgeWidth = 0) + COORDS="%d,%d,%d,%d"%(xLeftOffset+offsetA-LRectHeight, yZero+40+Zorder*(LRectWidth+3),\ + xLeftOffset+offsetA,yZero+40+Zorder*(LRectWidth+3)+LRectWidth) + HREF="javascript:showDatabase3('%s','%s','%s','');" % (showLocusForm,fd.RISet+"Geno", Lname) + Areas=HT.Area(shape='rect',coords=COORDS,href=HREF, title="Locus : " + Lname) + gifmap.areas.append(Areas) + ##piddle bug + if j == 0: + canvas.drawLine(startPosX,yZero,startPosX,yZero+40, color=lineColor) + startPosX += (self.ChrLengthDistList[j]+self.GraphInterval)*plotXScale + + canvas.drawLine(xLeftOffset, yZero, xLeftOffset+plotWidth, yZero, color=pid.black, width=X_AXIS_THICKNESS) # Draw the X axis itself + + def getColorForMarker(self, chrCount,flag):# no change is needed + chrColorDict={} + for i in range(chrCount): + if flag==1: # display blue and lightblue intercross + chrColorDict[i]=pid.black + elif flag==0: + if (i%2==0): + chrColorDict[i]=pid.blue + else: + chrColorDict[i]=pid.lightblue + else:#display different color for different chr + if i in [0,8,16]: + chrColorDict[i]=pid.black + elif i in [1,9,17]: + chrColorDict[i]=pid.red + elif i in [2,10,18]: + chrColorDict[i]=pid.lightgreen + elif i in [3,11,19]: + chrColorDict[i]=pid.blue + elif i in [4,12]: + chrColorDict[i]=pid.lightblue + elif i in [5,13]: + chrColorDict[i]=pid.hotpink + elif i in [6,14]: + chrColorDict[i]=pid.gold + elif i in [7,15]: + chrColorDict[i]=pid.grey + + return chrColorDict + + + def drawProbeSetPosition(self, canvas, plotXScale, offset= (40, 120, 80, 10), zoom = 1, startMb = None, endMb = None): + if len(self.traitList) != 1: + return + + xLeftOffset, xRightOffset, yTopOffset, yBottomOffset = offset + plotWidth = canvas.size[0] - xLeftOffset - xRightOffset + plotHeight = canvas.size[1] - yTopOffset - yBottomOffset + yZero = canvas.size[1] - yBottomOffset + fontZoom = zoom + if zoom == 2: + fontZoom = 1.5 + + try: + Chr = self.traitList[0].chr # self.traitListChr =self.traitList[0].chr=_vals need to change to chrList and mbList + Mb = self.traitList[0].mb # self.traitListMb =self.traitList[0].mb=_vals + except: + return + + if self.plotScale == 'physic': + if self.selectedChr > -1: + if self.genotype[0].name != Chr or Mb < self.startMb or Mb > self.endMb: + return + else: + locPixel = xLeftOffset + (Mb-self.startMb)*plotXScale + else: + locPixel = xLeftOffset + for i, _chr in enumerate(self.genotype): + if _chr.name != Chr: + locPixel += (self.ChrLengthDistList[i] + self.GraphInterval)*plotXScale + else: + locPixel += Mb*plotXScale + break + else: + if self.selectedChr > -1: + if self.genotype[0].name != Chr: + return + else: + for i, _locus in enumerate(self.genotype[0]): + #the trait's position is on the left of the first genotype + if i==0 and _locus.Mb >= Mb: + locPixel=-1 + break + + #the trait's position is between two traits + if i > 0 and self.genotype[0][i-1].Mb < Mb and _locus.Mb >= Mb: + locPixel = xLeftOffset + plotXScale*(self.genotype[0][i-1].cM+(_locus.cM-self.genotype[0][i-1].cM)*(Mb -self.genotype[0][i-1].Mb)/(_locus.Mb-self.genotype[0][i-1].Mb)) + break + + #the trait's position is on the right of the last genotype + if i==len(self.genotype[0]) and Mb>=_locus.Mb: + locPixel = -1 + else: + locPixel = xLeftOffset + for i, _chr in enumerate(self.genotype): + if _chr.name != Chr: + locPixel += (self.ChrLengthDistList[i] + self.GraphInterval)*plotXScale + else: + locPixel += (Mb*(_chr[-1].cM-_chr[0].cM)/self.ChrLengthCMList[i])*plotXScale + break + if locPixel >= 0: + traitPixel = ((locPixel, yZero), (locPixel-6, yZero+12), (locPixel+6, yZero+12)) + canvas.drawPolygon(traitPixel, edgeColor=pid.black, fillColor=self.TRANSCRIPT_LOCATION_COLOR, closed=1) + + if self.legendChecked: + startPosY = 15 + nCol = 2 + smallLabelFont = pid.Font(ttf="trebuc", size=12, bold=1) + leftOffset = xLeftOffset+(nCol-1)*200 + canvas.drawPolygon(((leftOffset+6, startPosY-6), (leftOffset, startPosY+6), (leftOffset+12, startPosY+6)), edgeColor=pid.black, fillColor=self.TRANSCRIPT_LOCATION_COLOR, closed=1) + canvas.drawString("Sequence Site", (leftOffset+15), (startPosY+5), smallLabelFont, self.TOP_RIGHT_INFO_COLOR) + + # build dict based on plink result, key is chr, value is list of [snp,BP,pValue] + def getPlinkResultDict(self,outputFileName='',thresholdPvalue=-1,ChrOrderIdNameDict={}): + + ChrList =self.ChrList + plinkResultDict={} + + plinkResultfp = open("%s%s.qassoc"% (webqtlConfig.TMPDIR, outputFileName), "rb") + + headerLine=plinkResultfp.readline()# read header line + line = plinkResultfp.readline() + + valueList=[] # initialize value list, this list will include snp, bp and pvalue info + pValueList=[] + count=0 + + while line: + #convert line from str to list + lineList=self.buildLineList(line=line) + + # only keep the records whose chromosome name is in db + if ChrOrderIdNameDict.has_key(int(lineList[0])) and lineList[-1] and lineList[-1].strip()!='NA': + + chrName=ChrOrderIdNameDict[int(lineList[0])] + snp = lineList[1] + BP = lineList[2] + pValue = float(lineList[-1]) + pValueList.append(pValue) + + if plinkResultDict.has_key(chrName): + valueList=plinkResultDict[chrName] + + # pvalue range is [0,1] + if thresholdPvalue >=0 and thresholdPvalue<=1: + if pValue < thresholdPvalue: + valueList.append((snp,BP,pValue)) + count+=1 + + plinkResultDict[chrName]=valueList + valueList=[] + else: + if thresholdPvalue>=0 and thresholdPvalue<=1: + if pValue < thresholdPvalue: + valueList.append((snp,BP,pValue)) + count+=1 + + if valueList: + plinkResultDict[chrName]=valueList + + valueList=[] + + + line =plinkResultfp.readline() + else: + line=plinkResultfp.readline() + + if pValueList: + minPvalue= min(pValueList) + else: + minPvalue=0 + + return count,minPvalue,plinkResultDict + + + ###################################################### + # input: line: str,one line read from file + # function: convert line from str to list; + # output: lineList list + ####################################################### + def buildLineList(self,line=None): + + lineList = string.split(string.strip(line),' ')# irregular number of whitespaces between columns + lineList =[ item for item in lineList if item <>''] + lineList = map(string.strip, lineList) + + return lineList + + #added by NL: automatically generate pheno txt file for PLINK based on strainList passed from dataEditing page + def genPhenoTxtFileForPlink(self,phenoFileName='', RISetName='', probesetName='', valueDict={}): + pedFileStrainList=self.getStrainNameFromPedFile(RISetName=RISetName) + outputFile = open("%s%s.txt"%(webqtlConfig.TMPDIR,phenoFileName),"wb") + headerLine = 'FID\tIID\t%s\n'%probesetName + outputFile.write(headerLine) + + newValueList=[] + + #if valueDict does not include some strain, value will be set to -9999 as missing value + for item in pedFileStrainList: + try: + value=valueDict[item] + value=str(value).replace('value=','') + value=value.strip() + except: + value=-9999 + + newValueList.append(value) + + + newLine='' + for i, strain in enumerate(pedFileStrainList): + j=i+1 + value=newValueList[i] + newLine+='%s\t%s\t%s\n'%(strain, strain, value) + + if j%1000==0: + outputFile.write(newLine) + newLine='' + + if newLine: + outputFile.write(newLine) + + outputFile.close() + + # get strain name from ped file in order + def getStrainNameFromPedFile(self, RISetName=''): + pedFileopen= open("%splink/%s.ped"%(webqtlConfig.HTMLPATH, RISetName),"r") + line =pedFileopen.readline() + strainNameList=[] + + while line: + lineList=string.split(string.strip(line),'\t') + lineList=map(string.strip,lineList) + + strainName=lineList[0] + strainNameList.append(strainName) + + line =pedFileopen.readline() + + return strainNameList + + ################################################################ + # Generate Chr list, Chr OrderId and Retrieve Length Information + ################################################################ + def getChrNameOrderIdLength(self,RISet=''): + + try: + query = """ + Select + Chr_Length.Name,Chr_Length.OrderId,Length from Chr_Length, InbredSet + where + Chr_Length.SpeciesId = InbredSet.SpeciesId AND + InbredSet.Name = '%s' + Order by OrderId + """ % (RISet) + self.cursor.execute(query) + + results =self.cursor.fetchall() + ChrList=[] + ChrLengthMbList=[] + ChrNameOrderIdDict={} + ChrOrderIdNameDict={} + + for item in results: + ChrList.append(item[0]) + ChrNameOrderIdDict[item[0]]=item[1] # key is chr name, value is orderId + ChrOrderIdNameDict[item[1]]=item[0] # key is orderId, value is chr name + ChrLengthMbList.append(item[2]) + + except: + ChrList=[] + ChrNameOrderIdDict={} + ChrLengthMbList=[] + + return ChrList,ChrNameOrderIdDict,ChrOrderIdNameDict,ChrLengthMbList diff --git a/wqflask/wqflask/marker_regression/__init__.py b/wqflask/wqflask/marker_regression/__init__.py new file mode 100644 index 00000000..e69de29b --- /dev/null +++ b/wqflask/wqflask/marker_regression/__init__.py diff --git a/wqflask/wqflask/marker_regression/marker_regression.py b/wqflask/wqflask/marker_regression/marker_regression.py new file mode 100755 index 00000000..ed01a3fa --- /dev/null +++ b/wqflask/wqflask/marker_regression/marker_regression.py @@ -0,0 +1,1648 @@ +from __future__ import absolute_import, print_function, division + +from base.trait import GeneralTrait +from base import data_set #import create_dataset + +from pprint import pformat as pf + +import time +import string +import math +#from math import * +#import piddle +import sys +import os +import httplib +import urllib + +from htmlgen import HTMLgen2 as HT +from utility import Plot +from wqflask.interval_analyst import GeneUtil +from base.trait import GeneralTrait +from base.data_set import create_dataset +from base.templatePage import templatePage +from utility import webqtlUtil +from base import webqtlConfig +from dbFunction import webqtlDatabaseFunction +from base.GeneralObject import GeneralObject + +import reaper +import cPickle +from utility.THCell import THCell +from utility.TDCell import TDCell + + +class MarkerRegression(object): + + #def __init__(self, start_vars): + # + # print("[mike] Now start_vars is:", pf(start_vars)) + # + # self.dataset = data_set.create_dataset(start_vars['dataset_name']) + # self.this_trait = GeneralTrait(dataset=self.dataset.name, + # name=start_vars['trait_id'], + # cellid=None) + # + # print("self.this_trait is: ", pf(self.this_trait)) + # print("self.dataset is: ", pf(self.dataset)) + + def __init__(self, start_vars): + #templatePage.__init__(self, fd) + + #if not self.openMysql(): + # return + + self.dataset = create_dataset(start_vars['dataset_name']) + + #self.initializeParameters(start_vars) + + #filename= webqtlUtil.genRandStr("Itvl_") + #ChrList,ChrNameOrderIdDict,ChrOrderIdNameDict,ChrtLengthMbList= self.getChrNameOrderIdLength(RISet=fd.RISet) + + if False: # For PLINK + + traitInfoList = string.split(string.strip(fd.identification),':') + probesetName = string.strip(traitInfoList[-1]) + plinkOutputFileName= webqtlUtil.genRandStr("%s_%s_"%(fd.RISet,probesetName)) + + # get related values from fd.allTraitData; the format of 'allTraitValueDict'is {strainName1: value=-0.2...} + fd.readData() + allTraitValueDict = fd.allTraitData + + #automatically generate pheno txt file for PLINK + self.genPhenoTxtFileForPlink(phenoFileName=plinkOutputFileName,RISetName=fd.RISet,probesetName=probesetName, valueDict=allTraitValueDict) + # os.system full path is required for input and output files; specify missing value is -9999 + plink_command = '%splink/plink --noweb --ped %splink/%s.ped --no-fid --no-parents --no-sex --no-pheno --map %splink/%s.map --pheno %s/%s.txt --pheno-name %s --missing-phenotype -9999 --out %s%s --assoc ' % (webqtlConfig.HTMLPATH, webqtlConfig.HTMLPATH, fd.RISet, webqtlConfig.HTMLPATH, fd.RISet, webqtlConfig.TMPDIR, plinkOutputFileName, probesetName, webqtlConfig.TMPDIR, plinkOutputFileName) + + os.system(plink_command) + + if fd.identification: + heading2 = HT.Paragraph('Trait ID: %s' % fd.identification) + heading2.__setattr__("class","subtitle") + self.dict['title'] = '%s: Genome Association' % fd.identification + else: + heading2 = "" + self.dict['title'] = 'Genome Association' + + if fd.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 = "" + + heading = HT.Paragraph('Trait Data Entered for %s Set' % fd.RISet) + heading.__setattr__("class","title") + + # header info part:Trait Data Entered for HLC Set & Trait ID: + headerdiv = HT.TR(HT.TD(heading, heading2,heading3, width='45%',valign='top', align='left', bgColor='#eeeeee')) + + self.ChrList=ChrList # get chr name from '1' to 'X' + self.ChrLengthMbList = ChrLengthMbList + + # build plink result dict based on chr, key is chr name, value is in list type including Snpname, bp and pvalue info + plinkResultDict={} + count,minPvalue,plinkResultDict =self.getPlinkResultDict(outputFileName=plinkOutputFileName,thresholdPvalue=self.pValue,ChrOrderIdNameDict=ChrOrderIdNameDict) + + # if can not find results which are matched with assigned p-value, system info will show up + if count >0: + + #for genome association report table + reportTable="" + # sortable table object + resultstable,tblobj,bottomInfo = self.GenReportForPLINK(ChrNameOrderIdDict=ChrNameOrderIdDict, RISet=fd.RISet,plinkResultDict=plinkResultDict,thresholdPvalue=self.pValue,chrList=self.ChrList) + + # creat object for result table for sort function + objfile = open('%s.obj' % (webqtlConfig.TMPDIR+filename), 'wb') + cPickle.dump(tblobj, objfile) + objfile.close() + + sortby = ("Index", "up") + reportTable =HT.Div(webqtlUtil.genTableObj(tblobj=tblobj, file=filename, sortby=sortby, tableID = "sortable", addIndex = "0"), Id="sortable") + + descriptionTable = HT.TableLite(border=0, cellpadding=0, cellspacing=0) + descriptionTable.append(HT.TR(HT.TD(reportTable, colspan=3))) + descriptionTable.append(HT.TR(HT.TD(HT.BR(),HT.BR()))) + descriptionTable.append(bottomInfo) + + # get each chr's length + self.ChrLengthMbList = map(lambda x: x/1000000.0, self.ChrLengthMbList) # change unit from bp to mb + self.ChrLengthMbSum = reduce(lambda x, y:x+y, self.ChrLengthMbList, 0.0)# get total length of all chrs + if self.ChrLengthMbList: + self.GraphInterval = self.ChrLengthMbSum/(len(self.ChrLengthMbList)*12) #Empirical Mb interval + else: + self.GraphInterval = 1 + + # for human data, there's no CM value + self.ChrLengthCMList = [] + self.ChrLengthCMSum = 0 + + # begin: common part with human data + intCanvas = pid.PILCanvas(size=(self.graphWidth,self.graphHeight)) + gifmap = self.plotIntMappingForPLINK(fd, intCanvas, startMb = self.startMb, endMb = self.endMb, plinkResultDict=plinkResultDict) + + intCanvas.save(os.path.join(webqtlConfig.IMGDIR, filename), format='png') + intImg=HT.Image('/image/'+filename+'.png', border=0, usemap='#WebQTLImageMap') + + TD_LR = HT.TR(HT.TD(HT.Blockquote(gifmap,intImg, HT.P()), bgColor='#eeeeee', height = 200)) + self.dict['body'] = str(headerdiv)+str(TD_LR)+str(resultstable)+str(HT.TR(HT.TD(descriptionTable))) + + else: + heading = "Genome Association" + detail = ['There is no association with marker that meets this criteria. Please provide a less stringend threshold. The minimun p-value is %s.'%minPvalue] + self.error(heading=heading,detail=detail) + return + + else: # QTLreaper result + #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 = self.dataset.group.read_genotype_file() + print("[black]:", genotype) + + _strains, _vals, _vars, N = fd.informativeStrains(_genotype.prgy, weightedRegression) + + if fd.identification: + heading2 = HT.Paragraph('Trait ID: %s' % fd.identification) + heading2.__setattr__("class","subtitle") + self.dict['title'] = '%s: Genome Association' % fd.identification + else: + heading2 = "" + self.dict['title'] = 'Genome Association' + + if fd.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 = "Genome Association" + detail = ['Fewer than %d strain data were entered for %s data set. No mapping attempted.' % (webqtlConfig.KMININFORMATIVE, fd.RISet)] + self.error(heading=heading,detail=detail) + return + else: + heading = HT.Paragraph('Trait Data Entered for %s Set' % fd.RISet) + heading.__setattr__("class","title") + + datadiv = HT.TD(heading, heading2,heading3, width='45%',valign='top', align='left', bgColor='#eeeeee') + resultstable,tblobj,bottomInfo = self.GenReport(ChrNameOrderIdDict,fd, _genotype, _strains, _vals, _vars) + #resultstable = self.GenReport(fd, _genotype, _strains, _vals, _vars) + + # creat object for result table for sort function + objfile = open('%s.obj' % (webqtlConfig.TMPDIR+filename), 'wb') + cPickle.dump(tblobj, objfile) + objfile.close() + + sortby = ("Index", "up") + reportTable =HT.Div(webqtlUtil.genTableObj(tblobj=tblobj, file=filename, sortby=sortby, tableID = "sortable", addIndex = "0"), Id="sortable") + + descriptionTable = HT.TableLite(border=0, cellpadding=0, cellspacing=0) + descriptionTable.append(HT.TR(HT.TD(reportTable, colspan=3))) + descriptionTable.append(HT.TR(HT.TD(HT.BR(),HT.BR()))) + descriptionTable.append(bottomInfo) + + self.traitList=_vals + + ##########################plot####################### + + ################################################################ + # Generate Chr list and Retrieve Length Information + ################################################################ + self.genotype= _genotype + self.ChrList = [("All", -1)] + + for i, indChr in enumerate(self.genotype): + self.ChrList.append((indChr.name, i)) + + self.cursor.execute(""" + Select + Length from Chr_Length, InbredSet + where + Chr_Length.SpeciesId = InbredSet.SpeciesId AND + InbredSet.Name = '%s' AND + Chr_Length.Name in (%s) + Order by + OrderId + """ % (fd.RISet, string.join(map(lambda X: "'%s'" % X[0], self.ChrList[1:]), ", "))) + + self.ChrLengthMbList = self.cursor.fetchall() + self.ChrLengthMbList = map(lambda x: x[0]/1000000.0, self.ChrLengthMbList) + self.ChrLengthMbSum = reduce(lambda x, y:x+y, self.ChrLengthMbList, 0.0) + if self.ChrLengthMbList: + self.MbGraphInterval = self.ChrLengthMbSum/(len(self.ChrLengthMbList)*12) #Empirical Mb interval + else: + self.MbGraphInterval = 1 + + self.ChrLengthCMList = [] + for i, _chr in enumerate(self.genotype): + self.ChrLengthCMList.append(_chr[-1].cM - _chr[0].cM) + self.ChrLengthCMSum = reduce(lambda x, y:x+y, self.ChrLengthCMList, 0.0)# used for calculate plot scale + + self.GraphInterval = self.MbGraphInterval #Mb + + # begin: common part with human data + intCanvas = pid.PILCanvas(size=(self.graphWidth,self.graphHeight)) + gifmap = self.plotIntMapping(fd, intCanvas, startMb = self.startMb, endMb = self.endMb, showLocusForm= "") + filename= webqtlUtil.genRandStr("Itvl_") + intCanvas.save(os.path.join(webqtlConfig.IMGDIR, filename), format='png') + intImg=HT.Image('/image/'+filename+'.png', border=0, usemap='#WebQTLImageMap') + + ################################################################ + # footnote goes here + ################################################################ + btminfo = HT.Paragraph(Id="smallsize") #Small('More information about this graph is available here.') + + if (self.additiveChecked): + btminfo.append(HT.BR(), 'A positive additive coefficient (', HT.Font('green', color='green'), ' line) indicates that %s alleles increase trait values. In contrast, a negative additive coefficient (' % fd.ppolar, HT.Font('red', color='red'), ' line) indicates that %s alleles increase trait values.' % fd.mpolar) + + + TD_LR = HT.TR(HT.TD(HT.Blockquote(gifmap,intImg, HT.P()), bgColor='#eeeeee', height = 200)) + + self.dict['body'] = str(datadiv)+str(TD_LR)+str(resultstable)+str(HT.TR(HT.TD(descriptionTable))) + + # end: common part with human data + + + + + # add by NL 10-2-2011 + def initializeParameters(self, fd): + """ + Initializes all of the MarkerRegressionPage class parameters, + acquiring most values from the formdata (fd) + """ + ################################### + # manhattam plot parameters + ################################### + + self.graphHeight = 600 + self.graphWidth = 1280 + self.plotScale = 'physic' + self.selectedChr = -1 + self.GRAPH_BACK_DARK_COLOR = pid.HexColor(0xF1F1F9) + self.GRAPH_BACK_LIGHT_COLOR = pid.HexColor(0xFBFBFF) + self.LRS_COLOR = pid.HexColor(0x0000FF) + self.LRS_LOD ='LRS' + self.lrsMax = float(fd.formdata.getvalue('lrsMax', 0)) + self.startMb = fd.formdata.getvalue('startMb', "-1") + self.endMb = fd.formdata.getvalue('endMb', "-1") + self.mappingMethodId = fd.formdata.getvalue('mappingMethodId', "0") + self.permChecked=True + self.multipleInterval=False + self.SIGNIFICANT_WIDTH = 5 + self.SUGGESTIVE_WIDTH = 5 + self.SIGNIFICANT_COLOR = pid.HexColor(0xEBC7C7) + self.SUGGESTIVE_COLOR = pid.gainsboro + self.colorCollection = [self.LRS_COLOR] + self.additiveChecked= True + self.ADDITIVE_COLOR_POSITIVE = pid.green + self.legendChecked =False + self.pValue=float(fd.formdata.getvalue('pValue',-1)) + + # allow user to input p-value greater than 1, + # in this case, the value will be treated as -lgP value. so the input value needs to be transferred to power of 10 format + if self.pValue >1: + self.pValue =10**-(self.pValue) + + try: + self.startMb = float(self.startMb) + self.endMb = float(self.endMb) + if self.startMb > self.endMb: + temp = self.startMb + self.startMb = self.endMb + self.endMb = temp + #minimal distance 10bp + if self.endMb - self.startMb < 0.00001: + self.endMb = self.startMb + 0.00001 + except: + self.startMb = self.endMb = -1 + + def GenReportForPLINK(self, ChrNameOrderIdDict={},RISet='',plinkResultDict= {},thresholdPvalue=-1,chrList=[]): + + 'Create an HTML division which reports any loci which are significantly associated with the submitted trait data.' + ######################################### + # Genome Association 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':RISet+"Geno",'CellID':'_', \ + 'RISet':RISet, 'incparentsf1':'on'} + for key in hddn.keys(): + locusForm.append(HT.Input(name=key, value=hddn[key], type='hidden')) + + regressionHeading = HT.Paragraph('Genome Association Report') + regressionHeading.__setattr__("class","title") + + filename= webqtlUtil.genRandStr("GenomeAsscociation_") + fpText = open('%s.txt' % (webqtlConfig.TMPDIR+filename), 'wb') + fpText.write('The loci meet the criteria of P-Value <= %3.6f.\n'%thresholdPvalue) + pValueInfo =HT.Paragraph('The loci meet the criteria of P-Value <= %3.6f.\n'%thresholdPvalue) + + textUrl = HT.Href(text = 'Download', url= '/tmp/'+filename+'.txt', target = "_blank", Class='fs12 fwn') + bottomInfo = HT.TR(HT.TD(HT.Paragraph(textUrl, ' result in tab-delimited text format.', HT.BR(), HT.BR(),Class="fs12 fwn"), colspan=3)) + + tblobj={} # build dict for genTableObj function; keys include header and body + tblobj_header = [] # value of key 'header' + tblobj_body=[] # value of key 'body' + reportHeaderRow=[] # header row list for tblobj_header (html part) + headerList=['Index','SNP Name','Chr','Mb','-log(P)'] + headerStyle="fs14 fwb ffl b1 cw cbrb" # style of the header + cellColorStyle = "fs13 b1 fwn c222" # style of the cells + + if headerList: + for ncol, item in enumerate(headerList): + reportHeaderRow.append(THCell(HT.TD(item, Class=headerStyle, valign='bottom',nowrap='ON'),text=item, idx=ncol)) + #download file for table headers' names + fpText.write('SNP_Name\tChromosome\tMb\t-log(P)\n') + + tblobj_header.append(reportHeaderRow) + tblobj['header']=tblobj_header + + index=1 + for chr in chrList: + + if plinkResultDict.has_key(chr): + if chr in ChrNameOrderIdDict.keys(): + chrOrderId =ChrNameOrderIdDict[chr] + else: + chrOrderId=chr + + valueList=plinkResultDict[chr] + + for value in valueList: + reportBodyRow=[] # row list for tblobj_body (html part) + snpName=value[0] + bp=value[1] + mb=int(bp)/1000000.0 + + try: + pValue =float(value[2]) + except: + pValue =1 + formattedPvalue = -math.log10(pValue) + + formattedPvalue = webqtlUtil.SciFloat(formattedPvalue) + dbSnprs=snpName.replace('rs','') + SnpHref = HT.Href(text=snpName, url="http://www.ncbi.nlm.nih.gov/projects/SNP/snp_ref.cgi?rs=%s"%dbSnprs, target="_blank") + + selectCheck=HT.Input(type="checkbox", Class="checkbox", name="index",value=index, onClick="highlight(this)") + reportBodyRow.append(TDCell(HT.TD(str(index),selectCheck, align='right',Class=cellColorStyle,nowrap='ON'),str(index),index)) + reportBodyRow.append(TDCell(HT.TD(SnpHref, Class=cellColorStyle,nowrap='ON'),snpName, snpName)) + reportBodyRow.append(TDCell(HT.TD(chr, Class=cellColorStyle, align="center",nowrap='ON'),chr, chrOrderId)) + reportBodyRow.append(TDCell(HT.TD('%3.6f'%mb, Class=cellColorStyle, align="center",nowrap='ON'),mb, mb)) + reportBodyRow.append(TDCell(HT.TD(formattedPvalue, Class=cellColorStyle, align="center",nowrap='ON'),formattedPvalue, float(formattedPvalue))) + + fpText.write('%s\t%s\t%3.6f\t%s\n' % (snpName, str(chr), mb, formattedPvalue)) + index+=1 + + tblobj_body.append(reportBodyRow) + + tblobj['body']=tblobj_body + rv=HT.TR(HT.TD(regressionHeading,pValueInfo, locusForm, HT.P(), width='55%',valign='top', align='left',bgColor='#eeeeee')) + + return rv, tblobj,bottomInfo + + + def GenReport(self, ChrNameOrderIdDict,fd, _genotype, _strains, _vals, _vars= []): + 'Create an HTML division which reports any loci which are significantly associated with the submitted trait data.' + #calculate QTL for each trait + self.qtlresults = [] + if webqtlUtil.ListNotNull(_vars): + qtlresults = _genotype.regression(strains = _strains, trait = _vals, variance = _vars) + LRSArray = _genotype.permutation(strains = _strains, trait = _vals, variance = _vars, nperm=fd.nperm) + else: + qtlresults = _genotype.regression(strains = _strains, trait = _vals) + LRSArray = _genotype.permutation(strains = _strains, trait = _vals,nperm=fd.nperm) + + self.qtlresults.append(qtlresults) + + filename= webqtlUtil.genRandStr("GenomeAsscociation_") + + # set suggestive, significant and highly significant LRS + 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) + + self.significance =fd.significance + self.suggestive = fd.suggestive + self.highlysignificant = LRSArray[int(fd.nperm*0.99-1)] + _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:] + + ######################################### + # Permutation Graph + ######################################### + 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) + + permutationHeading = HT.Paragraph('Histogram of Permutation Test') + permutationHeading.__setattr__("class","title") + + permutation = HT.TableLite() + permutation.append(HT.TR(HT.TD(img))) + + + ######################################### + # Genome Association 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('Genome Association Report') + regressionHeading.__setattr__("class","title") + # report is the info part above report table + 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') + fpText.write('Suggestive LRS =%3.2f\n'%self.suggestive) + fpText.write('Significant LRS =%3.2f\n'%self.significance) + fpText.write('Highly Significant LRS =%3.2f\n'%self.highlysignificant) + LRSInfo =HT.Paragraph(' Suggestive LRS =%3.2f\n'%fd.suggestive, HT.BR(), ' Significant LRS =%3.2f\n'%fd.significance,HT.BR(),' Highly Significant LRS =%3.2f\n' % self.highlysignificant) + + textUrl = HT.Href(text = 'Download', url= '/tmp/'+filename+'.txt', target = "_blank", Class='fs12 fwn') + + bottomInfo = HT.TR(HT.TD(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"))) + + tblobj={} # build dict for genTableObj function; keys include header and body + tblobj_header = [] # value of key 'header' + tblobj_body=[] # value of key 'body' + reportHeaderRow=[] # header row list for tblobj_header (html part) + headerStyle="fs14 fwb ffl b1 cw cbrb" # style of the header + cellColorStyle = "fs13 b1 fwn c222" # style of the cells + + headerList=['Index','LRS','Chr','Mb','Locus','Additive Effect'] + for ncol, item in enumerate(headerList): + reportHeaderRow.append(THCell(HT.TD(item, Class=headerStyle, valign='bottom',nowrap='ON'),text=item, idx=ncol)) + + if fd.genotype.type == 'intercross': + ncol =len(headerList) + reportHeaderRow.append(THCell(HT.TD('Dominance Effect', Class=headerStyle, valign='bottom',nowrap='ON'),text='Dominance Effect', idx=ncol)) + + #download file for table headers' names + fpText.write('LRS\tChromosome\tMb\tLocus\tAdditive Effect\tDominance Effect\n') + + index=1 + for ii in filtered: + #add by NL 06-20-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 + + if LRS > fd.significance: + lrs = HT.TD(HT.Font('%3.3f*' % LRS, color='#FF0000'),Class=cellColorStyle) + else: + lrs = HT.TD('%3.3f' % LRS,Class=cellColorStyle) + + if ii.locus.chr in ChrNameOrderIdDict.keys(): + chrOrderId =ChrNameOrderIdDict[ii.locus.chr] + else: + chrOrderId=ii.locus.chr + + reportBodyRow=[] # row list for tblobj_body (html part) + selectCheck=HT.Input(type="checkbox", Class="checkbox", name="index",value=index, onClick="highlight(this)") + reportBodyRow.append(TDCell(HT.TD(str(index),selectCheck, align='right',Class=cellColorStyle,nowrap='ON'),str(index),index)) + reportBodyRow.append(TDCell(lrs,LRS, LRS)) + reportBodyRow.append(TDCell(HT.TD(ii.locus.chr, Class=cellColorStyle, align="center",nowrap='ON'),ii.locus.chr, chrOrderId)) + reportBodyRow.append(TDCell(HT.TD('%3.6f'%ii.locus.Mb, Class=cellColorStyle, align="center",nowrap='ON'),ii.locus.Mb, ii.locus.Mb)) + reportBodyRow.append(TDCell(HT.TD(HT.Href(text=ii.locus.name, url = "javascript:showTrait('%s','%s');" % (locusFormName, ii.locus.name), Class='normalsize'), Class=cellColorStyle, align="center",nowrap='ON'),ii.locus.name, ii.locus.name)) + reportBodyRow.append(TDCell(HT.TD('%3.3f' % ii.additive, Class=cellColorStyle, align="center",nowrap='ON'),ii.additive, ii.additive)) + reportBodyRow.append(TDCell(HT.TD('%3.3f' % ii.dominance, Class=cellColorStyle, align="center",nowrap='ON'),ii.dominance, ii.dominance)) + + fpText.write('%2.3f\t%s\t%3.6f\t%s\t%2.3f\t%2.3f\n' % (LRS, ii.locus.chr, ii.locus.Mb, ii.locus.name, ii.additive, ii.dominance)) + index+=1 + tblobj_body.append(reportBodyRow) + else: + #download file for table headers' names + fpText.write('LRS\tChromosome\tMb\tLocus\tAdditive Effect\n') + + index=1 + for ii in filtered: + #add by NL 06-20-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 + + if LRS > fd.significance: + lrs = HT.TD(HT.Font('%3.3f*' % LRS, color='#FF0000'),Class=cellColorStyle) + else: + lrs = HT.TD('%3.3f' % LRS,Class=cellColorStyle) + + if ii.locus.chr in ChrNameOrderIdDict.keys(): + chrOrderId =ChrNameOrderIdDict[ii.locus.chr] + else: + chrOrderId=ii.locus.chr + + reportBodyRow=[] # row list for tblobj_body (html part) + selectCheck=HT.Input(type="checkbox", Class="checkbox", name="index",value=index, onClick="highlight(this)") + reportBodyRow.append(TDCell(HT.TD(str(index),selectCheck, align='right',Class=cellColorStyle,nowrap='ON'),str(index),index)) + reportBodyRow.append(TDCell(lrs,LRS, LRS)) + reportBodyRow.append(TDCell(HT.TD(ii.locus.chr, Class=cellColorStyle, align="center",nowrap='ON'),ii.locus.chr, chrOrderId)) + reportBodyRow.append(TDCell(HT.TD('%3.6f'%ii.locus.Mb, Class=cellColorStyle, align="center",nowrap='ON'),ii.locus.Mb, ii.locus.Mb)) + reportBodyRow.append(TDCell(HT.TD(HT.Href(text=ii.locus.name, url = "javascript:showTrait('%s','%s');" % (locusFormName, ii.locus.name), Class='normalsize'), Class=cellColorStyle, align="center",nowrap='ON'),ii.locus.name, ii.locus.name)) + reportBodyRow.append(TDCell(HT.TD('%3.3f' % ii.additive, Class=cellColorStyle, align="center",nowrap='ON'),ii.additive, ii.additive)) + + fpText.write('%2.3f\t%s\t%3.6f\t%s\t%2.3f\n' % (LRS, ii.locus.chr, ii.locus.Mb, ii.locus.name, ii.additive)) + index+=1 + tblobj_body.append(reportBodyRow) + + tblobj_header.append(reportHeaderRow) + tblobj['header']=tblobj_header + tblobj['body']=tblobj_body + + rv=HT.TD(regressionHeading,LRSInfo,report, locusForm, HT.P(),width='55%',valign='top', align='left', bgColor='#eeeeee') + if fd.genotype.type == 'intercross': + bottomInfo.append(HT.BR(), HT.BR(), HT.Strong('Dominance Effect'),' is the difference between the mean trait value of cases heterozygous at a marker and the average mean for the two groups homozygous at this marker: e.g., BD - (BB+DD)/2]. A positive dominance effect indicates that the average phenotype of BD heterozygotes exceeds the mean of BB and DD homozygotes. No dominance deviation can be computed for a set of recombinant inbred strains or for a backcross.') + return rv,tblobj,bottomInfo + + return rv,tblobj,bottomInfo + + def plotIntMappingForPLINK(self, fd, canvas, offset= (80, 120, 20, 80), zoom = 1, startMb = None, endMb = None, showLocusForm = "",plinkResultDict={}): + #calculating margins + xLeftOffset, xRightOffset, yTopOffset, yBottomOffset = offset + + fontZoom = zoom + if zoom == 2: + fontZoom = 1.5 + + xLeftOffset = int(xLeftOffset*fontZoom) + xRightOffset = int(xRightOffset*fontZoom) + yBottomOffset = int(yBottomOffset*fontZoom) + + cWidth = canvas.size[0] + cHeight = canvas.size[1] + plotWidth = cWidth - xLeftOffset - xRightOffset + plotHeight = cHeight - yTopOffset - yBottomOffset + startPixelX = xLeftOffset + endPixelX = (xLeftOffset + plotWidth) + + #Drawing Area Height + drawAreaHeight = plotHeight + if self.plotScale == 'physic' and self.selectedChr > -1: # for single chr + drawAreaHeight -= self.ENSEMBL_BAND_HEIGHT + self.UCSC_BAND_HEIGHT+ self.WEBQTL_BAND_HEIGHT + 3*self.BAND_SPACING+ 10*zoom + if self.geneChecked: + drawAreaHeight -= self.NUM_GENE_ROWS*self.EACH_GENE_HEIGHT + 3*self.BAND_SPACING + 10*zoom + else: + if self.selectedChr > -1: + drawAreaHeight -= 20 + else:# for all chrs + drawAreaHeight -= 30 + + #Image map + gifmap = HT.Map(name='WebQTLImageMap') + + newoffset = (xLeftOffset, xRightOffset, yTopOffset, yBottomOffset) + # Draw the alternating-color background first and get plotXScale + plotXScale = self.drawGraphBackgroundForPLINK(canvas, gifmap, offset=newoffset, zoom= zoom, startMb=startMb, endMb = endMb,plinkResultDict=plinkResultDict) + + # Draw X axis + self.drawXAxisForPLINK(fd, canvas, drawAreaHeight, gifmap, plotXScale, showLocusForm, offset=newoffset, zoom= zoom, startMb=startMb, endMb = endMb) + # Draw manhattam plot + self.drawManhattanPlotForPLINK(canvas, drawAreaHeight, gifmap, plotXScale, offset=newoffset, zoom= zoom, startMb=startMb, endMb = endMb,plinkResultDict=plinkResultDict,thresholdPvalue=self.pValue) + + return gifmap + + + def plotIntMapping(self, fd, canvas, offset= (80, 120, 20, 80), zoom = 1, startMb = None, endMb = None, showLocusForm = ""): + #calculating margins + xLeftOffset, xRightOffset, yTopOffset, yBottomOffset = offset + + fontZoom = zoom + if zoom == 2: + fontZoom = 1.5 + + xLeftOffset = int(xLeftOffset*fontZoom) + xRightOffset = int(xRightOffset*fontZoom) + yBottomOffset = int(yBottomOffset*fontZoom) + + cWidth = canvas.size[0] + cHeight = canvas.size[1] + plotWidth = cWidth - xLeftOffset - xRightOffset + plotHeight = cHeight - yTopOffset - yBottomOffset + startPixelX = xLeftOffset + endPixelX = (xLeftOffset + plotWidth) + + #Drawing Area Height + drawAreaHeight = plotHeight + if self.plotScale == 'physic' and self.selectedChr > -1: # for single chr + drawAreaHeight -= self.ENSEMBL_BAND_HEIGHT + self.UCSC_BAND_HEIGHT+ self.WEBQTL_BAND_HEIGHT + 3*self.BAND_SPACING+ 10*zoom + if self.geneChecked: + drawAreaHeight -= self.NUM_GENE_ROWS*self.EACH_GENE_HEIGHT + 3*self.BAND_SPACING + 10*zoom + else:# for all chrs + if self.selectedChr > -1: + drawAreaHeight -= 20 + else: + drawAreaHeight -= 30 + + #Image map + gifmap = HT.Map(name='WebQTLImageMap') + + newoffset = (xLeftOffset, xRightOffset, yTopOffset, yBottomOffset) + # Draw the alternating-color background first and get plotXScale + plotXScale = self.drawGraphBackground(canvas, gifmap, offset=newoffset, zoom= zoom, startMb=startMb, endMb = endMb) + + # Draw X axis + self.drawXAxis(fd, canvas, drawAreaHeight, gifmap, plotXScale, showLocusForm, offset=newoffset, zoom= zoom, startMb=startMb, endMb = endMb) + # Draw QTL curve + self.drawQTL(canvas, drawAreaHeight, gifmap, plotXScale, offset=newoffset, zoom= zoom, startMb=startMb, endMb = endMb) + + #draw legend + if self.multipleInterval: + self.drawMultiTraitName(fd, canvas, gifmap, showLocusForm, offset=newoffset) + elif self.legendChecked: + self.drawLegendPanel(fd, canvas, offset=newoffset) + else: + pass + + #draw position, no need to use a separate function + if fd.genotype.Mbmap: + self.drawProbeSetPosition(canvas, plotXScale, offset=newoffset) + + return gifmap + + + # functions for manhattam plot of markers + def drawManhattanPlotForPLINK(self, canvas, drawAreaHeight, gifmap, plotXScale, offset= (40, 120, 80, 10), zoom = 1, startMb = None, endMb = None,plinkResultDict={},thresholdPvalue=-1): + + xLeftOffset, xRightOffset, yTopOffset, yBottomOffset = offset + plotWidth = canvas.size[0] - xLeftOffset - xRightOffset + plotHeight = canvas.size[1] - yTopOffset - yBottomOffset + fontZoom = zoom + if zoom == 2: + fontZoom = 1.5 + + # INTERCROSS = (self.genotype.type=="intercross") + INTERCROSS ='' #?????? + + ChrLengthDistList = self.ChrLengthMbList + drawRegionDistance = self.ChrLengthMbSum + GraphInterval=self.GraphInterval + pvalueHeightThresh = drawAreaHeight - 80 #ZS: Otherwise the plot gets very close to the chromosome labels + + #draw the pvalue scale + #We first determine whether or not we are using a sliding scale. + #If so, we need to compute the maximum pvalue value to determine where the max y-value should be, and call this pvalueMax. + #pvalueTop is then defined to be above the pvalueMax by enough to add one additional pvalueScale increment. + #if we are using a set-scale, then we set pvalueTop to be the user's value, and pvalueMax doesn't matter. + + # for human data we use p value instead of lrs + pValueList=[] + for key in plinkResultDict: + valueList = plinkResultDict[key] + for item in valueList: + pValue = item[-1] + pValueList.append(pValue) + + formattedPValueList=[] + for pValue in pValueList: + try: + pValue=float(pValue) + except: + pValue =1 + formattedpValue = -math.log10(pValue) + formattedPValueList.append(formattedpValue) + + #sliding scale + pvalueMax = max(formattedPValueList) + #pvalueMax =pvalueMax +1 + # no permutation result for plink func: GenReport() + pvalueMin = int(-math.log10(thresholdPvalue)) + + if pvalueMax> 100: + pvalueScale = 20.0 + elif pvalueMax > 20: + pvalueScale = 5.0 + elif pvalueMax > 7.5: + pvalueScale = 2.5 + else: + pvalueScale = 1.0 + + # the base line for x-axis is -log(thresholdPvalue) + pvalueAxisList = Plot.frange(pvalueMin, pvalueMax, pvalueScale) + #make sure the user's value appears on the y-axis + #ZS: There is no way to do this without making the position of the points not directly proportional to a given distance on the y-axis + #tempPvalueMax=round(pvalueMax) + tempPvalueMax = pvalueAxisList[len(pvalueAxisList)-1] + pvalueScale + pvalueAxisList.append(tempPvalueMax) + + #ZS: I don't understand this; the if statement will be true for any number that isn't exactly X.5. + #if abs(tempPvalueMax-pvalueMax) <0.5: + # tempPvalueMax=tempPvalueMax+1 + # pvalueAxisList.append(tempPvalueMax) + + #draw the "pvalue" string to the left of the axis + pvalueScaleFont=pid.Font(ttf="verdana", size=14*fontZoom, bold=0) + pvalueLODFont=pid.Font(ttf="verdana", size=14*zoom*1.5, bold=0) + yZero = yTopOffset + plotHeight + + #yAxis label display area + yAxis_label ='-log(P)' + canvas.drawString(yAxis_label, xLeftOffset - canvas.stringWidth("999.99", font=pvalueScaleFont) - 10*zoom, \ + yZero - 150, font=pvalueLODFont, color=pid.black, angle=90) + + for i,item in enumerate(pvalueAxisList): + ypvalue = yZero - (float(i)/float(len(pvalueAxisList) - 1)) * pvalueHeightThresh + canvas.drawLine(xLeftOffset, ypvalue, xLeftOffset - 4, ypvalue, color=self.LRS_COLOR, width=1*zoom) + scaleStr = "%2.1f" % item + #added by NL 6-24-2011:Y-axis scale display + canvas.drawString(scaleStr, xLeftOffset-4-canvas.stringWidth(scaleStr, font=pvalueScaleFont)-5, ypvalue+3, font=pvalueScaleFont, color=self.LRS_COLOR) + + ChrList=self.ChrList + startPosX = xLeftOffset + + for i, chr in enumerate(ChrList): + + if plinkResultDict.has_key(chr): + plinkresultList = plinkResultDict[chr] + + m = 0 + #add by NL 06-24-2011: for mahanttam plot + symbolFont = pid.Font(ttf="fnt_bs", size=5,bold=0) + # color for point in each chr + chrCount=len(ChrList) + chrColorDict =self.getColorForMarker(chrCount=chrCount,flag=1) + for j, item in enumerate(plinkresultList): + try : + mb=float(item[1])/1000000.0 + except: + mb=0 + + try : + pvalue =float(item[-1]) + except: + pvalue =1 + + try: + snpName = item[0] + except: + snpName='' + + formattedPvalue = -math.log10(pvalue) + + Xc = startPosX + (mb-startMb)*plotXScale + Yc = yZero - (formattedPvalue-pvalueMin)*pvalueHeightThresh/(tempPvalueMax - pvalueMin) + canvas.drawString("5", Xc-canvas.stringWidth("5",font=symbolFont)/2+1,Yc+2,color=chrColorDict[i], font=symbolFont) + m += 1 + + startPosX += (ChrLengthDistList[i]+GraphInterval)*plotXScale + + canvas.drawLine(xLeftOffset, yZero, xLeftOffset, yTopOffset, color=self.LRS_COLOR, width=1*zoom) #the blue line running up the y axis + + def drawQTL(self, canvas, drawAreaHeight, gifmap, plotXScale, offset= (40, 120, 80, 10), zoom = 1, startMb = None, endMb = None): + + xLeftOffset, xRightOffset, yTopOffset, yBottomOffset = offset + plotWidth = canvas.size[0] - xLeftOffset - xRightOffset + plotHeight = canvas.size[1] - yTopOffset - yBottomOffset + fontZoom = zoom + if zoom == 2: + fontZoom = 1.5 + + INTERCROSS = (self.genotype.type=="intercross") + + ChrLengthDistList = self.ChrLengthMbList + GraphInterval=self.GraphInterval + LRSHeightThresh = drawAreaHeight + AdditiveHeightThresh = drawAreaHeight/2 + DominanceHeightThresh = drawAreaHeight/2 + + #draw the LRS scale + #We first determine whether or not we are using a sliding scale. + #If so, we need to compute the maximum LRS value to determine where the max y-value should be, and call this LRSMax. + #LRSTop is then defined to be above the LRSMax by enough to add one additional LRSScale increment. + #if we are using a set-scale, then we set LRSTop to be the user's value, and LRSMax doesn't matter. + + if self.LRS_LOD == 'LOD': + lodm = self.LODFACTOR + else: + lodm = 1.0 + + if self.lrsMax <= 0: #sliding scale + LRSMax = max(map(max, self.qtlresults)).lrs + #genotype trait will give infinite LRS + LRSMax = min(LRSMax, webqtlConfig.MAXLRS) + LRSMax = max(self.significance, LRSMax) + else: + LRSMax = self.lrsMax*lodm + + if LRSMax/lodm > 100: + LRSScale = 20.0 + elif LRSMax/lodm > 20: + LRSScale = 5.0 + elif LRSMax/lodm > 7.5: + LRSScale = 2.5 + else: + LRSScale = 1.0 + + LRSAxisList = Plot.frange(LRSScale, LRSMax/lodm, LRSScale) + #make sure the user's value appears on the y-axis + #update by NL 6-21-2011: round the LOD value to 100 when LRSMax is equal to 460 + LRSAxisList.append(round(LRSMax/lodm)) + + #draw the "LRS" or "LOD" string to the left of the axis + LRSScaleFont=pid.Font(ttf="verdana", size=14*fontZoom, bold=0) + LRSLODFont=pid.Font(ttf="verdana", size=14*zoom*1.5, bold=0) + yZero = yTopOffset + plotHeight + + #yAxis label display area + canvas.drawString(self.LRS_LOD, xLeftOffset - canvas.stringWidth("999.99", font=LRSScaleFont) - 10*zoom, \ + yZero - 150, font=LRSLODFont, color=pid.black, angle=90) + + for item in LRSAxisList: + yLRS = yZero - (item*lodm/LRSMax) * LRSHeightThresh + canvas.drawLine(xLeftOffset, yLRS, xLeftOffset - 4, yLRS, color=self.LRS_COLOR, width=1*zoom) + scaleStr = "%2.1f" % item + #added by NL 6-24-2011:Y-axis scale display + canvas.drawString(scaleStr, xLeftOffset-4-canvas.stringWidth(scaleStr, font=LRSScaleFont)-5, yLRS+3, font=LRSScaleFont, color=self.LRS_COLOR) + + + #"Significant" and "Suggestive" Drawing Routine + # ======= Draw the thick lines for "Significant" and "Suggestive" ===== (crowell: I tried to make the SNPs draw over these lines, but piddle wouldn't have it...) + if self.permChecked and not self.multipleInterval: + significantY = yZero - self.significance*LRSHeightThresh/LRSMax + suggestiveY = yZero - self.suggestive*LRSHeightThresh/LRSMax + + + startPosX = xLeftOffset + for i, _chr in enumerate(self.genotype): + rightEdge = int(startPosX + self.ChrLengthDistList[i]*plotXScale - self.SUGGESTIVE_WIDTH/1.5) + #added by NL 6-24-2011:draw suggestive line (grey one) + canvas.drawLine(startPosX+self.SUGGESTIVE_WIDTH/1.5, suggestiveY, rightEdge, suggestiveY, color=self.SUGGESTIVE_COLOR, + width=self.SUGGESTIVE_WIDTH*zoom, clipX=(xLeftOffset, xLeftOffset + plotWidth-2)) + #added by NL 6-24-2011:draw significant line (pink one) + canvas.drawLine(startPosX+self.SUGGESTIVE_WIDTH/1.5, significantY, rightEdge, significantY, color=self.SIGNIFICANT_COLOR, + width=self.SIGNIFICANT_WIDTH*zoom, clipX=(xLeftOffset, xLeftOffset + plotWidth-2)) + sugg_coords = "%d, %d, %d, %d" % (startPosX, suggestiveY-2, rightEdge + 2*zoom, suggestiveY+2) + sig_coords = "%d, %d, %d, %d" % (startPosX, significantY-2, rightEdge + 2*zoom, significantY+2) + if self.LRS_LOD == 'LRS': + sugg_title = "Suggestive LRS = %0.2f" % self.suggestive + sig_title = "Significant LRS = %0.2f" % self.significance + else: + sugg_title = "Suggestive LOD = %0.2f" % (self.suggestive/4.61) + sig_title = "Significant LOD = %0.2f" % (self.significance/4.61) + Areas1 = HT.Area(shape='rect',coords=sugg_coords,title=sugg_title) + Areas2 = HT.Area(shape='rect',coords=sig_coords,title=sig_title) + gifmap.areas.append(Areas1) + gifmap.areas.append(Areas2) + + startPosX += (self.ChrLengthDistList[i]+self.GraphInterval)*plotXScale + + + if self.multipleInterval: + lrsEdgeWidth = 1 + else: + additiveMax = max(map(lambda X : abs(X.additive), self.qtlresults[0])) + if INTERCROSS: + dominanceMax = max(map(lambda X : abs(X.dominance), self.qtlresults[0])) + else: + dominanceMax = -1 + lrsEdgeWidth = 2 + for i, qtlresult in enumerate(self.qtlresults): + m = 0 + startPosX = xLeftOffset + thisLRSColor = self.colorCollection[i] + + #add by NL 06-24-2011: for mahanttam plot + symbolFont = pid.Font(ttf="fnt_bs", size=5,bold=0) + + for j, _chr in enumerate(self.genotype): + chrCount=len(self.genotype) + chrColorDict =self.getColorForMarker(chrCount=chrCount,flag=1) + LRSCoordXY = [] + AdditiveCoordXY = [] + DominanceCoordXY = [] + for k, _locus in enumerate(_chr): + if self.plotScale == 'physic': + Xc = startPosX + (_locus.Mb-startMb)*plotXScale + else: + Xc = startPosX + (_locus.cM-_chr[0].cM)*plotXScale + # updated by NL 06-18-2011: + # fix the over limit LRS graph issue since genotype trait may give infinite LRS; + # for any lrs is over than 460(LRS max in this system), it will be reset to 460 + if qtlresult[m].lrs> 460 or qtlresult[m].lrs=='inf': + Yc = yZero - webqtlConfig.MAXLRS*LRSHeightThresh/LRSMax + else: + Yc = yZero - qtlresult[m].lrs*LRSHeightThresh/LRSMax + + LRSCoordXY.append((Xc, Yc)) + #add by NL 06-24-2011: for mahanttam plot + #self.significance/4.61 consider chr and LOD + # significantY = yZero - self.significance*LRSHeightThresh/LRSMax + # if Yc >significantY: + # canvas.drawString(":", Xc-canvas.stringWidth(":",font=symbolFont)/2+1,Yc+2,color=pid.black, font=symbolFont) + # else: + # canvas.drawString(":", Xc-canvas.stringWidth(":",font=symbolFont)/2+1,Yc+2,color=pid.black, font=symbolFont) + + # add by NL 06-27-2011: eliminate imputed value when locus name is equal to '-' + if (qtlresult[m].locus.name) and (qtlresult[m].locus.name!=' - '): + canvas.drawString("5", Xc-canvas.stringWidth("5",font=symbolFont)/2+1,Yc+2,color=chrColorDict[j], font=symbolFont) + + if not self.multipleInterval and self.additiveChecked: + Yc = yZero - qtlresult[m].additive*AdditiveHeightThresh/additiveMax + AdditiveCoordXY.append((Xc, Yc)) + if not self.multipleInterval and INTERCROSS and self.additiveChecked: + Yc = yZero - qtlresult[m].dominance*DominanceHeightThresh/dominanceMax + DominanceCoordXY.append((Xc, Yc)) + m += 1 + + startPosX += (ChrLengthDistList[j]+GraphInterval)*plotXScale + + + ###draw additive scale + if not self.multipleInterval and self.additiveChecked: + additiveScaleFont=pid.Font(ttf="verdana",size=12*fontZoom,bold=0) + additiveScale = Plot.detScaleOld(0,additiveMax) + additiveStep = (additiveScale[1]-additiveScale[0])/additiveScale[2] + additiveAxisList = Plot.frange(0, additiveScale[1], additiveStep) + maxAdd = additiveScale[1] + addPlotScale = AdditiveHeightThresh/additiveMax + + additiveAxisList.append(additiveScale[1]) + for item in additiveAxisList: + additiveY = yZero - item*addPlotScale + canvas.drawLine(xLeftOffset + plotWidth,additiveY,xLeftOffset+4+ plotWidth,additiveY,color=self.ADDITIVE_COLOR_POSITIVE, width=1*zoom) + scaleStr = "%2.3f" % item + canvas.drawString(scaleStr,xLeftOffset + plotWidth +6,additiveY+5,font=additiveScaleFont,color=self.ADDITIVE_COLOR_POSITIVE) + + canvas.drawLine(xLeftOffset+plotWidth,additiveY,xLeftOffset+plotWidth,yZero,color=self.ADDITIVE_COLOR_POSITIVE, width=1*zoom) + + canvas.drawLine(xLeftOffset, yZero, xLeftOffset, yTopOffset, color=self.LRS_COLOR, width=1*zoom) #the blue line running up the y axis + + def drawGraphBackgroundForPLINK(self, canvas, gifmap, offset= (80, 120, 80, 50), zoom = 1, startMb = None, endMb = None,plinkResultDict={} ): + + xLeftOffset, xRightOffset, yTopOffset, yBottomOffset = offset + plotWidth = canvas.size[0] - xLeftOffset - xRightOffset + plotHeight = canvas.size[1] - yTopOffset - yBottomOffset + fontZoom = zoom + if zoom == 2: + fontZoom = 1.5 + + #calculate plot scale + #XZ: all of these global variables should be passed from function signiture + ChrLengthDistList = self.ChrLengthMbList + drawRegionDistance = self.ChrLengthMbSum + GraphInterval=self.GraphInterval + ChrList =self.ChrList + + #multiple chromosome view + plotXScale = plotWidth / ((len(ChrList)-1)*GraphInterval + drawRegionDistance) + + startPosX = xLeftOffset + chrLabelFont=pid.Font(ttf="verdana",size=24*fontZoom,bold=0) + + for i, _chr in enumerate(ChrList): + + if (i % 2 == 0): + theBackColor = self.GRAPH_BACK_DARK_COLOR + else: + theBackColor = self.GRAPH_BACK_LIGHT_COLOR + # NL:resize chr width for drawing + if float(ChrLengthDistList[i])<90: + ChrLengthDistList[i]=90 + #draw the shaded boxes and the sig/sug thick lines + canvas.drawRect(startPosX, yTopOffset, startPosX + ChrLengthDistList[i]*plotXScale, \ + yTopOffset+plotHeight, edgeColor=pid.gainsboro,fillColor=theBackColor) + + chrNameWidth = canvas.stringWidth(_chr, font=chrLabelFont) + chrStartPix = startPosX + (ChrLengthDistList[i]*plotXScale -chrNameWidth)/2 + chrEndPix = startPosX + (ChrLengthDistList[i]*plotXScale +chrNameWidth)/2 + + canvas.drawString(_chr, chrStartPix, yTopOffset +20,font = chrLabelFont,color=pid.dimgray) + COORDS = "%d,%d,%d,%d" %(chrStartPix, yTopOffset, chrEndPix,yTopOffset +20) + + #add by NL 09-03-2010 + HREF = "javascript:changeView(%d,%s);" % (i,ChrLengthDistList) + Areas = HT.Area(shape='rect',coords=COORDS,href=HREF) + gifmap.areas.append(Areas) + startPosX += (ChrLengthDistList[i]+GraphInterval)*plotXScale + + return plotXScale + + + def drawGraphBackground(self, canvas, gifmap, offset= (80, 120, 80, 50), zoom = 1, startMb = None, endMb = None): + ##conditions + ##multiple Chromosome view + ##single Chromosome Physical + ##single Chromosome Genetic + xLeftOffset, xRightOffset, yTopOffset, yBottomOffset = offset + plotWidth = canvas.size[0] - xLeftOffset - xRightOffset + plotHeight = canvas.size[1] - yTopOffset - yBottomOffset + fontZoom = zoom + if zoom == 2: + fontZoom = 1.5 + + #calculate plot scale + if self.plotScale != 'physic': + self.ChrLengthDistList = self.ChrLengthCMList + drawRegionDistance = self.ChrLengthCMSum + else: + self.ChrLengthDistList = self.ChrLengthMbList + drawRegionDistance = self.ChrLengthMbSum + + if self.selectedChr > -1: #single chromosome view + spacingAmt = plotWidth/13.5 + i = 0 + for startPix in Plot.frange(xLeftOffset, xLeftOffset+plotWidth, spacingAmt): + if (i % 2 == 0): + theBackColor = self.GRAPH_BACK_DARK_COLOR + else: + theBackColor = self.GRAPH_BACK_LIGHT_COLOR + i += 1 + canvas.drawRect(startPix, yTopOffset, min(startPix+spacingAmt, xLeftOffset+plotWidth), \ + yTopOffset+plotHeight, edgeColor=theBackColor, fillColor=theBackColor) + + drawRegionDistance = self.ChrLengthDistList[self.selectedChr] + self.ChrLengthDistList = [drawRegionDistance] + if self.plotScale == 'physic': + plotXScale = plotWidth / (endMb-startMb) + else: + plotXScale = plotWidth / drawRegionDistance + + else: #multiple chromosome view + plotXScale = plotWidth / ((len(self.genotype)-1)*self.GraphInterval + drawRegionDistance) + + startPosX = xLeftOffset + chrLabelFont=pid.Font(ttf="verdana",size=24*fontZoom,bold=0) + + for i, _chr in enumerate(self.genotype): + + if (i % 2 == 0): + theBackColor = self.GRAPH_BACK_DARK_COLOR + else: + theBackColor = self.GRAPH_BACK_LIGHT_COLOR + + #draw the shaded boxes and the sig/sug thick lines + canvas.drawRect(startPosX, yTopOffset, startPosX + self.ChrLengthDistList[i]*plotXScale, \ + yTopOffset+plotHeight, edgeColor=pid.gainsboro,fillColor=theBackColor) + + chrNameWidth = canvas.stringWidth(_chr.name, font=chrLabelFont) + chrStartPix = startPosX + (self.ChrLengthDistList[i]*plotXScale -chrNameWidth)/2 + chrEndPix = startPosX + (self.ChrLengthDistList[i]*plotXScale +chrNameWidth)/2 + + canvas.drawString(_chr.name, chrStartPix, yTopOffset +20,font = chrLabelFont,color=pid.dimgray) + COORDS = "%d,%d,%d,%d" %(chrStartPix, yTopOffset, chrEndPix,yTopOffset +20) + + #add by NL 09-03-2010 + HREF = "javascript:changeView(%d,%s);" % (i,self.ChrLengthMbList) + Areas = HT.Area(shape='rect',coords=COORDS,href=HREF) + gifmap.areas.append(Areas) + startPosX += (self.ChrLengthDistList[i]+self.GraphInterval)*plotXScale + + return plotXScale + + # XZ: The only difference of function drawXAxisForPLINK and function drawXAxis are the function name and the self.plotScale condition. + def drawXAxisForPLINK(self, fd, canvas, drawAreaHeight, gifmap, plotXScale, showLocusForm, offset= (40, 120, 80, 10), zoom = 1, startMb = None, endMb = None): + xLeftOffset, xRightOffset, yTopOffset, yBottomOffset = offset + plotWidth = canvas.size[0] - xLeftOffset - xRightOffset + plotHeight = canvas.size[1] - yTopOffset - yBottomOffset + yZero = canvas.size[1] - yBottomOffset + fontZoom = zoom + if zoom == 2: + fontZoom = 1.5 + + #Parameters + ChrLengthDistList = self.ChrLengthMbList + GraphInterval=self.GraphInterval + + NUM_MINOR_TICKS = 5 # Number of minor ticks between major ticks + X_MAJOR_TICK_THICKNESS = 2 + X_MINOR_TICK_THICKNESS = 1 + X_AXIS_THICKNESS = 1*zoom + + # ======= Alex: Draw the X-axis labels (megabase location) + MBLabelFont = pid.Font(ttf="verdana", size=12*fontZoom, bold=0) + xMajorTickHeight = 15 # How high the tick extends below the axis + xMinorTickHeight = 5*zoom + xAxisTickMarkColor = pid.black + xAxisLabelColor = pid.black + fontHeight = 12*fontZoom # How tall the font that we're using is + spacingFromLabelToAxis = 20 + spacingFromLineToLabel = 3 + + if self.plotScale == 'physic': + strYLoc = yZero + spacingFromLabelToAxis + canvas.fontHeight(MBLabelFont) + ###Physical single chromosome view + if self.selectedChr > -1: + graphMbWidth = endMb - startMb + XScale = Plot.detScale(startMb, endMb) + XStart, XEnd, XStep = XScale + if XStep < 8: + XStep *= 2 + spacingAmtX = spacingAmt = (XEnd-XStart)/XStep + + j = 0 + while abs(spacingAmtX -int(spacingAmtX)) >= spacingAmtX/100.0 and j < 6: + j += 1 + spacingAmtX *= 10 + + formatStr = '%%2.%df' % j + + for counter, _Mb in enumerate(Plot.frange(XStart, XEnd, spacingAmt / NUM_MINOR_TICKS)): + if _Mb < startMb or _Mb > endMb: + continue + Xc = xLeftOffset + plotXScale*(_Mb - startMb) + if counter % NUM_MINOR_TICKS == 0: # Draw a MAJOR mark, not just a minor tick mark + canvas.drawLine(Xc, yZero, Xc, yZero+xMajorTickHeight, color=xAxisTickMarkColor, width=X_MAJOR_TICK_THICKNESS) # Draw the MAJOR tick mark + labelStr = str(formatStr % _Mb) # What Mbase location to put on the label + strWidth = canvas.stringWidth(labelStr, font=MBLabelFont) + drawStringXc = (Xc - (strWidth / 2.0)) + canvas.drawString(labelStr, drawStringXc, strYLoc, font=MBLabelFont, color=xAxisLabelColor, angle=0) + else: + canvas.drawLine(Xc, yZero, Xc, yZero+xMinorTickHeight, color=xAxisTickMarkColor, width=X_MINOR_TICK_THICKNESS) # Draw the MINOR tick mark + # end else + + ###Physical genome wide view + else: + distScale = 0 + startPosX = xLeftOffset + for i, distLen in enumerate(ChrLengthDistList): + if distScale == 0: #universal scale in whole genome mapping + if distLen > 75: + distScale = 25 + elif distLen > 30: + distScale = 10 + else: + distScale = 5 + for tickdists in range(distScale, ceil(distLen), distScale): + canvas.drawLine(startPosX + tickdists*plotXScale, yZero, startPosX + tickdists*plotXScale, yZero + 7, color=pid.black, width=1*zoom) + canvas.drawString(str(tickdists), startPosX+tickdists*plotXScale, yZero + 10*zoom, color=pid.black, font=MBLabelFont, angle=270) + startPosX += (ChrLengthDistList[i]+GraphInterval)*plotXScale + + megabaseLabelFont = pid.Font(ttf="verdana", size=14*zoom*1.5, bold=0) + canvas.drawString("Megabases", xLeftOffset + (plotWidth -canvas.stringWidth("Megabases", font=megabaseLabelFont))/2, + strYLoc + canvas.fontHeight(MBLabelFont) + 5*zoom, font=megabaseLabelFont, color=pid.black) + pass + + canvas.drawLine(xLeftOffset, yZero, xLeftOffset+plotWidth, yZero, color=pid.black, width=X_AXIS_THICKNESS) # Draw the X axis itself + + def drawXAxis(self, fd, canvas, drawAreaHeight, gifmap, plotXScale, showLocusForm, offset= (40, 120, 80, 10), zoom = 1, startMb = None, endMb = None): + xLeftOffset, xRightOffset, yTopOffset, yBottomOffset = offset + plotWidth = canvas.size[0] - xLeftOffset - xRightOffset + plotHeight = canvas.size[1] - yTopOffset - yBottomOffset + yZero = canvas.size[1] - yBottomOffset + fontZoom = zoom + if zoom == 2: + fontZoom = 1.5 + + #Parameters + NUM_MINOR_TICKS = 5 # Number of minor ticks between major ticks + X_MAJOR_TICK_THICKNESS = 2 + X_MINOR_TICK_THICKNESS = 1 + X_AXIS_THICKNESS = 1*zoom + + # ======= Alex: Draw the X-axis labels (megabase location) + MBLabelFont = pid.Font(ttf="verdana", size=12*fontZoom, bold=0) + xMajorTickHeight = 15 # How high the tick extends below the axis + xMinorTickHeight = 5*zoom + xAxisTickMarkColor = pid.black + xAxisLabelColor = pid.black + fontHeight = 12*fontZoom # How tall the font that we're using is + spacingFromLabelToAxis = 20 + spacingFromLineToLabel = 3 + + if self.plotScale == 'physic': + strYLoc = yZero + spacingFromLabelToAxis + canvas.fontHeight(MBLabelFont) + ###Physical single chromosome view + if self.selectedChr > -1: + graphMbWidth = endMb - startMb + XScale = Plot.detScale(startMb, endMb) + XStart, XEnd, XStep = XScale + if XStep < 8: + XStep *= 2 + spacingAmtX = spacingAmt = (XEnd-XStart)/XStep + + j = 0 + while abs(spacingAmtX -int(spacingAmtX)) >= spacingAmtX/100.0 and j < 6: + j += 1 + spacingAmtX *= 10 + + formatStr = '%%2.%df' % j + + for counter, _Mb in enumerate(Plot.frange(XStart, XEnd, spacingAmt / NUM_MINOR_TICKS)): + if _Mb < startMb or _Mb > endMb: + continue + Xc = xLeftOffset + plotXScale*(_Mb - startMb) + if counter % NUM_MINOR_TICKS == 0: # Draw a MAJOR mark, not just a minor tick mark + canvas.drawLine(Xc, yZero, Xc, yZero+xMajorTickHeight, color=xAxisTickMarkColor, width=X_MAJOR_TICK_THICKNESS) # Draw the MAJOR tick mark + labelStr = str(formatStr % _Mb) # What Mbase location to put on the label + strWidth = canvas.stringWidth(labelStr, font=MBLabelFont) + drawStringXc = (Xc - (strWidth / 2.0)) + canvas.drawString(labelStr, drawStringXc, strYLoc, font=MBLabelFont, color=xAxisLabelColor, angle=0) + else: + canvas.drawLine(Xc, yZero, Xc, yZero+xMinorTickHeight, color=xAxisTickMarkColor, width=X_MINOR_TICK_THICKNESS) # Draw the MINOR tick mark + # end else + + ###Physical genome wide view + else: + distScale = 0 + startPosX = xLeftOffset + for i, distLen in enumerate(self.ChrLengthDistList): + if distScale == 0: #universal scale in whole genome mapping + if distLen > 75: + distScale = 25 + elif distLen > 30: + distScale = 10 + else: + distScale = 5 + for tickdists in range(distScale, ceil(distLen), distScale): + canvas.drawLine(startPosX + tickdists*plotXScale, yZero, startPosX + tickdists*plotXScale, yZero + 7, color=pid.black, width=1*zoom) + canvas.drawString(str(tickdists), startPosX+tickdists*plotXScale, yZero + 10*zoom, color=pid.black, font=MBLabelFont, angle=270) + startPosX += (self.ChrLengthDistList[i]+self.GraphInterval)*plotXScale + + megabaseLabelFont = pid.Font(ttf="verdana", size=14*zoom*1.5, bold=0) + canvas.drawString("Megabases", xLeftOffset + (plotWidth -canvas.stringWidth("Megabases", font=megabaseLabelFont))/2, + strYLoc + canvas.fontHeight(MBLabelFont) + 5*zoom, font=megabaseLabelFont, color=pid.black) + pass + else: + ChrAInfo = [] + preLpos = -1 + distinctCount = 0.0 + if len(self.genotype) > 1: + for i, _chr in enumerate(self.genotype): + thisChr = [] + Locus0CM = _chr[0].cM + nLoci = len(_chr) + if nLoci <= 8: + for _locus in _chr: + if _locus.name != ' - ': + if _locus.cM != preLpos: + distinctCount += 1 + preLpos = _locus.cM + thisChr.append([_locus.name, _locus.cM-Locus0CM]) + else: + for j in (0, nLoci/4, nLoci/2, nLoci*3/4, -1): + while _chr[j].name == ' - ': + j += 1 + if _chr[j].cM != preLpos: + distinctCount += 1 + preLpos = _chr[j].cM + thisChr.append([_chr[j].name, _chr[j].cM-Locus0CM]) + ChrAInfo.append(thisChr) + else: + for i, _chr in enumerate(self.genotype): + thisChr = [] + Locus0CM = _chr[0].cM + for _locus in _chr: + if _locus.name != ' - ': + if _locus.cM != preLpos: + distinctCount += 1 + preLpos = _locus.cM + thisChr.append([_locus.name, _locus.cM-Locus0CM]) + ChrAInfo.append(thisChr) + + stepA = (plotWidth+0.0)/distinctCount + + LRectWidth = 10 + LRectHeight = 3 + offsetA = -stepA + lineColor = pid.lightblue + startPosX = xLeftOffset + for j, ChrInfo in enumerate(ChrAInfo): + preLpos = -1 + for i, item in enumerate(ChrInfo): + Lname,Lpos = item + if Lpos != preLpos: + offsetA += stepA + differ = 1 + else: + differ = 0 + preLpos = Lpos + Lpos *= plotXScale + if self.selectedChr > -1: + Zorder = i % 5 + else: + Zorder = 0 + if differ: + canvas.drawLine(startPosX+Lpos,yZero,xLeftOffset+offsetA,\ + yZero+25, color=lineColor) + canvas.drawLine(xLeftOffset+offsetA,yZero+25,xLeftOffset+offsetA,\ + yZero+40+Zorder*(LRectWidth+3),color=lineColor) + rectColor = pid.orange + else: + canvas.drawLine(xLeftOffset+offsetA, yZero+40+Zorder*(LRectWidth+3)-3,\ + xLeftOffset+offsetA, yZero+40+Zorder*(LRectWidth+3),color=lineColor) + rectColor = pid.deeppink + canvas.drawRect(xLeftOffset+offsetA, yZero+40+Zorder*(LRectWidth+3),\ + xLeftOffset+offsetA-LRectHeight,yZero+40+Zorder*(LRectWidth+3)+LRectWidth,\ + edgeColor=rectColor,fillColor=rectColor,edgeWidth = 0) + COORDS="%d,%d,%d,%d"%(xLeftOffset+offsetA-LRectHeight, yZero+40+Zorder*(LRectWidth+3),\ + xLeftOffset+offsetA,yZero+40+Zorder*(LRectWidth+3)+LRectWidth) + HREF="javascript:showDatabase3('%s','%s','%s','');" % (showLocusForm,fd.RISet+"Geno", Lname) + Areas=HT.Area(shape='rect',coords=COORDS,href=HREF, title="Locus : " + Lname) + gifmap.areas.append(Areas) + ##piddle bug + if j == 0: + canvas.drawLine(startPosX,yZero,startPosX,yZero+40, color=lineColor) + startPosX += (self.ChrLengthDistList[j]+self.GraphInterval)*plotXScale + + canvas.drawLine(xLeftOffset, yZero, xLeftOffset+plotWidth, yZero, color=pid.black, width=X_AXIS_THICKNESS) # Draw the X axis itself + + def getColorForMarker(self, chrCount,flag):# no change is needed + chrColorDict={} + for i in range(chrCount): + if flag==1: # display blue and lightblue intercross + chrColorDict[i]=pid.black + elif flag==0: + if (i%2==0): + chrColorDict[i]=pid.blue + else: + chrColorDict[i]=pid.lightblue + else:#display different color for different chr + if i in [0,8,16]: + chrColorDict[i]=pid.black + elif i in [1,9,17]: + chrColorDict[i]=pid.red + elif i in [2,10,18]: + chrColorDict[i]=pid.lightgreen + elif i in [3,11,19]: + chrColorDict[i]=pid.blue + elif i in [4,12]: + chrColorDict[i]=pid.lightblue + elif i in [5,13]: + chrColorDict[i]=pid.hotpink + elif i in [6,14]: + chrColorDict[i]=pid.gold + elif i in [7,15]: + chrColorDict[i]=pid.grey + + return chrColorDict + + + def drawProbeSetPosition(self, canvas, plotXScale, offset= (40, 120, 80, 10), zoom = 1, startMb = None, endMb = None): + if len(self.traitList) != 1: + return + + xLeftOffset, xRightOffset, yTopOffset, yBottomOffset = offset + plotWidth = canvas.size[0] - xLeftOffset - xRightOffset + plotHeight = canvas.size[1] - yTopOffset - yBottomOffset + yZero = canvas.size[1] - yBottomOffset + fontZoom = zoom + if zoom == 2: + fontZoom = 1.5 + + try: + Chr = self.traitList[0].chr # self.traitListChr =self.traitList[0].chr=_vals need to change to chrList and mbList + Mb = self.traitList[0].mb # self.traitListMb =self.traitList[0].mb=_vals + except: + return + + if self.plotScale == 'physic': + if self.selectedChr > -1: + if self.genotype[0].name != Chr or Mb < self.startMb or Mb > self.endMb: + return + else: + locPixel = xLeftOffset + (Mb-self.startMb)*plotXScale + else: + locPixel = xLeftOffset + for i, _chr in enumerate(self.genotype): + if _chr.name != Chr: + locPixel += (self.ChrLengthDistList[i] + self.GraphInterval)*plotXScale + else: + locPixel += Mb*plotXScale + break + else: + if self.selectedChr > -1: + if self.genotype[0].name != Chr: + return + else: + for i, _locus in enumerate(self.genotype[0]): + #the trait's position is on the left of the first genotype + if i==0 and _locus.Mb >= Mb: + locPixel=-1 + break + + #the trait's position is between two traits + if i > 0 and self.genotype[0][i-1].Mb < Mb and _locus.Mb >= Mb: + locPixel = xLeftOffset + plotXScale*(self.genotype[0][i-1].cM+(_locus.cM-self.genotype[0][i-1].cM)*(Mb -self.genotype[0][i-1].Mb)/(_locus.Mb-self.genotype[0][i-1].Mb)) + break + + #the trait's position is on the right of the last genotype + if i==len(self.genotype[0]) and Mb>=_locus.Mb: + locPixel = -1 + else: + locPixel = xLeftOffset + for i, _chr in enumerate(self.genotype): + if _chr.name != Chr: + locPixel += (self.ChrLengthDistList[i] + self.GraphInterval)*plotXScale + else: + locPixel += (Mb*(_chr[-1].cM-_chr[0].cM)/self.ChrLengthCMList[i])*plotXScale + break + if locPixel >= 0: + traitPixel = ((locPixel, yZero), (locPixel-6, yZero+12), (locPixel+6, yZero+12)) + canvas.drawPolygon(traitPixel, edgeColor=pid.black, fillColor=self.TRANSCRIPT_LOCATION_COLOR, closed=1) + + if self.legendChecked: + startPosY = 15 + nCol = 2 + smallLabelFont = pid.Font(ttf="trebuc", size=12, bold=1) + leftOffset = xLeftOffset+(nCol-1)*200 + canvas.drawPolygon(((leftOffset+6, startPosY-6), (leftOffset, startPosY+6), (leftOffset+12, startPosY+6)), edgeColor=pid.black, fillColor=self.TRANSCRIPT_LOCATION_COLOR, closed=1) + canvas.drawString("Sequence Site", (leftOffset+15), (startPosY+5), smallLabelFont, self.TOP_RIGHT_INFO_COLOR) + + # build dict based on plink result, key is chr, value is list of [snp,BP,pValue] + def getPlinkResultDict(self,outputFileName='',thresholdPvalue=-1,ChrOrderIdNameDict={}): + + ChrList =self.ChrList + plinkResultDict={} + + plinkResultfp = open("%s%s.qassoc"% (webqtlConfig.TMPDIR, outputFileName), "rb") + + headerLine=plinkResultfp.readline()# read header line + line = plinkResultfp.readline() + + valueList=[] # initialize value list, this list will include snp, bp and pvalue info + pValueList=[] + count=0 + + while line: + #convert line from str to list + lineList=self.buildLineList(line=line) + + # only keep the records whose chromosome name is in db + if ChrOrderIdNameDict.has_key(int(lineList[0])) and lineList[-1] and lineList[-1].strip()!='NA': + + chrName=ChrOrderIdNameDict[int(lineList[0])] + snp = lineList[1] + BP = lineList[2] + pValue = float(lineList[-1]) + pValueList.append(pValue) + + if plinkResultDict.has_key(chrName): + valueList=plinkResultDict[chrName] + + # pvalue range is [0,1] + if thresholdPvalue >=0 and thresholdPvalue<=1: + if pValue < thresholdPvalue: + valueList.append((snp,BP,pValue)) + count+=1 + + plinkResultDict[chrName]=valueList + valueList=[] + else: + if thresholdPvalue>=0 and thresholdPvalue<=1: + if pValue < thresholdPvalue: + valueList.append((snp,BP,pValue)) + count+=1 + + if valueList: + plinkResultDict[chrName]=valueList + + valueList=[] + + + line =plinkResultfp.readline() + else: + line=plinkResultfp.readline() + + if pValueList: + minPvalue= min(pValueList) + else: + minPvalue=0 + + return count,minPvalue,plinkResultDict + + + ###################################################### + # input: line: str,one line read from file + # function: convert line from str to list; + # output: lineList list + ####################################################### + def buildLineList(self,line=None): + + lineList = string.split(string.strip(line),' ')# irregular number of whitespaces between columns + lineList =[ item for item in lineList if item <>''] + lineList = map(string.strip, lineList) + + return lineList + + #added by NL: automatically generate pheno txt file for PLINK based on strainList passed from dataEditing page + def genPhenoTxtFileForPlink(self,phenoFileName='', RISetName='', probesetName='', valueDict={}): + pedFileStrainList=self.getStrainNameFromPedFile(RISetName=RISetName) + outputFile = open("%s%s.txt"%(webqtlConfig.TMPDIR,phenoFileName),"wb") + headerLine = 'FID\tIID\t%s\n'%probesetName + outputFile.write(headerLine) + + newValueList=[] + + #if valueDict does not include some strain, value will be set to -9999 as missing value + for item in pedFileStrainList: + try: + value=valueDict[item] + value=str(value).replace('value=','') + value=value.strip() + except: + value=-9999 + + newValueList.append(value) + + + newLine='' + for i, strain in enumerate(pedFileStrainList): + j=i+1 + value=newValueList[i] + newLine+='%s\t%s\t%s\n'%(strain, strain, value) + + if j%1000==0: + outputFile.write(newLine) + newLine='' + + if newLine: + outputFile.write(newLine) + + outputFile.close() + + # get strain name from ped file in order + def getStrainNameFromPedFile(self, RISetName=''): + pedFileopen= open("%splink/%s.ped"%(webqtlConfig.HTMLPATH, RISetName),"r") + line =pedFileopen.readline() + strainNameList=[] + + while line: + lineList=string.split(string.strip(line),'\t') + lineList=map(string.strip,lineList) + + strainName=lineList[0] + strainNameList.append(strainName) + + line =pedFileopen.readline() + + return strainNameList + + ################################################################# + ## Generate Chr list, Chr OrderId and Retrieve Length Information + ################################################################# + #def getChrNameOrderIdLength(self,RISet=''): + # + # try: + # query = """ + # Select + # Chr_Length.Name,Chr_Length.OrderId,Length from Chr_Length, InbredSet + # where + # Chr_Length.SpeciesId = InbredSet.SpeciesId AND + # InbredSet.Name = '%s' + # Order by OrderId + # """ % (RISet) + # self.cursor.execute(query) + # + # results =self.cursor.fetchall() + # ChrList=[] + # ChrLengthMbList=[] + # ChrNameOrderIdDict={} + # ChrOrderIdNameDict={} + # + # for item in results: + # ChrList.append(item[0]) + # ChrNameOrderIdDict[item[0]]=item[1] # key is chr name, value is orderId + # ChrOrderIdNameDict[item[1]]=item[0] # key is orderId, value is chr name + # ChrLengthMbList.append(item[2]) + # + # except: + # ChrList=[] + # ChrNameOrderIdDict={} + # ChrLengthMbList=[] + # + # return ChrList,ChrNameOrderIdDict,ChrOrderIdNameDict,ChrLengthMbList |