path: root/web/webqtl/correlation/correlationFunction.py

# 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 NL 2011/03/23


import math
import rpy2.robjects
import pp
import string

from utility import webqtlUtil
from base.webqtlTrait import webqtlTrait
from dbFunction import webqtlDatabaseFunction



#XZ: The input 'controls' is String. It contains the full name of control traits.
#XZ: The input variable 'strainlst' is List. It contains the strain names of primary trait.
#XZ: The returned tcstrains is the list of list [[],[]...]. So are tcvals and tcvars. The last returned parameter is list of numbers.
#XZ, 03/29/2010: For each returned control trait, there is no None value in it.
def controlStrains(controls, strainlst):

    controls = controls.split(',')

    cvals = {}
    for oneTraitName in controls:
        oneTrait = webqtlTrait(fullname=oneTraitName, cursor=webqtlDatabaseFunction.getCursor() )
        oneTrait.retrieveData()
        cvals[oneTraitName] = oneTrait.data

    tcstrains = []
    tcvals = []
    tcvars = []

    for oneTraitName in controls:
        strains = []
        vals = []
        vars = []

        for _strain in strainlst:
            if cvals[oneTraitName].has_key(_strain):
                _val = cvals[oneTraitName][_strain].val
                if _val != None:
                    strains.append(_strain)
                    vals.append(_val)
                    vars.append(None)

        tcstrains.append(strains)
        tcvals.append(vals)
        tcvars.append(vars)

    return tcstrains, tcvals, tcvars, [len(x) for x in tcstrains]



#XZ, 03/29/2010: After execution of functon "controlStrains" and "fixStrains", primary trait and control traits have the same strains and in the same order. There is no 'None' value in them.
def fixStrains(_strains,_controlstrains,_vals,_controlvals,_vars,_controlvars):
    """Corrects strains, vals, and vars so that all contrain only those strains common
    to the reference trait and all control traits."""

    def dictify(strains,vals,vars):
        subdict = {}
        for i in xrange(len(strains)):
            subdict[strains[i]] = (vals[i],vars[i])
        return subdict

    #XZ: The 'dicts' is a list of dictionary. The first element is the dictionary of reference trait. The rest elements are for control traits.
    dicts = []
    dicts.append(dictify(_strains,_vals,_vars))

    nCstrains = len(_controlstrains)
    for i in xrange(nCstrains):
        dicts.append(dictify(_controlstrains[i],_controlvals[i],_controlvars[i]))

    _newstrains = []
    _vals = []
    _vars = []
    _controlvals = [[] for x in xrange(nCstrains)]
    _controlvars = [[] for x in xrange(nCstrains)]

    for strain in _strains:
        inall = True
        for d in dicts:
            if strain not in d: 
                inall = False
                break
        if inall:
            _newstrains.append(strain)
            _vals.append(dicts[0][strain][0])
            _vars.append(dicts[0][strain][1])
            for i in xrange(nCstrains):
                _controlvals[i].append(dicts[i+1][strain][0])
                _controlvars[i].append(dicts[i+1][strain][1])

    return _newstrains,  _vals, _controlvals, _vars, _controlvars


#XZ, 6/15/2010: If there is no identical control traits, the returned list is empty.
#else, the returned list has two elements of control trait name.
def findIdenticalControlTraits ( controlVals, controlNames ):
    nameOfIdenticalTraits = []

    controlTraitNumber = len(controlVals)

    if controlTraitNumber > 1:

        #XZ: reset the precision of values and convert to string type
        for oneTraitVal in controlVals:
            for oneStrainVal in oneTraitVal:
                oneStrainVal = '%.3f' % oneStrainVal                

        for i, oneTraitVal in enumerate( controlVals ):
            for j in range(i+1, controlTraitNumber):
                if oneTraitVal == controlVals[j]:
                    nameOfIdenticalTraits.append(controlNames[i])
                    nameOfIdenticalTraits.append(controlNames[j])

    return nameOfIdenticalTraits

