Changes in [1798:f1e73901a1d3:1799:2f952a765716] in orange-bioinformatics


Ignore:
Location:
_bioinformatics
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • _bioinformatics/obiGeneSetSig.py

    r1775 r1786  
    11from __future__ import absolute_import 
    22 
     3import random 
    34import math 
    45from collections import defaultdict 
     
    731732        return at 
    732733 
     734def estimate_linear_fit(data, i): 
     735    """ 
     736    Chen et 2008 write about t-score of the least square fit in the 
     737    original article, but here we can use just the t-test, because the 
     738    t-scores obtained are exactly the same (the authors used the 
     739    the more general version as they supported continous outcomes). 
     740    """ 
     741     
     742    """ 
     743    #implementation that could also support continous outcomes 
     744    x = numpy.array([ ex[i].value for ex in data ]) 
     745    y = numpy.array([ int(ex[-1]) for ex in data ]) #here classes are 0 and 1 
     746    ret = scipy.stats.linregress(x,y) 
     747    b = ret[0] 
     748    seb = ret[4] 
     749    return b/seb 
     750    """ 
     751    """ 
     752    #per mathword.wolfram.com - results are the same 
     753    c = numpy.cov(x,y) 
     754    n = len(x) 
     755    b = c[0,1]/c[0,0] #the same result as from li 
     756    s = math.sqrt((c[1,1]*n - b*c[0,1]*n)/(n-2)) 
     757    seb = s/math.sqrt(c[0,0]*n) 
     758    return b/seb 
     759    """ 
     760    return obiExpression.MA_t_test()(i, data) 
     761 
     762 
     763class SPCA(ParametrizedTransformation): 
     764    """ Per Chen et al. Supervised principal component analysis for 
     765    gene set enrichment of microarray data with continuous or survival 
     766    outcomes. Bioinformatics, 2008.  """ 
     767 
     768    def __init__(self, **kwargs): 
     769        self.threshold = kwargs.pop("threshold", None) 
     770        self.top = kwargs.pop("top", None) 
     771        self.atleast = kwargs.pop("atleast", 0) 
     772        super(SPCA, self).__init__(**kwargs) 
     773 
     774    def _get_par(self, datao): 
     775        scores = [ estimate_linear_fit(datao, a) for a in range(len(datao.domain.attributes)) ] 
     776        scores = list(enumerate(map(abs, scores))) 
     777        select = None 
     778        if self.threshold is not None: 
     779            select = [ i for i,s in scores if s > self.threshold ] 
     780        elif self.top is not None: 
     781            select = nth(sorted(scores, key=lambda x: -x[1])[:self.top], 0) 
     782 
     783        if len(select) < self.atleast: 
     784            select = nth(sorted(scores, key=lambda x: -x[1])[:self.atleast], 0) 
     785 
     786        if select == None: 
     787            somethingWrongWithSelection 
     788 
     789        doms = Orange.data.Domain([ datao.domain.attributes[i] for i in select ], datao.domain.class_var) 
     790        datas = Orange.data.Table(doms, datao) 
     791 
     792        if len(select) == 0: 
     793            return select, None 
     794        else: 
     795            return set(select), pca(datas) 
     796 
     797    def _use_par(self, arr, constructt): 
     798        select, constructt = constructt 
     799 
     800        if len(select) == 0: 
     801            return 0. 
     802 
     803        arr = [ arr[i].value for i in range(len(arr.domain.attributes))  
     804            if i in select ] 
     805 
     806        evals, evect, xmean = constructt 
     807 
     808        arr = arr - xmean # same input transformation - a row in a matrix 
     809        ev0 = evect[0] #this is a row in a matrix - do a dot product 
     810        a = numpy.dot(arr, ev0) 
     811 
     812        return a 
     813 
     814def _shuffleClass(data, rand): 
     815    """ Destructive! """ 
     816    locations = range(len(data)) 
     817    rand.shuffle(locations) 
     818    attribute = -1 
     819    l = [None]*len(data) 
     820    for i in range(len(data)): 
     821        l[locations[i]] = data[i][attribute] 
     822    for i in range(len(data)): 
     823        data[i][attribute] = l[i] 
     824 
     825class SPCA_ttperm(SPCA): 
     826    """ Set threshold with a permutation test. """ 
     827 
     828    def __init__(self, **kwargs): 
     829        self.pval = kwargs.pop("pval", 0.01) #target p-value 
     830        self.perm = kwargs.pop("perm", 100) #number of class permutation 
     831        self.sperm = kwargs.pop("sperm", 100) #sampled attributes per permutation 
     832        super(SPCA_ttperm, self).__init__(**kwargs) 
     833 
     834    def __call__(self, data, *args, **kwargs): 
     835        joined = [] 
     836        rand = random.Random(0) 
     837        nat = len(data.domain.attributes) 
     838 
     839        datap = Orange.data.Table(data.domain, data) 
     840 
     841        for p in range(self.perm): 
     842            _shuffleClass(datap, rand) 
     843            if self.sperm is not None: 
     844                ti = rand.sample(xrange(nat), self.sperm) 
     845            else: 
     846                ti = range(nat) 
     847            joined.extend([ obiExpression.MA_t_test()(i, datap)  
     848                for i in ti ]) 
     849 
     850        joined = map(abs, joined) 
     851        joined.sort(reverse=True) 
     852 
     853        t = joined[int(self.pval*len(joined))] 
     854 
     855        self.threshold = t 
     856        return super(SPCA_ttperm, self).__call__(data, *args, **kwargs) 
     857     
     858 
    733859if __name__ == "__main__": 
    734860 
    735861    data = Orange.data.Table("iris") 
     862 
    736863    gsets = obiGeneSets.collections({ 
    737864        #"ALL": ['sepal length', 'sepal width', 'petal length', 'petal width'], 
     
    761888        print '\n'.join([ a + ": " +str(b) for a,b in ol]) 
    762889 
    763     ass = Assess(data, matcher=matcher, gene_sets=gsets, class_values=choosen_cv, min_part=0.0) 
     890    ass = SPCA(data, matcher=matcher, gene_sets=gsets, class_values=choosen_cv, min_part=0.0, top=0) 
    764891    ar = to_old_dic(ass.domain, data[:5]) 
    765892    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  
Note: See TracChangeset for help on using the changeset viewer.