# multitrait.py # a tool to analyze the correlations between several different traits and the traits # in a given dataset # # Parameters: # correlation -- either "pearson" or "spearman" depending on which ones we want to use # # filename -- an input file containing the traits to analyze # # progress -- if set, this parameter outputs a static progress page # and uses a META redirect to trigger the real computation # # targetDatabaseType: # one of "ProbeSet", "Publish", "Genotype" depending on the type of database # we will use for the analysis # # targetDatabaseId: # the id (*Freeze.Id in the database) of the particular database we will analyze # # threshold -- a float between 0 and 1 to determine which coefficents we wil l consider # # firstRun -- either 0 or 1 # whether to automatically pick reasonable defaults for the other three parameters # # outputType -- either "html" or "text" # # Author: Stephen Pitts # June 15, 2004 #Xiaodong changed the dependancy structure import copy import sys import cgi import os import os.path import math import time import numarray import tempfile import string import cgitb #all tracebacks come out as HTMLified CGI,useful when we have a random crash in the middle from base import templatePage from base.webqtlTrait import webqtlTrait from utility import webqtlUtil from base import webqtlConfig import trait import correlation import htmlModule cgitb.enable() # where this program's data files are RootDir = webqtlConfig.IMGDIR # XZ, 09/10/2008: add module name 'webqtlConfig.' RootDirURL = "/image/" # XZ, 09/10/2008: This parameter is not used in this module tempfile.tempdir = RootDir tempfile.template = "multitrait" # MultitraitException: used if something goes wrong # maybe in the future we should make exceptions more granular class MultitraitException(Exception): def __init__(self, message): self.message = message def __repr__(self): return "MultitraitException: %s" % self.message # buildParamDict: Cursor -> ParamDict # to process and validate CGI arguments # see the comment at the top of this file for valid cgi # parameters def buildParamDict(cursor, fd): params = {} fs = fd.formdata #cgi.FieldStorage() params["progress"] = fs.getfirst("progress", "0") params["filename"] = fs.getfirst("filename", "") if params["filename"] == "": raise MultitraitException("Required parameter filename missing.") params["targetDatabase"] = fs.getfirst("targetDatabase", "U74Av2RMA_Raw_ProbeSet_March04") params["firstRun"] = webqtlUtil.safeInt(fs.getfirst("firstRun", "0"),0) params["threshold"] = webqtlUtil.safeFloat(fs.getfirst("threshold", "0.5"), 0.5) params["subsetSize"] = webqtlUtil.safeInt(fs.getfirst("subsetSize", "10"), 10) if params["subsetSize"] < -1: params["subsetSize"] = -1 params["correlation"] = fs.getfirst("correlation", "pearson") params["subsetCount"] = webqtlUtil.safeInt(fs.getfirst("subsetCount", 10), 10) if params["subsetCount"] < -1: params["subsetCount"] = -1 #params["outputType"] = fs.getfirst("outputType", "html") #if params["outputType"] not in ("html", "text"): # params["outputType"] = "html" if params["correlation"] not in ("pearson", "spearman"): params["correlation"] = "pearson" params["correlationName"] = params["correlation"].capitalize() # one of two cases: # 1) We have just come from a submit, so there are a bunch of display* # but no displaySets. Thus, the code down there converts the display* # to displaySets so the GET request doesn't get too long # 2) We have just been redirected from a progress page which already has # a converted displaySets for us. displaySets = webqtlUtil.safeInt(fs.getfirst("displaySets","0"), 0) if displaySets == 0: for key in fs.keys(): if key[:7] == "display": #print "Hit display key %s
" % key try: whichSet = int(key[7:]) # prevent malicious attacks whichSet = min(whichSet, 512) displaySets += pow(2, whichSet) except ValueError: pass params["displaySets"] = displaySets #print "In the beginning, display sets was %s: %s
" % (displaySets, # str(binaryDecompose(displaySets))) # if we are just gonna display a progress page, then there's no # reason to look up detailed database information #if params["progress"] == "1": # return params a,b = trait.dbNameToTypeId(cursor, params["targetDatabase"]) # XZ, 09/10/2008: add module name params["targetDatabaseType"] = a params["targetDatabaseId"] = b params["targetDatabaseName"] = params["targetDatabase"] return params # readInputFile: DB cursor -> string -> string, (arrayof Trait) def readInputFile(cursor, filename): """ To read an input file with n lines in the following format ,, and retrieve and populate traits with appropriate data from the database Also, for our purposes. we store the database type and database id in fields attached to the trait instances. We use this information to generate Javascript popups with trait information. In addition, we read the strain of mice that the traits are from so we can only show those databases to correlate against. """ handle = open(filename) line = handle.readline() inbredSetName = line.strip() line = handle.readline() traits = [] # XZ, 09/10/2008: In this while loop block, I changed the original variable name 'trait' to 'oneTrait' while line != "": line = line.strip() dbType, dbId, tName = line.split(",") if dbType == "ProbeSet": oneTrait = trait.queryProbeSetTraitByName(cursor, tName) # XZ, 09/10/2008: add module name oneTrait.populateDataId(cursor, dbId) oneTrait.dbName = trait.dbTypeIdToName(cursor, dbType, dbId) # XZ, 09/10/2008: add module name elif dbType == "Geno": speciesId = trait.getSpeciesIdByDbTypeId(cursor, dbType, dbId) oneTrait = trait.queryGenotypeTraitByName(cursor, speciesId, tName) # XZ, 09/10/2008: add module name oneTrait.populateDataId(cursor, dbId) oneTrait.dbName = trait.dbTypeIdToName(cursor, dbType, dbId) # XZ, 09/10/2008: add module name elif dbType == "Publish": oneTrait = trait.queryPublishTraitByName(cursor, dbId, tName) # XZ, 09/10/2008: add module name oneTrait.populateDataId(cursor, dbId) oneTrait.dbName = trait.dbTypeIdToName(cursor, dbType, dbId) # XZ, 09/10/2008: add module name elif dbType == "Temp": oneTrait = trait.queryTempTraitByName(cursor, tName) # XZ, 09/10/2008: add module name oneTrait.populateDataId(cursor, dbId) oneTrait.dbName = "Temp" oneTrait.populateStrainData(cursor) traits.append(oneTrait) line = handle.readline() return inbredSetName, traits # loadDatabase: Cursor -> ParamDict -> arrayof Trait def loadDatabase(cursor, p): """ To load a set of traits as specified by the targetDatabaseId and targetDatabaseType parameters Cursor should be a fastCursor from the webqtl library (i.e. a MySQLdb SSCursor). Calling populateStrainData 20,000 or so times on a ProbeSet is really inefficent, so I wrote an optimized queryPopulatedProbeSetTraits in the trait module that uses a join to get all of the rows in bulk, store the resultset on the server, and do all sorts of nice buffering. It's about two or three times faster. """ if p["targetDatabaseType"] == "ProbeSet": # XZ, 09/10/2008: add module name dbTraits = trait.queryPopulatedProbeSetTraits(cursor, p["targetDatabaseId"]) elif p["targetDatabaseType"] == "Publish": # XZ, 09/10/2008: add module name dbTraits = trait.queryPublishTraits(cursor, p["targetDatabaseId"]) psd = trait.PublishTrait.populateStrainData elif p["targetDatabaseType"] == "Geno": # XZ, 09/10/2008: add module name dbTraits = trait.queryGenotypeTraits(cursor, p["targetDatabaseId"]) psd = trait.GenotypeTrait.populateStrainData else: print "Unknown target database type %s" % p["targetDatabaseType"] if p["targetDatabaseType"] != "ProbeSet": map(psd, dbTraits, [cursor]*len(dbTraits)) return dbTraits def runProbeSetCorrelations(cursor, p, traits): """ To run the correlations between the traits and the database. This function computes a correlation coefficent between each trait and every entry in the database, and partitions the database into a disjoint array of arrays which it returns. The length of the return array is 2^n, where n is the length of the trait array. Which constitutent element a of the return array a given trait ends up in is determined by the following formula i = i_02^0 + ... + i_(n-1)2^(n-1) where i_0 is 1 if corr(a,trait 0) >= threshold and 0 otherwise Since most of the several thousand database traits will end up with i=0, we don't return them, so the first element of the return array will be empty. A particular element of subarray j of the return array contains a 2-tuple (trait,kvalues). The variable trait is obviously the particular database trait that matches the user traits l_1, ..., l_m to which subarray j corresponds. kvalues is a list of the correlation values linking trait to l_1, ..., l_m, so the length of kvalues is the number of 1s in the binary representation of j (there must be a better way to describe this length). The return array is an array of 2-tuples. The first element of each tuple is the index of the particular subarray, and the second element is the subarray itself. The array is sorted in descending order by the number of 1's in the binary representation of the index so the first few subarrays are the ones that correspond to the largest sets. Each subarray is then sorted by the average of the magnitude of the individual correlation values. """ kMin = p["threshold"] traitArrays = {} # TODO: Add Spearman support freezeId = p["targetDatabaseId"] if p["correlation"] == "pearson": correlations = correlation.calcProbeSetPearsonMatrix(cursor, freezeId, traits) #XZ, 09/10/2008: add module name else: correlations = correlation.calcProbeSetSpearmanMatrix(freezeId, traits) #XZ, 09/10/2008: add module name # now we test all of the correlations in bulk test = numarray.absolute(correlations) test = numarray.greater_equal(test, kMin) test = test.astype(numarray.Int8) #print test db = trait.queryProbeSetTraits(cursor, freezeId) #XZ, 09/10/2008: add module name for i in range(len(db)): cIndex = 0 prods = [] for j in range(len(traits)): if test[i,j] == 1: cIndex += pow(2, j) prods.append(correlations[i,j]) if cIndex != 0: if not traitArrays.has_key(cIndex): traitArrays[cIndex] = [] traitArrays[cIndex].append((db[i], prods)) # sort each inner list of traitArrays # so the matched traits appear in descending order by the # average magnitude of the correlation def customCmp(traitPair, traitPair2): magAvg1 = numarray.average(map(abs, traitPair[1])) magAvg2 = numarray.average(map(abs, traitPair2[1])) # invert the sign to get descending order return -cmp(magAvg1, magAvg2) for traitArray in traitArrays.values(): traitArray.sort(customCmp) # sort the outer list of traitArrays traitArrays2 = [] i = 0 for key in traitArrays.keys(): a = traitArrays[key] if len(a) > 0: traitArrays2.append((key,a,len(binaryDecompose(key)), len(a))) # we sort by the number of 1's in the binary output # and then by the size of the list, both in descending order def customCmp2(aL,bL): a = -cmp(aL[2], bL[2]) if a == 0: return -cmp(aL[3], bL[3]) else: return a traitArrays2.sort(customCmp2) return traitArrays2 def runCorrelations(p, strainCount, traits, db): """ To run the correlations between the traits and the database. This function computes a correlation coefficent between each trait and every entry in the database, and partitions the database into a disjoint array of arrays which it returns. The length of the return array is 2^n, where n is the length of the trait array. Which constitutent element a of the return array a given trait ends up in is determined by the following formula i = i_02^0 + ... + i_(n-1)2^(n-1) where i_0 is 1 if corr(a,trait 0) >= threshold and 0 otherwise Since most of the several thousand database traits will end up with i=0, we don't return them, so the first element of the return array will be empty. A particular element of subarray j of the return array contains a 2-tuple (trait,kvalues). The variable trait is obviously the particular database trait that matches the user traits l_1, ..., l_m to which subarray j corresponds. kvalues is a list of the correlation values linking trait to l_1, ..., l_m, so the length of kvalues is the number of 1s in the binary representation of j (there must be a better way to describe this length). The return array is an array of 2-tuples. The first element of each tuple is the index of the particular subarray, and the second element is the subarray itself. The array is sorted in descending order by the number of 1's in the binary representation of the index so the first few subarrays are the ones that correspond to the largest sets. Each subarray is then sorted by the average of the magnitude of the individual correlation values. """ kMin = p["threshold"] traitArrays = {} # TODO: Add Spearman support if p["correlation"] == "pearson": correlations = correlation.calcPearsonMatrix(db, traits, strainCount) #XZ, 09/10/2008: add module name else: correlations = correlation.calcSpearmanMatrix(db, traits, strainCount) #XZ, 09/10/2008: add module name # now we test all of the correlations in bulk test = numarray.absolute(correlations) test = numarray.greater_equal(test, kMin) test = test.astype(numarray.Int8) #print test for i in range(len(db)): cIndex = 0 prods = [] for j in range(len(traits)): if test[i,j] == 1: cIndex += pow(2, j) prods.append(correlations[i,j]) if cIndex != 0: if not traitArrays.has_key(cIndex): traitArrays[cIndex] = [] traitArrays[cIndex].append((db[i], prods)) # sort each inner list of traitArrays # so the matched traits appear in descending order by the # average magnitude of the correlation def customCmp(traitPair, traitPair2): magAvg1 = numarray.average(map(abs, traitPair[1])) magAvg2 = numarray.average(map(abs, traitPair2[1])) # invert the sign to get descending order return -cmp(magAvg1, magAvg2) for traitArray in traitArrays.values(): traitArray.sort(customCmp) # sort the outer list of traitArrays traitArrays2 = [] i = 0 for key in traitArrays.keys(): a = traitArrays[key] if len(a) > 0: traitArrays2.append((key,a,len(binaryDecompose(key)), len(a))) # we sort by the number of 1's in the binary output # and then by the size of the list, both in descending order def customCmp2(aL,bL): a = -cmp(aL[2], bL[2]) if a == 0: return -cmp(aL[3], bL[3]) else: return a traitArrays2.sort(customCmp2) return traitArrays2 # XZ, 09/09/2008: In multiple trait correlation result page, # XZ, 09/09/2008: click "Download a text version of the above results in CSV format" # TraitCorrelationText: a class to display trait correlations # as textual output class TraitCorrelationText: # build a text shell to describe the given trait correlations # this method sets self.output; use str(self) to actually # get the text page # # traits is a list of traits and traitArray is a # list of 3-tuples: index, traits', garbage # where index is a binary-encoded description of which subset of # traits the list traits' matches # # traits' is a list of 3-tuples as well: trait, correlations, garbage # where trait is a particular trait and correlations is a list of float # correlations (matching traits above) def __init__(self, p, traits, traitArray): output = "Correlation Comparison\n" output += "from WebQTL and the University of Tennessee Health Science Center\n" output += "initiated at " + time.asctime(time.gmtime()) + " UTC\n\n" output += self.showOptionPanel(p) output += self.showSelectedTraits(traits) output += self.showSummaryCorrelationResults(p, traits, traitArray) output += self.showDetailedCorrelationResults(p, traits, traitArray) self.output = output # showOptionPanel: ParamDict -> string # to display the options used to run this correlation def showOptionPanel(self, params): output = "Correlation Comparison Options:\n" output += "Target database,%s\n" % params["targetDatabase"] output += "Correlation type,%s\n" % params["correlationName"] output += "Threshold,%f\n" % params["threshold"] #output += "Subsets to Show,%d\n" % params["subsetCount"] #output += "Traits to Show Per Subset,%d\n\n" % params["subsetSize"] return output # showSelectedTraits: (listof Trait) -> string # to display the traits compared with the database # note: we can't use tabular output because the traits could be of # different types and produce different fields def showSelectedTraits(self, traits): output = "Selected Traits:\n" for trait in traits: output += '"' + trait.longName() + '"' + "\n" output += "\n" return output # showSummaryCorrelationResults: ParamDict -> (listof Trait) -> # TraitArray -> string # see comment for __init__ for a description of TraitArray # # to show a summary (sets and sizes) of the correlation results # as well as an X to indicate whether they will be included # in the detailed output def showSummaryCorrelationResults(self, p, traits, traitArray): output = "Correlation Comparison Summary:\n" #if p["subsetCount"] != -1: # ourSubsetCount = min(p["subsetCount"], len(traitArray)) #else: ourSubsetCount = len(traitArray) displayDecomposition = binaryDecompose(p["displaySets"]) for j in range(ourSubsetCount): i = traitArray[j][0] traitSubarray = traitArray[j][1] if len(traitSubarray) == 0: continue targetTraits = decomposeIndex(traits, i) traitDesc = string.join(map(trait.Trait.shortName, targetTraits), # XZ, 09/10/2008: add module name ", ") if j in displayDecomposition: checked = "X" else: checked = "" output += '"%s","%s","%d"\n' % (checked, traitDesc, len(traitSubarray)) output += "\n" return output # showDetailedCorrelationResults: ParamDict -> (listof Trait) -> # TraitArray -> string # # to show a detailed list of the correlation results; that is, # to completely enumerate each subset of traitArray using the # filtering parameters in p def showDetailedCorrelationResults(self, p, traits, traitArray): output = "Correlation Comparison Details:\n" displayDecomposition = binaryDecompose(p["displaySets"]) displayDecomposition.sort() def formatCorr(c): return "%.4f" % c for j in displayDecomposition: i = traitArray[j][0] traitSubarray = traitArray[j][1] if len(traitSubarray) == 0: continue targetTraits = decomposeIndex(traits, i) extraColumnHeaders = map(trait.Trait.shortName, targetTraits) # XZ, 09/10/2008: add module name traitDesc = string.join(extraColumnHeaders, ", ") #if(p["subsetSize"] != -1 and len(traitSubarray) > p["subsetSize"]): # traitDesc += ",(showing top %s of %s)" % (p["subsetSize"], # len(traitSubarray)) # traitSubarray = traitSubarray[0:p["subsetSize"]] output += "%s\n" % traitDesc output += traitSubarray[0][0].csvHeader([], extraColumnHeaders) output += "\n" for oneTrait, corr in traitSubarray:#XZ, 09/10/2008: change original variable name 'trait' to 'oneTrait' corr = map(formatCorr, corr) output += oneTrait.csvRow([], corr) + "\n" output += "\n" return output # __str__ : string # to return self.output as the string representation of this page # self.output is built in __init__ def __str__(self): return self.output # TraitCorrelationPage: a class to display trait correlations # for now this is just one HTML file, so we don't even write it # to a temporary file somewhere class TraitCorrelationPage(templatePage.templatePage): """ Using the templatePage class, we build an HTML shell for the core data here: the trait correlation lists. The way templatePage works, we build the page in pieces in the __init__ method and later on use the inherited write method to render the page. """ def __init__(self, fd, p, cursor, traits, traitArray, inbredSetName, txtFilename): templatePage.templatePage.__init__(self, fd) self.dict["title"] = "Correlation Comparison" self.dict["basehref"] = "" # NL: deleted js1 content part, since it has not been used in this project self.dict["js1"] = "" self.dict["js2"] = "" body = "