#XZ, 6/15/2010: If there is no identical control traits, the returned list is empty.
#else, the returned list has two elements of control trait name.
#primaryVal is of list type. It contains value of primary trait.
#primaryName is of string type.
#controlVals is of list type. Each element is list too. Each element contain value of one control trait.
#controlNames is of list type.
def findIdenticalTraits (primaryVal, primaryName, controlVals, controlNames ):
    nameOfIdenticalTraits = []

    #XZ: reset the precision of values and convert to string type
    for oneStrainVal in primaryVal:
        oneStrainVal = '%.3f' % oneStrainVal

    for oneTraitVal in controlVals:
        for oneStrainVal in oneTraitVal:
            oneStrainVal = '%.3f' % oneStrainVal

    controlTraitNumber = len(controlVals)

    if controlTraitNumber > 1:
        for i, oneTraitVal in enumerate( controlVals ):
            for j in range(i+1, controlTraitNumber):
                if oneTraitVal == controlVals[j]:
                    nameOfIdenticalTraits.append(controlNames[i])
                    nameOfIdenticalTraits.append(controlNames[j])
                    break

    if len(nameOfIdenticalTraits) == 0:
        for i, oneTraitVal in enumerate( controlVals ):
            if primaryVal == oneTraitVal:
                nameOfIdenticalTraits.append(primaryName)
                nameOfIdenticalTraits.append(controlNames[i])
                break

    return nameOfIdenticalTraits



#XZ, 03/29/2010: The strains in primaryVal, controlVals, targetVals must be of the same number and in same order.
#XZ: No value in primaryVal and controlVals could be None.

