diff options
Diffstat (limited to 'wqflask/base')
| -rwxr-xr-x | wqflask/base/data_set.py | 212 | ||||
| -rwxr-xr-x | wqflask/base/trait.py | 101 | ||||
| -rwxr-xr-x | wqflask/base/webqtlConfig.py | 2 |
3 files changed, 184 insertions, 131 deletions
diff --git a/wqflask/base/data_set.py b/wqflask/base/data_set.py index 9b0a3dcc..07fe9cd9 100755 --- a/wqflask/base/data_set.py +++ b/wqflask/base/data_set.py @@ -38,6 +38,7 @@ from base import species from dbFunction import webqtlDatabaseFunction from utility import webqtlUtil from utility.benchmark import Bench +from wqflask.my_pylmm.pyLMM import chunks from MySQLdb import escape_string as escape from pprint import pformat as pf @@ -46,7 +47,7 @@ from pprint import pformat as pf DS_NAME_MAP = {} def create_dataset(dataset_name): - print("dataset_name:", dataset_name) + #print("dataset_name:", dataset_name) query = """ SELECT DBType.Name @@ -68,10 +69,17 @@ def create_dataset(dataset_name): dataset_class = globals()[dataset_ob] return dataset_class(dataset_name) +def create_in_clause(items): + """Create an in clause for mysql""" + in_clause = ', '.join("'{}'".format(x) for x in mescape(*items)) + in_clause = '( {} )'.format(in_clause) + return in_clause + + def mescape(*items): """Multiple escape""" - escaped = [escape(item) for item in items] - print("escaped is:", escaped) + escaped = [escape(str(item)) for item in items] + #print("escaped is:", escaped) return escaped @@ -82,8 +90,8 @@ class Markers(object): self.markers = json.load(json_data_fh) def add_pvalues(self, p_values): - print("length of self.markers:", len(self.markers)) - print("length of p_values:", len(p_values)) + #print("length of self.markers:", len(self.markers)) + #print("length of p_values:", len(p_values)) # THIS IS only needed for the case when we are limiting the number of p-values calculated if len(self.markers) < len(p_values): @@ -153,7 +161,7 @@ class DatasetGroup(object): self.f1list = None self.parlist = None self.get_f1_parent_strains() - print("parents/f1s: {}:{}".format(self.parlist, self.f1list)) + #print("parents/f1s: {}:{}".format(self.parlist, self.f1list)) self.species = webqtlDatabaseFunction.retrieve_species(self.name) @@ -162,7 +170,7 @@ class DatasetGroup(object): def get_markers(self): - print("self.species is:", self.species) + #print("self.species is:", self.species) if self.species == "human": marker_class = HumanMarkers else: @@ -235,6 +243,7 @@ class DataSet(object): self.retrieve_other_names() self.group = DatasetGroup(self) # sets self.group and self.group_id and gets genotype + self.group.read_genotype_file() self.species = species.TheSpecies(self) @@ -284,14 +293,14 @@ class DataSet(object): self.name, self.name, self.name)) - print("query_args are:", query_args) + #print("query_args are:", query_args) - print(""" - SELECT Id, Name, FullName, ShortName - FROM %s - WHERE public > %s AND - (Name = '%s' OR FullName = '%s' OR ShortName = '%s') - """ % (query_args)) + #print(""" + # SELECT Id, Name, FullName, ShortName + # FROM %s + # WHERE public > %s AND + # (Name = '%s' OR FullName = '%s' OR ShortName = '%s') + # """ % (query_args)) self.id, self.name, self.fullname, self.shortname = g.db.execute(""" SELECT Id, Name, FullName, ShortName @@ -615,34 +624,33 @@ class MrnaAssayDataSet(DataSet): and ProbeSetFreezeId = {} """.format(escape(str(self.id))) results = g.db.execute(query).fetchall() - print("After get_trait_list query") + #print("After get_trait_list query") trait_data = {} for trait in results: print("Retrieving sample_data for ", trait[0]) trait_data[trait[0]] = self.retrieve_sample_data(trait[0]) - print("After retrieve_sample_data") + #print("After retrieve_sample_data") return trait_data def get_trait_data(self): - sample_ids = [] - for sample in self.group.samplelist: - query = """ - SELECT Strain.Id FROM Strain, Species - WHERE Strain.Name = '{}' - and Strain.SpeciesId=Species.Id - and Species.name = '{}' - """.format(*mescape(sample, self.group.species)) - this_id = g.db.execute(query).fetchone()[0] - sample_ids.append('%d' % this_id) - print("sample_ids size: ", len(sample_ids)) + self.samplelist = self.group.samplelist + self.group.parlist + self.group.f1list + query = """ + SELECT Strain.Name, Strain.Id FROM Strain, Species + WHERE Strain.Name IN {} + and Strain.SpeciesId=Species.Id + and Species.name = '{}' + """.format(create_in_clause(self.samplelist), *mescape(self.group.species)) + results = dict(g.db.execute(query).fetchall()) + sample_ids = [results[item] for item in self.samplelist] # MySQL limits the number of tables that can be used in a join to 61, # so we break the sample ids into smaller chunks - chunk_count = 50 - n = len(sample_ids) / chunk_count - if len(sample_ids) % chunk_count: - n += 1 - print("n: ", n) + # Postgres doesn't have that limit, so we can get rid of this after we transition + chunk_size = 50 + number_chunks = int(math.ceil(len(sample_ids) / chunk_size)) + trait_sample_data = [] + for sample_ids_step in chunks.divide_into_chunks(sample_ids, number_chunks): + #XZ, 09/24/2008: build one temporary table that only contains the records associated with the input GeneId #tempTable = None #if GeneId and db.type == "ProbeSet": @@ -656,24 +664,21 @@ class MrnaAssayDataSet(DataSet): # TissueProbeSetFreezeId=tissueProbeSetFreezeId, # method=method, # returnNumber=returnNumber) - trait_sample_data = [] - for step in range(int(n)): - temp = [] - sample_ids_step = sample_ids[step*chunk_count:min(len(sample_ids), (step+1)*chunk_count)] - for item in sample_ids_step: - temp.append('T%s.value' % item) + + temp = ['T%s.value' % item for item in sample_ids_step] query = "SELECT {}.Name,".format(escape(self.type)) data_start_pos = 1 query += string.join(temp, ', ') query += ' FROM ({}, {}XRef, {}Freeze) '.format(*mescape(self.type, self.type, self.type)) - #XZ, 03/04/2009: Xiaodong changed Data to %sData and changed parameters from %(item,item, db.type,item,item) to %(db.type, item,item, db.type,item,item) + for item in sample_ids_step: query += """ left join {}Data as T{} on T{}.Id = {}XRef.DataId and T{}.StrainId={}\n """.format(*mescape(self.type, item, item, self.type, item, item)) + query += """ WHERE {}XRef.{}FreezeId = {}Freeze.Id and {}Freeze.Name = '{}' @@ -681,23 +686,24 @@ class MrnaAssayDataSet(DataSet): order by {}.Id """.format(*mescape(self.type, self.type, self.type, self.type, self.name, self.type, self.type, self.type, self.type)) - print("query: ", query) results = g.db.execute(query).fetchall() trait_sample_data.append(results) - + trait_count = len(trait_sample_data[0]) self.trait_data = collections.defaultdict(list) + # put all of the separate data together into a dictionary where the keys are # trait names and values are lists of sample values - for j in range(trait_count): - trait_name = trait_sample_data[0][j][0] - for i in range(int(n)): - self.trait_data[trait_name] += trait_sample_data[i][j][data_start_pos:] - + for trait_counter in range(trait_count): + trait_name = trait_sample_data[0][trait_counter][0] + for chunk_counter in range(int(number_chunks)): + self.trait_data[trait_name] += ( + trait_sample_data[chunk_counter][trait_counter][data_start_pos:]) + def get_trait_info(self, trait_list=None, species=''): - # Note: setting trait_list to [] is probably not a great idea. + # Note: setting trait_list to [] is probably not a great idea. if not trait_list: trait_list = [] @@ -706,9 +712,7 @@ class MrnaAssayDataSet(DataSet): if not this_trait.haveinfo: this_trait.retrieveInfo(QTL=1) - if this_trait.symbol: - pass - else: + if not this_trait.symbol: this_trait.symbol = "N/A" #XZ, 12/08/2008: description @@ -716,62 +720,56 @@ class MrnaAssayDataSet(DataSet): description_string = str(this_trait.description).strip() target_string = str(this_trait.probe_target_description).strip() - description_display = '' - if len(description_string) > 1 and description_string != 'None': description_display = description_string else: description_display = this_trait.symbol - if len(description_display) > 1 and description_display != 'N/A' and len(target_string) > 1 and target_string != 'None': + if (len(description_display) > 1 and description_display != 'N/A' and + len(target_string) > 1 and target_string != 'None'): description_display = description_display + '; ' + target_string.strip() # Save it for the jinja2 template this_trait.description_display = description_display - #print(" xxxxdd [%s]: %s" % (type(this_trait.description_display), description_display)) #XZ: trait_location_value is used for sorting trait_location_repr = 'N/A' trait_location_value = 1000000 if this_trait.chr and this_trait.mb: - print("this_trait.chr is:", this_trait.chr) - print("this_trait.mb is:", this_trait.mb) - try: - trait_location_value = float(this_trait.chr)*1000 + float(this_trait.mb) - except: - if this_trait.chr.upper() == 'X': - trait_location_value = 20*1000 + this_trait.mb - else: - trait_location_value = ord(str(this_trait.chr).upper()[0])*1000 + this_trait.mb - - this_trait.location_repr = 'Chr %s: %.4f Mb' % (this_trait.chr, float(this_trait.mb) ) + #Checks if the chromosome number can be cast to an int (i.e. isn't "X" or "Y") + #This is so we can convert the location to a number used for sorting + trait_location_value = self.convert_location_to_value(this_trait.chr, this_trait.mb) + #try: + # trait_location_value = int(this_trait.chr)*1000 + this_trait.mb + #except ValueError: + # if this_trait.chr.upper() == 'X': + # trait_location_value = 20*1000 + this_trait.mb + # else: + # trait_location_value = (ord(str(this_trait.chr).upper()[0])*1000 + + # this_trait.mb) + + #ZS: Put this in function currently called "convert_location_to_value" + this_trait.location_repr = 'Chr %s: %.4f Mb' % (this_trait.chr, + float(this_trait.mb)) this_trait.location_value = trait_location_value - #this_trait.trait_location_value = trait_location_value - #XZ, 01/12/08: This SQL query is much faster. + #Get mean expression value query = ( -"""select ProbeSetXRef.mean from ProbeSetXRef, ProbeSet - where ProbeSetXRef.ProbeSetFreezeId = %s and - ProbeSet.Id = ProbeSetXRef.ProbeSetId and - ProbeSet.Name = '%s' + """select ProbeSetXRef.mean from ProbeSetXRef, ProbeSet + where ProbeSetXRef.ProbeSetFreezeId = %s and + ProbeSet.Id = ProbeSetXRef.ProbeSetId and + ProbeSet.Name = '%s' """ % (escape(str(this_trait.dataset.id)), escape(this_trait.name))) - print("query is:", pf(query)) + #print("query is:", pf(query)) result = g.db.execute(query).fetchone() + + mean = result[0] if result else 0 - if result: - if result[0]: - mean = result[0] - else: - mean=0 - else: - mean = 0 - - #XZ, 06/05/2009: It is neccessary to turn on nowrap - this_trait.mean = repr = "%2.3f" % mean + this_trait.mean = "%2.3f" % mean #LRS and its location this_trait.LRS_score_repr = 'N/A' @@ -790,23 +788,39 @@ class MrnaAssayDataSet(DataSet): result = self.cursor.fetchone() if result: - if result[0] and result[1]: - LRS_Chr = result[0] - LRS_Mb = result[1] - - #XZ: LRS_location_value is used for sorting - try: - LRS_location_value = int(LRS_Chr)*1000 + float(LRS_Mb) - except: - if LRS_Chr.upper() == 'X': - LRS_location_value = 20*1000 + float(LRS_Mb) - else: - LRS_location_value = ord(str(LRS_chr).upper()[0])*1000 + float(LRS_Mb) + #if result[0] and result[1]: + # lrs_chr = result[0] + # lrs_mb = result[1] + lrs_chr, lrs_mb = result + #XZ: LRS_location_value is used for sorting + lrs_location_value = self.convert_location_to_value(lrs_chr, lrs_mb) + + #try: + # lrs_location_value = int(lrs_chr)*1000 + float(lrs_mb) + #except: + # if lrs_chr.upper() == 'X': + # lrs_location_value = 20*1000 + float(lrs_mb) + # else: + # lrs_location_value = (ord(str(LRS_chr).upper()[0])*1000 + + # float(lrs_mb)) + + this_trait.LRS_score_repr = '%3.1f' % this_trait.lrs + this_trait.LRS_score_value = this_trait.lrs + this_trait.LRS_location_repr = 'Chr %s: %.4f Mb' % (lrs_chr, float(lrs_mb)) + + + def convert_location_to_value(self, chromosome, mb): + try: + location_value = int(chromosome)*1000 + float(mb) + except ValueError: + if chromosome.upper() == 'X': + location_value = 20*1000 + float(mb) + else: + location_value = (ord(str(chromosome).upper()[0])*1000 + + float(mb)) + + return location_value - this_trait.LRS_score_repr = LRS_score_repr = '%3.1f' % this_trait.lrs - this_trait.LRS_score_value = LRS_score_value = this_trait.lrs - this_trait.LRS_location_repr = LRS_location_repr = 'Chr %s: %.4f Mb' % (LRS_Chr, float(LRS_Mb) ) - def get_sequence(self): query = """ SELECT @@ -912,7 +926,7 @@ class TempDataSet(DataSet): def geno_mrna_confidentiality(ob): dataset_table = ob.type + "Freeze" - print("dataset_table [%s]: %s" % (type(dataset_table), dataset_table)) + #print("dataset_table [%s]: %s" % (type(dataset_table), dataset_table)) query = '''SELECT Id, Name, FullName, confidentiality, AuthorisedUsers FROM %s WHERE Name = %%s''' % (dataset_table) diff --git a/wqflask/base/trait.py b/wqflask/base/trait.py index dde8b8d8..db76ddea 100755 --- a/wqflask/base/trait.py +++ b/wqflask/base/trait.py @@ -1,6 +1,8 @@ from __future__ import absolute_import, division, print_function import string +import resource + from htmlgen import HTMLgen2 as HT @@ -15,22 +17,38 @@ from pprint import pformat as pf from flask import Flask, g -class GeneralTrait: +def print_mem(stage=""): + mem = resource.getrusage(resource.RUSAGE_SELF).ru_maxrss + print("{}: {}".format(stage, mem/1024)) + +class GeneralTrait(object): """ Trait class defines a trait in webqtl, can be either Microarray, Published phenotype, genotype, or user input trait """ - def __init__(self, **kw): - print("in GeneralTrait") - self.dataset = kw.get('dataset') # database name + def __init__(self, get_qtl_info=False, **kw): + # xor assertion + assert bool(kw.get('dataset')) != bool(kw.get('dataset_name')), "Needs dataset ob. or name"; + if kw.get('dataset_name'): + self.dataset = create_dataset(kw.get('dataset_name')) + else: + self.dataset = kw.get('dataset') self.name = kw.get('name') # Trait ID, ProbeSet ID, Published ID, etc. self.cellid = kw.get('cellid') self.identification = kw.get('identification', 'un-named trait') self.haveinfo = kw.get('haveinfo', False) self.sequence = kw.get('sequence') # Blat sequence, available for ProbeSet self.data = kw.get('data', {}) + + # Sets defaultst + self.locus = None + self.lrs = None + self.pvalue = None + self.mean = None + self.num_overlap = None + if kw.get('fullname'): name2 = value.split("::") @@ -39,13 +57,12 @@ class GeneralTrait: # self.cellid is set to None above elif len(name2) == 3: self.dataset, self.name, self.cellid = name2 - - self.dataset = create_dataset(self.dataset) # Todo: These two lines are necessary most of the time, but perhaps not all of the time # So we could add a simple if statement to short-circuit this if necessary - self.retrieve_info() + self.retrieve_info(get_qtl_info=get_qtl_info) self.retrieve_sample_data() + def get_name(self): @@ -78,7 +95,7 @@ class GeneralTrait: #desc = self.handle_pca(desc) stringy = desc return stringy - + def display_name(self): @@ -229,7 +246,7 @@ class GeneralTrait: #def items(self): # return self.__dict__.items() - def retrieve_info(self, QTL=False): + def retrieve_info(self, get_qtl_info=False): assert self.dataset, "Dataset doesn't exist" if self.dataset.type == 'Publish': query = """ @@ -251,7 +268,7 @@ class GeneralTrait: PublishXRef.InbredSetId = PublishFreeze.InbredSetId AND PublishFreeze.Id = %s """ % (self.name, self.dataset.id) - traitInfo = g.db.execute(query).fetchone() + trait_info = g.db.execute(query).fetchone() #XZ, 05/08/2009: Xiaodong add this block to use ProbeSet.Id to find the probeset instead of just using ProbeSet.Name #XZ, 05/08/2009: to avoid the problem of same probeset name from different platforms. elif self.dataset.type == 'ProbeSet': @@ -268,8 +285,8 @@ class GeneralTrait: """ % (escape(display_fields_string), escape(self.dataset.name), escape(self.name)) - traitInfo = g.db.execute(query).fetchone() - print("traitInfo is: ", pf(traitInfo)) + trait_info = g.db.execute(query).fetchone() + #print("trait_info is: ", pf(trait_info)) #XZ, 05/08/2009: We also should use Geno.Id to find marker instead of just using Geno.Name # to avoid the problem of same marker name from different species. elif self.dataset.type == 'Geno': @@ -286,23 +303,24 @@ class GeneralTrait: """ % (escape(display_fields_string), escape(self.dataset.name), escape(self.name)) - traitInfo = g.db.execute(query).fetchone() - print("traitInfo is: ", pf(traitInfo)) + trait_info = g.db.execute(query).fetchone() + #print("trait_info is: ", pf(trait_info)) else: #Temp type query = """SELECT %s FROM %s WHERE Name = %s """ % (string.join(self.dataset.display_fields,','), self.dataset.type, self.name) - traitInfo = g.db.execute(query).fetchone() + trait_info = g.db.execute(query).fetchone() #self.cursor.execute(query) - #traitInfo = self.cursor.fetchone() - if traitInfo: + #trait_info = self.cursor.fetchone() + if trait_info: self.haveinfo = True #XZ: assign SQL query result to trait attributes. for i, field in enumerate(self.dataset.display_fields): - setattr(self, field, str(traitInfo[i])) + print(" mike: {} -> {} - {}".format(field, type(trait_info[i]), trait_info[i])) + setattr(self, field, trait_info[i]) if self.dataset.type == 'Publish': self.confidential = 0 @@ -310,6 +328,10 @@ class GeneralTrait: self.confidential = 1 self.homologeneid = None + + print("self.geneid is:", self.geneid) + print(" type:", type(self.geneid)) + print("self.dataset.group.name is:", self.dataset.group.name) if self.dataset.type == 'ProbeSet' and self.dataset.group and self.geneid: #XZ, 05/26/2010: From time to time, this query get error message because some geneid values in database are not number. #XZ: So I have to test if geneid is number before execute the query. @@ -321,6 +343,8 @@ class GeneralTrait: # geneidIsNumber = False #if geneidIsNumber: + + query = """ SELECT HomologeneId @@ -332,6 +356,7 @@ class GeneralTrait: InbredSet.SpeciesId = Species.Id AND Species.TaxonomyId = Homologene.TaxonomyId """ % (escape(str(self.geneid)), escape(self.dataset.group.name)) + print("-> query is:", query) result = g.db.execute(query).fetchone() #else: # result = None @@ -339,26 +364,40 @@ class GeneralTrait: if result: self.homologeneid = result[0] - if QTL: + if get_qtl_info: if self.dataset.type == 'ProbeSet' and not self.cellid: - traitQTL = g.db.execute(""" + query = """ SELECT ProbeSetXRef.Locus, ProbeSetXRef.LRS, ProbeSetXRef.pValue, ProbeSetXRef.mean FROM ProbeSetXRef, ProbeSet WHERE ProbeSetXRef.ProbeSetId = ProbeSet.Id AND - ProbeSet.Name = "%s" AND - ProbeSetXRef.ProbeSetFreezeId =%s - """, (self.name, self.dataset.id)).fetchone() + ProbeSet.Name = "{}" AND + ProbeSetXRef.ProbeSetFreezeId ={} + """.format(self.name, self.dataset.id) + trait_qtl = g.db.execute(query).fetchone() #self.cursor.execute(query) - #traitQTL = self.cursor.fetchone() - if traitQTL: - self.locus, self.lrs, self.pvalue, self.mean = traitQTL + #trait_qtl = self.cursor.fetchone() + if trait_qtl: + self.locus, self.lrs, self.pvalue, self.mean = trait_qtl + if self.locus: + query = """ + select Geno.Chr, Geno.Mb from Geno, Species + where Species.Name = '{}' and + Geno.Name = '{}' and + Geno.SpeciesId = Species.Id + """.format(self.dataset.group.species, self.locus) + print("query is:", query) + result = g.db.execute(query).fetchone() + self.locus_chr = result[0] + self.locus_mb = result[1] else: - self.locus = self.lrs = self.pvalue = self.mean = "" + self.locus = self.locus_chr = self.locus_mb = self.lrs = self.pvalue = self.mean = "" + + if self.dataset.type == 'Publish': - traitQTL = g.db.execute(""" + trait_qtl = g.db.execute(""" SELECT PublishXRef.Locus, PublishXRef.LRS FROM @@ -369,9 +408,9 @@ class GeneralTrait: PublishFreeze.Id =%s """, (self.name, self.dataset.id)).fetchone() #self.cursor.execute(query) - #traitQTL = self.cursor.fetchone() - if traitQTL: - self.locus, self.lrs = traitQTL + #trait_qtl = self.cursor.fetchone() + if trait_qtl: + self.locus, self.lrs = trait_qtl else: self.locus = self.lrs = "" else: diff --git a/wqflask/base/webqtlConfig.py b/wqflask/base/webqtlConfig.py index 49afb631..a811c3cd 100755 --- a/wqflask/base/webqtlConfig.py +++ b/wqflask/base/webqtlConfig.py @@ -52,7 +52,7 @@ ENSEMBLETRANSCRIPT_URL="http://useast.ensembl.org/Mus_musculus/Lucene/Details?sp SECUREDIR = GNROOT + 'secure/' COMMON_LIB = GNROOT + 'support/admin' HTMLPATH = GNROOT + 'web/' -PYLMM_PATH = HTMLPATH + 'plink/' +PYLMM_PATH = '/home/zas1024/' SNP_PATH = '/mnt/xvdf1/snps/' IMGDIR = HTMLPATH +'image/' IMAGESPATH = HTMLPATH + 'images/' |