Correlation Comparison

" body += "

Run at %s UTC

" % time.asctime(time.gmtime()) body += """

The correlation comparison tool identifies intersecting sets of traits that are correlated with your selections at a specified threshold. A correlation comparison involves the following steps:

  1. Correlate: Choose a Target Database, a Correlation Type, and a Correlation Threshold. For your initial correlation, leave Number of Subsets to Show and Traits to Show per Subset at their default values of 10. Using the Correlation Options panel, you can adjust the Correlation Threshold, Number of Subsets to Show, and Traits to Show per Subset.

  2. Add to Collection: You can use the check boxes in the Correlation Comparison Details panel and the buttons at the bottom of the page to add these results to your selections page for further analysis in WebQTL.

  3. Filter: Using the Correlation Comparison Summary panel, choose which subsets you would like to display for export. Note that if you change the parameters in the Correlation Options panel, you will need to re-apply your filter.

  4. Export: Once you are satisfied with your report, use the export link at the bottom of the page to save the report as a comma-separated (CSV) text file which you can then import into Excel or another tool. Note: the exported report will list all subsets in the summary view and only those traits in the subsets you have selected in the Filter step.

""" # body += """ #

The correlation # comparison tool identifies the intersecting sets of traits that are # correlated with your selections. A correlation comparison involves # the following steps:

