about summary refs log tree commit diff
path: root/web/webqtl/textUI/cmdInterval.py
diff options
context:
space:
mode:
authorroot2012-05-08 18:39:56 -0500
committerroot2012-05-08 18:39:56 -0500
commitea46f42ee640928b92947bfb204c41a482d80937 (patch)
tree9b27a4eb852d12539b543c3efee9d2a47ef470f3 /web/webqtl/textUI/cmdInterval.py
parent056b5253fc3857b0444382aa39944f6344dc1ceb (diff)
downloadgenenetwork2-ea46f42ee640928b92947bfb204c41a482d80937.tar.gz
Add all the source codes into the github.
Diffstat (limited to 'web/webqtl/textUI/cmdInterval.py')
-rwxr-xr-xweb/webqtl/textUI/cmdInterval.py174
1 files changed, 174 insertions, 0 deletions
diff --git a/web/webqtl/textUI/cmdInterval.py b/web/webqtl/textUI/cmdInterval.py
new file mode 100755
index 00000000..0b97c7c3
--- /dev/null
+++ b/web/webqtl/textUI/cmdInterval.py
@@ -0,0 +1,174 @@
+# Copyright (C) University of Tennessee Health Science Center, Memphis, TN.
+#
+# This program is free software: you can redistribute it and/or modify it
+# under the terms of the GNU Affero General Public License
+# as published by the Free Software Foundation, either version 3 of the
+# License, or (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+# See the GNU Affero General Public License for more details.
+#
+# This program is available from Source Forge: at GeneNetwork Project
+# (sourceforge.net/projects/genenetwork/).
+#
+# Contact Drs. Robert W. Williams and Xiaodong Zhou (2010)
+# at rwilliams@uthsc.edu and xzhou15@uthsc.edu
+#
+#
+#
+# This module is used by GeneNetwork project (www.genenetwork.org)
+#
+# Created by GeneNetwork Core Team 2010/08/10
+#
+# Last updated by GeneNetwork Core Team 2010/10/20
+
+import string
+import os
+
+import reaper
+
+from base import webqtlConfig
+from cmdClass import cmdClass
+		
+#########################################
+#      Interval Mapping Class
+#########################################
+class cmdInterval(cmdClass):
+
+	def __init__(self,fd=None):
+
+		cmdClass.__init__(self,fd)
+
+		if not webqtlConfig.TEXTUI:
+			self.contents.append("Please send your request to http://robot.genenetwork.org")
+			return
+
+
+		self.example = '###Example : <a href="%s%s?cmd=%s&probeset=100001_at&probe=136415&db=bra03-03Mas5&sort=pos&return=100&chr=12">%s%s?cmd=%s&probeset=100001_at&probe=136415&db=bra03-03Mas5&sort=pos&return=100&chr=12</a>' % (webqtlConfig.CGIDIR, webqtlConfig.SCRIPTFILE, self.cmdID, webqtlConfig.CGIDIR, webqtlConfig.SCRIPTFILE, self.cmdID)
+		if self.accessError:
+			return
+		self.sort = None
+		self.step = 0.01
+		self.peak = 1
+		self.chr = None
+		self.sort = None
+		self.returnnumber = 20
+		if self.error:
+			self.error = 1
+			self.contents.append(self.example)
+			return
+		else:
+			try:
+				self.sort = self.data.getvalue('sort')
+				if string.lower(self.sort) == 'pos':
+					self.sort = 'pos'
+				else:
+					self.sort = 'lrs'
+			except:
+				self.sort = None
+			
+			try:
+				self.returnnumber = int(self.data.getvalue('return'))
+			except:
+				self.returnnumber = 20
+			try:	
+				self.chr = self.data.getvalue('chr')
+			except:
+				self.chr = None
+				
+			self.readDB()
+
+	def  readDB(self):
+		prefix, dbId = self.getDBId(self.database)
+		if not prefix or not dbId or (self.probe and string.lower(self.probe) in ("all","mm","pm")):
+			self.contents.append("###Error: source trait doesn't exist or SELECT more than one trait.")
+			self.contents.append(self.example)
+			self.contents.append(self.accessCode)
+			return
+		RISet = self.getRISet(prefix, dbId)
+		traitdata, heads = self.getTraitData(prefix, dbId, self.probeset, self.probe)
+		if not traitdata:
+			self.contents.append("###Error: source trait doesn't exist or SELECT more than one trait.")
+			self.contents.append(self.example)
+			self.contents.append(self.accessCode)
+			return
+			
+		dataset0 = reaper.Dataset()
+		dataset0.read(os.path.join(webqtlConfig.GENODIR, RISet + '.geno'))
+		strainList = list(dataset0.prgy)
+		dataset = dataset0.addinterval()
+		if self.chr != None:
+			for _chr in dataset:
+				if string.lower(_chr.name) ==  string.lower(self.chr):
+					dataset.chromosome = [_chr]
+					break
+		
+		strains = []
+		trait = []
+		_prgy = dataset.prgy
+		for item in traitdata:
+			if item[0] in _prgy:
+				strains.append(item[0])	
+				trait.append(item[1])
+
+		qtlscan = dataset.regression(strains, trait)
+		LRS = dataset.permutation(strains, trait)
+		nperm = len(LRS)
+		
+		#print inter1[0]
+		returnPeak = []
+		nqtl = len(qtlscan)
+		if self.peak:
+			for i in range(nqtl):
+				if i == 0 or qtlscan[i].locus.chr != qtlscan[i-1].locus.chr:
+					if qtlscan[i].lrs < qtlscan[i+1].lrs:
+						continue
+				elif i == nqtl-1 or qtlscan[i].locus.chr != qtlscan[i+1].locus.chr:
+					if qtlscan[i].lrs < qtlscan[i-1].lrs:
+						continue
+				else:
+					if qtlscan[i].lrs < qtlscan[i+1].lrs or qtlscan[i].lrs < qtlscan[i-1].lrs:
+						continue
+				returnPeak.append(qtlscan[i])
+		else:
+			returnPeak = qtlscan[:]
+			
+		if returnPeak:
+			self.contents.append("Locus\tLRS\tChr\tAdditive\tp-value\tcM") 
+			qtlresult = []
+			for item in returnPeak:
+				p_value = reaper.pvalue(item.lrs,LRS)
+				qtlresult.append((item.locus.name,item.lrs,item.locus.chr,item.additive,p_value, item.locus.cM))
+			if self.sort == 'lrs':
+				qtlresult.sort(self.cmpLRS2)
+			for item in qtlresult:
+				self.contents.append("%s\t%2.4f\t%s\t%2.4f\t%1.4f\t%s" % item)
+		else:
+			self.contents.append("###Error: Error occurs while regression.")
+			return
+	
+	def cmpPValue(self,A,B):
+		try:
+			if A[-1] > B[-1]:
+				return 1
+			elif A[-1] ==  B[-1]:
+				return 0
+			else:
+				return -1	
+		except:
+			return 0
+
+	def cmpLRS2(self,A,B):
+		try:
+			if A[1] < B[1]:
+				return 1
+			elif A[1] == B[1]:
+				return 0
+			else:
+				return -1	
+		except:
+			return 0
+
+