Changeset 1774:d0585dc25172 in orange-bioinformatics for _bioinformatics/obiGeneSetSig.py


Ignore:
Timestamp:
05/07/13 16:48:33 (12 months ago)
Author:
markotoplak
Branch:
default
Message:

obiGeneSetSig: added another method based on likelihood ratios by Su et al (2009).

File:
1 edited

Legend:

Unmodified
Added
Removed
  • _bioinformatics/obiGeneSetSig.py

    r1772 r1774  
    2828 
    2929    def _match_instance(self, instance, geneset, takegenes=None): 
     30        """ 
     31        Return 
     32            - a gene matcher with the instance as a target 
     33            - { name: attribute indices } of an instance 
     34            - genes names on the data set that were matched by the gene set 
     35 
     36        If takegenes is a list of indices, use only genes from 
     37        the gene set with specified indices. 
     38        """ 
    3039        nm, name_ind = self._mat_ni(instance) 
    3140        genes = [ nm.umatch(gene) for gene in geneset ] 
     
    126135            return 0 
    127136 
     137def estimate_gaussian_per_class(data, i, a=None, b=None): 
     138    cv = data.domain.class_var 
     139 
     140    if a == None: a = cv.values[0] 
     141    if b == None: b = cv.values[1] 
     142 
     143    def avWCVal(value): 
     144        return [ex[i].value for ex in data if ex[-1].value == value and not ex[i].isSpecial() ] 
     145 
     146    list1 = avWCVal(a) 
     147    list2 = avWCVal(b) 
     148 
     149    mi1 = mi2 = st1 = st2 = None 
     150 
     151    try: 
     152        mi1 = statc.mean(list1) 
     153        st1 = statc.std(list1) 
     154    except: 
     155        pass 
     156 
     157    try: 
     158        mi2 = statc.mean(list2) 
     159        st2 = statc.std(list2) 
     160    except: 
     161        pass 
     162     
     163    return mi1, st1, mi2, st2 
     164 
    128165class AT_edelmanParametricLearner(object): 
    129166    """ 
     
    140177 
    141178    def __call__(self, i, data): 
    142         cv = data.domain.classVar 
    143         #print data.domain 
     179        cv = data.domain.class_var 
    144180 
    145181        if self.a == None: self.a = cv.values[0] 
    146182        if self.b == None: self.b = cv.values[1] 
    147183 
    148         def avWCVal(value): 
    149             return [ex[i].value for ex in data if ex[-1].value == value and not ex[i].isSpecial() ] 
    150  
    151         list1 = avWCVal(self.a) 
    152         list2 = avWCVal(self.b) 
    153  
    154         mi1 = mi2 = st1 = st2 = None 
    155  
    156         try: 
    157             mi1 = statc.mean(list1) 
    158             st1 = statc.std(list1) 
    159         except: 
    160             pass 
    161      
    162         try: 
    163             mi2 = statc.mean(list2) 
    164             st2 = statc.std(list2) 
    165         except: 
    166             pass 
     184        mi1, st1, mi2, st2 = estimate_gaussian_per_class(data, i, a=self.a, b=self.b) 
    167185 
    168186        return AT_edelmanParametric(mi1=mi1, mi2=mi2, st1=st1, st2=st2) 
     
    336354 
    337355def mat_ni(data, matcher): 
     356    """ Return (in a tuple): 
     357        - a gene matcher that matches to the attribute names of data 
     358        - a dictionary attribute names -> indices in the data set. 
     359    """ 
    338360    nm = matcher([at.name for at in data.domain.attributes]) 
    339361    name_ind = dict((n.name,i) for i,n in enumerate(data.domain.attributes)) 
     
    616638        return at 
    617639 
     640def compute_llr(data, inds, cache): 
     641    """ 
     642    Compute CORG for this geneset specified with gene inds 
     643    in the example table. Output is the list of gene inds 
     644    in CORG. 
     645    """ 
     646    def gaussianc(data, at, cache=None): 
     647        """ Cached attribute  tscore calculation """ 
     648        if cache != None and at in cache: return cache[at] 
     649        ma = estimate_gaussian_per_class(data, at) 
     650        if cache != None: cache[at] = ma 
     651        return ma 
     652 
     653    gf = [ gaussianc(data, at, cache) for at in inds ] 
     654    return gf 
     655 
     656""" To avoid scipy overhead """ 
     657from math import pi 
     658_norm_pdf_C = math.sqrt(2*pi) 
     659_norm_pdf_logC = math.log(_norm_pdf_C) 
     660def _norm_logpdf(x, mi, std): 
     661    return -(x-mi)**2 / (2.0*std**2) - _norm_pdf_logC - math.log(std) 
     662 
     663def _llrlogratio(v, mi1, std1, mi2, std2): 
     664    if mi1 == None or std1 == None or mi2 == None or std2 == None: 
     665        return 0. #problem with estimation 
     666    #lpdf1 = scipy.stats.norm.logpdf(v, mi1, std1) 
     667    #lpdf2 = scipy.stats.norm.logpdf(v, mi2, std2) 
     668    lpdf1 = _norm_logpdf(v, mi1, std1) #avoids scipy's overhead 
     669    lpdf2 = _norm_logpdf(v, mi2, std2) 
     670    r = lpdf1 - lpdf2 
     671    return r 
     672 
     673class LLR(ParametrizedTransformation): 
     674    """  
     675    From Su et al: Accurate and Reliable Cancer Classification Base 
     676    on Probalistic Inference of Pathway Activity. Plos One, 2009. 
     677    """ 
     678 
     679    def __init__(self, **kwargs): 
     680        self.normalize = kwargs.pop("normalize", False) #normalize final results 
     681        super(LLR, self).__init__(**kwargs) 
     682 
     683    def __call__(self, *args, **kwargs): 
     684        self._gauss_cache = {} #caching of gaussian estimates 
     685        self._normalizec = {} 
     686        return super(LLR, self).__call__(*args, **kwargs) 
     687 
     688    def build_feature(self, data, gs): 
     689 
     690        at = Orange.feature.Continuous(name=str(gs)) 
     691        geneset = list(gs.genes) 
     692 
     693        nm, name_ind, genes, takegenes, to_geneset = self._match_data(data, geneset, odic=True) 
     694 
     695        gsi = [ name_ind[g] for g in genes ] 
     696        gausse = compute_llr(data, gsi, self._gauss_cache) 
     697        genes_gs = [ to_geneset[g] for g in genes ] 
     698 
     699        if self.normalize: # per (3) in the paper 
     700            #compute log ratios for all samples and genes from this gene set 
     701            for i, gene_gs, g in zip(gsi, genes_gs, gausse): 
     702                if gene_gs not in self._normalizec: #skip if computed already 
     703                    r = [ _llrlogratio(ex[i].value, *g) for ex in data ] 
     704                    self._normalizec[gene_gs] = (statc.mean(r), statc.std(r)) 
     705 
     706        def t(ex, w, genes_gs=genes_gs, gausse=gausse, normalizec=self._normalizec): 
     707            nm2, name_ind2, genes2 = self._match_instance(ex, genes_gs, None) 
     708            gsvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] 
     709 
     710            vals = [ _llrlogratio(v, *g) if v != "?" else 0.0 for v,g in zip(gsvalues, gausse) ] 
     711 
     712            if len(normalizec): #normalize according to (3) 
     713                vals2 = [] 
     714                for v,g in zip(vals, genes_gs): 
     715                    m,s = self._normalizec[g] 
     716                    vals2.append((v-m)/s) 
     717                vals = vals2 
     718             
     719            return sum(vals) 
     720 
     721      
     722        at.get_value_from = t 
     723        return at 
     724 
    618725if __name__ == "__main__": 
    619726 
Note: See TracChangeset for help on using the changeset viewer.