def determinePartialsByR (primaryVal, controlVals, targetVals, targetNames, method='p'):

    def compute_partial ( primaryVal, controlVals, targetVals, targetNames, method ):

        rpy2.robjects.r("""
pcor.test <- function(x,y,z,use="mat",method="p",na.rm=T){
        # The partial correlation coefficient between x and y given z
        #
        # pcor.test is free and comes with ABSOLUTELY NO WARRANTY.
        #
        # x and y should be vectors
        #
        # z can be either a vector or a matrix
        #
        # use: There are two methods to calculate the partial correlation coefficient.
        #        One is by using variance-covariance matrix ("mat") and the other is by using recursive formula ("rec").
        #        Default is "mat".
        #
        # method: There are three ways to calculate the correlation coefficient,
        #           which are Pearson's ("p"), Spearman's ("s"), and Kendall's ("k") methods.
        #           The last two methods which are Spearman's and Kendall's coefficient are based on the non-parametric analysis.
        #           Default is "p".
        #
        # na.rm: If na.rm is T, then all the missing samples are deleted from the whole dataset, which is (x,y,z).
        #        If not, the missing samples will be removed just when the correlation coefficient is calculated.
        #          However, the number of samples for the p-value is the number of samples after removing
        #          all the missing samples from the whole dataset.
        #          Default is "T".

        x <- c(x)
        y <- c(y)
        z <- as.data.frame(z)

        if(use == "mat"){
                p.use <- "Var-Cov matrix"
                pcor = pcor.mat(x,y,z,method=method,na.rm=na.rm)
        }else if(use == "rec"){
                p.use <- "Recursive formula"
                pcor = pcor.rec(x,y,z,method=method,na.rm=na.rm)
        }else{
                stop("use should be either rec or mat!\n")
        }

        # print the method
        if(gregexpr("p",method)[[1]][1] == 1){
                p.method <- "Pearson"
        }else if(gregexpr("s",method)[[1]][1] == 1){
                p.method <- "Spearman"
        }else if(gregexpr("k",method)[[1]][1] == 1){
                p.method <- "Kendall"
        }else{
                stop("method should be pearson or spearman or kendall!\n")
        }

        # sample number
        n <- dim(na.omit(data.frame(x,y,z)))[1]

        # given variables' number
        gn <- dim(z)[2]

        # p-value
        if(p.method == "Kendall"){
                statistic <- pcor/sqrt(2*(2*(n-gn)+5)/(9*(n-gn)*(n-1-gn)))
                p.value <- 2*pnorm(-abs(statistic))

        }else{
                statistic <- pcor*sqrt((n-2-gn)/(1-pcor^2))
                p.value <- 2*pnorm(-abs(statistic))
        }

        data.frame(estimate=pcor,p.value=p.value,statistic=statistic,n=n,gn=gn,Method=p.method,Use=p.use)
}

# By using var-cov matrix
pcor.mat <- function(x,y,z,method="p",na.rm=T){

        x <- c(x)
        y <- c(y)
        z <- as.data.frame(z)

        if(dim(z)[2] == 0){
                stop("There should be given data\n")
        }

        data <- data.frame(x,y,z)

        if(na.rm == T){
                data = na.omit(data)
        }

        xdata <- na.omit(data.frame(data[,c(1,2)]))
        Sxx <- cov(xdata,xdata,m=method)

        xzdata <- na.omit(data)
        xdata <- data.frame(xzdata[,c(1,2)])
        zdata <- data.frame(xzdata[,-c(1,2)])
        Sxz <- cov(xdata,zdata,m=method)

        zdata <- na.omit(data.frame(data[,-c(1,2)]))
        Szz <- cov(zdata,zdata,m=method)

        # is Szz positive definite?
        zz.ev <- eigen(Szz)$values
        if(min(zz.ev)[1]<0){
                stop("\'Szz\' is not positive definite!\n")
        }

        # partial correlation
        Sxx.z <- Sxx - Sxz %*% solve(Szz) %*% t(Sxz)

        rxx.z <- cov2cor(Sxx.z)[1,2]

        rxx.z
}

# By using recursive formula
pcor.rec <- function(x,y,z,method="p",na.rm=T){
        #

        x <- c(x)
        y <- c(y)
        z <- as.data.frame(z)

        if(dim(z)[2] == 0){
                stop("There should be given data\n")
        }

        data <- data.frame(x,y,z)

        if(na.rm == T){
                data = na.omit(data)
        }

        # recursive formula
        if(dim(z)[2] == 1){
                tdata <- na.omit(data.frame(data[,1],data[,2]))
                rxy <- cor(tdata[,1],tdata[,2],m=method)

                tdata <- na.omit(data.frame(data[,1],data[,-c(1,2)]))
                rxz <- cor(tdata[,1],tdata[,2],m=method)

                tdata <- na.omit(data.frame(data[,2],data[,-c(1,2)]))
                ryz <- cor(tdata[,1],tdata[,2],m=method)

                rxy.z <- (rxy - rxz*ryz)/( sqrt(1-rxz^2)*sqrt(1-ryz^2) )

                return(rxy.z)
        }else{
                x <- c(data[,1])
                y <- c(data[,2])
                z0 <- c(data[,3])
                zc <- as.data.frame(data[,-c(1,2,3)])

                rxy.zc <- pcor.rec(x,y,zc,method=method,na.rm=na.rm)
                rxz0.zc <- pcor.rec(x,z0,zc,method=method,na.rm=na.rm)
                ryz0.zc <- pcor.rec(y,z0,zc,method=method,na.rm=na.rm)

                rxy.z <- (rxy.zc - rxz0.zc*ryz0.zc)/( sqrt(1-rxz0.zc^2)*sqrt(1-ryz0.zc^2) )
                return(rxy.z)
        }
}
""")

        R_pcorr_function = rpy2.robjects.r['pcor.test']
        R_corr_test = rpy2.robjects.r['cor.test']

        primary = rpy2.robjects.FloatVector(range(len(primaryVal)))
        for i in range(len(primaryVal)):
            primary[i] = primaryVal[i]

        control = rpy2.robjects.r.matrix(rpy2.robjects.FloatVector( range(len(controlVals)*len(controlVals[0])) ), ncol=len(controlVals))
        for i in range(len(controlVals)):
            for j in range(len(controlVals[0])):
                control[i*len(controlVals[0]) + j] = controlVals[i][j]

        allcorrelations = []

        for targetIndex, oneTargetVals in enumerate(targetVals):

            this_primary = None
            this_control = None
            this_target = None

            if None in oneTargetVals:

                goodIndex = []
                for i in range(len(oneTargetVals)):
                    if oneTargetVals[i] != None:
                        goodIndex.append(i)

                this_primary = rpy2.robjects.FloatVector(range(len(goodIndex)))
                for i in range(len(goodIndex)):
                    this_primary[i] = primaryVal[goodIndex[i]]

                this_control = rpy2.robjects.r.matrix(rpy2.robjects.FloatVector( range(len(controlVals)*len(goodIndex)) ), ncol=len(controlVals))
                for i in range(len(controlVals)):
                    for j in range(len(goodIndex)):
                        this_control[i*len(goodIndex) + j] = controlVals[i][goodIndex[j]]

                this_target = rpy2.robjects.FloatVector(range(len(goodIndex)))
                for i in range(len(goodIndex)):
                    this_target[i] = oneTargetVals[goodIndex[i]]

            else:
                this_primary = primary
                this_control = control
                this_target = rpy2.robjects.FloatVector(range(len(oneTargetVals)))
                for i in range(len(oneTargetVals)):
                    this_target[i] = oneTargetVals[i]

            one_name = targetNames[targetIndex]
            one_N = len(this_primary)

            #calculate partial correlation
            one_pc_coefficient = 'NA'
            one_pc_p = 1

            try:
                if method == 's':
                    result = R_pcorr_function(this_primary, this_target, this_control, method='s')
                else:
                    result = R_pcorr_function(this_primary, this_target, this_control)

                #XZ: In very few cases, the returned coefficient is nan.
                #XZ: One way to detect nan is to compare the number to itself. NaN is always != NaN
                if result[0][0] == result[0][0]:
                    one_pc_coefficient = result[0][0]
                    #XZ: when the coefficient value is 1 (primary trait and target trait are the same),
                    #XZ: occationally, the returned p value is nan instead of 0.
                    if result[1][0] == result[1][0]:
                        one_pc_p = result[1][0]
                    elif abs(one_pc_coefficient - 1) < 0.0000001:
                        one_pc_p = 0
            except:
                pass

            #calculate zero order correlation
            one_corr_coefficient = 0
            one_corr_p = 1

            try:
                if method == 's':
                    R_result = R_corr_test(this_primary, this_target, method='spearman')
                else:
                    R_result = R_corr_test(this_primary, this_target)

                one_corr_coefficient = R_result[3][0]
                one_corr_p = R_result[2][0]
            except:
                pass

            traitinfo = [ one_name, one_N, one_pc_coefficient, one_pc_p, one_corr_coefficient, one_corr_p ]

            allcorrelations.append(traitinfo)

        return allcorrelations
    #End of function compute_partial


    allcorrelations = []

    target_trait_number = len(targetVals)

    if target_trait_number < 1000:
        allcorrelations = compute_partial ( primaryVal, controlVals, targetVals, targetNames, method )
    else:
        step = 1000
        job_number = math.ceil( float(target_trait_number)/step )

        job_targetVals_lists = []
        job_targetNames_lists = []

        for job_index in range( int(job_number) ):
            starti = job_index*step
            endi = min((job_index+1)*step, target_trait_number)

            one_job_targetVals_list = []
            one_job_targetNames_list = []

            for i in range( starti, endi ):
                one_job_targetVals_list.append( targetVals[i] )
                one_job_targetNames_list.append( targetNames[i] )

            job_targetVals_lists.append( one_job_targetVals_list )
            job_targetNames_lists.append( one_job_targetNames_list )

        ppservers = ()
        # Creates jobserver with automatically detected number of workers
        job_server = pp.Server(ppservers=ppservers)

        jobs = []
        results = []

        for i, one_job_targetVals_list in enumerate( job_targetVals_lists ):
            one_job_targetNames_list = job_targetNames_lists[i]
            #pay attention to modules from outside
            jobs.append( job_server.submit(func=compute_partial, args=( primaryVal, controlVals, one_job_targetVals_list, one_job_targetNames_list, method), depfuncs=(), modules=("rpy2.robjects",)) )

        for one_job in jobs:
            one_result = one_job()
            results.append( one_result )

        for one_result in results:
            for one_traitinfo in one_result:
                allcorrelations.append( one_traitinfo )

    return allcorrelations



