Changes in [1810:ad758982f7eb:1811:d89e4d232162] in orange-bioinformatics


Ignore:
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • _bioinformatics/obiGeneSetSig.py

    r1775 r1803  
    11from __future__ import absolute_import 
    22 
     3import random 
    34import math 
    45from collections import defaultdict 
     
    165166     
    166167    if common_if_extreme and extreme(): 
    167         print "extreme", st1, st2, 
    168168        st1 = st2 = statc.std(list1 + list2) 
    169         print "new", st1, st2 
    170169 
    171170    return mi1, st1, mi2, st2 
     
    623622    """ 
    624623 
    625     def __call__(self, *args, **kwargs): 
     624    def build_features(self, *args, **kwargs): 
    626625        self.tscorecache = {} #reset a cache 
    627         return super(CORGs, self).__call__(*args, **kwargs) 
     626        return super(CORGs, self).build_features(*args, **kwargs) 
    628627 
    629628    def build_feature(self, data, gs): 
     
    670669 
    671670def _llrlogratio(v, mi1, std1, mi2, std2): 
    672     if mi1 == None or std1 == None or mi2 == None or std2 == None: 
     671    if mi1 == None or std1 == None or mi2 == None or std2 == None or std1 == 0 or std2 == 0: 
    673672        return 0. #problem with estimation 
    674673    #lpdf1 = scipy.stats.norm.logpdf(v, mi1, std1) 
     
    686685 
    687686    def __init__(self, **kwargs): 
    688         self.normalize = kwargs.pop("normalize", False) #normalize final results 
     687        self.normalize = kwargs.pop("normalize", True) #normalize final results 
    689688        super(LLR, self).__init__(**kwargs) 
    690689 
    691     def __call__(self, *args, **kwargs): 
     690    def build_features(self, *args, **kwargs): 
    692691        self._gauss_cache = {} #caching of gaussian estimates 
    693692        self._normalizec = {} 
    694         return super(LLR, self).__call__(*args, **kwargs) 
     693        return super(LLR, self).build_features(*args, **kwargs) 
    695694 
    696695    def build_feature(self, data, gs): 
     
    721720                vals2 = [] 
    722721                for v,g in zip(vals, genes_gs): 
    723                     m,s = self._normalizec[g] 
    724                     vals2.append((v-m)/s) 
     722                    m,s = normalizec[g] 
     723                    if s == 0: #disregard attributes without differences 
     724                        vals2.append(0.) 
     725                    else: 
     726                        vals2.append((v-m)/s) 
    725727                vals = vals2 
    726728             
    727729            return sum(vals) 
    728  
    729730      
    730731        at.get_value_from = t 
    731732        return at 
    732733 
     734class LLR_slow(ParametrizedTransformation): 
     735    """ Slow and rough implementation of LLR (testing correctness).""" 
     736 
     737    def _get_par(self, datao): 
     738        gaussiane = [ estimate_gaussian_per_class(datao, at, common_if_extreme=True) for at in range(len(datao.domain.attributes)) ] 
     739        normalizec = [] 
     740        for i,g in zip(range(len(datao.domain.attributes)), gaussiane): 
     741            r = [ _llrlogratio(ex[i].value, *g) for ex in datao ] 
     742            normalizec.append((statc.mean(r), statc.std(r))) 
     743        return gaussiane, normalizec 
     744 
     745    def _use_par(self, arr, constructt): 
     746        gaussiane, normalizec = constructt 
     747        arr = [ arr[i].value for i in range(len(arr.domain.attributes)) ] 
     748        return sum ( (_llrlogratio(v, *g)-m)/s for v,g,n in zip(arr, gaussiane, normalizec)) 
     749 
     750 
     751def estimate_linear_fit(data, i): 
     752    """ 
     753    Chen et 2008 write about t-score of the least square fit in the 
     754    original article, but here we can use just the t-test, because the 
     755    t-scores obtained are exactly the same (the authors used the 
     756    the more general version as they supported continous outcomes). 
     757    """ 
     758     
     759    """ 
     760    #implementation that could also support continous outcomes 
     761    x = numpy.array([ ex[i].value for ex in data ]) 
     762    y = numpy.array([ int(ex[-1]) for ex in data ]) #here classes are 0 and 1 
     763    ret = scipy.stats.linregress(x,y) 
     764    b = ret[0] 
     765    seb = ret[4] 
     766    return b/seb 
     767    """ 
     768    """ 
     769    #per mathword.wolfram.com - results are the same 
     770    c = numpy.cov(x,y) 
     771    n = len(x) 
     772    b = c[0,1]/c[0,0] #the same result as from li 
     773    s = math.sqrt((c[1,1]*n - b*c[0,1]*n)/(n-2)) 
     774    seb = s/math.sqrt(c[0,0]*n) 
     775    return b/seb 
     776    """ 
     777    return obiExpression.MA_t_test()(i, data) 
     778 
     779 
     780class SPCA(ParametrizedTransformation): 
     781    """ Per Chen et al. Supervised principal component analysis for 
     782    gene set enrichment of microarray data with continuous or survival 
     783    outcomes. Bioinformatics, 2008.  """ 
     784 
     785    def __init__(self, **kwargs): 
     786        self.threshold = kwargs.pop("threshold", None) 
     787        self.top = kwargs.pop("top", None) 
     788        self.atleast = kwargs.pop("atleast", 0) 
     789        super(SPCA, self).__init__(**kwargs) 
     790 
     791    def _get_par(self, datao): 
     792        scores = [ estimate_linear_fit(datao, a) for a in range(len(datao.domain.attributes)) ] 
     793        scores = list(enumerate(map(abs, scores))) 
     794        select = None 
     795        if self.threshold is not None: 
     796            select = [ i for i,s in scores if s > self.threshold ] 
     797        elif self.top is not None: 
     798            select = nth(sorted(scores, key=lambda x: -x[1])[:self.top], 0) 
     799 
     800        if len(select) < self.atleast: 
     801            select = nth(sorted(scores, key=lambda x: -x[1])[:self.atleast], 0) 
     802 
     803        if select == None: 
     804            somethingWrongWithSelection 
     805 
     806        doms = Orange.data.Domain([ datao.domain.attributes[i] for i in select ], datao.domain.class_var) 
     807        datas = Orange.data.Table(doms, datao) 
     808 
     809        if len(select) == 0: 
     810            return select, None 
     811        else: 
     812            return set(select), pca(datas) 
     813 
     814    def _use_par(self, arr, constructt): 
     815        select, constructt = constructt 
     816 
     817        if len(select) == 0: 
     818            return 0. 
     819 
     820        arr = [ arr[i].value for i in range(len(arr.domain.attributes))  
     821            if i in select ] 
     822 
     823        evals, evect, xmean = constructt 
     824 
     825        arr = arr - xmean # same input transformation - a row in a matrix 
     826        ev0 = evect[0] #this is a row in a matrix - do a dot product 
     827        a = numpy.dot(arr, ev0) 
     828 
     829        return a 
     830 
     831def _shuffleClass(data, rand): 
     832    """ Destructive! """ 
     833    locations = range(len(data)) 
     834    rand.shuffle(locations) 
     835    attribute = -1 
     836    l = [None]*len(data) 
     837    for i in range(len(data)): 
     838        l[locations[i]] = data[i][attribute] 
     839    for i in range(len(data)): 
     840        data[i][attribute] = l[i] 
     841 
     842class SPCA_ttperm(SPCA): 
     843    """ Set threshold with a permutation test. """ 
     844 
     845    def __init__(self, **kwargs): 
     846        self.pval = kwargs.pop("pval", 0.01) #target p-value 
     847        self.perm = kwargs.pop("perm", 100) #number of class permutation 
     848        self.sperm = kwargs.pop("sperm", 100) #sampled attributes per permutation 
     849        super(SPCA_ttperm, self).__init__(**kwargs) 
     850 
     851    def build_features(self, data, *args, **kwargs): 
     852        joined = [] 
     853        rand = random.Random(0) 
     854        nat = len(data.domain.attributes) 
     855 
     856        datap = Orange.data.Table(data.domain, data) 
     857 
     858        for p in range(self.perm): 
     859            _shuffleClass(datap, rand) 
     860            if self.sperm is not None: 
     861                ti = rand.sample(xrange(nat), self.sperm) 
     862            else: 
     863                ti = range(nat) 
     864            joined.extend([ obiExpression.MA_t_test()(i, datap)  
     865                for i in ti ]) 
     866 
     867        joined = map(abs, joined) 
     868        joined.sort(reverse=True) 
     869 
     870        t = joined[int(self.pval*len(joined))] 
     871 
     872        self.threshold = t 
     873        return super(SPCA_ttperm, self).build_features(data, *args, **kwargs) 
     874     
     875 
    733876if __name__ == "__main__": 
    734877 
    735878    data = Orange.data.Table("iris") 
     879 
    736880    gsets = obiGeneSets.collections({ 
    737881        #"ALL": ['sepal length', 'sepal width', 'petal length', 'petal width'], 
     
    761905        print '\n'.join([ a + ": " +str(b) for a,b in ol]) 
    762906 
    763     ass = Assess(data, matcher=matcher, gene_sets=gsets, class_values=choosen_cv, min_part=0.0) 
     907    ass = LLR(data, matcher=matcher, gene_sets=gsets, class_values=choosen_cv, min_part=0.0, normalize=True) 
     908    #ass = LLR_slow(data, matcher=matcher, gene_sets=gsets, class_values=choosen_cv, min_part=0.0) 
    764909    ar = to_old_dic(ass.domain, data[:5]) 
    765910    pp2(ar) 
  • _bioinformatics/obiOMIM.py

    r1797 r1799  
    88 
    99class disease(object): 
    10         """ A class representing a disease in the OMIM database 
    11         """ 
    12         regex = re.compile(r'(?P<name>.*?),? (?P<id>[0-9]{3,6} )?(?P<m1>\([123?]\) )?(?P<m2>\([123?]\) )? *$') 
    13         __slots__ = ["name", "id", "mapping"] 
    14         def __init__(self, morbidmap_line): 
    15             string = morbidmap_line.split("|", 1)[0] 
    16             match = self.regex.match(string) 
    17     #        print string 
    18     #        print match.groups() 
    19             self.name, self.id, self.mapping = [s.strip() if s else s for s in match.groups()[:3]] 
    20             if match.group("m2"): 
    21                 self.mapping += " " + match.group("m2").strip() 
    22                                                                                      
     10    """ A class representing a disease in the OMIM database 
     11    """ 
     12    regex = re.compile(r'(?P<name>.*?),? (?P<id>[0-9]{3,6} )?(?P<m1>\([123?]\) )?(?P<m2>\([123?]\) )? *$') 
     13    __slots__ = ["name", "id", "mapping"] 
     14    def __init__(self, morbidmap_line): 
     15        string = morbidmap_line.split("|", 1)[0] 
     16        match = self.regex.match(string) 
     17#        print string 
     18#        print match.groups() 
     19        self.name, self.id, self.mapping = [s.strip() if s else s for s in match.groups()[:3]] 
     20        if match.group("m2"): 
     21            self.mapping += " " + match.group("m2").strip() 
     22                                                                                 
    2323class OMIM(object): 
    2424    VERSION = 1 
     
    3333 
    3434        self.load(filename) 
    35                                                                                                                                                                                      
     35     
    3636    @classmethod 
    3737    def download_from_NCBI(cls, file=None): 
     
    4141        shutil.copyfileobj(stream, file, length=10) 
    4242        file.close() 
    43                                                                                                                                                                                                                                                          
     43 
    4444    @classmethod 
    4545    def get_instance(cls): 
     
    5050        instance.__dict__ = cls._shared_dict 
    5151        return instance  
    52                                                                                                                                                                                                                                                                                                                                      
     52     
    5353    def load(self, filename): 
    5454        file = open(filename, "rb") 
    5555        lines = file.read().splitlines() 
    5656        self._disease_dict = dict([(disease(line), line) for line in lines if line]) 
    57                                                                                                                                                                                                                                                                                                                                                                      
     57     
    5858    def diseases(self): 
    5959        return self._disease_dict.keys() 
    60                                                                                                                                                                                                                                                                                                                                                                                      
     60     
    6161    def genes(self): 
    6262        return sorted(set(reduce(list.__add__, [self.disease_genes(disease) for disease in self.diseases()], []))) 
    63                                                                                                                                                                                                                                                                                                                                                                                                  
     63     
    6464    def disease_genes(self, disease): 
    6565        return self._disease_dict[disease].split("|")[1].split(", ") 
    66                                                                                                                                                                                                                                                                                                                                                                                                              
     66     
    6767    def gene_diseases(self): 
    6868        d = defaultdict(set) 
     
    7171                d[gene].add(disease) 
    7272        return d 
    73                                                                                                                                                                                                                                                                                                                                                                                                                                                                              
     73 
    7474def diseases(): 
    7575    """ Return all disease objects 
    7676    """ 
    7777    return OMIM.get_instance().diseases() 
    78                                                                                                                                                                                                                                                                                                                                                                                                                                                                                              
     78 
    7979def genes(): 
    8080    """ Return a set of all genes referenced in OMIM  
  • setup.py

    r1733 r1800  
    1616DOCUMENTATION_NAME = 'Orange Bioinformatics' 
    1717 
    18 VERSION = '2.5a7' 
     18VERSION = '2.5a8' 
    1919 
    2020DESCRIPTION = 'Orange Bioinformatics add-on for Orange data mining software package.' 
Note: See TracChangeset for help on using the changeset viewer.