Changeset 414:d7764b1f5cbc in orange-bioinformatics


Ignore:
Timestamp:
07/07/08 22:30:57 (6 years ago)
Author:
markotoplak
Branch:
default
Convert:
33142236615f18cb3b4acd6083a8f0f0a0a9cbe4
Message:

GSEA remove suffix of geneset files.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • obiGsea.py

    r411 r414  
    4141    return lambda d: [ meas(i,d) for i in range(len(d.domain.attributes)) ] 
    4242 
    43 def lvar (inlist): 
    44     n = inlist.size 
    45     mn = numpy.mean(inlist) 
    46     deviations = inlist - mn 
    47     return numpy.dot(deviations,deviations)/float(n-1) 
    48  
    49 def stdev(l): 
    50     return math.sqrt(lvar(l)) 
    51  
    52 def stdevm(l): 
    53     m = numpy.mean(l) 
    54     std = stdev(l) 
    55     #print std, 0.2*abs(1.0 if m == 0 else m) 
    56     #return minmally 2*|mi|, where mi=0 is adjusted to mi=1 
    57     return max(std, 0.2*abs(1.0 if m == 0 else m)) 
    58  
    59 class MA_signalToNoise2: 
    60     """ 
    61     Returns signal to noise measurement: difference of means of two classes 
    62     divided by the sum of standard deviations for both classes.  
    63     """ 
    64  
    65     def __init__(self, a=None, b=None): 
    66         """ 
    67         a and b are choosen class values. 
    68         """ 
    69         self.a = a 
    70         self.b = b 
    71  
    72     def __call__(self, i, data): 
    73         cv = data.domain.classVar 
    74         #print data.domain 
    75  
    76         datao = data 
    77  
    78         #for faster computation. to save dragging many attributes along 
    79         dom2 = orange.Domain([data.domain.attributes[i]], data.domain.classVar) 
    80         data = orange.ExampleTable(dom2, data) 
    81         i = 0 
    82  
    83         if self.a == None: self.a = 0 
    84         if self.b == None: self.b = 1 
    85  
    86         a,c = data.toNumpyMA("A/C") 
    87         a = a[:,0] 
    88         #print a.data 
    89         masked = numpy.where(a.mask==False) 
    90         a = a.data 
    91         c = c.data 
    92         a = a[masked] 
    93         c = c[masked] 
    94  
    95         def avWCValNP(value): 
    96             return  a[numpy.where(c == value)] 
    97  
    98         exa = avWCValNP(0) 
    99         exb = avWCValNP(1) 
    100  
    101         try: 
    102             rval = (numpy.mean(exa)-numpy.mean(exb))/(stdevm(exa)+stdevm(exb)) 
    103             return rval 
    104         except: 
    105             #return some "middle" value 
    106             return 0 
    107  
    10843def orderedPointersCorr(lcor): 
    10944    """ 
     
    11752    return ordered 
    11853 
    119 def enrichmentScoreRankedOld(subset, lcor, ordered, p=1.0): 
     54def enrichmentScoreRanked(subset, lcor, ordered, p=1.0, rev2=None): 
    12055    """ 
    12156    Input data and subset.  
     
    12661 
    12762    Returns enrichment score on given data. 
    128     """ 
    129  
    130     #preglej ponovitve 
    131     #if len(subset) != len(set(subset)): print "CUDNO", subset  
    132  
    133     subset = set(subset) 
    134  
    135     #add if gene is not in the subset 
    136     notInA = -(1. / (len(lcor)-len(subset))) 
    137     #base for addition if gene is in the subset 
    138     cors = [ abs(lcor[i])**p for i in subset ] 
    139     sumcors = sum(cors) 
    140  
    141     #this should not happen 
    142     if sumcors == 0.0: 
    143         return (0.0, None) 
    144      
    145     inAb = 1./sumcors 
    146  
    147     ess = [0.0] 
    148      
    149     for i in ordered: 
    150         ess.append(ess[-1] + \ 
    151             (inAb*abs(lcor[i]**p) if i in subset else notInA) 
    152         ) 
    153  
    154     maxEs = max(ess) 
    155     minEs = min(ess) 
    156      
    157     return (maxEs if abs(maxEs) > abs(minEs) else minEs, ess[1:]) 
    158  
    159  
    160 def enrichmentScoreRanked(subset, lcor, ordered, p=1.0, rev2=None): 
    161     """ 
    162     Input data and subset.  
    163      
    164     subset: list of attribute indices of the input data belonging 
    165         to the same set. 
    166     lcor: correlations with class for each attribute in a list.  
    167  
    168     Returns enrichment score on given data. 
    169     """ 
     63 
     64    This implementation efficiently handles "sparse" genesets (that 
     65    cover only a small subset of all genes in the dataset). 
     66    """ 
     67 
     68    #print lcor 
    17069 
    17170    subset = set(subset) 
     
    243142    print "REAL", (maxEs if abs(maxEs) > abs(minEs) else minEs, ess[1:]) 
    244143 
    245     aaaaa = Aaaa 
    246144    """ 
    247145    return (maxSum if abs(maxSum) > abs(minSum) else minSum, []) 
     
    301199    or negative portion of the distribution corresponding to the sign  
    302200    of the observed ES(S). 
    303  
    304     WHAT DOES THAT MEAN? 
    305201    """ 
    306202     
     
    311207        else:  
    312208            return float(len([ a for a in esnull if a >= es ]))/ \ 
    313                 len([ a for a in esnull if a > 0]) 
     209                len([ a for a in esnull if a >= 0]) 
    314210    except: 
    315211        return 1.0 
     
    346242  
    347243    lcor = rankingf(data) 
     244    #print lcor 
     245 
    348246    ordered = orderedPointersCorr(lcor) 
    349247 
     
    436334def gseaSignificance(enrichmentScores, enrichmentNulls): 
    437335 
     336    #print enrichmentScores 
     337 
    438338    enrichmentPVals = [] 
    439339    nEnrichmentScores = [] 
     
    443343        es = enrichmentScores[i] 
    444344        enrNull = enrichmentNulls[i] 
     345        #print es, enrNull 
    445346 
    446347        enrichmentPVals.append(gseapval(es, enrNull)) 
     
    894795        fn = os.path.join(path, file) 
    895796        name = info.get(file, file) 
     797        if name == file: 
     798            #remove suffix if no info 
     799            if name.lower()[-4:] == ".gmt" or name.lower()[-4:] == ".pck": 
     800                name = name[:-4] 
     801 
    896802        out.append( (name, fn) ) 
    897803 
     
    941847""" 
    942848 
     849 
     850 
     851def etForAttribute(datal,a): 
     852    tables = len(datal) 
     853 
     854    def getAttrVals(data, attr): 
     855        dom2 = orange.Domain([data.domain[attr]], False) 
     856        dataa = orange.ExampleTable(dom2, data) 
     857        return [ a[0].native() for a in dataa ] 
     858 
     859    domainl = [] 
     860    valuesl = [] 
     861 
     862    for id, data in enumerate(datal): 
     863        v = getAttrVals(data,a) 
     864        valuesl.append(v) 
     865        domainl.append(orange.FloatVariable(name=("v"+str(id)))) 
     866 
     867    classvals = getAttrVals(data, datal[0].domain.classVar) 
     868    valuesl += [ classvals ] 
     869 
     870    dom = orange.Domain(domainl, datal[0].domain.classVar) 
     871    examples = [ list(a) for a in zip(*valuesl) ] 
     872 
     873    datat = orange.ExampleTable(dom, examples) 
     874 
     875    return datat 
     876 
     877 
     878def evaluateEtWith(fn): 
     879    """ 
     880    fn - evaluates example talbe 
     881    """ 
     882 
     883    def newf(datal): 
     884        res = [] 
     885        for a in datal[0].domain.attributes: 
     886            et = etForAttribute(datal, a) 
     887            res.append(fn(et)) 
     888        return res 
     889 
     890    return newf 
     891 
     892 
    943893if  __name__=="__main__": 
    944894 
     
    976926     
    977927    #data = orange.ExampleTable("leukemiaGSEA.tab") 
    978     data = orange.ExampleTable("demo.tab") 
     928    #data = orange.ExampleTable("demo.tab") 
     929    data = orange.ExampleTable("sterolTalkHepa.tab") 
     930 
     931    print data.domain.classVar.values 
     932 
    979933    print "loaded data" 
    980934 
     
    983937        #print "done" 
    984938 
    985     def etForAttribute(datal,a): 
    986         tables = len(datal) 
    987  
    988         def getAttrVals(data, attr): 
    989             dom2 = orange.Domain([data.domain[attr]], False) 
    990             dataa = orange.ExampleTable(dom2, data) 
    991             return [ a[0].native() for a in dataa ] 
    992  
    993         domainl = [] 
    994         valuesl = [] 
    995  
    996         for id, data in enumerate(datal): 
    997             v = getAttrVals(data,a) 
    998             valuesl.append(v) 
    999             domainl.append(orange.FloatVariable(name=("v"+str(id)))) 
    1000  
    1001         classvals = getAttrVals(data, datal[0].domain.classVar) 
    1002         valuesl += [ classvals ] 
    1003  
    1004         dom = orange.Domain(domainl, datal[0].domain.classVar) 
    1005         examples = [ list(a) for a in zip(*valuesl) ] 
    1006  
    1007         datat = orange.ExampleTable(dom, examples) 
    1008  
    1009         return datat 
    1010  
    1011  
    1012     def evaluateWith(fn): 
    1013         """ 
    1014         fn - evaluates example talbe 
    1015         """ 
    1016  
    1017         def newf(datal): 
    1018             res = [] 
    1019             for a in datal[0].domain.attributes: 
    1020                 et = etForAttribute(datal, a) 
    1021                 res.append(fn(et)) 
    1022             return res 
    1023  
    1024         return newf 
    1025      
    1026939    def give1(data): 
    1027940        return 1 
     
    1032945 
    1033946    def knn(data): 
    1034         learner = orange.kNNLearner(k=int(len(data))) 
     947        learner = orange.kNNLearner(k=int(sqrt(len(data)))) 
    1035948        res = orngTest.crossValidation([ learner ], data) 
    1036         return orngStat.AUC(res)[0] 
     949        return orngStat.AUC(res)[0]**3 
    1037950 
    1038951 
    1039952    gen1 = collections(["c2.all.v2.5.symbols.gmt"], default=False) 
     953    print gen1.keys()[:100] 
     954 
     955    import sys; sys.exit(0) 
    1040956 
    1041957    #data = orange.ExampleTable(data.domain, data[:1]) 
    1042958 
    1043     #res2 = runGSEA(data, n=5, geneSets=gen1, permutation="class", callback=novi, atLeast=3) 
    1044     res2 = runGSEA([data], n=100, geneSets=gen1, permutation="class", callback=novi, atLeast=3, rankingf=evaluateWith(knn)) 
    1045      
    1046     print '\n'.join([ str(a) + ": " +str(b[2]) for a,b in sorted(res2.items(), key=lambda x: x[1][2])]) 
    1047  
    1048  
    1049  
    1050  
     959    #res2 = runGSEA(data, n=100, geneSets=gen1, permutation="class", callback=novi, atLeast=3) 
     960    res2 = runGSEA([data, data], n=100, geneSets=gen1, permutation="class", callback=novi, atLeast=3, rankingf=evaluateEtWith(knn)) 
     961     
     962    print '\n'.join([ str(a) + ": " +str(b[:3]) for a,b in sorted(res2.items(), key=lambda x: x[1][2])]) 
     963 
     964 
     965 
     966 
Note: See TracChangeset for help on using the changeset viewer.