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:
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.
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.
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.
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:
#
# 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.
# 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.
# 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.
# 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.
# 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.
#
#
# """
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
'''
return output
# showSelectedTraits: listof Trait -> string
# to show a list of the selected traits
def showSelectedTraits(self, traits, p, inbredSetName):
output = '''
"
return output
# showSummaryCorrelationResults
# show just the number of traits in each subarray
def showSummaryCorrelationResults(self, p, traits, traitArray):
output = '''
'''
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 += '''
''' % 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']