#XZ, April 30, 2010: The input primaryTrait and targetTrait are instance of webqtlTrait
#XZ: The primaryTrait and targetTrait should have executed retrieveData function
def calZeroOrderCorr (primaryTrait, targetTrait, method='pearson'):

    #primaryTrait.retrieveData()

    #there is no None value in primary_val
    primary_strain, primary_val, primary_var = primaryTrait.exportInformative()

    #targetTrait.retrieveData()

    #there might be None value in target_val
    target_val = targetTrait.exportData(primary_strain, type="val")

    R_primary = rpy2.robjects.FloatVector(range(len(primary_val)))
    for i in range(len(primary_val)):
        R_primary[i] = primary_val[i]

    N = len(target_val)

    if None in target_val:
        goodIndex = []
        for i in range(len(target_val)):
            if target_val[i] != None:
                goodIndex.append(i)

        N = len(goodIndex)

        R_primary = rpy2.robjects.FloatVector(range(len(goodIndex)))
        for i in range(len(goodIndex)):
            R_primary[i] = primary_val[goodIndex[i]]

        R_target = rpy2.robjects.FloatVector(range(len(goodIndex)))
        for i in range(len(goodIndex)):
            R_target[i] = target_val[goodIndex[i]]

    else:
        R_target = rpy2.robjects.FloatVector(range(len(target_val)))
        for i in range(len(target_val)):
            R_target[i] = target_val[i]

    R_corr_test = rpy2.robjects.r['cor.test']

    if method == 'spearman':
        R_result = R_corr_test(R_primary, R_target, method='spearman')
    else:
        R_result = R_corr_test(R_primary, R_target)

    corr_result = []
    corr_result.append( R_result[3][0] )
    corr_result.append( N )
    corr_result.append( R_result[2][0] )
    
    return corr_result
	
