source: orange-bioinformatics/_bioinformatics/obimiRNA.py @ 1636:10d234fdadb9

Revision 1636:10d234fdadb9, 14.6 KB checked in by mitar, 2 years ago (diff)

Restructuring because we will not be using namespaces.

Line 
1from __future__ import absolute_import, division
2
3from collections import defaultdict
4import math, os, random, re, urllib
5
6from Orange.orng import orngServerFiles as osf
7import statc
8
9from . import obiGene as ge, obiGO as go, obiKEGG as kg, obiProb as op, obiTaxonomy
10
11mirnafile = osf.localpath_download('miRNA','miRNA.txt')
12premirnafile = osf.localpath_download('miRNA','premiRNA.txt')
13
14################################################################################################################
15################################################################################################################
16
17def __build_lib(filename, labels=True, MATtoPRE=True, ACCtoID=True, clust=False):
18    """
19    build_lib() function takes as input a filename
20    and gives as output some variables there will be used in
21    the module.
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
46
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
163
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'}
167
168toTaxo = dict(zip(fromTaxo.values(), fromTaxo.keys()))
169
170
171class miRNAException(Exception):
172    pass
173
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:
190            raise miRNAException("ids() Error: Check the input value.")
191
192           
193       
194class mat_miRNA:
195    pass
196
197class pre_miRNA:
198    pass
199 
200     
201def get_info(objectID,type='mat'):
202        """
203        get_info() function takes a miRNA identifier as input
204        and returns a miRNA object.
205        """
206        if type == 'mat':
207            objectID = re.sub('mir','miR',objectID)
208            if objectID in IDs:
209                attr = [line.rstrip() for line in open(mirnafile).readlines()][0].split('\t')
210           
211                to_return = mat_miRNA()
212                setattr(to_return, attr[0], objectID)
213           
214                for n,a in enumerate(attr[1:]):
215                    setattr(to_return, a, miRNA_lib[objectID][n])
216           
217                return to_return
218            else:
219                raise miRNAException("get_info() Error: Check the input value.")
220           
221        elif type == 'pre':
222            objectID = re.sub('miR','mir',objectID)
223            if objectID in preIDs:
224                attr = [line.rstrip() for line in open(premirnafile).readlines()][0].split('\t')
225           
226                to_return = pre_miRNA()
227                setattr(to_return, attr[0], objectID)
228           
229                for n,a in enumerate(attr[1:]):
230                    setattr(to_return, a, premiRNA_lib[objectID][n])
231                           
232                return to_return
233            else:
234                raise miRNAException("get_info() Error: Check the input value.")
235        else:
236           raise miRNAException("get_info() Error: Check type value.")
237
238                   
239def cluster(clusterID, type='name'):
240    """
241    cluster() function take a cluster identifier or a premiRNA
242    and return the list of premiRNAs clustered together."
243    """
244    if type=='name':
245        if clusterID in clusters:
246            return clusters[clusterID]
247        else:
248            raise miRNAException("cluster() Error: ClusterID not found in premiRNA names.")
249   
250    elif type=='num':
251        if clusterID in num_toClusters:
252            return num_toClusters[clusterID]
253        else:
254            raise miRNAException("cluster() Error: ClusterID not found in clusters' list.")
255    else:
256        raise miRNAException("cluster() Error: Check the input value.")
257
258
259def fromACC_toID(accession):
260    """
261    fromACC_toID() takes a miRNA accession number
262    and returns a miRNA id.
263    """
264    if accession in ACCtoID:
265        return ACCtoID[accession]
266    if accession in preACCtoID:
267        return preACCtoID[accession]
268    else:
269        print "Accession not found."
270        return False
271   
272
273def __reverseDict(old_dict):
274    """
275    switch dictionary: keys <--> values;
276    """
277    new_dict = {}
278    for k,valuesList in old_dict.items():
279        for val in valuesList:
280            if val != [] and not(val in new_dict):
281                new_dict[val]=[]
282            if not(k in new_dict[val]):
283                new_dict[val].append(k)
284    return new_dict
285
286
287def get_geneMirnaLib(org=None):
288    """
289    build dictionary gene:[miRNAs]
290    """
291    mirnaGenes={}
292    for m in ids(org):
293        mirnaGenes[m] = get_info(m).targets.split(',')
294       
295    return __reverseDict(mirnaGenes)
296
297
298def get_GO(mirna_list, annotations, enrichment=False, pval=0.1, goSwitch=True):
299    """
300    get_GO() takes as input a list of miRNAs of the organism for which the annotations are defined.
301    If goSwitch is False, get_GO() returns a dictionary that has miRNAs as keys and GO IDs as values;
302    in the other case it returns a dictionary with GO IDs as keys and miRNAs as values.
303    """
304   
305    from . import obiGene
306    genematcher = obiGene.matcher([obiGene.GMGO(annotations.taxid)] + \
307        ([obiGene.GMDicty()] if annotations.taxid == "352472"  else []))
308    genematcher.set_targets(annotations.geneNames)
309   
310    mirna_list = list(set(mirna_list))
311    mirAnnotations = {}
312   
313    for m in mirna_list:
314        genes = get_info(m).targets.split(',')
315        genes = filter(None, map(genematcher.umatch, genes))
316       
317        if enrichment==False:
318            mirna_ann=[]
319            for gene in genes:
320                gene_annotations = annotations.geneAnnotations[gene]
321                for ga in gene_annotations:
322                    mirna_ann.append(ga)
323            mirAnnotations[m] = list(set([an.GO_ID for an in mirna_ann]))
324        elif enrichment==True:
325            if len(genes):
326                resP = annotations.GetEnrichedTerms(genes,aspect='P')
327                resC = annotations.GetEnrichedTerms(genes,aspect='C')
328                resF = annotations.GetEnrichedTerms(genes,aspect='F')
329                res = {}
330                res.update(resP)
331                res.update(resC)
332                res.update(resF)
333                tups = [(pVal,go_id) for go_id, (ge,pVal,ref) in res.items()]
334                tups.sort()           
335                p_correct = op.FDR([p for p,go_id in tups])           
336                mirAnnotations[m] = [tups[i][1] for i, p in enumerate(p_correct) if p < pval]
337            else:
338                mirAnnotations[m]=[]
339
340    if goSwitch:
341        return __reverseDict(mirAnnotations)
342    else:
343        return mirAnnotations
344 
345
346
347def filter_GO(mirna_goid, annotations, treshold=0.04, reverse=True):   
348    """
349    filter_GO() takes as input a dictionary like {mirna:[list of GO_IDs]} and
350    remove the most common GO IDs in each list using the TF-IDF criterion.
351    """       
352    uniqGO = list(set(reduce(lambda x,y: x+y, mirna_goid.values())))       
353   
354    goIDF = {}   
355    for n,go in enumerate(uniqGO):       
356        goIDF[go] = np.log(len(mirna_goid)/len(filter(lambda x: go in x, mirna_goid.values())))       
357
358    data = []
359    new_dict={}
360    for m,goList in mirna_goid.items():
361        TF_IDF ={}
362        genes = get_info(m).targets.split(',')
363        mirnaAnnotations = [annotations.GetAnnotatedTerms(g).keys() for g in genes]
364        for go in goList:
365            TF = len(filter(lambda x: go in x, mirnaAnnotations))/len(genes)
366            TF_IDF[go] = TF*goIDF[go]
367        if treshold:
368            t = treshold
369        else:
370            t=0
371        new_dict[m] = filter(lambda x: TF_IDF[x] > t, goList)   
372        #data.append(TF_IDF.values())
373       
374#    if treshold:
375#        data = sorted(reduce(lambda x,y: x+y, data))   
376#        tresholds = filter(lambda x: x[1]>treshold[0] and x[1]<treshold[1], [(d,statc.percentileofscore(data, d)) for d in list(set(data))])
377#        t = np.mean([t[0] for t in tresholds])
378#        print t
379#    else:
380#        t=0
381#   
382#    new_dict={}     
383#    for m,goList in mirna_goid.items():       
384#        new_dict[m] = filter(lambda x: TF_IDF[x] > t, goList)
385   
386    if reverse:
387        return __reverseDict(new_dict)
388    else:
389        return new_dict
390
391
392
393def get_pathways(mirna_list, organism='hsa', enrichment=False, pVal=0.1, pathSwitch=True):
394    """
395    get_pathways() takes as input a list of miRNAs and returns a dictionary that has miRNAs as keys
396    and pathways IDs as values; if the switch is set on True,
397    it returns a dictionary with pathways IDs as keys and miRNAs as values.
398    """
399    gmkegg = ge.GMKEGG(organism)
400    org = kg.KEGGOrganism(organism)
401    gmkegg.set_targets(org.get_genes())     
402   
403    genes = list(set(reduce(lambda x,y: x+y,[get_info(m).targets.split(',') for m in mirna_list])))
404    keggNames = dict([(g,gmkegg.umatch(g)) for g in genes if gmkegg.umatch(g)])
405       
406    mirnaPathways = {}
407    for m in mirna_list:
408        kegg_genes = [keggNames[g] for g in get_info(m).targets.split(',') if g in keggNames]
409        if enrichment:
410            mirnaPathways[m] = [path_id for path_id,(geneList,p,geneNum) in org.get_enriched_pathways_by_genes(kegg_genes).items() if p < pVal]
411        else:
412            paths = filter(None,[list(org.get_pathways_by_genes([k])) for k in kegg_genes])                   
413            if paths:
414                mirnaPathways[m] = list(set(reduce(lambda x,y: x+y,paths)))
415            else:
416                mirnaPathways[m] = []
417   
418    if pathSwitch:
419        return __reverseDict(mirnaPathways)
420    else:
421        return mirnaPathways
422       
423
424def removeOldMirnas(mirna_list, getOnlyMature=False):
425    """
426    removeOldMirnas() takes a list of miRNAs as input and
427    divides them in two lists, accordin if they're still present
428    on miRBase or not.
429    """
430    old = []
431    for m in mirna_list:
432        try:
433            mat = get_info(m)
434        except Exception:
435            if getOnlyMature:
436                old.append(m)
437            else:               
438                try:
439                    pre = get_info(m,type='pre')
440                except Exception:
441                    old.append(m)
442    return [old, filter(lambda x: not(x in old), mirna_list)]
443   
444#######################
445
446
447
448
449
450
451
452
453
454
455
456     
Note: See TracBrowser for help on using the repository browser.