#
    #
  1. Correlate: Choose a Target Database, a Correlation Type, and a Correlation Threshold. # For the initial correlation, leave Subsets to Show and Traits to Show per Subset # at their default values of 10.

  2. #
  3. Refine Correlation: Using the Correlation Options panel, # adjust the Correlation Threshold, Subsets to Show, and Traits to # Show per Subset until you have a reasonable number of traits.

  4. #
  5. Filter: Using the Correlation Comparison Summary panel, choose which subsets you would # like to see. Note that if you change the parameters in the Correlation Options panel, you will # loose the filter you have selected.

  6. #
  7. Export: Once you are satisfied with your report, use the export # link at the bottom of the page to save the report as a comma-separated (CSV) text file which # you can then import into Excel or another tool. Note: the exported report # will show all subsets in the summary view and all traits in each subset you have # selected in the Filter step. #

  8. Shopping Cart: In addition, you can use the # check boxes in the Correlation Comparison Details panel and the # buttons at the bottom of the page to add the traits you have found to the shopping cart.

    #
  9. #
# """ body += self.showOptionPanel(p, cursor, inbredSetName) body += self.showSelectedTraits(traits, p, inbredSetName) if p["firstRun"] == 0: body += self.showCorrelationResults(p, inbredSetName, traits, traitArray) exportParams = copy.copy(p) exportParams["outputType"] = "text" body += ('''