#####################################################################################
#Input: primaryValue(list): one list of expression values of one probeSet, 
#       targetValue(list): one list of expression values of one probeSet, 
#		method(string): indicate correlation method ('pearson' or 'spearman')
#Output: corr_result(list): first item is Correlation Value, second item is tissue number,
#                           third item is PValue
#Function: get correlation value,Tissue quantity ,p value result by using R;
#Note : This function is special case since both primaryValue and targetValue are from
#the same dataset. So the length of these two parameters is the same. They are pairs.
#Also, in the datatable TissueProbeSetData, all Tissue values are loaded based on 
#the same tissue order
#####################################################################################	

def calZeroOrderCorrForTiss (primaryValue=[], targetValue=[], method='pearson'):

    R_primary = rpy2.robjects.FloatVector(range(len(primaryValue)))
    N = len(primaryValue)
    for i in range(len(primaryValue)):
        R_primary[i] = primaryValue[i]

    R_target = rpy2.robjects.FloatVector(range(len(targetValue)))
    for i in range(len(targetValue)):
        R_target[i]=targetValue[i]

    R_corr_test = rpy2.robjects.r['cor.test']
    if method =='spearman':
        R_result = R_corr_test(R_primary, R_target, method='spearman')
    else:
        R_result = R_corr_test(R_primary, R_target)
	
    corr_result =[]
    corr_result.append( R_result[3][0])
    corr_result.append( N )
    corr_result.append( R_result[2][0])

    return corr_result




