Committing a bunch of changes related to integrating GEMMA and

adding the correlation matrix page

1 files changed, 97 insertions, 0 deletions

diff --git a/web/webqtl/markerRegression/MarkerRegressionPage.py b/web/webqtl/markerRegression/MarkerRegressionPage.py
index 7f830b4b..43988d29 100755
--- a/web/webqtl/markerRegression/MarkerRegressionPage.py
+++ b/web/webqtl/markerRegression/MarkerRegressionPage.py
@@ -151,6 +151,100 @@ class MarkerRegressionPage(templatePage):
 				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:
@@ -1459,6 +1553,9 @@ class MarkerRegressionPage(templatePage):
 			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={}):