Changes in [1778:a94aec22f398:1779:736a0fc7fccf] in orange-bioinformatics


Ignore:
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • _bioinformatics/obiGeneSetSig.py

    r1728 r1775  
    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 ] 
     
    7483            # while the transformed data table uses cross-validation 
    7584            # internally 
    76             folds = 5 
    77             cvi = Orange.data.sample.SubsetIndicesCV(data, folds) 
     85            if self.cv == True: 
     86                cvi = Orange.data.sample.SubsetIndicesCV(data, 5) 
     87            elif self.cv != False: 
     88                cvi = self.cv(data) 
    7889            data_cv = [ [] for _ in range(len(data)) ] 
    79             for f in range(folds): 
     90            for f in set(cvi): 
    8091                learn = data.select(cvi, f, negate=True) 
    8192                test = data.select(cvi, f) 
     
    124135            return 0 
    125136 
     137def estimate_gaussian_per_class(data, i, a=None, b=None, common_if_extreme=False): 
     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    def extreme(): 
     164        return st1 == 0 or st2 == 0 
     165     
     166    if common_if_extreme and extreme(): 
     167        print "extreme", st1, st2, 
     168        st1 = st2 = statc.std(list1 + list2) 
     169        print "new", st1, st2 
     170 
     171    return mi1, st1, mi2, st2 
     172 
    126173class AT_edelmanParametricLearner(object): 
    127174    """ 
     
    138185 
    139186    def __call__(self, i, data): 
    140         cv = data.domain.classVar 
    141         #print data.domain 
     187        cv = data.domain.class_var 
    142188 
    143189        if self.a == None: self.a = cv.values[0] 
    144190        if self.b == None: self.b = cv.values[1] 
    145191 
    146         def avWCVal(value): 
    147             return [ex[i].value for ex in data if ex[-1].value == value and not ex[i].isSpecial() ] 
    148  
    149         list1 = avWCVal(self.a) 
    150         list2 = avWCVal(self.b) 
    151  
    152         mi1 = mi2 = st1 = st2 = None 
    153  
    154         try: 
    155             mi1 = statc.mean(list1) 
    156             st1 = statc.std(list1) 
    157         except: 
    158             pass 
    159      
    160         try: 
    161             mi2 = statc.mean(list2) 
    162             st2 = statc.std(list2) 
    163         except: 
    164             pass 
     192        mi1, st1, mi2, st2 = estimate_gaussian_per_class(data, i, a=self.a, b=self.b) 
    165193 
    166194        return AT_edelmanParametric(mi1=mi1, mi2=mi2, st1=st1, st2=st2) 
     
    289317        return attributes 
    290318    
    291 def setSig_example_geneset(ex, data, no_unknowns): 
     319def setSig_example_geneset(ex, data, no_unknowns, check_same=False): 
    292320    """ Gets learning data and example with the same domain, both 
    293321    containing only genes from the gene set. """ 
     
    298326        vals1 = ex1.native(0)[:-1] 
    299327        vals2 = ex2.native(0)[:-1] 
     328 
     329        if check_same and vals1 == vals2: 
     330            return 10 #they are the same 
    300331 
    301332        #leaves undefined elements out 
     
    324355    #create distances to all learning data - save or other class 
    325356    for c in data: 
    326         distances[classValueMap[c[-1].value]].append(pearson(c, ex)) 
     357        p = pearson(c, ex) 
     358        if p != 10: 
     359             distances[classValueMap[c[-1].value]].append(pearson(c, ex)) 
    327360 
    328361    return ttest(distances[0], distances[1]) 
    329362 
    330363def mat_ni(data, matcher): 
     364    """ Return (in a tuple): 
     365        - a gene matcher that matches to the attribute names of data 
     366        - a dictionary attribute names -> indices in the data set. 
     367    """ 
    331368    nm = matcher([at.name for at in data.domain.attributes]) 
    332369    name_ind = dict((n.name,i) for i,n in enumerate(data.domain.attributes)) 
     
    358395    def __init__(self, **kwargs): 
    359396        self.no_unknowns = kwargs.pop("no_unknowns", False) 
     397        self.check_same = kwargs.pop("check_same", False) 
    360398        super(SetSig, self).__init__(**kwargs) 
    361399 
     
    378416            example = Orange.data.Instance(domain, exvalues) 
    379417 
    380             return setSig_example_geneset(example, datao, self.no_unknowns) #only this one is setsig specific 
     418            return setSig_example_geneset(example, datao, self.no_unknowns, check_same=self.check_same) #only this one is setsig specific 
    381419      
    382420        at.get_value_from = t 
     
    412450 
    413451            return self._use_par(ex, constructt) 
    414              
     452         
    415453        at.get_value_from = t 
     454        at.dbg = constructt #for debugging 
     455         
    416456        return at 
    417457 
     
    569609    for a in range(2, len(sortedinds)+1): 
    570610        tg = S(sortedinds[:a]) 
    571         if tg > g: 
     611        if tg > g: #improvement 
    572612            g = tg 
    573613            bg = a 
     
    575615            break 
    576616         
    577     return sortedinds[:a] 
     617    return sortedinds[:bg] #FIXED - one too many was taken 
    578618 
    579619class CORGs(ParametrizedTransformation): 
     
    606646        return at 
    607647 
     648def compute_llr(data, inds, cache): 
     649    """ 
     650    Compute CORG for this geneset specified with gene inds 
     651    in the example table. Output is the list of gene inds 
     652    in CORG. 
     653    """ 
     654    def gaussianc(data, at, cache=None): 
     655        """ Cached attribute  tscore calculation """ 
     656        if cache != None and at in cache: return cache[at] 
     657        ma = estimate_gaussian_per_class(data, at, common_if_extreme=True) 
     658        if cache != None: cache[at] = ma 
     659        return ma 
     660 
     661    gf = [ gaussianc(data, at, cache) for at in inds ] 
     662    return gf 
     663 
     664""" To avoid scipy overhead """ 
     665from math import pi 
     666_norm_pdf_C = math.sqrt(2*pi) 
     667_norm_pdf_logC = math.log(_norm_pdf_C) 
     668def _norm_logpdf(x, mi, std): 
     669    return -(x-mi)**2 / (2.0*std**2) - _norm_pdf_logC - math.log(std) 
     670 
     671def _llrlogratio(v, mi1, std1, mi2, std2): 
     672    if mi1 == None or std1 == None or mi2 == None or std2 == None: 
     673        return 0. #problem with estimation 
     674    #lpdf1 = scipy.stats.norm.logpdf(v, mi1, std1) 
     675    #lpdf2 = scipy.stats.norm.logpdf(v, mi2, std2) 
     676    lpdf1 = _norm_logpdf(v, mi1, std1) #avoids scipy's overhead 
     677    lpdf2 = _norm_logpdf(v, mi2, std2) 
     678    r = lpdf1 - lpdf2 
     679    return r 
     680 
     681class LLR(ParametrizedTransformation): 
     682    """  
     683    From Su et al: Accurate and Reliable Cancer Classification Base 
     684    on Probalistic Inference of Pathway Activity. Plos One, 2009. 
     685    """ 
     686 
     687    def __init__(self, **kwargs): 
     688        self.normalize = kwargs.pop("normalize", False) #normalize final results 
     689        super(LLR, self).__init__(**kwargs) 
     690 
     691    def __call__(self, *args, **kwargs): 
     692        self._gauss_cache = {} #caching of gaussian estimates 
     693        self._normalizec = {} 
     694        return super(LLR, self).__call__(*args, **kwargs) 
     695 
     696    def build_feature(self, data, gs): 
     697 
     698        at = Orange.feature.Continuous(name=str(gs)) 
     699        geneset = list(gs.genes) 
     700 
     701        nm, name_ind, genes, takegenes, to_geneset = self._match_data(data, geneset, odic=True) 
     702 
     703        gsi = [ name_ind[g] for g in genes ] 
     704        gausse = compute_llr(data, gsi, self._gauss_cache) 
     705        genes_gs = [ to_geneset[g] for g in genes ] 
     706 
     707        if self.normalize: # per (3) in the paper 
     708            #compute log ratios for all samples and genes from this gene set 
     709            for i, gene_gs, g in zip(gsi, genes_gs, gausse): 
     710                if gene_gs not in self._normalizec: #skip if computed already 
     711                    r = [ _llrlogratio(ex[i].value, *g) for ex in data ] 
     712                    self._normalizec[gene_gs] = (statc.mean(r), statc.std(r)) 
     713 
     714        def t(ex, w, genes_gs=genes_gs, gausse=gausse, normalizec=self._normalizec): 
     715            nm2, name_ind2, genes2 = self._match_instance(ex, genes_gs, None) 
     716            gsvalues = [ vou(ex, gn, name_ind2) for gn in genes2 ] 
     717 
     718            vals = [ _llrlogratio(v, *g) if v != "?" else 0.0 for v,g in zip(gsvalues, gausse) ] 
     719 
     720            if len(normalizec): #normalize according to (3) 
     721                vals2 = [] 
     722                for v,g in zip(vals, genes_gs): 
     723                    m,s = self._normalizec[g] 
     724                    vals2.append((v-m)/s) 
     725                vals = vals2 
     726             
     727            return sum(vals) 
     728 
     729      
     730        at.get_value_from = t 
     731        return at 
     732 
    608733if __name__ == "__main__": 
    609734 
  • server_update/updatemiRNA.py

    r1722 r1773  
    2626        sendMail('Uncorrect format of miRBase data-file.')         
    2727        return False 
     28 
     29IDpat = re.compile('ID\s*(\S*)\s*standard;') 
     30ACpat = re.compile('AC\s*(\S*);') 
     31RXpat = re.compile('RX\s*PUBMED;\s(\d*).') 
     32FT1pat = re.compile('FT\s*miRNA\s*(\d{1,}\.\.\d{1,})') 
     33FT2pat = re.compile('FT\s*/accession="(MIMAT[0-9]*)"') 
     34FT3pat = re.compile('FT\s*/product="(\S*)"') 
     35SQpat = re.compile('SQ\s*(.*other;)') 
     36seqpat = re.compile('\s*([a-z\s]*)\s*\d*') 
     37 
    2838 
    2939     
     
    97107                 
    98108                if r[0:2] == 'ID': 
    99                     preID = str(re.findall('ID\s*(\S*)\s*standard;',r)[0]) 
     109                    preID = str(IDpat.findall(r)[0]) 
    100110                    print preID 
    101111                         
    102112                elif r[0:2] == 'AC': 
    103                     preACC = str(re.findall('AC\s*(\S*);',r)[0]) 
     113                    preACC = str(ACpat.findall(r)[0]) 
    104114                    web_addr = 'http://www.mirbase.org/cgi-bin/mirna_entry.pl?acc=%s' % preACC 
    105115                         
    106                 elif r[0:2] == 'RX' and not(re.findall('RX\s*PUBMED;\s(\d*).',r)==[]): 
    107                     pubIDs.append(str(re.findall('RX\s*PUBMED;\s(\d*).',r)[0])) 
    108                              
    109                 elif r[0:2]=='FT' and not(re.findall('FT\s*miRNA\s*(\d{1,}\.\.\d{1,})',r)==[]): 
    110                     loc_mat = str(re.findall('FT\s*miRNA\s*(\d{1,}\.\.\d{1,})',r)[0]) 
     116                elif r[0:2] == 'RX' and not(RXpat.findall(r)==[]): 
     117                    pubIDs.append(str(RXpat.findall(r)[0])) 
     118 
     119                elif r[0:2]=='FT' and not(FT1pat.findall(r)==[]): 
     120                    loc_mat = str(FT1pat.findall(r)[0]) 
    111121                         
    112122                    if not(loc_mat==[]): 
    113123                         my_locs.append(loc_mat) 
    114124                 
    115                 elif r[0:2]=='FT' and not(re.findall('FT\s*/accession="(MIMAT[0-9]*)"', r)==[]): 
    116                      mat_acc = str(re.findall('FT\s*/accession="(MIMAT[0-9]*)"', r)[0]) 
     125                elif r[0:2]=='FT' and not(FT2pat.findall(r)==[]): 
     126                     mat_acc = str(FT2pat.findall(r)[0]) 
    117127                         
    118128                     if matACCs == '': 
     
    124134                         my_accs.append(mat_acc)     
    125135                                 
    126                 elif r[0:2]=='FT' and not(re.findall('FT\s*/product="(\S*)"', r)==[]): 
    127                      mat_id = str(re.findall('FT\s*/product="(\S*)"', r)[0]) 
     136                elif r[0:2]=='FT' and not(FT3pat.findall(r)==[]): 
     137                     mat_id = str(FT3pat.findall(r)[0]) 
    128138                         
    129139                     if matIDs == '': 
     
    137147                elif r[0:2]=='SQ': 
    138148             
    139                      preSQ_INFO = str(re.findall('SQ\s*(.*other;)', r)[0]) 
     149                     preSQ_INFO = str(SQpat.findall(r)[0]) 
    140150                     seq = 'on' 
    141151             
    142152                elif r[0:2]=='  ' and seq == 'on': 
    143                      preSQ.append(str(re.findall('\s*([a-z\s]*)\s*\d*',r)[0]).replace(' ','')) 
     153                     preSQ.append(str(seqpat.findall(r)[0]).replace(' ','')) 
    144154                      
    145155            ### cluster search 
     
    274284                 
    275285                if type_file == 'mat': 
    276                     serverFiles.upload("miRNA", label, filename, title="miRNA: %s mature form" % org, tags=["tag1", "tag2"]) 
     286                    serverFiles.upload("miRNA", label, filename, title="miRNA: %s mature form" % org, tags=["miRNA"]) 
    277287                    serverFiles.unprotect("miRNA", label) 
    278288                    print '%s mat uploaded' % org 
     
    282292                     
    283293                elif type_file == 'pre': 
    284                     serverFiles.upload("miRNA", label, filename, title="miRNA: %s pre-form" % org, tags=["tag1", "tag2"]) 
     294                    serverFiles.upload("miRNA", label, filename, title="miRNA: %s pre-form" % org, tags=["miRNA"]) 
    285295                    serverFiles.unprotect("miRNA", label) 
    286296                    print '%s pre uploaded' % org 
     
    292302                    print 'Check the label.' 
    293303     
    294     serverFiles.upload("miRNA", "miRNA.txt", miRNA_path) 
     304    serverFiles.upload("miRNA", "miRNA.txt", miRNA_path, title="miRNA: miRNA library", tags=["miRNA"] ) 
    295305    serverFiles.unprotect("miRNA", "miRNA.txt") 
    296306    print '\nmiRNA.txt uploaded' 
    297307     
    298     serverFiles.upload("miRNA", "premiRNA.txt", premiRNA_path) 
     308    serverFiles.upload("miRNA", "premiRNA.txt", premiRNA_path, title="miRNA: pre-form library", tags=["miRNA"]) 
    299309    serverFiles.unprotect("miRNA", "premiRNA.txt") 
    300310    print 'premiRNA.txt uploaded\n' 
Note: See TracChangeset for help on using the changeset viewer.