def batchCalTissueCorr(primaryTraitValue=[], SymbolValueDict={}, method='pearson'):

    def cal_tissue_corr(primaryTraitValue, oneSymbolValueDict, method ):

        oneSymbolCorrDict = {}
        oneSymbolPvalueDict = {}

        R_corr_test = rpy2.robjects.r['cor.test']

        R_primary = rpy2.robjects.FloatVector(range(len(primaryTraitValue)))

        for i in range(len(primaryTraitValue)):
            R_primary[i] = primaryTraitValue[i]

        for (oneTraitSymbol, oneTraitValue) in oneSymbolValueDict.iteritems():
            R_target = rpy2.robjects.FloatVector(range(len(oneTraitValue)))
            for i in range(len(oneTraitValue)):
                R_target[i] = oneTraitValue[i]

            if method =='spearman':
                R_result = R_corr_test(R_primary, R_target, method='spearman')
            else:
                R_result = R_corr_test(R_primary, R_target)

            oneSymbolCorrDict[oneTraitSymbol] = R_result[3][0]
            oneSymbolPvalueDict[oneTraitSymbol] = R_result[2][0]

        return(oneSymbolCorrDict, oneSymbolPvalueDict)



    symbolCorrDict = {}
    symbolPvalueDict = {}

    items_number = len(SymbolValueDict)

    if items_number <= 1000:
        symbolCorrDict, symbolPvalueDict = cal_tissue_corr(primaryTraitValue, SymbolValueDict, method)
    else:
        items_list = SymbolValueDict.items()

        step = 1000
        job_number = math.ceil( float(items_number)/step )

        job_oneSymbolValueDict_list = []

        for job_index in range( int(job_number) ):
            starti = job_index*step
            endi = min((job_index+1)*step, items_number)

            oneSymbolValueDict = {}

            for i in range( starti, endi ):
                one_item = items_list[i]
                one_symbol = one_item[0]
                one_value = one_item[1]
                oneSymbolValueDict[one_symbol] = one_value

            job_oneSymbolValueDict_list.append( oneSymbolValueDict )


        ppservers = ()
        # Creates jobserver with automatically detected number of workers
        job_server = pp.Server(ppservers=ppservers)

        jobs = []
        results = []

        for i, oneSymbolValueDict in enumerate( job_oneSymbolValueDict_list ):

            #pay attention to modules from outside
            jobs.append( job_server.submit(func=cal_tissue_corr, args=(primaryTraitValue, oneSymbolValueDict, method), depfuncs=(), modules=("rpy2.robjects",)) )

        for one_job in jobs:
            one_result = one_job()
            results.append( one_result )

        for one_result in results:
            oneSymbolCorrDict, oneSymbolPvalueDict = one_result
            symbolCorrDict.update( oneSymbolCorrDict )
            symbolPvalueDict.update( oneSymbolPvalueDict )

    return (symbolCorrDict, symbolPvalueDict)

###########################################################################
#Input: cursor, GeneNameLst (list), TissueProbeSetFreezeId
#output: geneIdDict,dataIdDict,ChrDict,MbDict,descDict,pTargetDescDict (Dict)      
#function: get multi dicts for short and long label functions, and for getSymbolValuePairDict and 
# getGeneSymbolTissueValueDict to build dict to get CorrPvArray
#Note: If there are multiple probesets for one gene, select the one with highest mean.
###########################################################################	
def getTissueProbeSetXRefInfo(cursor=None,GeneNameLst=[],TissueProbeSetFreezeId=0):
	Symbols =""	
	symbolList =[]
	geneIdDict ={}
	dataIdDict = {}
	ChrDict = {}
	MbDict = {}
	descDict = {}
	pTargetDescDict = {}

	count = len(GeneNameLst)

	# Added by NL 01/06/2011
	# Note that:inner join is necessary in this query to get distinct record in one symbol group with highest mean value
	# Duo to the limit size of TissueProbeSetFreezeId table in DB, performance of inner join is acceptable.
	if count==0:
		query='''		
				select t.Symbol,t.GeneId, t.DataId,t.Chr, t.Mb,t.description,t.Probe_Target_Description 
				from (
					select Symbol, max(Mean) as maxmean 
					from TissueProbeSetXRef 
					where TissueProbeSetFreezeId=%s and Symbol!='' and Symbol Is Not Null group by Symbol) 
				as x inner join TissueProbeSetXRef as t on t.Symbol = x.Symbol and t.Mean = x.maxmean;				
			'''%TissueProbeSetFreezeId
				
	else:
		for i, item in enumerate(GeneNameLst):
		
			if i == count-1:
				Symbols += "'%s'" %item
			else:
				Symbols += "'%s'," %item
				
		Symbols = "("+ Symbols+")"
		query='''
				select t.Symbol,t.GeneId, t.DataId,t.Chr, t.Mb,t.description,t.Probe_Target_Description 
				from (
					select Symbol, max(Mean) as maxmean 
					from TissueProbeSetXRef 
					where TissueProbeSetFreezeId=%s and Symbol in %s group by Symbol) 
				as x inner join TissueProbeSetXRef as t on t.Symbol = x.Symbol and t.Mean = x.maxmean;		
			'''% (TissueProbeSetFreezeId,Symbols)

	try: 
	
		cursor.execute(query)
		results =cursor.fetchall()
		resultCount = len(results)
		# Key in all dicts is the lower-cased symbol
		for i, item in enumerate(results):	
			symbol = item[0]
			symbolList.append(symbol)
			
			key =symbol.lower()
			geneIdDict[key]=item[1]
			dataIdDict[key]=item[2]			
			ChrDict[key]=item[3]
			MbDict[key]=item[4]
			descDict[key]=item[5]
			pTargetDescDict[key]=item[6]

	except:
		symbolList = None
		geneIdDict=None
		dataIdDict=None
		ChrDict=None
		MbDict=None
		descDict=None
		pTargetDescDict=None
		
	return symbolList,geneIdDict,dataIdDict,ChrDict,MbDict,descDict,pTargetDescDict
		
