path: root/web/webqtl/markerRegression/MarkerRegressionPage.py

# 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
			
		if self.mappingMethodId == '5': # For GEMMA
			
			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
			gemma_command = './%sgemma/gemma -bfile %s -k %s.kin -lmm 1 -o %s_output' % (webqtlConfig.HTMLPATH, fd.RISet, fd.RISet, fd.RISet)
			#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(gemma_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:]
		
		
		
		#########################################
		#      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('&nbsp;&nbsp;&nbsp;&nbsp;Suggestive LRS =%3.2f\n'%fd.suggestive, HT.BR(), '&nbsp;&nbsp;&nbsp;&nbsp;Significant LRS =%3.2f\n'%fd.significance,HT.BR(),'&nbsp;&nbsp;&nbsp;&nbsp;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