Export these results

Download a text version of the above results in CSV format. This text version differs from the version you see on this page in two ways. First, the summary view shows all subsets. Second, the details view shows all traits in the subsets that you have selected.

''' % txtFilename) body += "" self.dict["body"] = body # showOptionPanel: ParamDict -> Cursor -> String -> String # to build an option panel for the multitrait correlation # we expect the database list to be a list of 2-tuples # the first element of each tuple is the short name # and the second element of the tuple is the long name def showOptionPanel(self, params, cursor, inbredSetName): output = '''

Correlation Options

" corrSelected = ["",""] if params["correlation"] == "pearson": corrSelected[0] = "SELECTED" else: corrSelected[1] = "SELECTED" output += (''' ''' % (corrSelected[0], corrSelected[1])) output += ('' % params["threshold"]) output += ('' % params["subsetCount"]) output += ('' % params["subsetSize"]) # a cosmetic change to hopefully make this form a bit easier to use # if params["firstRun"] == 1: # applyName = "Correlate" # else: # applyName = "Refine Correlation" output += '''
Target Database: ''' % (webqtlConfig.CGIDIR, webqtlConfig.SCRIPTFILE, params["filename"]) output += htmlModule.genDatabaseMenu(db = cursor, public=0, RISetgp = inbredSetName, selectname="targetDatabase", selected=params["targetDatabase"]) output += "
Correlation Method:
Correlation Threshold:
Subsets to Show (-1 to show all subsets):
Traits to Show per Subset (-1 to show all traits):
''' return output # showSelectedTraits: listof Trait -> string # to show a list of the selected traits def showSelectedTraits(self, traits, p, inbredSetName): output = '''
''' % (webqtlConfig.CGIDIR, webqtlConfig.SCRIPTFILE, inbredSetName) output += "

Selected Traits

" output += '' flip = 1 colors = ["FFFFFF", "cccccc"] for trait in traits: # we take advantage of the secret dbName attribute that # loadDatabase fills in descriptionString = trait.genHTML() if trait.db.type == 'Publish' and trait.confidential: descriptionString = trait.genHTML(privilege=self.privilege, userName=self.userName, authorized_users=trait.authorized_users) output += ''' ''' % (colors[flip], trait.db.name, trait.db.name, trait.db.name, trait.name, descriptionString) flip = not flip output += "
DatabaseTrait
%s %s
" return output # showSummaryCorrelationResults # show just the number of traits in each subarray def showSummaryCorrelationResults(self, p, traits, traitArray): output = '''
''' % (webqtlConfig.CGIDIR, webqtlConfig.SCRIPTFILE, p["filename"], p["correlation"], p["threshold"], p["subsetCount"], p["subsetSize"], p["targetDatabase"]) output += ''' ''' # figure out a scale for the summary graph # for now we set max = 300 pixels wide if p["subsetCount"] != -1: ourSubsetCount = min(p["subsetCount"], len(traitArray)) else: ourSubsetCount = len(traitArray) screenWidth = 600 lengths = [] for j in range(ourSubsetCount): lengths.append(len(traitArray[j][1])) maxLength = max(lengths) displayDecomposition = binaryDecompose(p["displaySets"]) flip = 0 colors = ["FFFFFF", "cccccc"] for j in range(ourSubsetCount): i = traitArray[j][0] traitSubarray = traitArray[j][1] if len(traitSubarray) == 0: continue targetTraits = decomposeIndex(traits, i) traitDesc = string.join(map(webqtlTrait.displayName, targetTraits), ", ") if j in displayDecomposition: checked = "CHECKED" else: checked = "" barWidth = (len(traitSubarray) * screenWidth) / maxLength output += ('''''' % (colors[flip], j, checked, traitDesc, len(traitSubarray), barWidth)) flip = not flip output += '''
Trait Subsets Intersecting Set Size
%s %s
''' return output # showDetailedCorrelationResults # actually show the traits in each subarray def showDetailedCorrelationResults(self, p, inbredSetName, traits, traitArray): output = "

Correlation Comparison Details

" # the hidden form below powers all of the JavaScript links, # the shopping cart links, and the correlation plot links output += '''
''' % (webqtlConfig.CGIDIR, webqtlConfig.SCRIPTFILE, p["targetDatabase"], inbredSetName) displayDecomposition = binaryDecompose(p["displaySets"]) # necessary to ensure that subset order is the same in the # summary and the detailed view displayDecomposition.sort() # here's a trick: the first trait we show must have the widest row because it correlates # with the largest set of input traits firstSubset = traitArray[displayDecomposition[0]] firstTrait = firstSubset[1][0][0] extraColumnCount = firstSubset[2] totalColumnCount = 1 + len(firstTrait.row()) + extraColumnCount output += "\n" for j in displayDecomposition: i = traitArray[j][0] traitSubarray = traitArray[j][1] # we don't display trait combinations for which there are # no correlations if len(traitSubarray) == 0: continue # generate a description of the traits that this particular array # matches highly targetTraits = decomposeIndex(traits, i) extraColumnHeaders = map(webqtlTrait.displayName, targetTraits) traitDesc = string.join(extraColumnHeaders, ", ") # massage extraColumnHeaders so that they can be wrapped for i in range(len(extraColumnHeaders)): ech = extraColumnHeaders[i] ech = ech.replace("-", " ") ech = ech.replace("_", " ") extraColumnHeaders[i] = ech # pad extraColumnHeaders if we have less columns than the max paddingNeeded = extraColumnCount - len(extraColumnHeaders) if paddingNeeded > 0: extraColumnHeaders.extend(paddingNeeded * [" "]) # we limit the output to the top ones if(p["subsetSize"] != -1 and len(traitSubarray) > p["subsetSize"]): traitDesc += " (showing top %s of %s)" % (p["subsetSize"], len(traitSubarray)) traitSubarray = traitSubarray[0:p["subsetSize"]] # combine that description with actual database traits themselves # and the correlation values output += '' % (totalColumnCount, traitDesc) #output += '

%s

\n

%s

\n'% traitDesc # we assume that every trait in traitSubarray is the same type # of trait flip = 0 colors = ["FFFFFF", "cccccc"] output += traitSubarray[0][0].tableRowHeader([" "], extraColumnHeaders, colors[0]) for traitPair in traitSubarray: corr = [] traitPair[0].dbName = p['targetDatabase'] trait = traitPair[0] for i in range(len(traitPair[1])): corrValue = traitPair[1][i] corrPlotLink = (''' %.2f ''' % (p["targetDatabaseName"], trait.name, targetTraits[i].db.name, targetTraits[i].name, "0", corrValue)) corr.append(corrPlotLink) corr.extend(paddingNeeded * [" "]) checkbox = ('' % (p["targetDatabaseName"], trait.name)) flip = not flip output += traitPair[0].tableRow([checkbox], corr, colors[flip]) #output += "
" i += 1 output += ' ' % totalColumnCount output += "" # print form buttons if there were checkboxes above output += '''
''' % inbredSetName return output # showCorrelationResults: ParamDict -> listof Trait -> tupleof (int,arrayof trait) -> String # to build an output display for the multitrait correlation results def showCorrelationResults(self, p, inbredSetName, traits, traitArray): output = '''

