source: orange-bioinformatics/_bioinformatics/obimiRNA.py @ 1834:91f840200ae3

Revision 1834:91f840200ae3, 14.7 KB checked in by Flashpoint <vid.flashpoint@…>, 10 months ago (diff)

Changed the exception value from KEGG.OrganismNotFound to obiTaxonomy.UnknownSpeciesIdentifier

RevLine 
[1632]1from __future__ import absolute_import, division
[1018]2
[1424]3from collections import defaultdict
[1632]4import math, os, random, re, urllib
[1424]5
[1632]6from Orange.orng import orngServerFiles as osf
7import statc
8
[1716]9from . import obiGene as ge, obiGO as go, obiKEGG as kg, obiProb as op, obiTaxonomy
[1018]10
11mirnafile = osf.localpath_download('miRNA','miRNA.txt')
12premirnafile = osf.localpath_download('miRNA','premiRNA.txt')
13
[1031]14################################################################################################################
15################################################################################################################
[1018]16
[1024]17def __build_lib(filename, labels=True, MATtoPRE=True, ACCtoID=True, clust=False):
[1023]18    """
[1047]19    build_lib() function takes as input a filename
20    and gives as output some variables there will be used in
21    the module.
[1023]22    """
23    content = [l.rstrip() for l in open(filename).readlines()][1:]
24    to_return = []
25   
26    ids = [l.split('\t')[0] for l in content]
27    to_return.append(ids)
28   
29    if labels: 
30        to_return.append(list(set(i.split('-')[0] for i in ids)))
31   
32    elements = [line.rstrip().split('\t') for line in content] 
33    to_return.append(dict((elem[0],elem[1:]) for elem in elements))
34   
35    if MATtoPRE:
36        to_return.append(dict([(e[0],e[3]) for e in elements]))
37   
38    if ACCtoID:
39        to_return.append(dict([(e[1],e[0]) for e in elements]))
40       
41    if clust:
42        to_return.append(dict([(e[0],e[5]) for e in elements]))
43       
44   
45    return to_return
[1018]46
[1424]47def open_microcosm(org="mus_musculus", version="v5"):
48    """ Open the miRna targets from the EBI microcosm site.
49    """
50    import urllib2, zipfile
51    import StringIO
52    stream = urllib2.urlopen("ftp://ftp.ebi.ac.uk/pub/databases/microcosm/{version}/arch.{version}.txt.{org}.zip".format(org=org, version=version))
53    contents = stream.read()
54    stream = StringIO.StringIO(contents)
55    return stream
56   
57def parse_targets_microcosm_v5(file, max_pvalue=None, min_score=None):
58    """ Parse miRna targets (version 'v5.txt' format file from microcosm)
59    """
60   
61    ids = set() #mirna ids
62    labels = set() # org codes (i.e hsa, mmu ..)
63    #["mirna accession", "seq", "mirna id", "targets"]
64    mirna_lib = defaultdict(lambda: ["", "", "", set()])
65    mat_to_pre = {}
66    if isinstance(file, basestring):
67        file = open(file, "rb")
68    import csv
69    reader = csv.reader(file, delimiter="\t")
70    for line in reader:
71        if not line or line[0].startswith("##"):
72            continue
73       
74        seq = line[1]
75        org = seq.split("-", 1)[0] #SEQ
76        ids.add(seq)
77        labels.add(org)
78       
79        mirna_lib[seq] #In case some mirna dont have any targets
80       
81        if max_pvalue is not None:
82            if float(line[10]) > max_pvalue:
83                continue
84        if min_score is not None:
85            if float(line[9]) < min_score:
86                continue
87           
88        mirna_lib[seq][-1].add(line[11]) #TRANSCRIPT_ID
89       
90    for key, value in mirna_lib.iteritems():
91        value[-1] = ",".join(sorted(value[-1]) or ["None"])
92       
93    ids = sorted(ids)
94    labels = sorted(labels)
95    mirna_lib = dict(mirna_lib)
96    return ids, labels, mirna_lib, {}, {}
97
98
99def load_miRNA_microCosm(org="mus_musculus", max_pvalue=None, min_score=None):
100    """ Load miRNA's from microcosm into the global scope (currently
101    only Mus musculus is supported)
102   
103    """
104    global IDs, LABELS, miRNA_lib, mat_toPre, ACCtoID
105    global preIDs, premiRNA_lib, preACCtoID, clusters
106    global num_toClusters,  clusters_toNum
107   
108    file = osf.localpath_download("miRNA", "v5.txt.{org}".format(org=org))
109    [IDs, LABELS, miRNA_lib, mat_toPre, ACCtoID] = parse_targets_microcosm_v5(file,
110                                max_pvalue=max_pvalue, min_score=min_score)
111    [preIDs, premiRNA_lib, preACCtoID, clusters] = [], {}, {}, {}
112    num_toClusters, clusters_toNum = {}, {}
113
114   
115load_miRNA = load_miRNA_microCosm
116
117
118def load_miRNA_TargetScan():
119    """ This loads miRNAs from miRBase and targets from TargetScan.
120    Will also load pre-miRNAs
121    """
122    global IDs, LABELS, miRNA_lib, mat_toPre, ACCtoID
123    global preIDs, premiRNA_lib, preACCtoID, clusters
124    global num_toClusters,  clusters_toNum
125   
126    [IDs, LABELS, miRNA_lib, mat_toPre, ACCtoID] = __build_lib(mirnafile, 1,1,1,0)
127    [preIDs, premiRNA_lib, preACCtoID, clusters] = __build_lib(premirnafile,0,0,1,1)
128   
129    num_toClusters = {}
130    clusters_toNum = {}
131    n=0
132    for k,v in clusters.items():
133        if v !='None':
134            g = v.split(',')
135            g.append(k)       
136            group = sorted(g)
137            if not(group in num_toClusters.values()):
138                 num_toClusters[n] = group
139                 for e in group:
140                     clusters_toNum[e]=n
141                 n += 1
142     
143### library:
144load_miRNA()
145
146#[IDs, LABELS, miRNA_lib, mat_toPre, ACCtoID] = __build_lib(mirnafile, 1,1,1,0)
147#[preIDs, premiRNA_lib, preACCtoID, clusters] = __build_lib(premirnafile,0,0,1,1)
148#
149#num_toClusters = {}
150#clusters_toNum = {}
151#n=0
152#for k,v in clusters.items():
153#    if v !='None':
154#        g = v.split(',')
155#        g.append(k)       
156#        group = sorted(g)
157#        if not(group in num_toClusters.values()):
158#             num_toClusters[n] = group
159#             for e in group:
160#                 clusters_toNum[e]=n
161#             n += 1
162
[1018]163
[1023]164fromTaxo = {3702:'ath', 9913:'bta', 6239:'cel', 3055:'cre', 7955:'dre',\
165             352472:'ddi', 7227:'dme', 9606:'hsa', 10090:'mmu', 4530:'osa',\
166              10116:'rno', 8355:'xla', 4577:'zma'}
[1018]167
[1424]168toTaxo = dict(zip(fromTaxo.values(), fromTaxo.keys()))
[1041]169
170
[1023]171class miRNAException(Exception):
172    pass
[1018]173
[1023]174def ids(taxid=None):
175    """
176    ids() functions takes an organism identifier (for human can be 9606 or hsa)
177    and returns all the miRNAs related to that organism. If there
178    is no argument, it returns all the miRNAs in the library.
179    """
180   
181    if not(taxid):
182        return IDs
183    else:
184        taxid = fromTaxo.get(taxid, taxid)
185       
186        if (taxid in toTaxo):
187            return [e for e in IDs if e.split('-')[0]==taxid]
188               
189        else:
[1716]190            from . import obiKEGG
[1834]191            raise obiTaxonomy.UnknownSpeciesIdentifier(taxid)
[1018]192       
193class mat_miRNA:
194    pass
195
196class pre_miRNA:
197    pass
198 
199     
[1023]200def get_info(objectID,type='mat'):
201        """
202        get_info() function takes a miRNA identifier as input
203        and returns a miRNA object.
204        """
205        if type == 'mat':
[1031]206            objectID = re.sub('mir','miR',objectID)
[1023]207            if objectID in IDs:
208                attr = [line.rstrip() for line in open(mirnafile).readlines()][0].split('\t')
[1018]209           
[1023]210                to_return = mat_miRNA()
211                setattr(to_return, attr[0], objectID)
[1018]212           
[1023]213                for n,a in enumerate(attr[1:]):
214                    setattr(to_return, a, miRNA_lib[objectID][n])
[1018]215           
[1023]216                return to_return
217            else:
218                raise miRNAException("get_info() Error: Check the input value.")
219           
220        elif type == 'pre':
[1080]221            objectID = re.sub('miR','mir',objectID)
[1023]222            if objectID in preIDs:
223                attr = [line.rstrip() for line in open(premirnafile).readlines()][0].split('\t')
224           
225                to_return = pre_miRNA()
226                setattr(to_return, attr[0], objectID)
227           
228                for n,a in enumerate(attr[1:]):
229                    setattr(to_return, a, premiRNA_lib[objectID][n])
230                           
231                return to_return
232            else:
233                raise miRNAException("get_info() Error: Check the input value.")
[1018]234        else:
[1023]235           raise miRNAException("get_info() Error: Check type value.")
236
237                   
238def cluster(clusterID, type='name'):
239    """
240    cluster() function take a cluster identifier or a premiRNA
241    and return the list of premiRNAs clustered together."
242    """
243    if type=='name':
244        if clusterID in clusters:
245            return clusters[clusterID]
246        else:
247            raise miRNAException("cluster() Error: ClusterID not found in premiRNA names.")
[1018]248   
[1023]249    elif type=='num':
250        if clusterID in num_toClusters:
251            return num_toClusters[clusterID]
252        else:
253            raise miRNAException("cluster() Error: ClusterID not found in clusters' list.")
254    else:
255        raise miRNAException("cluster() Error: Check the input value.")
[1018]256
257
[1031]258def fromACC_toID(accession):
259    """
[1047]260    fromACC_toID() takes a miRNA accession number
261    and returns a miRNA id.
[1031]262    """
263    if accession in ACCtoID:
264        return ACCtoID[accession]
[1077]265    if accession in preACCtoID:
266        return preACCtoID[accession]
[1031]267    else:
268        print "Accession not found."
269        return False
[1047]270   
[1077]271
272def __reverseDict(old_dict):
[1047]273    """
[1077]274    switch dictionary: keys <--> values;
[1047]275    """
[1077]276    new_dict = {}
277    for k,valuesList in old_dict.items():
278        for val in valuesList:
279            if val != [] and not(val in new_dict):
280                new_dict[val]=[]
281            if not(k in new_dict[val]):
282                new_dict[val].append(k)
283    return new_dict
284
285
286def get_geneMirnaLib(org=None):
287    """
288    build dictionary gene:[miRNAs]
289    """
290    mirnaGenes={}
291    for m in ids(org):
292        mirnaGenes[m] = get_info(m).targets.split(',')
293       
294    return __reverseDict(mirnaGenes)
[1041]295
296
[1047]297def get_GO(mirna_list, annotations, enrichment=False, pval=0.1, goSwitch=True):
[1041]298    """
[1077]299    get_GO() takes as input a list of miRNAs of the organism for which the annotations are defined.
[1047]300    If goSwitch is False, get_GO() returns a dictionary that has miRNAs as keys and GO IDs as values;
301    in the other case it returns a dictionary with GO IDs as keys and miRNAs as values.
[1041]302    """
[1047]303   
[1632]304    from . import obiGene
[1047]305    genematcher = obiGene.matcher([obiGene.GMGO(annotations.taxid)] + \
306        ([obiGene.GMDicty()] if annotations.taxid == "352472"  else []))
307    genematcher.set_targets(annotations.geneNames)
308   
309    mirna_list = list(set(mirna_list))
[1041]310    mirAnnotations = {}
[1031]311   
[1047]312    for m in mirna_list:
[1041]313        genes = get_info(m).targets.split(',')
[1047]314        genes = filter(None, map(genematcher.umatch, genes))
[1041]315       
[1047]316        if enrichment==False:
[1041]317            mirna_ann=[]
318            for gene in genes:
319                gene_annotations = annotations.geneAnnotations[gene]
[1047]320                for ga in gene_annotations:
321                    mirna_ann.append(ga)
[1077]322            mirAnnotations[m] = list(set([an.GO_ID for an in mirna_ann]))
[1047]323        elif enrichment==True:
[1058]324            if len(genes):
325                resP = annotations.GetEnrichedTerms(genes,aspect='P')
326                resC = annotations.GetEnrichedTerms(genes,aspect='C')
327                resF = annotations.GetEnrichedTerms(genes,aspect='F')
328                res = {}
329                res.update(resP)
330                res.update(resC)
331                res.update(resF)
332                tups = [(pVal,go_id) for go_id, (ge,pVal,ref) in res.items()]
333                tups.sort()           
[1077]334                p_correct = op.FDR([p for p,go_id in tups])           
[1058]335                mirAnnotations[m] = [tups[i][1] for i, p in enumerate(p_correct) if p < pval]
[1077]336            else:
337                mirAnnotations[m]=[]
[1047]338
339    if goSwitch:
[1077]340        return __reverseDict(mirAnnotations)
[1047]341    else:
342        return mirAnnotations
343 
[1077]344
345
[1100]346def filter_GO(mirna_goid, annotations, treshold=0.04, reverse=True):   
[1077]347    """
[1093]348    filter_GO() takes as input a dictionary like {mirna:[list of GO_IDs]} and
[1077]349    remove the most common GO IDs in each list using the TF-IDF criterion.
[1100]350    """       
[1077]351    uniqGO = list(set(reduce(lambda x,y: x+y, mirna_goid.values())))       
352   
353    goIDF = {}   
354    for n,go in enumerate(uniqGO):       
355        goIDF[go] = np.log(len(mirna_goid)/len(filter(lambda x: go in x, mirna_goid.values())))       
356
[1100]357    data = []
[1077]358    new_dict={}
359    for m,goList in mirna_goid.items():
360        TF_IDF ={}
361        genes = get_info(m).targets.split(',')
362        mirnaAnnotations = [annotations.GetAnnotatedTerms(g).keys() for g in genes]
363        for go in goList:
364            TF = len(filter(lambda x: go in x, mirnaAnnotations))/len(genes)
365            TF_IDF[go] = TF*goIDF[go]
366        if treshold:
[1100]367            t = treshold
[1077]368        else:
[1100]369            t=0
370        new_dict[m] = filter(lambda x: TF_IDF[x] > t, goList)   
371        #data.append(TF_IDF.values())
372       
373#    if treshold:
374#        data = sorted(reduce(lambda x,y: x+y, data))   
375#        tresholds = filter(lambda x: x[1]>treshold[0] and x[1]<treshold[1], [(d,statc.percentileofscore(data, d)) for d in list(set(data))])
376#        t = np.mean([t[0] for t in tresholds])
377#        print t
378#    else:
379#        t=0
380#   
381#    new_dict={}     
382#    for m,goList in mirna_goid.items():       
383#        new_dict[m] = filter(lambda x: TF_IDF[x] > t, goList)
[1077]384   
385    if reverse:
386        return __reverseDict(new_dict)
387    else:
388        return new_dict
389
390
391
392def get_pathways(mirna_list, organism='hsa', enrichment=False, pVal=0.1, pathSwitch=True):
393    """
394    get_pathways() takes as input a list of miRNAs and returns a dictionary that has miRNAs as keys
395    and pathways IDs as values; if the switch is set on True,
396    it returns a dictionary with pathways IDs as keys and miRNAs as values.
397    """
398    gmkegg = ge.GMKEGG(organism)
399    org = kg.KEGGOrganism(organism)
400    gmkegg.set_targets(org.get_genes())     
401   
402    genes = list(set(reduce(lambda x,y: x+y,[get_info(m).targets.split(',') for m in mirna_list])))
403    keggNames = dict([(g,gmkegg.umatch(g)) for g in genes if gmkegg.umatch(g)])
404       
405    mirnaPathways = {}
406    for m in mirna_list:
407        kegg_genes = [keggNames[g] for g in get_info(m).targets.split(',') if g in keggNames]
408        if enrichment:
[1715]409            mirnaPathways[m] = [path_id for path_id,(geneList,p,geneNum) in org.get_enriched_pathways(kegg_genes).items() if p < pVal]
[1077]410        else:
411            paths = filter(None,[list(org.get_pathways_by_genes([k])) for k in kegg_genes])                   
412            if paths:
413                mirnaPathways[m] = list(set(reduce(lambda x,y: x+y,paths)))
414            else:
415                mirnaPathways[m] = []
416   
417    if pathSwitch:
418        return __reverseDict(mirnaPathways)
419    else:
420        return mirnaPathways
421       
422
[1100]423def removeOldMirnas(mirna_list, getOnlyMature=False):
[1102]424    """
425    removeOldMirnas() takes a list of miRNAs as input and
426    divides them in two lists, accordin if they're still present
427    on miRBase or not.
428    """
[1100]429    old = []
430    for m in mirna_list:
431        try:
432            mat = get_info(m)
433        except Exception:
434            if getOnlyMature:
435                old.append(m)
436            else:               
437                try:
438                    pre = get_info(m,type='pre')
439                except Exception:
440                    old.append(m)
441    return [old, filter(lambda x: not(x in old), mirna_list)]
[1041]442   
[1023]443#######################
[1018]444
[1041]445
446
447
448
449
450
451
452
453
454
[1047]455     
Note: See TracBrowser for help on using the repository browser.