###########################################################################
#Input: cursor, symbolList (list), dataIdDict(Dict)
#output: symbolValuepairDict (dictionary):one dictionary of Symbol and Value Pair,
#        key is symbol, value is one list of expression values of one probeSet;
#function: get one dictionary whose key is gene symbol and value is tissue expression data (list type).
#Attention! All keys are lower case!
###########################################################################	
def getSymbolValuePairDict(cursor=None,symbolList=None,dataIdDict={}):		
	symbolList = map(string.lower, symbolList)
	symbolValuepairDict={}
	valueList=[]

	for key in symbolList:
		if dataIdDict.has_key(key):
			DataId = dataIdDict[key]
			
			valueQuery = "select value from TissueProbeSetData where Id=%s" % DataId
			try :
				cursor.execute(valueQuery)
				valueResults = cursor.fetchall()
				for item in valueResults:
					item =item[0]
					valueList.append(item)	
				symbolValuepairDict[key] = valueList
				valueList=[]
			except:
				symbolValuepairDict[key] = None
			
	return symbolValuepairDict


########################################################################################################
#input: cursor, symbolList (list), dataIdDict(Dict): key is symbol
#output: SymbolValuePairDict(dictionary):one dictionary of Symbol and Value Pair.
#        key is symbol, value is one list of expression values of one probeSet.
#function: wrapper function for getSymbolValuePairDict function
#          build gene symbol list if necessary, cut it into small lists if necessary,
#          then call getSymbolValuePairDict function and merge the results.
########################################################################################################

def getGeneSymbolTissueValueDict(cursor=None,symbolList=None,dataIdDict={}):
        limitNum=1000
        count = len(symbolList)

        SymbolValuePairDict = {}

        if count !=0 and count <=limitNum:
                SymbolValuePairDict = getSymbolValuePairDict(cursor=cursor,symbolList=symbolList,dataIdDict=dataIdDict)

        elif count >limitNum:
                SymbolValuePairDict={}
                n = count/limitNum
                start =0
                stop =0

                for i in range(n):
                        stop =limitNum*(i+1)
                        gList1 = symbolList[start:stop]
                        PairDict1 = getSymbolValuePairDict(cursor=cursor,symbolList=gList1,dataIdDict=dataIdDict)
                        start =limitNum*(i+1)

                        SymbolValuePairDict.update(PairDict1)

                if stop < count:
                        stop = count
                        gList2 = symbolList[start:stop]
                        PairDict2 = getSymbolValuePairDict(cursor=cursor,symbolList=gList2,dataIdDict=dataIdDict)
                        SymbolValuePairDict.update(PairDict2)

        return SymbolValuePairDict

########################################################################################################
#input: cursor, GeneNameLst (list), TissueProbeSetFreezeId(int)
#output: SymbolValuePairDict(dictionary):one dictionary of Symbol and Value Pair.
#        key is symbol, value is one list of expression values of one probeSet.
#function: wrapper function of getGeneSymbolTissueValueDict function
#          for CorrelationPage.py
########################################################################################################

def getGeneSymbolTissueValueDictForTrait(cursor=None,GeneNameLst=[],TissueProbeSetFreezeId=0):
	SymbolValuePairDict={}
	symbolList,geneIdDict,dataIdDict,ChrDict,MbDict,descDict,pTargetDescDict = getTissueProbeSetXRefInfo(cursor=cursor,GeneNameLst=GeneNameLst,TissueProbeSetFreezeId=TissueProbeSetFreezeId)
	if symbolList:
		SymbolValuePairDict = getGeneSymbolTissueValueDict(cursor=cursor,symbolList=symbolList,dataIdDict=dataIdDict)
	return SymbolValuePairDict
	
