This repository has been archived by the owner on Nov 2, 2021. It is now read-only.
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgene_associations.py
1626 lines (1491 loc) · 83.1 KB
/
gene_associations.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
###gene_associations
#Copyright 2005-2008 J. David Gladstone Institutes, San Francisco California
#Author Nathan Salomonis - [email protected]
#Permission is hereby granted, free of charge, to any person obtaining a copy
#of this software and associated documentation files (the "Software"), to deal
#in the Software without restriction, including without limitation the rights
#to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
#copies of the Software, and to permit persons to whom the Software is furnished
#to do so, subject to the following conditions:
#THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
#INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A
#PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT
#HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
#OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
#SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
"""This module contains fairly generic methods for reading different GO-Elite gene relationship
files, user input and denominator files, and summarizing quantitative and gene annotation data."""
import sys, string
import os.path
import unique
import math
import statistics
import gene_associations
import OBO_import
import export
import re
###### File Import Functions ######
def filepath(filename):
fn = unique.filepath(filename)
return fn
def read_directory(sub_dir):
dir_list = unique.read_directory(sub_dir); dir_list2 = []
###Code to prevent folder names from being included
for entry in dir_list:
if entry[-4:] == ".txt" or entry[-4:] == ".csv" or ".owl" in entry or ".gpml" in entry or ".xml" in entry or ".gmt" in entry: dir_list2.append(entry)
return dir_list2
###### Classes ######
class GrabFiles:
def setdirectory(self,value):
self.data = value
def display(self):
print self.data
def searchdirectory(self,search_term):
#self is an instance while self.data is the value of the instance
all_matching,file_dir,file = getDirectoryFiles(self.data,str(search_term))
#if len(file)<1: print self.data,search_term,'not found', filepath(self.data)
return file_dir,file
def getAllFiles(self,search_term):
#self is an instance while self.data is the value of the instance
all_matching,file_dir,file = getDirectoryFiles(self.data,str(search_term))
#if len(file)<1: print search_term,'not found'
return all_matching
def getMatchingFolders(self,search_term):
dir_list = unique.read_directory(self.data); matching = ''
root_dir = filepath('')
for folder in dir_list:
if search_term in folder and '.' not in folder:
matching = filepath(self.data[1:]+'/'+folder)
if root_dir not in matching:
matching = root_dir + matching ### python bug????
return matching
def getDirectoryFiles(import_dir, search_term):
exact_file = ''; exact_file_dir=''; all_matching=[]
dir_list = read_directory(import_dir) #send a sub_directory to a function to identify all files in a directory
for data in dir_list: #loop through each file in the directory to output results
present = verifyFile(import_dir+'/'+data) ### Check to see if the file is present in formatting the full file_dir
if present == True:
affy_data_dir = import_dir+'/'+data
else: affy_data_dir = import_dir[1:]+'/'+data
if search_term in affy_data_dir:
if 'version.txt' not in affy_data_dir: exact_file_dir = affy_data_dir; exact_file = data; all_matching.append(exact_file_dir)
return all_matching, exact_file_dir,exact_file
def verifyFile(filename):
present = False
try:
fn=filepath(filename)
for line in open(fn,'rU').xreadlines(): present = True; break
except Exception: present = False
return present
class GeneRelationships:
def __init__(self,uid,gene,terms,uid_system,gene_system):
self._uid = uid; self._gene = gene; self._terms = terms
self._uid_system = uid_system; self._gene_system = gene_system
def UniqueID(self): return self._uid
def GeneID(self): return self._gene
def GOIDs(self): return self._terms
def GOIDInts(self):
goid_list = self._terms; goid_list_int = []
for goid in goid_list: goid_list_int.append(int(goid))
return goid_list_int
def GOIDStrings(self):
goid_list = self._terms; goid_list_str = []
for goid in goid_list:
if ':' not in goid: goid_list_str.append('GO:'+goid)
return goid_list_str
def MAPPs(self):
mapp_list = self._terms
return mapp_list
def UIDSystem(self): return self.uid_system
def GeneIDSystem(self): return self._gene_system
def setUIDValues(self,uid_values):
self._uid_values = uid_values
def UIDValues(self): return self._uid_values
def Report(self):
output = self.UniqueID()
return output
def __repr__(self): return self.Report()
class GeneAnnotations:
def __init__(self,gene,symbol,name,gene_system):
self._gene = gene; self._symbol = symbol; self._name = name;
self._gene_system = gene_system
def GeneID(self): return self._gene
def Symbol(self): return self._symbol
def SymbolLower(self): return string.lower(self._symbol)
def Description(self): return self._name
def GeneIDSystem(self): return self._gene_system
def Report(self):
output = self.GeneID()+'|'+self.Symbol()
return output
def __repr__(self): return self.Report()
class RedundantRelationships:
def __init__(self,redundant_names,inverse_names):
self._redundant_names = redundant_names; self._inverse_names = inverse_names
def RedundantNames(self):
redundant_with_terms = string.join(self._redundant_names,'|')
return redundant_with_terms
def InverseNames(self):
redundant_with_terms = string.join(self._inverse_names,'|')
return redundant_with_terms
def Report(self):
output = self.RedundantNames()+'|'+self.InverseNames()
return output
def __repr__(self): return self.Report()
def cleanUpLine(line):
data = string.replace(line,'\n','')
data = string.replace(data,'\c','')
data = string.replace(data,'\r','')
data = string.replace(data,'"','')
return data
def importGeneric(filename):
fn=filepath(filename); key_db = {}
for line in open(fn,'rU').xreadlines():
data = cleanUpLine(line)
t = string.split(data,'\t')
key_db[t[0]] = t[1:]
return key_db
###### Begin Gene-Relationship Parsing ######
def importGeneData(species_code,mod):
if 'export' in mod: export_status,mod = mod
else: export_status = 'no'
###Pass in species_code, mod
program_type,database_dir = unique.whatProgramIsThis()
gene_import_dir = '/'+database_dir+'/'+species_code+'/gene'
g = GrabFiles(); g.setdirectory(gene_import_dir)
filedir,filename = g.searchdirectory(mod) ###Identify gene files corresponding to a particular MOD
system_type = string.replace(filename,'.txt','')
fn=filepath(filedir); gene_annotations={}; x=0
for line in open(fn,'rU').xreadlines():
data = cleanUpLine(line)
t = string.split(data,'\t')
if x == 0: x = 1
else:
gene = t[0]; symbol = t[1]; name = t[2]
s = GeneAnnotations(gene,symbol,name,system_type)
gene_annotations[gene] = s
if export_status == 'export':
export_dir = database_dir+'/'+species_code+'/uid-gene/'+mod+'-Symbol'+'.txt'
data = export.ExportFile(export_dir)
data.write('GeneID\tSymbol\n')
for gene in gene_annotations:
s = gene_annotations[gene]
if len(s.Symbol())>0:
data.write(gene+'\t'+s.Symbol()+'\n')
data.close()
else:
return gene_annotations
def buildNestedAssociations(species_code):
print "Nested gene association, not yet build for",species_code
print "Proceeding to build nested associations for this species\n"
export_databases = 'no'; genmapp_mod = 'Ensembl' ### default
import GO_Elite; system_codes,source_types,mod_types = GO_Elite.getSourceData()
OBO_import.buildNestedOntologyAssociations(species_code,export_databases,mod_types,genmapp_mod)
def importGeneToOntologyData(species_code,mod,gotype,ontology_type):
program_type,database_dir = unique.whatProgramIsThis()
if gotype == 'nested':
geneGO_import_dir = '/'+database_dir+'/'+species_code+'/nested'
if ontology_type == 'GeneOntology': ontology_type = 'GO.'
mod += '_to_Nested'
elif gotype == 'null':
geneGO_import_dir = '/'+database_dir+'/'+species_code+'/gene-go'
gg = GrabFiles(); gg.setdirectory(geneGO_import_dir)
try: filedir,file = gg.searchdirectory(mod+'-'+ontology_type) ###Identify gene files corresponding to a particular MOD
except Exception:
if gotype == 'nested':
buildNestedAssociations(species_code)
filedir,file = gg.searchdirectory(mod+'-'+ontology_type)
global gene_to_go; x=0
fn=filepath(filedir); gene_to_go={}
for line in open(fn,'rU').xreadlines():
if gotype == 'nested' and x==0: x = 1
else:
if 'Ontology' not in line and 'ID' not in line: ### Header included in input from AP or BioMart
data = cleanUpLine(line)
t = string.split(data,'\t')
if len(t)>1:
try: gene = t[0]; goid = t[1]
except IndexError: print [line],t, 'Ontology-gene relationship not valid...please clean up',filedir;sys.exit()
goid = string.replace(goid,'GO:','') ### Included in input from AP
### If supplied without proceeding zeros
if len(goid)<7 and len(goid)>0:
extended_goid=goid
while len(extended_goid)< 7: extended_goid = '0'+extended_goid
goid = extended_goid
if len(gene)>0 and len(goid)>0:
if ':' not in goid:
goid = 'GO:'+goid
try: gene_to_go[gene].append(goid)
except KeyError: gene_to_go[gene]= [goid]
return gene_to_go
def importGeneMAPPData(species_code,mod):
program_type,database_dir = unique.whatProgramIsThis()
geneMAPP_import_dir = '/'+database_dir+'/'+species_code+'/gene-mapp'
gm = GrabFiles(); gm.setdirectory(geneMAPP_import_dir)
filedir,file = gm.searchdirectory(mod) ### Identify gene files corresponding to a particular MOD
global gene_to_mapp; x = 0
fn=filepath(filedir); gene_to_mapp={}
for line in open(fn,'rU').xreadlines():
data = cleanUpLine(line)
if x==0: x=1
else:
t = string.split(data,'\t')
gene = t[0]; mapp = t[2]
try: gene_to_mapp[gene].append(mapp)
except KeyError: gene_to_mapp[gene]= [mapp]
return gene_to_mapp
def exportCustomPathwayMappings(gene_to_custom,mod,system_codes,custom_sets_folder):
if '.txt' in custom_sets_folder:
export_dir = custom_sets_folder
else:
export_dir = custom_sets_folder+'/CustomGeneSets/custom_gene_set.txt'
#print 'Exporting:',export_dir
try: data = export.ExportFile(export_dir)
except Exception: data = export.ExportFile(export_dir[1:])
for system_code in system_codes:
if system_codes[system_code] == mod: mod_code = system_code
data.write('GeneID\t'+mod+'\tPathway\n'); relationships=0
gene_to_custom = eliminate_redundant_dict_values(gene_to_custom)
for gene in gene_to_custom:
for pathway in gene_to_custom[gene]:
values = string.join([gene,mod_code,pathway],'\t')+'\n'
data.write(values); relationships+=1
data.close()
#print relationships,'Custom pathway-to-ID relationships exported...'
def importGeneCustomData(species_code,system_codes,custom_sets_folder,mod):
print 'Importing custom pathway relationships...'
#print 'Trying to import text data'
gene_to_custom = importTextCustomData(species_code,system_codes,custom_sets_folder,mod)
#print len(gene_to_custom)
#print 'Trying to import gmt data'
gmt_data = parseGMT(custom_sets_folder)
#print 'Trying to import gpml data'
gpml_data,pathway_db = parseGPML(custom_sets_folder)
#print 'Trying to import biopax data'
biopax_data = parseBioPax(custom_sets_folder)
#print 'Unifying gene systems for biopax'
gene_to_BioPax = unifyGeneSystems(biopax_data,species_code,mod); #print len(gene_to_BioPax)
#print 'Unifying gene systems for WP'
gene_to_WP = unifyGeneSystems(gpml_data,species_code,mod); #print len(gene_to_WP)
#print 'Unifying gene systems for gmt'
gene_to_GMT = unifyGeneSystems(gmt_data,species_code,mod); #print len(gene_to_WP)
#print 'Combine WP-biopax'
gene_to_xml = combineDBs(gene_to_WP,gene_to_BioPax)
#print 'Combine xml-text'
gene_to_custom = combineDBs(gene_to_xml,gene_to_custom)
#print 'Combine gmt-other'
gene_to_custom = combineDBs(gene_to_GMT,gene_to_custom)
exportCustomPathwayMappings(gene_to_custom,mod,system_codes,custom_sets_folder)
"""
### Combine WikiPathway associations with the custom
try: gene_to_mapp = importGeneMAPPData(species_code,mod)
except Exception: gene_to_mapp = {}
for gene in gene_to_mapp:
for mapp in gene_to_mapp[gene]:
try: gene_to_custom[gene].append(mapp)
except KeyError: gene_to_custom[gene]= [mapp]"""
return gene_to_custom
def importTextCustomData(species_code,system_codes,custom_sets_folder,mod):
program_type,database_dir = unique.whatProgramIsThis()
gm = GrabFiles(); gm.setdirectory(custom_sets_folder); system = None
filedirs = gm.getAllFiles('.txt')
global gene_to_custom; gene_to_custom={}
for filedir in filedirs:
try:
file = string.split(filedir,'/')[-1]
print "Reading custom gene set",filedir
fn=filepath(filedir); x = 1
for line in open(fn,'rU').xreadlines():
data = cleanUpLine(line)
if x==0: x=1
else:
x+=1
t = string.split(data,'\t')
try:
gene = t[0]; mapp = t[2]; system_code = t[1]
if system_code in system_codes: system = system_codes[system_code]
else:
if x == 3: print system_code, "is not a recognized system code. Skipping import of",file; break
except Exception:
if len(t)>0:
gene = t[0]
if len(t)==1: ### Hence, no system code is provided and only one gene-set is indicated per file
source_data = predictIDSource(t[0],system_codes)
if len(source_data)>0: system = source_data; mapp = file
else:
if x == 3: print file, 'is not propperly formatted (skipping import of relationships)'; break
elif len(t)==2:
if t[1] in system_codes: system = system_codes[t[1]]; mapp = file ### Hence, system code is provided by only one gene-set is indicated per file
else:
source_data = predictIDSource(t[0],system_codes)
if len(source_data)>0: system = source_data; mapp = t[1]
else:
if x == 3: print file, 'is not propperly formatted (skipping import of relationships)'; break
else: continue ### Skip line
try: gene_to_custom[gene].append(mapp)
except KeyError: gene_to_custom[gene]= [mapp]
#print [system, mod, len(gene_to_custom)]
### If the system code is not the MOD - Convert to the MOD
if (system != mod) and (system != None):
gene_to_custom2={}
mod_source = 'hide',mod+'-'+system+'.txt'
try: gene_to_source_id = getGeneToUid(species_code,mod_source)
except Exception: print mod_source,'relationships not found. Skipping import of',file; break
source_to_gene = OBO_import.swapKeyValues(gene_to_source_id)
if system == 'Symbol': source_to_gene = lowerAllIDs(source_to_gene)
for source_id in gene_to_custom:
original_source_id = source_id ### necessary when Symbol
if system == 'Symbol': source_id = string.lower(source_id)
if source_id in source_to_gene:
for gene in source_to_gene[source_id]:
gene_to_custom2[gene] = gene_to_custom[original_source_id]
gene_to_custom = gene_to_custom2
except Exception:
print file, 'not formatted propperly!'
return gene_to_custom
def grabFileRelationships(filename):
filename = string.replace(filename,'.txt','')
system1,system2 = string.split(filename,'-')
return system1,system2
def importUidGeneSimple(species_code,mod_source):
program_type,database_dir = unique.whatProgramIsThis()
geneUID_import_dir = '/'+database_dir+'/'+species_code+'/uid-gene'
ug = GrabFiles(); ug.setdirectory(geneUID_import_dir)
filedir,file = ug.searchdirectory(mod_source) ### Identify gene files corresponding to a particular MOD
gene_to_uid={}; x = 0
fn=filepath(filedir)
for line in open(fn,'rU').xreadlines():
data = cleanUpLine(line)
if x==0: x=1
else:
t = string.split(data,'\t')
gene = t[0]; uid = t[1]
try: gene_to_uid[gene].append(uid)
except KeyError: gene_to_uid[gene]= [uid]
return gene_to_uid
def augmentEnsemblGO(species_code):
ontology_type = 'GeneOntology'
entrez_ens = importUidGeneSimple(species_code,'EntrezGene-Ensembl')
try: ens_go=importGeneToOntologyData(species_code,'Ensembl','null',ontology_type)
except Exception: ens_go = {}
try: entrez_go=importGeneToOntologyData(species_code,'Entrez','null',ontology_type)
except Exception: entrez_go = {}
ens_go_translated = {}
for entrez in entrez_go:
if entrez in entrez_ens:
if len(entrez_ens[entrez])<3: ### Limit bad associations
#print entrez,entrez_ens[entrez];kill
for ens in entrez_ens[entrez]: ens_go_translated[ens] = entrez_go[entrez]
### Add these relationships to the original
for ens in ens_go_translated:
try: ens_go[ens] = unique.unique(ens_go[ens] + ens_go_translated[ens])
except Exception: ens_go[ens] = ens_go_translated[ens]
program_type,database_dir = unique.whatProgramIsThis()
export_dir = database_dir+'/'+species_code+'/gene-go/Ensembl-GeneOntology.txt'
print 'Exporting augmented:',export_dir
try: data = export.ExportFile(export_dir)
except Exception: data = export.ExportFile(export_dir[1:])
data.write('GeneID'+'\tGO ID\n')
for gene in ens_go:
for goid in ens_go[gene]:
values = string.join([gene,goid],'\t')+'\n'
data.write(values)
data.close()
### Build Nested
export_databases = 'no'; genmapp_mod = 'Ensembl'
full_path_db,path_id_to_goid,null = OBO_import.buildNestedOntologyAssociations(species_code,export_databases,['Ensembl'],genmapp_mod)
def importOntologyUIDGeneData(species_code,mod_source,gene_to_go,denominator_source_ids):
program_type,database_dir = unique.whatProgramIsThis()
geneUID_import_dir = '/'+database_dir+'/'+species_code+'/uid-gene'
ug = GrabFiles(); ug.setdirectory(geneUID_import_dir)
#print "begining to parse",mod_source, 'from',geneGO_import_dir
filedir,filename = ug.searchdirectory(mod_source) ###Identify gene files corresponding to a particular MOD
fn=filepath(filedir); uid_to_go={}; count=0; x=0
uid_system,gene_system = grabFileRelationships(filename)
for line in open(fn,'rU').xreadlines(): count+=1
original_increment = int(count/10); increment = original_increment
for line in open(fn,'rU').xreadlines():
data = cleanUpLine(line); x+=1
if program_type == 'GO-Elite':
if x == increment: increment+=original_increment; print '*',
t = string.split(data,'\t')
uid = t[1]; gene = t[0]
try:
if len(denominator_source_ids)>0:
null=denominator_source_ids[uid] ### requires that the source ID be in the list of analyzed denominators
goid = gene_to_go[gene]
y = GeneRelationships(uid,gene,goid,uid_system,gene_system)
try: uid_to_go[uid].append(y)
except KeyError: uid_to_go[uid] = [y]
except Exception: null=[]
try:
if len(denominator_source_ids)>0:
null=denominator_source_ids[uid] ### requires that the source ID be in the list of analyzed denominators
except Exception: null=[]
return uid_to_go,uid_system
def importGeneSetUIDGeneData(species_code,mod_source,gene_to_mapp,denominator_source_ids):
program_type,database_dir = unique.whatProgramIsThis()
geneUID_import_dir = '/'+database_dir+'/'+species_code+'/uid-gene'
ug = GrabFiles(); ug.setdirectory(geneUID_import_dir)
#print "begining to parse",mod_source, 'from',geneGO_import_dir
filedir,filename = ug.searchdirectory(mod_source) ###Identify gene files corresponding to a particular MOD
fn=filepath(filedir); uid_to_mapp={}; count=0; x=0
uid_system,gene_system = grabFileRelationships(filename)
for line in open(fn,'rU').xreadlines(): count+=1
original_increment = int(count/10); increment = original_increment
for line in open(fn,'rU').xreadlines():
data = cleanUpLine(line); x+=1
if program_type == 'GO-Elite':
if x == increment: increment+=original_increment; print '*',
t = string.split(data,'\t')
uid = t[1]; gene = t[0]
try:
if len(denominator_source_ids)>0:
null=denominator_source_ids[uid] ### requires that the source ID be in the list of analyzed denominators
mapp_name = gene_to_mapp[gene]
y = GeneRelationships(uid,gene,mapp_name,uid_system,gene_system)
try: uid_to_mapp[uid].append(y)
except KeyError: uid_to_mapp[uid] = [y]
except Exception: null=[]
return uid_to_mapp,uid_system
def eliminate_redundant_dict_values(database):
db1={}
for key in database: list = unique.unique(database[key]); list.sort(); db1[key] = list
return db1
def getGeneToUid(species_code,mod_source):
if 'hide' in mod_source: show_progress, mod_source = mod_source
else: show_progress = 'yes'
program_type,database_dir = unique.whatProgramIsThis()
import_dir = '/'+database_dir+'/'+species_code+'/uid-gene'
ug = GrabFiles(); ug.setdirectory(import_dir)
filedir,filename = ug.searchdirectory(mod_source) ###Identify gene files corresponding to a particular MOD
fn=filepath(filedir); gene_to_uid={}; count = 0; x=0
uid_system,gene_system = grabFileRelationships(filename)
for line in open(fn,'r').xreadlines(): count+=1
original_increment = int(count/10); increment = original_increment
for line in open(fn,'rU').xreadlines():
data = cleanUpLine(line); x+=1
if x == increment and show_progress == 'yes':
increment+=original_increment; print '*',
t = string.split(data,'\t')
uid = t[1]; gene = t[0]
try: gene_to_uid[gene].append(uid)
except KeyError: gene_to_uid[gene] = [uid]
gene_to_uid = eliminate_redundant_dict_values(gene_to_uid)
return gene_to_uid
def getGeneToUidNoExon(species_code,mod_source):
gene_to_uid = getGeneToUid(species_code,mod_source)
try:
probeset_db = simpleExonImporter(species_code); filtered_db={}
for gene in gene_to_uid:
for uid in gene_to_uid[gene]:
try: probeset_db[uid]
except KeyError: ### Only inlcude probesets not in the exon database
try: filtered_db[gene].append(uid)
except KeyError: filtered_db[gene] = [uid]
return filtered_db
except Exception: return gene_to_uid
def simpleExonImporter(species_code):
filename = 'AltDatabase/'+species_code+'/exon/'+species_code+'_Ensembl_probesets.txt'
fn=filepath(filename); probeset_db={}
for line in open(fn,'rU').xreadlines():
data = cleanUpLine(line)
data = string.split(data,'\t')
probeset_db[data[0]]=[]
return probeset_db
def predictIDSource(id,system_codes):
au = ['A','B','C','D','E','F','G','H','I','J','K','L','M','N','O','P','Q','R','S','T','U','V','W','X','Y','Z']
nm = ['1','2','3','4','5','6','7','8','9']
affy_suffix = '_at'; ensembl_prefix = 'ENS'; source_data = ''; id_type = ''
if len(id)>3:
if affy_suffix == id[-3:]: id_type = 'Affymetrix'
elif ensembl_prefix == id[:3]: id_type = 'Ensembl'
elif id[2] == '.': id_type = 'Unigene'
elif (id[0] in au and id[1] in nm) or '_' in id: id_type = 'UniProt'
else:
try:
int_val = int(id)
if len(id) == 7: id_type = 'Affymetrix' ###All newer Affymetrix transcript_cluster_ids and probesets (can be mistaken for EntrezGene IDs)
else: id_type = 'EntrezGene'
except ValueError: null = []
###If the user changes the names of the above id_types, we need to verify that the id_type is in the database
if len(id_type)>0:
for code in system_codes:
if system_codes[code] == id_type: source_data = id_type
return source_data
def addNewCustomSystem(filedir,system,save_option,species_code):
gene_annotations={}
if 'update' in save_option:
print "Importing existing",species_code,system,"relationships."
gene_annotations = importGeneData(species_code,mod)
fn=filepath(filedir); length3=0; lengthnot3=0
for line in open(fn,'rU').xreadlines():
data = cleanUpLine(line)
t = string.split(data,'\t')
if len(t[0])>0:
try:
### Allow the user to only include one column here
gene = t[0]
try: symbol = t[1]
except Exception: symbol = gene
try: name = ''
except Exception: name = ''
s = GeneAnnotations(gene,symbol,name,'')
gene_annotations[gene] = s
except Exception:
print 'Unexpected error in the following line'; print t, filedir;kill
if len(gene_annotations)>0:
print 'Writing new annotation file:',species_code,system
export_dir = 'Databases/'+species_code+'/gene/'+system+'.txt'
data = export.ExportFile(export_dir)
data.write('ID\tSymbol\tDescription\n')
for gene in gene_annotations:
s = gene_annotations[gene]
data.write(string.join([gene,s.Symbol(),s.Description()],'\t')+'\n')
data.close()
return 'exported'
else: return 'not-exported'
def importGeneSetsIntoDatabase(source_file,species_code,mod):
source_filename = export.findFilename(source_file)
export.deleteFolder('BuildDBs/temp') ### Delete any previous data
destination_dir = filepath('BuildDBs/temp/'+source_filename)
export.customFileCopy(source_file,destination_dir)
print destination_dir
custom_sets_folder = export.findParentDir(destination_dir)
import GO_Elite; system_codes,source_types,mod_types = GO_Elite.getSourceData()
gene_to_custom = importGeneCustomData(species_code,system_codes,custom_sets_folder,mod)
return gene_to_custom
def addNewCustomRelationships(filedir,relationship_file,save_option,species_code):
relationship_filename = export.findFilename(relationship_file)
mod,data_type = string.split(relationship_filename,'-')
mod_id_to_related={}
if 'update' in save_option or '.txt' not in filedir:
if 'gene-mapp' in relationship_file:
if '.txt' not in filedir: ### Then these are probably gmt, gpml or owl
mod_id_to_related = importGeneSetsIntoDatabase(filedir,species_code,mod)
else:
mod_id_to_related = importGeneMAPPData(species_code,mod)
elif 'gene-go' in relationship_file:
mod_id_to_related=importGeneToOntologyData(species_code,mod,'null',data_type)
else: mod_id_to_related=importUidGeneSimple(species_code,mod+'-'+data_type+'.txt')
fn=filepath(filedir); length2=0; lengthnot2=0
if '.txt' in fn:
for line in open(fn,'rU').xreadlines():
data = cleanUpLine(line)
t = string.split(data,'\t')
data_length = len(t)
### Check the length of the input file
if data_length==2 or data_length==3: #This can occur if users load a custom gene set file created by GO-Elite
length2+=1
if len(t[0])>0 and len(t[-1])>0:
if ',' in t[0]: keys = string.split(t[0],',') ### These can be present in the MGI phenotype ontology gene association files
else: keys = [t[0]]
for key in keys:
try: mod_id_to_related[key].append(t[-1])
except KeyError: mod_id_to_related[key] = [t[-1]]
else: lengthnot2+=1
else: length2+=1
if length2>lengthnot2:
###Ensure that the input file is almost completely 2 columns
print 'Writing new relationship file:',species_code,relationship_file
if 'Databases' in relationship_file:
export_dir = relationship_file ### Sometimes we just include the entire path for clarification
if 'mapp' in export_dir: data_type = 'MAPP'
else: data_type = 'Ontology'
elif data_type == 'MAPP': export_dir = 'Databases/'+species_code+'/gene-mapp/'+relationship_file+'.txt'
elif data_type == 'Ontology':
export_dir = 'Databases/'+species_code+'/gene-go/'+relationship_file+'.txt'
else: export_dir = 'Databases/'+species_code+'/uid-gene/'+relationship_file+'.txt'
if data_type == 'Ontology':
### To trigger a rebuild of the ontology nested, must deleted the existing nested for this ontology
nested_path = filepath(string.replace(relationship_file,'gene-go','nested'))
nested_path = string.replace(nested_path,'GeneOntology','GO')
obo=string.split(nested_path,'-')[-1]
nested_path = string.replace(nested_path,'-'+obo,'_to_Nested-'+obo) ### Will break if there is another '-' in the name
try: os.remove(nested_path)
except Exception: null=[]
data = export.ExportFile(export_dir)
data.write('ID\t\tRelated ID\n')
for mod_id in mod_id_to_related:
for related_id in mod_id_to_related[mod_id]:
if data_type == 'MAPP': ###Pathway files have 3 rather than 2 columns
data.write(string.join([mod_id,'',related_id],'\t')+'\n')
else: data.write(string.join([mod_id,related_id],'\t')+'\n')
data.close()
return 'exported'
else: return 'not-exported'
def importUIDsForMAPPFinderQuery(filedir,system_codes,return_uid_values):
fn=filepath(filedir); x=0; uid_list = {}; uid_value_db = {}; source_data = ''; value_len_list = []
system_code_found = 'no'; first_uid=''; first_line = 'no lines in file'; underscore = False
source_data_db={}; error=''; bad_alphanumerics=0
for line in open(fn,'rU').xreadlines():
data = cleanUpLine(line)
t = string.split(data,'\t')
tnull = string.join(t,''); ### Indicates whether the line is blank - even if tabs are present
alphanumeric = string.join(re.findall(r"\w",data))
if x==0:
x=1; value_headers = t[1:]
else:
if len(tnull)>0:
if len(alphanumeric)==0: bad_alphanumerics+=1
first_uid = t[0]; first_line = t
x+=1; source_data_read = ''
uid = t[0]
uid = string.replace(uid,'---','') ###occurs for Affymetrix gene ID input
uids = string.split(uid,' /// ') ###occurs for Affymetrix gene IDs where there are multiple associations
for uid in uids:
if len(uid) > 0 and uid != ' ':
if uid[0]==' ': uid = uid[1:]
if uid[-1]==' ': uid = uid[:-1]
all_alphanumeric = re.findall(r"\w",uid)
if len(all_alphanumeric)>0: uid_list[uid]=[]
if len(t)>1: ###Therefore there is either gene values or system data associated with imported primary IDs
try:
system = t[1]
system = string.replace(system,' ','')
if system in system_codes:
#if len(source_data)<1: source_data = system_codes[t[1]]; system_code_found = 'yes'
source_data = system_codes[system]; system_code_found = 'yes'; source_data_read = source_data
if len(t)==2: values = [] ###Therefore, there is system code data but no associated values
else: values = t[2:]
else: values = t[1:]
value_len_list.append(len(values))
if len(values)>0:
if len(uid) > 0 and uid != ' ':
all_alphanumeric = re.findall(r"\w",uid)
if len(all_alphanumeric)>0: uid_value_db[uid] = values
if '_' in uid: underscore=True
except Exception:
source_data_read =''
if len(source_data_read)<1:
for uid in uids:
source_data_read = predictIDSource(uid,system_codes)
if len(source_data_read)>0: source_data = source_data_read
if len(source_data_read)>0: source_data_db[source_data_read]=[]
first_uid = string.replace(first_uid,'Worksheet','!!!!')
filenames = string.split(filedir,'/'); filename = filenames[-1]
if '!!!!' in first_uid:
error = 'WARNING!!! There appears to be a formatting file issue with the file:\n"'+filename+'"\nPlease correct and re-run (should be tab-delimited text).'
elif len(source_data_db)>1:
error = 'WARNING!!! There is more than one gene system (e.g., Ensembl and EntrezGene) in the file:\n"'+filename+'"\nPlease correct and re-run.'
elif source_data == '':
error = 'WARNING!!! No System Code identified in:\n"'+filename+'"\nPlease provide in input text file and re-run.\n If System Code column present, the file format may be incorrect.'
try: error +='Possible system code: '+system+' not recognized.'
except Exception: None
elif x>0 and bad_alphanumerics>1:
error = 'WARNING!!! Invalid text file encoding found (invalid alphanumeric values). Please resave text file in a standard tab-delimited format'
if underscore and system == 'S':
error += '\nSwissProt IDs of the type P53_HUMAN (P04637) are not recognized as propper UniProt IDs. Consider using alternative compatible IDs (e.g., P04637).'
if return_uid_values == 'yes' and len(value_len_list)>0:
value_len_list.sort(); longest_value_list = value_len_list[-1]
for uid in uid_value_db:
values = uid_value_db[uid]
if len(values) != longest_value_list:
while len(values) != longest_value_list: values.append('')
if system_code_found == 'yes': value_headers = value_headers[1:]###Therfore, column #2 has system code data
return uid_list,uid_value_db,value_headers,error
elif return_uid_values == 'yes': return uid_list,uid_value_db,value_headers,error
else: return uid_list,source_data,error
def getAllDenominators(denom_search_dir,system_codes):
import mappfinder; denom_all_db={}
m = GrabFiles(); m.setdirectory(denom_search_dir)
if len(denom_search_dir)>0: denom_dir_list = mappfinder.readDirText(denom_search_dir)
else: denom_dir_list = []
for denominator_file in denom_dir_list:
gene_file_dir, gene_file = m.searchdirectory(denominator_file)
denom_gene_list,source_data_denom,error_message = importUIDsForMAPPFinderQuery(denom_search_dir+'/'+gene_file,system_codes,'no')
for i in denom_gene_list: denom_all_db[i]=[]
return denom_all_db
def grabNestedGeneToOntologyAssociations(species_code,mod,source_data,system_codes,denom_search_dir,ontology_type):
program_type,database_dir = unique.whatProgramIsThis()
global printout; printout = 'yes'
gene_annotations = importGeneData(species_code,mod)
### Filter source IDs for those in all user denominator files
denominator_source_ids = getAllDenominators(denom_search_dir,system_codes)
try: gene_to_go = importGeneToOntologyData(species_code,mod,'nested',ontology_type)
except Exception:
print "Warning...the MOD you have selected:",mod,ontology_type,"is missing the appropriate relationship files",
print "necessary to run GO-Elite. Either replace the missing files ("+database_dir+'/'+species_code+') or select a different MOD at runtime.'
print 'Exiting program.'; sys.exit()
if source_data != mod:
mod_source = mod+'-'+source_data+'.txt'
uid_to_go,uid_system = importOntologyUIDGeneData(species_code,mod_source,gene_to_go,denominator_source_ids)
gene_to_go = convertGeneDBtoGeneObjects(gene_to_go,mod)
for gene in gene_to_go: uid_to_go[gene] = gene_to_go[gene]
return uid_to_go, uid_system, gene_annotations
else:
gene_to_go = convertGeneDBtoGeneObjects(gene_to_go,mod)
return gene_to_go, mod, gene_annotations
def grabNestedGeneToPathwayAssociations(species_code,mod,source_data,system_codes,custom_sets_folder,denom_search_dir,ontology_type):
program_type,database_dir = unique.whatProgramIsThis()
global printout; printout = 'yes'
gene_annotations = importGeneData(species_code,mod)
### Filter source IDs for those in all user denominator files
denominator_source_ids = getAllDenominators(denom_search_dir,system_codes)
try:
if ontology_type == 'WikiPathways': ontology_type = 'MAPP'
gene_to_mapp = importGeneMAPPData(species_code,mod+'-'+ontology_type)
except Exception:
### If a custom collection, shouldn't be found in the gene-mapp directory
try: gene_to_mapp = importGeneCustomData(species_code,system_codes,custom_sets_folder,mod)
except Exception: gene_to_mapp = {}
if source_data != mod:
mod_source = mod+'-'+source_data+'.txt'
uid_to_mapp,uid_system = importGeneSetUIDGeneData(species_code,mod_source,gene_to_mapp,denominator_source_ids)
gene_to_mapp = convertGeneDBtoGeneObjects(gene_to_mapp,mod)
### Combine the gene and UID data in order to later filter out gene's linked to an arrayid that are not linked to the pathway (common for some pathways)
for gene in gene_to_mapp: uid_to_mapp[gene] = gene_to_mapp[gene]
return uid_to_mapp, uid_system, gene_annotations
else:
gene_to_mapp = convertGeneDBtoGeneObjects(gene_to_mapp,mod)
return gene_to_mapp, mod, gene_annotations
def convertGeneDBtoGeneObjects(gene_to_pathway,mod):
gene_to_pathway_objects={}
for gene in gene_to_pathway:
pathways = gene_to_pathway[gene]
y = GeneRelationships(gene,gene,pathways,mod,mod) ###Note: in the analagous function for unique_ids, pathways is called pathway but is still a list (correct this at some point)
try: gene_to_pathway_objects[gene].append(y)
except KeyError: gene_to_pathway_objects[gene] = [y]
return gene_to_pathway_objects
def matchInputIDsToGOEliteTerms(gene_file_dir,go_elite_output_dir,system_codes,mappfinder_file,nested_collapsed_go_tree,uid_to_go,gene_annotations,full_go_name_db,uid_system,prev_combined_associations,prev_combined_gene_ranking):
global combined_associations ###This variable exists as a global under GO-Elite main as well, but needs to be re-initialized here as well
global combined_gene_ranking; global value_headers
combined_associations = prev_combined_associations
combined_gene_ranking = prev_combined_gene_ranking
try: uid_list,uid_value_db,value_headers,error = importUIDsForMAPPFinderQuery(gene_file_dir,system_codes,'yes')
except Exception: uid_list,uid_value_db,value_headers,error = importUIDsForMAPPFinderQuery(gene_file_dir[1:],system_codes,'yes')
filename_list = string.split(mappfinder_file,"/")
filename = filename_list[-1]
related_goid_db={}
for uid in uid_list:
try:
for y in uid_to_go[uid]: ### Different instances for each gene associated with the uid (e.g. multiple probeset to gene associations)
try: related_goids = y.GOIDs()
except Exception: print uid, y;kill
###If there are numerical values associated with the primary unique ID, store these with the GeneRelationships instance data
if uid in uid_value_db:
uid_values = uid_value_db[uid]
y.setUIDValues(uid_values)
else: y.setUIDValues([])
for goid in related_goids:
if goid in nested_collapsed_go_tree:
try: related_goid_db[goid].append(y)
except KeyError: related_goid_db[goid] = [y] ###store probesets linked to each goid
except Exception: null=[]
###Reduce related_goid_db to unique GO to gene relationships (rather than unique ids)
for goid in related_goid_db:
gene_db={}
for y in related_goid_db[goid]:
if y.GeneID() in uid_to_go: ### Ensures that genes not in the GO term, but linked to an array id are excluded
try: gene_db[y.GeneID()].append(y)
except KeyError: gene_db[y.GeneID()]=[y]
related_goid_db[goid] = gene_db
go_elite_gene_associations={}; gene_uid_data={}
for parent_goid in nested_collapsed_go_tree:
go_associated = [] ###Add all genes belonging to the parent
if parent_goid in related_goid_db:
gene_db = related_goid_db[parent_goid]
for gene in gene_db: ###Too complicated to store this information in the next list: Create a database a call later
gene_uid_data[gene] = gene_db[gene]
go_associated.append(gene)
go_associated = unique.unique(go_associated)
go_elite_gene_associations[parent_goid] = go_associated
return exportGeneGOAssociations(filename,go_elite_output_dir,go_elite_gene_associations,gene_uid_data,gene_annotations,full_go_name_db,uid_system)
def matchInputIDsToMAPPEliteTerms(gene_file_dir,go_elite_output_dir,system_codes,mappfinder_file,uid_to_mapp,zscore_mapp,gene_annotations,uid_system,prev_combined_associations,prev_combined_gene_ranking):
global combined_associations ###This variable exists as a global under GO-Elite main as well, but needs to be re-initialized here as well
global combined_gene_ranking; global value_headers
combined_associations = prev_combined_associations
combined_gene_ranking = prev_combined_gene_ranking
try: uid_list,uid_value_db,value_headers,error = importUIDsForMAPPFinderQuery(gene_file_dir,system_codes,'yes')
except Exception: uid_list,uid_value_db,value_headers,error = importUIDsForMAPPFinderQuery(gene_file_dir[1:],system_codes,'yes')
filename_list = string.split(mappfinder_file,"/")
filename = filename_list[-1]
related_mapp_db={}
for uid in uid_list:
if uid in uid_to_mapp: ###Find probesets in input list with MAPP annotataions
for y in uid_to_mapp[uid]:
related_mapps = y.MAPPs()
###If there are numerical values associated with the primary unique ID, store these with the GeneRelationships instance data
if uid in uid_value_db:
uid_values = uid_value_db[uid]
y.setUIDValues(uid_values)
else: y.setUIDValues([])
for mapp in related_mapps:
if mapp in zscore_mapp: ###Thus this MAPP is a GO-elite Local term
try: related_mapp_db[mapp].append(y)
except KeyError: related_mapp_db[mapp] = [y] ###store probesets linked to each MAPP
###Reduce related_mapp_db to unique GO to gene relationships (rather than unique ids)
for mapp in related_mapp_db:
gene_db={}
for y in related_mapp_db[mapp]:
if y.GeneID() in uid_to_mapp: ### Ensures that genes not on the pathway, but linked to an array id are excluded
try: gene_db[y.GeneID()].append(y)
except KeyError: gene_db[y.GeneID()]=[y]
related_mapp_db[mapp] = gene_db
go_elite_gene_associations={}; gene_uid_data={}
for mapp in zscore_mapp: ###Redundant with above, but preserving the structure of the analagous GO processing function
mapp_associated = []
if mapp in related_mapp_db:
gene_db = related_mapp_db[mapp]
for gene in gene_db: ###Too complicated to store this information in the next list: Create a database a call later
gene_uid_data[gene] = gene_db[gene]
mapp_associated.append(gene)
mapp_associated = unique.unique(mapp_associated)
go_elite_gene_associations[mapp] = mapp_associated
full_mapp_name_db={}
return exportGeneGOAssociations(filename,go_elite_output_dir,go_elite_gene_associations,gene_uid_data,gene_annotations,full_mapp_name_db,uid_system)
def identifyRedundantPathways(go_elite_gene_associations):
for goid_or_mapp1 in go_elite_gene_associations:
for goid_or_mapp2 in go_elite_gene_associations:
if goid_or_mapp1 != goid_or_mapp2:
genes1 = go_elite_gene_associations[goid_or_mapp1]; total_genes1 = len(genes1)
genes2 = go_elite_gene_associations[goid_or_mapp2]; total_genes1 = len(genes2)
common_genes1 = []; common_genes2 = []
for gene in genes1:
if gene in genes2: common_genes1.append(gene)
for gene in genes2:
if gene in gene1: common_genes2.append(gene)
###Can get recipricol hits (unique_common_genes1 = unique_common_genes2), or hits in just one
def exportGeneGOAssociations(filename,go_elite_output_dir,go_elite_gene_associations,gene_uid_data,gene_annotations,full_go_name_db,uid_system):
gene_annotation_filename = filename[0:-4]+ '_' +uid_system+"-gene-associations.txt"
gene_ranking_filename = filename[0:-4]+ '_' +uid_system+"-Gene-Rankings.txt"
original_gene_annotation_filename = gene_annotation_filename
original_gene_ranking_filename = gene_ranking_filename
if len(go_elite_output_dir) == 0:
gene_annotation_filename = "output/gene_associations/" + gene_annotation_filename
gene_ranking_filename = "output/gene_associations/" + gene_ranking_filename
new_dir = 'output/gene_associations'
else:
gene_annotation_filename = go_elite_output_dir+'/gene_associations/'+gene_annotation_filename
gene_ranking_filename = go_elite_output_dir+'/gene_associations/'+gene_ranking_filename
new_dir = go_elite_output_dir
go_gene_annotation_db={}; go_values_db={}
data = export.ExportFile(gene_annotation_filename)
title = ['GeneID','Symbol','Description','UID-str','parent GOID','GO Name']+value_headers
title = string.join(title,'\t')
combined_associations[original_gene_annotation_filename]=[''] ###create a new key (filename) and append all entries to it
combined_associations[original_gene_annotation_filename].append(title)
data.write(title+'\n'); gene_centric_results={}; gene_to_goid={}
for parent_goid in go_elite_gene_associations:
try: go_name = full_go_name_db[parent_goid]
except KeyError: go_name = parent_goid
for gene in go_elite_gene_associations[parent_goid]:
uid_data = gene_uid_data[gene]
uid_list=[]; uid_value_db={}; indexes=[]; uid_value_list = []; uid_mean_value_list = []
for y in uid_data:
uid_list.append(y.UniqueID())
index = 0
for value in y.UIDValues():
try: uid_value_db[index].append(value)
except KeyError: uid_value_db[index] = [value]
indexes.append(index); index+=1
indexes = unique.unique(indexes); indexes.sort()
for index in indexes:
value_str = string.join(uid_value_db[index],'|')
try: value_mean = statistics.avg(uid_value_db[index])
except ValueError: value_mean = 0
uid_value_list.append(value_str); uid_mean_value_list.append(value_mean)
uid_str = string.join(uid_list,'|')
try:
s = gene_annotations[gene]
info = [s.GeneID(),s.Symbol(),s.Description(),uid_str,str(parent_goid),go_name]+uid_value_list
except KeyError:
s = GeneAnnotations(gene,'','','')
info = [gene,'','',uid_str,str(parent_goid),go_name]
###Record these associations to include in the main GO-elite output results file
try: go_gene_annotation_db[parent_goid].append(s)
except KeyError: go_gene_annotation_db[parent_goid] = [s]
try: go_values_db[parent_goid].append(uid_mean_value_list)
except KeyError: go_values_db[parent_goid] = [uid_mean_value_list]
try: gene_centric_results[gene,uid_str].append(go_name)
except KeyError: gene_centric_results[gene,uid_str] = [go_name]
try: gene_to_goid[gene].append(parent_goid)
except KeyError: gene_to_goid[gene] = [parent_goid]
try: info = string.join(info,'\t')
except TypeError: print info;kill
combined_associations[original_gene_annotation_filename].append(info)
data.write(info+'\n')
data.close()
#print 'Nested gene associations (re-derived) written for GO-Elite results to:\n',gene_annotation_filename
data = export.ExportFile(gene_ranking_filename); sorted_gene_ranking=[]
title = ['GeneID','Symbol','Description','#GO-Terms Associated','Percent GO-Terms Associated','UID-str','GO Names']
title = string.join(title,'\t')
combined_gene_ranking[original_gene_ranking_filename]=['']
combined_gene_ranking[original_gene_ranking_filename].append(title)
data.write(title+'\n')
for (gene,uid_str) in gene_centric_results:
go_terms = gene_centric_results[(gene,uid_str)]
num_go_terms = len(go_terms); total_go_elite_terms = len(go_elite_gene_associations); percent = float(num_go_terms)/total_go_elite_terms
go_name_str = string.join(go_terms,' // ')
try:
s = gene_annotations[gene]
info = [s.GeneID(),s.Symbol(),s.Description(),str(num_go_terms),str(percent),uid_str,go_name_str]
except KeyError:
info = [gene,'','',str(num_go_terms),str(percent),uid_str,go_name_str]
info = string.join(info,'\t')
sorted_gene_ranking.append((num_go_terms,info))
sorted_gene_ranking.sort(); sorted_gene_ranking.reverse()
for (rank,info) in sorted_gene_ranking:
combined_gene_ranking[original_gene_ranking_filename].append(info)