source: orange-bioinformatics/obimiRNA.py @ 1424:56f4abf21046

Revision 1424:56f4abf21046, 14.7 KB checked in by ales_erjavec <ales.erjavec@…>, 3 years ago (diff)

Added support for loading miRNA targets from microCosm files.

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