Correlation Comparison Summary

%s correlations were computed for each of the selected traits with each trait in the %s database. Subsets of database traits for which correlations were higher than %s or lower than -%s are shown below based on which traits they correlated highly with. The top %s subsets, ranked by the number of input traits that they correspond with, are shown, and at most %s traits in each subset are shown.

''' % (p["correlationName"], p["targetDatabase"], p["targetDatabaseName"], p["threshold"], p["threshold"], p["subsetCount"], p["subsetSize"]) totalTraits = 0 for j in range(len(traitArray)): totalTraits += len(traitArray[j][1]) if totalTraits == 0: output += """

No shared corrrelates were found with your given traits at this threshold. You may wish to lower the correlation threshold or choose different traits.

""" else: output += self.showSummaryCorrelationResults(p, traits, traitArray) output += self.showDetailedCorrelationResults(p, inbredSetName, traits, traitArray) return output # decomposeIndex: (listof Trait) -> Int -> # (listof Trait) # to use i to partition T into a sublist # each bit in i controls the inclusion or exclusion of a trait def decomposeIndex(traits, i): targetTraits = [] for j in range(len(traits)): # look, mom, a bitwise and! # expression below tests whether the jth bit is # set in i # see runCorrelation for how we decompose the # array index if (i & pow(2,j)) == pow(2,j): targetTraits.append(traits[j]) return targetTraits # binaryDecompose: int -> (listof int) # to decompose a number into its constituent powers of 2 # returns a list of the exponents a_1...a_n such that the input m # is m = 2^a_1 + ... + 2^a_n def binaryDecompose(n): if n == 0: return [] # we start with the highest power of 2 <= this number # and work our way down, subtracting powers of 2 start = long(math.floor(math.log(n)/math.log(2))) exponents = [] while start >= 0: if n >= long(math.pow(2, start)): n -= math.pow(2,start) exponents.append(start) start -= 1 return exponents # powerOf : int -> int -> boolean # to determine whether m is a power of n; # more precisely, whether there exists z in Z s.t. # n^z = m def powerOf(m, n): trialZ = math.floor(math.log(m)/math.log(n)) return pow(n,trialZ) == m class compCorrPage(templatePage.templatePage): def __init__(self,fd): templatePage.templatePage.__init__(self, fd) if not self.openMysql(): return cursor = self.cursor params = buildParamDict(cursor, fd) # get the input data inbredSetName, traits = readInputFile(cursor, RootDir + params["filename"]) # and what we are comparing the data to dbTraits = [] if params["targetDatabaseType"] != "ProbeSet": dbTraits = loadDatabase(cursor, params) # run the comparison itself strainCount = trait.queryStrainCount(cursor) # XZ, 09/10/2008: add module name if params["targetDatabaseType"] == "ProbeSet": results = runProbeSetCorrelations(cursor, params, traits) else: results = runCorrelations(params, strainCount, traits, dbTraits) # try to be smart about what to output: # we want to limit the number of traits shown, at least initially # and since traitArray is already sorted with most interesting # subsets first, we simply pick up the first 500 or so traits # that we find if params["displaySets"] == 0: selectedTraits = 0 for j in range(len(results)): #print "Scanning subarray %d" % j if selectedTraits <= 200: params["displaySets"] += pow(2, j) selectedTraits += len(results[j][1]) traitList = [] for oneTrait in traits: # XZ, 09/10/2008: change the original variable name 'trait' to 'oneTrait' traitName = oneTrait.dbName+'::'+oneTrait.name # XZ, 09/10/2008: change the original variable name 'trait' to 'oneTrait' aTrait = webqtlTrait(cursor=self.cursor, fullname=traitName) traitList.append(aTrait) # and generate some output txtOutputFilename = tempfile.mktemp() txtOutputHandle = open(txtOutputFilename, "w") txtOutput = TraitCorrelationText(params, traits, results) txtOutputHandle.write(str(txtOutput)) txtOutputHandle.close() txtOutputFilename = os.path.split(txtOutputFilename)[1] self.dict['body'] = TraitCorrelationPage(fd, params, cursor, traitList, results, inbredSetName, txtOutputFilename).dict['body']