# 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