-
Notifications
You must be signed in to change notification settings - Fork 29
/
LineageProfilerIterate.py
executable file
·4821 lines (4335 loc) · 224 KB
/
LineageProfilerIterate.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
### Based on code from AltAnalyze's LineageProfiler (http://altanalyze.org)
#Author Nathan Salomonis - nsalomonis@gmail.com
#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 script iterates the LineageProfiler algorithm (correlation based classification method) to identify sample types relative to one
of two references given one or more gene models. The main function is runLineageProfiler.
The program performs the following actions:
1) Import a tab-delimited reference expression file with three columns (ID, biological group 1, group 2) and a header row (biological group names)
2) Import a tab-delimited expression file with gene IDs (column 1), sample names (row 1) and normalized expression values (e.g., delta CT values)
3) (optional - import existing models) Import a tab-delimited file with comma delimited gene-models for analysis
4) (optional - find new models) Identify all possible combinations of gene models for a supplied model size variable (e.g., --s 7)
5) Iterate through any supplied or identified gene models to obtain predictions for novel or known sample types
6) Export prediction results for all analyzed models to the folder CellClassification.
7) (optional) Print the top 20 scores and models for all possible model combinations of size --s
"""
import sys, string, shutil
import math
import os.path
import copy
import time
import getopt
try: import scipy
except Exception: pass
import traceback
import warnings
import random
import collections
try: import unique ### Not required (used in AltAnalyze)
except Exception: None
try: import export ### Not required (used in AltAnalyze)
except Exception: None
command_args = string.join(sys.argv,' ')
if len(sys.argv[1:])>1 and '-' in command_args and '--GUI' not in command_args:
runningCommandLine = True
else:
runningCommandLine = False
#from stats_scripts import salstat_stats; reload(salstat_stats)
try:
from scipy import stats
use_scipy = True
except Exception:
use_scipy = False ### scipy is not required but is used as a faster implementation of Fisher Exact Test when present
def filepath(filename):
try: fn = unique.filepath(filename)
except Exception: fn = filename
return fn
def exportFile(filename):
try: export_data = export.ExportFile(filename)
except Exception: export_data = open(filename,'w')
return export_data
def makeUnique(item):
db1={}; list1=[]; k=0
for i in item:
try: db1[i]=[]
except TypeError: db1[tuple(i)]=[]; k=1
for i in db1:
if k==0: list1.append(i)
else: list1.append(list(i))
list1.sort()
return list1
def cleanUpLine(line):
line = string.replace(line,'\n','')
line = string.replace(line,'\c','')
data = string.replace(line,'\r','')
data = string.replace(data,'"','')
return data
def verifyFile(filename):
status = True
try:
fn=filepath(filename)
for line in open(fn,'rU').xreadlines(): status = True;break
except Exception: status = False
return status
def getHeader(filename):
header=[]
fn=filepath(filename)
for line in open(fn,'rU').xreadlines():
data = cleanUpLine(line)
header = string.split(data,'\t')[1:];break
if header[0] == 'row_clusters-flat': ### remove this cluster annotation column
header = header[1:]
return header
def getGroupsFromExpFile(filename):
""" Import cells and clusters from an Expression Heatmap """
cluster_db = collections.OrderedDict()
fn=filepath(filename)
count=0
for line in open(fn,'rU').xreadlines():
data = cleanUpLine(line)
t = string.split(data,'\t')[1:]
if count==0: header = t
else: clusters = t
count+=1
if count == 2: ### Exit on the 3rd row
break
if header[0] == 'row_clusters-flat': ### remove this cluster annotation column
header = header[1:]
clusters = clusters[1:]
index = 0
for cell in header:
cluster_db[cell] = clusters[index]
index+=1
return cluster_db
def getDataMatrix(filename):
header=[]
matrix={}
fn=filepath(filename)
for line in open(fn,'rU').xreadlines():
data = cleanUpLine(line)
if len(header)==0:
header = string.split(data,'\t')[1:]
else:
t = string.split(data,'\t')
values = map(float,t[1:])
matrix[t[0]]=values
return header,matrix
def returnLargeGlobalVars():
### Prints all large global variables retained in memory (taking up space)
all = [var for var in globals() if (var[:2], var[-2:]) != ("__", "__")]
for var in all:
try:
if len(globals()[var])>1:
print var, len(globals()[var])
except Exception: null=[]
def clearObjectsFromMemory(db_to_clear):
db_keys={}
for key in db_to_clear: db_keys[key]=[]
for key in db_keys:
try: del db_to_clear[key]
except Exception:
try:
for i in key: del i ### For lists of tuples
except Exception: del key ### For plain lists
def int_check(value):
val_float = float(value)
val_int = int(value)
if val_float == val_int:
integer_check = 'yes'
if val_float != val_int:
integer_check = 'no'
return integer_check
def IQR(array):
k1 = 75
k2 = 25
array.sort()
n = len(array)
value1 = float((n*k1)/100)
value2 = float((n*k2)/100)
if int_check(value1) == 'no':
k1_val = int(value1) + 1
if int_check(value1) == 'yes':
k1_val = int(value1)
if int_check(value2) == 'no':
k2_val = int(value2) + 1
if int_check(value2) == 'yes':
k2_val = int(value2)
try: median_val = scipy.median(array)
except Exception: median_val = Median(array)
upper75th = array[k1_val]
lower25th = array[k2_val]
int_qrt_range = upper75th - lower25th
T1 = lower25th-(1.5*int_qrt_range)
T2 = upper75th+(1.5*int_qrt_range)
return lower25th,median_val,upper75th,int_qrt_range,T1,T2
class IQRData:
def __init__(self,maxz,minz,medz,iq1,iq3):
self.maxz = maxz; self.minz = minz
self.medz = medz; self.iq1 = iq1
self.iq3 = iq3
def Max(self): return self.maxz
def Min(self): return self.minz
def Medium(self): return self.medz
def IQ1(self): return self.iq1
def IQ3(self): return self.iq3
def SummaryValues(self):
vals = string.join([str(self.IQ1()),str(self.Min()),str(self.Medium()),str(self.Max()),str(self.IQ3())],'\t')
return vals
def importGeneModels(filename):
x=0
geneModels={}; thresholds=None
#filename = None ### Override file import with default reference data (hard-coded)
if filename != None:
fn=filepath(filename)
fileRead = open(fn,'rU').xreadlines()
else:
fileRead = defaultGeneModels()
for line in fileRead:
try:
data = cleanUpLine(line)
t = string.split(data,'\t')
except Exception:
t = line
genes = t[0]
genes = string.replace(genes,"'",'')
genes = string.replace(genes,' ',',')
genes = string.split(genes,',')
if t>1:
try: thresholds = map(float,t[1:])
except Exception: thresholds = None
try:
if len(thresholds)==0: thresholds = None
except Exception: pass
models=[]
for gene in genes:
if len(gene)>0:
models.append(gene)
if len(models)>0:
geneModels[tuple(models)] = thresholds
return geneModels
def convertClusterOrderToGroupLabels(filename,refCells=None):
""" Lookup for cluster number to assigned cell-type name (used for cellHarmony centroid alignment)"""
clusterNumber = 0
fn=filepath(filename)
prior_label=None
groupNumber_label_db={}
cell_to_label_db={}
for line in open(fn,'rU').xreadlines():
data = cleanUpLine(line)
t = string.split(data,'\t')
label = t[-1]
cell = t[0]
if refCells != None:
if cell not in refCells: ### skip this cell/header/extra notation
alt_cell = cell+'.Reference'
if alt_cell not in refCells:
continue
else:
ICGS_cluster_number = refCells[alt_cell] ### Use the original cluster number if avaialble
clusterNumber = ICGS_cluster_number
cell_to_label_db[alt_cell]=label
else:
ICGS_cluster_number = refCells[cell] ### Use the original cluster number if avaialble
clusterNumber = ICGS_cluster_number
cell_to_label_db[cell]=label
else:
try: ICGS_cluster_number = t[-2]
except: ICGS_cluster_number = t[-1]
if label != prior_label:
clusterNumber+=1
if str(ICGS_cluster_number) not in groupNumber_label_db:
if label == ICGS_cluster_number:
ICGS_cluster_number = clusterNumber ### If the two columns are the same in the labels file
groupNumber_label_db[str(ICGS_cluster_number)]=label
prior_label = label
#for i in groupNumber_label_db: print [i,groupNumber_label_db[i]]
return groupNumber_label_db, cell_to_label_db
######### Below code deals is specific to this module #########
def runLineageProfiler(species,array_type,exp_input,exp_output,
codingtype,compendium_platform,modelSize=None,customMarkers=False,
geneModels=False,permute=False,useMulti=False,fl=None,label_file=None):
""" This code differs from LineageProfiler.py in that it is able to iterate through the LineageProfiler functions with distinct geneModels
that are either supplied by the user or discovered from all possible combinations. """
#global inputType
global exp_output_file; exp_output_file = exp_output; global targetPlatform
global tissues; global sample_headers; global collapse_override
global analysis_type; global coding_type; coding_type = codingtype
global tissue_to_gene; tissue_to_gene = {}; global platform; global cutoff
global customMarkerFile; global delim; global keyed_by; global pearson_list
global Permute; Permute=permute; global useMultiRef; useMultiRef = useMulti
pearson_list={}
cellHarmony=True
#global tissue_specific_db
try: returnCentroids = fl.ReturnCentroids()
except Exception: returnCentroids = False
print 'platform:',array_type
collapse_override = True
if 'ICGS' in customMarkers or 'MarkerGene' in customMarkers or returnCentroids:
""" When performing cellHarmony, build an ICGS expression reference with log2 TPM values rather than fold """
print 'Converting ICGS folds to ICGS expression values as a reference first...'
try: customMarkers = convertICGSClustersToExpression(customMarkers,exp_input,returnCentroids=returnCentroids,species=species,fl=fl)
except:
print traceback.format_exc()
print "Using the supplied reference file only (not importing raw expression)...Proceeding without differential expression analyses..."
pass
if label_file != None and label_file != '':
refCells = None
reference_exp_file = string.replace(customMarkers,'-centroid.txt','.txt')
status = verifyFile(reference_exp_file) ### Organized cells in the final cluster order
if status==True:
refCells = getGroupsFromExpFile(reference_exp_file)
label_file,cell_to_label_db = convertClusterOrderToGroupLabels(label_file,refCells=refCells)
customMarkerFile = customMarkers
if geneModels == False or geneModels == None: geneModels = []
else:
geneModels = importGeneModels(geneModels)
exp_input = string.replace(exp_input,'\\','/')
exp_input = string.replace(exp_input,'//','/')
exp_output = string.replace(exp_output,'\\','/')
exp_output = string.replace(exp_output,'//','/')
delim = "/"
print '\nRunning cellHarmony analysis on',string.split(exp_input,delim)[-1][:-4]
global correlate_by_order; correlate_by_order = 'no'
global rho_threshold; rho_threshold = -1
global correlate_to_tissue_specific; correlate_to_tissue_specific = 'no'
cutoff = 0.01
global value_type
platform = array_type
if 'stats.' in exp_input:
value_type = 'calls'
else:
value_type = 'expression'
tissue_specific_db={}; sample_headers=[]; tissues=[]
if len(array_type)==2:
### When a user-supplied expression is provided (no ExpressionOutput files provided - importGeneIDTranslations)
array_type, vendor = array_type
if ':' in vendor:
vendor = string.split(vendor,':')[1]
if array_type == 'RNASeq':
vendor = 'Symbol'
platform = array_type
else: vendor = 'Symbol'
if 'RawSplice' in exp_input or 'FullDatasets' in exp_input or coding_type == 'AltExon':
analysis_type = 'AltExon'
if platform != compendium_platform: ### If the input IDs are not Affymetrix Exon 1.0 ST probesets, then translate to the appropriate system
translate_to_genearray = 'no'
targetPlatform = compendium_platform
translation_db = importExonIDTranslations(array_type,species,translate_to_genearray)
keyed_by = 'translation'
else: translation_db=[]; keyed_by = 'primaryID'; targetPlatform = compendium_platform
else:
try:
### Get arrayID to Ensembl associations
if vendor != 'Not needed':
### When no ExpressionOutput files provided (user supplied matrix)
translation_db = importVendorToEnsemblTranslations(species,vendor,exp_input)
else:
translation_db = importGeneIDTranslations(exp_output)
keyed_by = 'translation'
targetPlatform = compendium_platform
analysis_type = 'geneLevel'
except Exception:
translation_db=[]; keyed_by = 'primaryID'; targetPlatform = compendium_platform; analysis_type = 'geneLevel'
targetPlatform = compendium_platform ### Overides above
try: importTissueSpecificProfiles(species,tissue_specific_db)
except Exception:
try:
try:
targetPlatform = 'exon'
importTissueSpecificProfiles(species,tissue_specific_db)
except Exception:
try:
targetPlatform = 'gene'
importTissueSpecificProfiles(species,tissue_specific_db)
except Exception:
targetPlatform = "3'array"
importTissueSpecificProfiles(species,tissue_specific_db)
except Exception:
print 'No compatible compendiums present...'
print traceback.format_exc()
forceError
gene_expression_db, sample_headers = importGeneExpressionValuesSimple(exp_input,translation_db)
### Make sure each sample ID is unique
samples_added=[]; i=1
### Organize sample expression, with the same gene order as the tissue expression set
for s in sample_headers:
if s in samples_added:
samples_added.append(s+"."+str(i)) ### Ensure only unique sample IDs exist
i+=1
else:
samples_added.append(s)
sample_headers = samples_added
pruneTissueSpecific=False
for gene in tissue_specific_db:
if gene not in gene_expression_db:
pruneTissueSpecific = True
break
if pruneTissueSpecific:
tissue_specific_db2={}
for gene in gene_expression_db:
if gene in tissue_specific_db:
tissue_specific_db2[gene] = tissue_specific_db[gene]
elif gene in translation_db:
altGeneID = translation_db[gene]
if altGeneID in tissue_specific_db:
tissue_specific_db2[gene] = tissue_specific_db[altGeneID]
tissue_specific_db = tissue_specific_db2
all_marker_genes=[]
for gene in tissue_specific_db:
all_marker_genes.append(gene)
#print [modelSize]
if len(geneModels)>0:
allPossibleClassifiers = geneModels
elif modelSize == None or modelSize == 'optimize' or modelSize == 'no':
allPossibleClassifiers={}
allPossibleClassifiers[tuple(all_marker_genes)]=None
else:
### A specific model size has been specified (e.g., find all 10-gene models)
allPossibleClassifiers = getRandomSets(all_marker_genes,modelSize)
num=1
all_models=[]
if len(allPossibleClassifiers)<16:
print 'Using:'
for model in allPossibleClassifiers:
print 'model',num, 'with',len(model),'genes' #model
num+=1
all_models+=model
#all_models = unique.unique(all_models)
#print len(all_models);sys.exit()
### This is the main analysis function
print 'Number of reference samples to compare to:',len(tissues)
if len(tissues)<20:
print tissues
if modelSize != 'optimize':
hit_list, hits, fails, prognostic_class_db,sample_diff_z, evaluate_size, prognostic_class1_db, prognostic_class2_db = iterateLineageProfiler(exp_input,
tissue_specific_db, allPossibleClassifiers,translation_db,compendium_platform,modelSize,species,gene_expression_db, sample_headers)
else:
summary_hit_list=[]
try: evaluate_size = len(allPossibleClassifiers[0])
except Exception:
### Occurs when custom models loaded
for i in allPossibleClassifiers:
evaluate_size = len(i); break
hit_list, hits, fails, prognostic_class_db,sample_diff_z, evaluate_size, prognostic_class1_db, prognostic_class2_db = iterateLineageProfiler(exp_input,
tissue_specific_db, allPossibleClassifiers,translation_db,compendium_platform,None,species,gene_expression_db, sample_headers)
while evaluate_size > 4:
hit_list.sort()
top_model = hit_list[-1][-1]
top_model_score = hit_list[-1][0]
"""
try: ### Used for evaluation only - gives the same top models
second_model = hit_list[-2][-1]
second_model_score = hit_list[-2][0]
if second_model_score == top_model_score:
top_model = second_model_score ### Try this
print 'selecting secondary'
except Exception: None
"""
allPossibleClassifiers = [hit_list[-1][-1]]
hit_list, hits, fails, prognostic_class_db,sample_diff_z, evaluate_size, prognostic_class1_db, prognostic_class2_db = iterateLineageProfiler(exp_input,
tissue_specific_db, allPossibleClassifiers,translation_db,compendium_platform,modelSize,species,gene_expression_db, sample_headers)
summary_hit_list+=hit_list
hit_list = summary_hit_list
root_dir = string.join(string.split(exp_output_file,'/')[:-1],'/')+'/'
dataset_name = string.replace(string.split(exp_input,'/')[-1][:-4],'exp.','')
output_classification_file = root_dir+'CellClassification/'+dataset_name+'-CellClassification.txt'
try: os.mkdir(root_dir+'CellClassification')
except Exception: None
export_summary = exportFile(output_classification_file)
models = []
for i in allPossibleClassifiers:
i = string.replace(str(i),"'",'')[1:-1]
models.append(i)
### If multiple class-headers with a common phenotype (different source), combine into a single report
class_list=[]
processed=0
for h in tissues:
if ':' in h and collapse_override==False:
try:
phenotype, source = string.split(h,':')
processed+=1
if phenotype not in class_list:
class_list.append(phenotype)
except Exception: pass
if len(class_list)==2 and len(tissues) == processed and collapse_override==False and cellHarmony==False: ### Ensures all reference headers have : in them
tissue_list = class_list
collapse = True
else:
tissue_list = tissues
collapse = False
print ''
class_headers = map(lambda x: x+' Predicted Hits',tissue_list)
export_header = ['Samples']+class_headers+['Composite Classification Score','Combined Correlation DiffScore','Predicted Class','Max-Rho']
if label_file != None:
export_header.append('CentroidLabel')
headers = string.join(export_header,'\t')+'\n'
export_summary.write(headers)
sorted_results=[] ### sort the results
try: numberOfModels = len(allPossibleClassifiers)
except Exception: numberOfModels = 1
accuracy=[]
ar=[]
non=[]
no_intermediate_accuracy=[]
verboseReport = False
for sample in prognostic_class_db:
if len(tissues)==2:
class1_score = prognostic_class1_db[sample]
class2_score = prognostic_class2_db[sample]
zscore_distribution = map(lambda x: str(x[0]), sample_diff_z[sample])
pearson_max_values = map(lambda x: str(x[1]), sample_diff_z[sample])
dist_list=[]
for i in zscore_distribution:
try: dist_list.append(float(i))
except Exception: None ### Occurs for 'NA's
#try: median_score = scipy.median(dist_list)
#except Exception: median_score = Median(dist_list)
try: sum_score = sum(dist_list) #scipy.median
except Exception: sum_score = sum(dist_list)
correlations=[]
for i in pearson_max_values:
try: correlations.append(float(i))
except Exception: None ### Occurs for 'NA's
median_correlation = scipy.median(correlations)
if median_correlation<0.8 and verboseReport:
print 'Sample: %s has a low median model Pearson correlation coefficient (%s)' % (sample,str(median_correlation))
if verboseReport==False:
print '.',
class_db = prognostic_class_db[sample]
class_scores=[]; class_scores_str=[]; class_scores_refs=[]; collapsed_pheno_scores={}
for tissue in tissues:
if collapse and collapse_override==False:
phenotype,source = string.split(tissue,':')
try: collapsed_pheno_scores[phenotype]+=class_db[tissue]
except Exception: collapsed_pheno_scores[phenotype]=class_db[tissue]
else:
class_scores_str.append(str(class_db[tissue]))
class_scores.append(class_db[tissue])
class_scores_refs.append((class_db[tissue],tissue))
if collapse and collapse_override == False:
for phenotype in tissue_list:
class_scores.append(collapsed_pheno_scores[phenotype]) ### Collapse the scores and report in the original phenotype order
class_scores_str.append(str(collapsed_pheno_scores[phenotype]))
class_scores_refs.append((collapsed_pheno_scores[phenotype],phenotype))
"""
for tissue in tissues:
class_scores_str.append(str(class_db[tissue]))
class_scores.append(class_db[tissue])
class_scores_refs.append((class_db[tissue],tissue))
"""
overall_prog_score = str(max(class_scores)-min(class_scores))
if len(tissues)==2 and cellHarmony==False:
class_scores_str = [str(class1_score),str(class2_score)] ### range of positive and negative scores for a two-class test
if class1_score == 0 and class2_score == 0:
call = 'Intermediate Risk '+ tissues[0]
elif class1_score == numberOfModels:
call = 'High Risk '+ tissues[0]
elif class2_score == numberOfModels:
call = 'Low Risk '+ tissues[0]
elif class1_score == 0:
call = 'Intermediate Risk '+ tissues[0]
elif class2_score == 0:
call = 'Intermediate Risk '+ tissues[0]
else:
call = 'Itermediate Risk '+ tissues[0]
overall_prog_score = str(class1_score-class2_score)
else:
class_scores_refs.sort()
call=class_scores_refs[-1][1] ### This is the reference with the max reported score
if call == tissue_list[-1]: ### Usually the classifier of interest should be listed first in the reference file, not second
overall_prog_score = str(float(overall_prog_score)*-1)
sum_score = sum_score*-1
values = [sample]+class_scores_str+[overall_prog_score,str(sum_score),call]
if label_file != None:
sampleLabel=''
call2 = string.replace(call,'.Reference','')
if call in label_file:
groupLabel = label_file[call] ### This is the unique label for the cluster number
elif call in cell_to_label_db:
groupLabel = cell_to_label_db[call] ### This is the unique label for the cluster number
else:
print [call]
for c in label_file:
print [c]
kill
values = string.join(values+zscore_distribution[:-1]+[str(max(correlations)),groupLabel],'\t')+'\n' ### Export line for cellHarmony classification results
else:
values = string.join(values+zscore_distribution[:-1]+[str(max(correlations))],'\t')+'\n' ### Export line for cellHarmony classification results
if ':' in sample:
sample = string.split(sample,':')[0]
if ':' in call:
call = string.split(call,':')[0]
if call==sample:
accuracy.append(float(1))
if float(overall_prog_score) > 10 or float(overall_prog_score) < -10:
no_intermediate_accuracy.append(float(1))
if 'non' in call: non.append(float(1))
else: ar.append(float(1))
else:
accuracy.append(float(0))
if float(overall_prog_score) > 10 or float(overall_prog_score) < -10:
no_intermediate_accuracy.append(float(0))
if 'non' in call: non.append(float(0))
else: ar.append(float(0))
sorted_results.append([float(overall_prog_score),sum_score,values])
sample_diff_z[sample] = dist_list
if verboseReport:
print len(no_intermediate_accuracy)
print no_intermediate_accuracy
print 'Overall Acuracy:',Average(accuracy)*100
print 'Sensititivity:', sum(ar), len(ar)
print 'Specificity:', sum(non), len(non)
print str(Average(accuracy)*100)+'\t'+str(Average(ar)*100)+'\t'+str(Average(non)*100)+'\t'+str(Average(no_intermediate_accuracy)*100)+'\t'+str(sum(ar))+'\t'+str(len(ar))+'\t'+str(sum(non))+'\t'+str(len(non))
else:
print '\nClassification analysis completed...'
sorted_results.sort()
sorted_results.reverse()
for i in sorted_results:
export_summary.write(i[-1])
export_summary.close()
print '\nResults file written to:',root_dir+'CellClassification/'+dataset_name+'-CellClassification.txt','\n'
hit_list.sort(); hit_list.reverse()
top_hit_list=[]
top_hit_db={}
hits_db={}; fails_db={}
### Only look at the max correlation for each sample
max_pearson_list=[]
for sample in pearson_list:
pearson_list[sample].sort()
for rho in pearson_list[sample][-2:]: ### get the top two correlations
max_pearson_list.append(rho)
avg_pearson_rho = Average(max_pearson_list)
maxPearson = max(max_pearson_list)
try:
for i in sample_diff_z:
zscore_distribution = sample_diff_z[i]
maxz = max(zscore_distribution); minz = min(zscore_distribution)
sample_diff_z[i] = string.join(map(str,zscore_distribution),'\t')
try:
lower25th,medz,upper75th,int_qrt_range,T1,T2 = IQR(zscore_distribution)
if float(maxz)>float(T2): maxz = T2
if float(minz) < float(T1): minz = T1
#iqr = IQRData(maxz,minz,medz,lower25th,upper75th)
#sample_diff_z[i] = iqr
except Exception:
pass
for i in hits:
try: hits_db[i]+=1
except Exception: hits_db[i]=1
for i in fails:
try: fails_db[i]+=1
except Exception: fails_db[i]=1
for i in fails_db:
if i not in hits:
try:
#print i+'\t'+'0\t'+str(fails_db[i])+'\t'+ sample_diff_z[i]
None
except Exception:
#print i
None
except Exception:
pass
exportModelScores = True
if modelSize != False:
#print 'Returning all model overall scores'
hits=[]
for i in hits_db:
hits.append([hits_db[i],i])
hits.sort()
hits.reverse()
for i in hits:
if i[1] in fails_db: fail = fails_db[i[1]]
else: fail = 0
try:
#print i[1]+'\t'+str(i[0])+'\t'+str(fail)+'\t'+sample_diff_z[i[1]]
None
except Exception:
#print i[1]
None
if modelSize == 'optimize': threshold = 80
else: threshold = 0
#print 'threshold:',threshold
for i in hit_list:
if i[0]>threshold:
top_hit_list.append(i[-1])
top_hit_db[tuple(i[-1])]=i[0]
if len(geneModels) > 0 and exportModelScores==False:
for i in hit_list:
#print i[:5],i[-1],i[-2] ### print all
pass
else:
"""
print 'Returning all over 90'
for i in hit_list:
if i[0]>85:
print i[:5],i[-1],i[-2] ### print all
sys.exit()"""
#print 'Top hits'
output_model_file = root_dir+'CellClassification/'+dataset_name+'-ModelScores.txt'
export_summary = exportFile(output_model_file)
print 'Exporting top-scoring models to:',output_model_file
title = 'Classification-Rate\tClass1-Hits\tClass1-Total\tClass2-Hits\tClass2-Total\tModel\tModel-Gene-Number\n'
export_summary.write(title)
for i in hit_list: #hit_list[:500]
overall_scores=[]
for x in i[:5]: overall_scores.append(str(x))
model = string.replace(str(i[-1])[1:-1],"'",'')
values = string.join(overall_scores+[model]+[str(i[-2])],'\t')+'\n'
export_summary.write(values)
export_summary.close()
"""
try:
if hit_list[0][0] == hit_list[20][0]:
for i in hit_list[20:]:
if hit_list[0][0] == i[0]:
print i[:5],i[-1],i[-2]
else: sys.exit()
except Exception: None ### Occurs if less than 20 entries here
"""
print 'Average Pearson correlation coefficient:', avg_pearson_rho
if avg_pearson_rho<0.9 and verboseReport:
print '\n\nWARNING!!!!!!!!!'
print '\tThe average Pearson correlation coefficient for all example models is less than 0.9.'
print '\tYour data may not be comparable to the provided reference (quality control may be needed).\n\n'
elif verboseReport:
print 'No unusual warning.\n'
reference_exp_file = customMarkers
query_exp_file = exp_input
classification_file = output_classification_file
harmonizeClassifiedSamples(species,reference_exp_file, query_exp_file, classification_file,fl=fl)
return top_hit_db
def iterateLineageProfiler(exp_input,tissue_specific_db,allPossibleClassifiers,translation_db,compendium_platform,
modelSize,species,gene_expression_db,sampleHeaders):
classifyBasedOnRho=True
hit_list=[]
### Iterate through LineageProfiler for all gene models (allPossibleClassifiers)
times = 1; k=1000; l=1000; hits=[]; fails=[]; f=0; s=0; sample_diff_z={}; prognostic_class1_db={}; prognostic_class2_db={}
prognostic_class_db={}
begin_time = time.time()
try: evaluate_size=len(allPossibleClassifiers[0]) ### Number of reference markers to evaluate
except Exception:
for i in allPossibleClassifiers: evaluate_size = len(i); break
if modelSize=='optimize':
evaluate_size -= 1
allPossibleClassifiers = getRandomSets(allPossibleClassifiers[0],evaluate_size)
### Determine if we should collapse the entries or not based on common phenotype references
class_list=[]; processed=0; alternate_class={}
for h in tissues:
if ':' in h and 'ENS' not in h:
try:
phenotype, source = string.split(h,':'); processed+=1
if phenotype not in class_list: class_list.append(phenotype)
except Exception: pass
try:
alternate_class[class_list[0]] = class_list[1]
alternate_class[class_list[1]] = class_list[0]
except Exception: pass
if len(class_list)==2 and len(tissues) == processed and collapse_override==False: ### Ensures all reference headers have : in them
tissue_list = class_list; collapse = True
else: tissue_list = tissues; collapse = False
cellHarmony=True
mean_percent_positive=[]
for classifiers in allPossibleClassifiers:
try: thresholds = allPossibleClassifiers[classifiers]
except Exception: thresholds = None
tissue_to_gene={}; expession_subset=[]; sample_headers=[]; classifier_specific_db={}
for gene in classifiers:
try: classifier_specific_db[gene] = tissue_specific_db[gene]
except Exception: None
#print len(gene_expression_db), len(classifier_specific_db), len(expession_subset), len(translation_db)
expession_subset = filterGeneExpressionValues(gene_expression_db,classifier_specific_db,translation_db,expession_subset)
### If the incorrect gene system was indicated re-run with generic parameters
if len(expession_subset)==0:
translation_db=[]; keyed_by = 'primaryID'; targetPlatform = compendium_platform; analysis_type = 'geneLevel'
tissue_specific_db={}
importTissueSpecificProfiles(species,tissue_specific_db)
expession_subset = filterGeneExpressionValues(gene_expression_db,tissue_specific_db,translation_db,expession_subset)
if len(sample_diff_z)==0: ### Do this for the first model examine only
for h in sampleHeaders:
sample_diff_z[h]=[] ### Create this before any data is added, since some models will exclude data for some samples (missing dCT values)
if len(expession_subset)!=len(classifiers): f+=1
#if modelSize=='optimize': print len(expession_subset), len(classifiers);sys.exit()
if (len(expession_subset) != len(classifiers)) and modelSize=='optimize':
print "Please provide a reference set of equal length or smaller to the input analysis set"; kill
#print len(expession_subset), len(classifiers);sys.exit()
if len(expession_subset)==len(classifiers): ### Sometimes a gene or two are missing from one set
s+=1
#print classifiers,'\t',
zscore_output_dir,tissue_scores,sampleHeaders = analyzeTissueSpecificExpressionPatterns(tissue_specific_db,expession_subset,sampleHeaders)
#except Exception: print len(classifier_specific_db), classifiers; error
headers = list(tissue_scores['headers']); del tissue_scores['headers']
if times == k:
end_time = time.time()
print int(end_time-begin_time),'seconds'
k+=l
times+=1; index=0; positive=0; positive_score_diff=0
sample_number = (len(headers)-1)
population1_denom=0; population1_pos=0; population2_pos=0; population2_denom=0
diff_positive=[]; diff_negative=[]
while index < sample_number:
### The scores are now a tuple of (Z-score,original_pearson_rho)
scores = map(lambda x: tissue_scores[x][index], tissue_scores)
zscores_only = map(lambda x: tissue_scores[x][index][0], tissue_scores)
scores_copy = list(scores); scores_copy.sort()
max_pearson_model = scores_copy[-1][1] ### e.g., tumor rho (this is the one we want to QC on later)
min_pearson_model = scores_copy[-2][1] ### e.g., non-tumor rho
diff_rho = max_pearson_model - min_pearson_model
diff_z = (scores_copy[-1][0]-scores_copy[-2][0])*100 ### Diff between the top two scores (z-scores are the first item)
if classifyBasedOnRho == True:
diff_z = diff_rho*10
positive_class=None
j=0
for tissue in tissue_scores:
if ':' in tissue and 'ENS' not in tissue:
group_name = string.split(tissue,':')[0]
else:
group_name = tissue
if scores[j][0] == max(zscores_only):
hit_score = 1; positive_class = tissue
else: hit_score = 0
if len(tissues)>2 or cellHarmony==True:
if group_name+':' in headers[index+1] and hit_score==1:
g = string.split(headers[index+1],':')[0]+':'
if g in group_name+':': ### reciprocol of above
positive+=1
try:
class_db = prognostic_class_db[headers[index+1]]
try: class_db[tissue]+=hit_score
except Exception: class_db[tissue]=hit_score
except Exception:
class_db={}
class_db[tissue]=hit_score
prognostic_class_db[headers[index+1]] = class_db
j+=1
if collapse and collapse_override==False:
phenotype, source = string.split(positive_class,':')
baseline_positive_class = alternate_class[phenotype]+':'+source
denom_score,denom_rho = tissue_scores[baseline_positive_class][index]
old_diff = diff_z
diff_z = scores_copy[-1][0]-denom_score ### Diff between the top two scores of the SAME SOURCE
diff_rho = (max_pearson_model - denom_rho)
if classifyBasedOnRho == True:
diff_z = diff_rho*10
#print headers[index+1], scores_copy[-1]
if len(tissues)==2:
if ':' in headers[index+1]:
pheno = string.split(headers[index+1],':')[0]
else:
pheno = None
diff_z = tissue_scores[tissues[0]][index][0]-tissue_scores[tissues[-1]][index][0] ### z-scores are the first item and pearson-rho is the second
diff_rho = (tissue_scores[tissues[0]][index][1]-tissue_scores[tissues[-1]][index][1])
if classifyBasedOnRho == True:
diff_z = diff_rho*10
if thresholds == None:
threshold1 = 0; threshold2 = 0
else:
threshold1, threshold2 = thresholds ### emperically derived cutoffs provided by the user for each model (e.g., mean+2SD of diff_rho)
if headers[index+1] not in prognostic_class1_db:
prognostic_class1_db[headers[index+1]]=0 ### Create a default value for each sample
if headers[index+1] not in prognostic_class2_db:
prognostic_class2_db[headers[index+1]]=0 ### Create a default value for each sample
if diff_z>threshold1:
prognostic_class1_db[headers[index+1]]+=1
elif diff_z<threshold2:
prognostic_class2_db[headers[index+1]]+=1
if diff_z>0 and (tissues[0] == pheno):
positive+=1; positive_score_diff+=abs(diff_z)
population1_pos+=1; diff_positive.append(abs(diff_z))
hits.append(headers[index+1]) ### see which are correctly classified
elif diff_z<0 and (tissues[-1] == pheno):
positive+=1; positive_score_diff+=abs(diff_z)
population2_pos+=1; diff_positive.append(abs(diff_z))
hits.append(headers[index+1]) ### see which are correctly classified
elif diff_z>0 and (tissues[-1] == pheno): ### Incorrectly classified
diff_negative.append(abs(diff_z))
fails.append(headers[index+1])
elif diff_z<0 and (tissues[0] == pheno): ### Incorrectly classified
#print headers[index+1]
diff_negative.append(abs(diff_z))
fails.append(headers[index+1])
if (tissues[0] == pheno):
population1_denom+=1
else:
population2_denom+=1
sample_diff_z[headers[index+1]].append((diff_z,max([max_pearson_model,min_pearson_model]))) ### Added pearson max here
index+=1
try: percent_positive = (float(positive)/float(index))*100
except ZeroDivisionError:
print 'WARNING!!!! No matching genes. Make sure your gene IDs are the same ID type as your reference.'
forceNoMatchingGeneError
mean_percent_positive.append(percent_positive)
if len(tissues)==2 and cellHarmony==False:
try:
pos = float(population1_pos)/float(population1_denom)
neg = float(population2_pos)/float(population2_denom)
#percent_positive = (pos+neg)/2
except Exception: pos = 0; neg = 0
hit_list.append([percent_positive,population1_pos, population1_denom,population2_pos,population2_denom,[Average(diff_positive),Average(diff_negative)],positive_score_diff,len(classifiers),classifiers])
else:
hit_list.append([percent_positive,len(classifiers),classifiers])
for sample in sample_diff_z:
if len(sample_diff_z[sample]) != (times-1): ### Occurs when there is missing data for a sample from the analyzed model
sample_diff_z[sample].append(('NA','NA')) ### add a null result
#print Average(mean_percent_positive), '\tAverage'
return hit_list, hits, fails, prognostic_class_db, sample_diff_z, evaluate_size, prognostic_class1_db, prognostic_class2_db
def factorial(n):
### Code from http://docs.python.org/lib/module-doctest.html
if not n >= 0:
raise ValueError("n must be >= 0")
if math.floor(n) != n:
raise ValueError("n must be exact integer")
if n+1 == n: # catch a value like 1e300
raise OverflowError("n too large")
result = 1
factor = 2
while factor <= n:
result *= factor
factor += 1
return result
def choose(n,x):
"""Equation represents the number of ways in which x objects can be selected from a total of n objects without regard to order."""
#(n x) = n!/(x!(n-x)!)
f = factorial
result = f(n)/(f(x)*f(n-x))
return result