########################################################################################################
#Input: cursor(cursor): MySQL connnection cursor; 
#       priGeneSymbolList(list): one list of gene symbol;
#       symbolValuepairDict(dictionary): one dictionary of Symbol and Value Pair,
#               key is symbol, value is one list of expression values of one probeSet;
#Output: corrArray(array): array of Correlation Value, 
#        pvArray(array): array of PValue;
#Function: build corrArray, pvArray for display by calling  calculation function:calZeroOrderCorrForTiss
########################################################################################################

def getCorrPvArray(cursor=None,priGeneSymbolList=[],symbolValuepairDict={}):
        # setting initial value for corrArray, pvArray equal to 0
        Num = len(priGeneSymbolList)

        corrArray = [([0] * (Num))[:] for i in range(Num)] 
        pvArray = [([0] * (Num))[:] for i in range(Num)]
        i = 0
        for pkey in priGeneSymbolList:
                j = 0
                pkey = pkey.strip().lower()# key in symbolValuepairDict is low case
                if symbolValuepairDict.has_key(pkey):
                        priValue = symbolValuepairDict[pkey]
                        for tkey in priGeneSymbolList:
                                tkey = tkey.strip().lower()# key in symbolValuepairDict is low case
                                if priValue and symbolValuepairDict.has_key(tkey):
                                        tarValue = symbolValuepairDict[tkey]

                                        if tarValue:
                                                if i>j:
                                                        # corrArray stores Pearson Correlation values
                                                        # pvArray stores Pearson P-Values
                                                        pcorr_result =calZeroOrderCorrForTiss(primaryValue=priValue,targetValue=tarValue)
                                                        corrArray[i][j] =pcorr_result[0]
                                                        pvArray[i][j] =pcorr_result[2]
                                                elif i<j:
                                                        # corrArray stores Spearman Correlation values
                                                        # pvArray stores Spearman P-Values
                                                        scorr_result =calZeroOrderCorrForTiss(primaryValue=priValue,targetValue=tarValue,method='spearman')
                                                        corrArray[i][j] =scorr_result[0]
                                                        pvArray[i][j] =scorr_result[2]
                                                else:
                                                        # on the diagonal line, correlation value is 1, P-Values is 0
                                                        corrArray[i][j] =1
                                                        pvArray[i][j] =0
                                                j+=1
                                        else:
                                                corrArray[i][j] = None
                                                pvArray[i][j] = None
                                                j+=1
                                else:
                                        corrArray[i][j] = None
                                        pvArray[i][j] = None
                                        j+=1
                else: 
                        corrArray[i][j] = None
                        pvArray[i][j] = None

                i+=1

        return corrArray, pvArray

########################################################################################################
#Input: cursor(cursor): MySQL connnection cursor; 
#       primaryTraitSymbol(string): one gene symbol;
#		TissueProbeSetFreezeId (int): Id of related TissueProbeSetFreeze
#       method: '0' default value, Pearson Correlation; '1', Spearman Correlation
#Output: symbolCorrDict(Dict): Dict of Correlation Value, key is symbol 
#        symbolPvalueDict(Dict): Dict of PValue,key is symbol ;
#Function: build symbolCorrDict, symbolPvalueDict for display by calling  calculation function:calZeroOrderCorrForTiss
########################################################################################################
def calculateCorrOfAllTissueTrait(cursor=None, primaryTraitSymbol=None, TissueProbeSetFreezeId=None,method='0'):

	symbolCorrDict = {}
	symbolPvalueDict = {}

	primaryTraitSymbolValueDict = getGeneSymbolTissueValueDictForTrait(cursor=cursor, GeneNameLst=[primaryTraitSymbol], TissueProbeSetFreezeId=TissueProbeSetFreezeId)
	primaryTraitValue = primaryTraitSymbolValueDict.values()[0]
	
	SymbolValueDict = getGeneSymbolTissueValueDictForTrait(cursor=cursor, GeneNameLst=[], TissueProbeSetFreezeId=TissueProbeSetFreezeId)
	
	if method =='1':
		symbolCorrDict, symbolPvalueDict = batchCalTissueCorr(primaryTraitValue,SymbolValueDict,method='spearman')
	else:
		symbolCorrDict, symbolPvalueDict = batchCalTissueCorr(primaryTraitValue,SymbolValueDict)


	return (symbolCorrDict, symbolPvalueDict)