diff --git a/ClearMap/Analysis/Label.py b/ClearMap/Analysis/Label.py index f8c2f38..5fa544a 100644 --- a/ClearMap/Analysis/Label.py +++ b/ClearMap/Analysis/Label.py @@ -55,7 +55,7 @@ """ -LabelRecord = collections.namedtuple('LabelRecord', 'id, name, acronym, color, parent, collapse'); +LabelRecord = collections.namedtuple('LabelRecord', 'id, name, acronym, color, structureOrder, parent, collapse'); """Structure of a label for a annotated region""" @@ -77,6 +77,7 @@ class LabelInfo(object): names = None; acronyms = None; colors = None; + structureOrder = None; parents = None; levels = None; collapse = None; @@ -92,7 +93,7 @@ def initialize(slf, annotationFile = DefaultAnnotationFile): reader = csv.reader(dfile); #skip header reader.next(); - labels = [LabelRecord._make((int(row[0]), row[1], row[2], [int(row[3]), int(row[4]), int(row[5])], _labelToInt(row[7]), _collapseToBool(row[9]))) for row in reader]; + labels = [LabelRecord._make((int(row[0]), row[1], row[2], [int(row[3]), int(row[4]), int(row[5])], _labelToInt(row[6]), _labelToInt(row[7]), _collapseToBool(row[9]))) for row in reader]; dfile.close(); @@ -100,6 +101,7 @@ def initialize(slf, annotationFile = DefaultAnnotationFile): slf.names = { x.id : x.name for x in labels}; slf.acronyms = { x.id : x.acronym for x in labels}; slf.colors = { x.id : x.color for x in labels}; + slf.structureOrders = { x.id : x.structureOrder for x in labels}; slf.parents = { x.id : x.parent for x in labels}; slf.collapse = { x.id : x.collapse for x in labels}; @@ -125,6 +127,9 @@ def acronym(slf, iid): def color(slf, iid): return slf.colors[iid]; + def structureOrder(slf, iid): + return slf.structureOrders[iid]; + def parent(slf, iid): return slf.parents[iid]; @@ -290,6 +295,10 @@ def labelToAcronym(label): def labelToColor(label): global Label; return [Label.color(x) for x in label]; + +def labelToStructureOrder(label): + global Label; + return [Label.structureOrder(x) for x in label]; diff --git a/ClearMap/Scripts/analyze_template.py b/ClearMap/Scripts/analyze_template.py index 3163f9b..6dad3ac 100644 --- a/ClearMap/Scripts/analyze_template.py +++ b/ClearMap/Scripts/analyze_template.py @@ -1,40 +1,49 @@ # -*- coding: utf-8 -*- """ -Example script to set up the parameters for the image processing pipeline +Example script to set up the parameters for analyzing the data generated in the +process_template scripts. There are 2 parts to this pipeline. One is the fast +and slow voxel analysis. The other is the fast and slow Region of Interest +statistical analysis. """ # import modules: - import ClearMap.Analysis.Statistics as stat -import ClearMap.Analysis.Label as lbl +from ClearMap.Analysis.Label import labelToName, labelToStructureOrder import ClearMap.IO.IO as io import ClearMap.Alignment.Resampling as rsp import numpy, os - +import ClearMap.Analysis.Tools.MultipleComparisonCorrection as FDR # Base Directory, usually where your experiment is saved: baseDirectory = '/home/yourname/experiment' +PathReg = '/home/mtllab/Documents/warping'; +AnnotationFile = os.path.join(PathReg, 'annotation_25_right_fullWD.tif'); - - -# Voxel-based statistics: -######################### - -#Load the data (heat maps generated previously ) -group1 = ['/home/yourname/experiment/sample1/cells_heatmap.tif', +group1vox = ['/home/yourname/experiment/sample1/cells_heatmap.tif', '/home/yourname/experiment/sample2/cells_heatmap.tif', '/home/yourname/experiment/sample3/cells_heatmap.tif', '/home/yourname/experiment/sample4/cells_heatmap.tif'] -group2 = ['/home/yourname/experiment/sample5/cells_heatmap.tif', +group2vox = ['/home/yourname/experiment/sample5/cells_heatmap.tif', '/home/yourname/experiment/sample6/cells_heatmap.tif', '/home/yourname/experiment/sample7/cells_heatmap.tif', '/home/yourname/experiment/sample8/cells_heatmap.tif'] - -g1 = stat.readDataGroup(group1); -g2 = stat.readDataGroup(group2); +group1_roi = ['/home/yourname/experiment/sample1/cells_transformed_to_Atlas.npy', + '/home/yourname/experiment/sample2/cells_transformed_to_Atlas.npy', + '/home/yourname/experiment/sample3/cells_transformed_to_Atlas.npy', + '/home/yourname/experiment/sample4/cells_transformed_to_Atlas.npy'] +group2_roi = ['/home/yourname/experiment/sample5/cells_transformed_to_Atlas.npy', + '/home/yourname/experiment/sample6/cells_transformed_to_Atlas.npy', + '/home/yourname/experiment/sample7/cells_transformed_to_Atlas.npy', + '/home/yourname/experiment/sample8/cells_transformed_to_Atlas.npy'] + +# Voxel-based statistics: +######################### +#Load the data (heat maps generated previously ) +g1 = stat.readDataGroup(group1vox); +g2 = stat.readDataGroup(group2vox); #Generated average and standard deviation maps ############################################## @@ -50,10 +59,6 @@ io.writeData(os.path.join(baseDirectory, 'group2_fast_mean.raw'), rsp.sagittalToCoronalData(g2a)); io.writeData(os.path.join(baseDirectory, 'group2_fast_std.raw'), rsp.sagittalToCoronalData(g2s)); - - - - #Generate the p-values map ########################## #pcutoff: only display pixels below this level of significance @@ -70,44 +75,15 @@ # Regions-based statistics: ########################### +group1i = [fn.replace('cells_transformed_to_Atlas', 'intensities') for fn in group1_roi]; +group2i = [fn.replace('cells_transformed_to_Atlas', 'intensities') for fn in group2_roi]; - -group1 = ['/home/yourname/experiment/sample1/cells_transformed_to_Atlas.npy', - '/home/yourname/experiment/sample2/cells_transformed_to_Atlas.npy', - '/home/yourname/experiment/sample3/cells_transformed_to_Atlas.npy', - '/home/yourname/experiment/sample4/cells_transformed_to_Atlas.npy', - '/home/yourname/experiment/sample5/cells_transformed_to_Atlas.npy'] - - -group2 = ['/home/yourname/experiment/sample6/cells_transformed_to_Atlas.npy', - '/home/yourname/experiment/sample7/cells_transformed_to_Atlas.npy', - '/home/yourname/experiment/sample8/cells_transformed_to_Atlas.npy', - '/home/yourname/experiment/sample9/cells_transformed_to_Atlas.npy', - '/home/yourname/experiment/sample10/cells_transformed_to_Atlas.npy'] - -group1i = [fn.replace('cells_transformed_to_Atlas', 'intensities') for fn in group1]; -group2i = [fn.replace('cells_transformed_to_Atlas', 'intensities') for fn in group2]; - - - -PathReg = '/home/mtllab/Documents/warping'; -AnnotationFile = os.path.join(PathReg, 'annotation_25_right_fullWD.tif'); - - - -#ids, pc1, pc1i = stat.countPointsGroupInRegions(group1, intensityGroup = group1i, withIds = True, labeledImage = lbl.DefaultLabeledImageFile, withCounts = True); -#pc2, pc2i = stat.countPointsGroupInRegions(group2, intensityGroup = group2i, withIds = False, labeledImage = lbl.DefaultLabeledImageFile, withCounts = True); - -ids, pc1, pc1i = stat.countPointsGroupInRegions(group1, intensityGroup = group1i, returnIds = True, labeledImage = AnnotationFile, returnCounts = True, collapse=True); -pc2, pc2i = stat.countPointsGroupInRegions(group2, intensityGroup = group2i, returnIds = False, labeledImage = AnnotationFile, returnCounts = True, collapse=True); - +ids, pc1, pc1i = stat.countPointsGroupInRegions(group1_roi, intensityGroup = group1i, returnIds = True, labeledImage = AnnotationFile, returnCounts = True, collapse=True); +pc2, pc2i = stat.countPointsGroupInRegions(group2_roi, intensityGroup = group2i, returnIds = False, labeledImage = AnnotationFile, returnCounts = True, collapse=True); pvals, psign = stat.tTestPointsInRegions(pc1, pc2, pcutoff = None, signed = True); pvalsi, psigni = stat.tTestPointsInRegions(pc1i, pc2i, pcutoff = None, signed = True, equal_var = True); -import ClearMap.Analysis.Tools.MultipleComparisonCorrection as FDR - - iid = pvalsi < 1; ids0 = ids[iid]; @@ -124,16 +100,16 @@ #make table - -dtypes = [('id','int64'),('mean1','f8'),('std1','f8'),('mean2','f8'),('std2','f8'),('pvalue', 'f8'),('qvalue', 'f8'),('psign', 'int64')]; -for i in range(len(group1)): +dtypes = [('id','int64'),('structureOrder', 'int64'),('mean1','f8'),('std1','f8'),('mean2','f8'),('std2','f8'),('pvalue', 'f8'),('qvalue', 'f8'),('psign', 'int64')]; +for i in range(len(group1_roi)): dtypes.append(('count1_%d' % i, 'f8')); -for i in range(len(group2)): +for i in range(len(group2_roi)): dtypes.append(('count2_%d' % i, 'f8')); dtypes.append(('name', 'a256')); table = numpy.zeros(ids0.shape, dtype = dtypes) table["id"] = ids0; +table["structureOrder"] = labelToStructureOrder(ids0); table["mean1"] = pc1i0.mean(axis = 1)/1000000; table["std1"] = pc1i0.std(axis = 1)/1000000; table["mean2"] = pc2i0.mean(axis = 1)/1000000; @@ -142,11 +118,11 @@ table["qvalue"] = qvalsi0; table["psign"] = psigni0; -for i in range(len(group1)): +for i in range(len(group1_roi)): table["count1_%d" % i] = pc10[:,i]; -for i in range(len(group2)): +for i in range(len(group2_roi)): table["count2_%d" % i] = pc20[:,i]; -table["name"] = lbl.labelToName(ids0); +table["name"] = labelToName(ids0); #sort by qvalue @@ -163,19 +139,18 @@ f.close(); ############################# - - #make table -dtypes = [('id','int64'),('mean1','f8'),('std1','f8'),('mean2','f8'),('std2','f8'),('pvalue', 'f8'),('qvalue', 'f8'),('psign', 'int64')]; -for i in range(len(group1)): +dtypes = [('id','int64'), ('structureOrder', 'int64'),('mean1','f8'),('std1','f8'),('mean2','f8'),('std2','f8'),('pvalue', 'f8'),('qvalue', 'f8'),('psign', 'int64')]; +for i in range(len(group1_roi)): dtypes.append(('count1_%d' % i, 'f8')); -for i in range(len(group2)): +for i in range(len(group2_roi)): dtypes.append(('count2_%d' % i, 'f8')); dtypes.append(('name', 'a256')); table = numpy.zeros(ids0.shape, dtype = dtypes) table["id"] = ids0; +table["structureOrder"] = labelToStructureOrder(ids0); table["mean1"] = pc10.mean(axis = 1); table["std1"] = pc10.std(axis = 1); table["mean2"] = pc20.mean(axis = 1); @@ -184,12 +159,11 @@ table["qvalue"] = qvals0; table["psign"] = psigni0; -for i in range(len(group1)): +for i in range(len(group1_roi)): table["count1_%d" % i] = pc10[:,i]; -for i in range(len(group2)): +for i in range(len(group2_roi)): table["count2_%d" % i] = pc20[:,i]; -table["name"] = lbl.labelToName(ids0); - +table["name"] = labelToName(ids0); #sort by qvalue ii = numpy.argsort(pvals0); @@ -206,49 +180,16 @@ ############################################################################# +group1i = [fn.replace('cells_transformed_to_Atlas', 'intensities') for fn in group1_roi]; +group2i = [fn.replace('cells_transformed_to_Atlas', 'intensities') for fn in group2_roi]; - -baseDirectory = '/home/mtllab/Documents/bicuculline' - - -group2 = [#'/home/mtllab/Documents/bicuculline/2/cells_heatmap.tif', - # '/home/mtllab/Documents/bicuculline/3/cells_transformed_to_Atlas.npy', - '/home/mtllab/Documents/bicuculline/7/cells_transformed_to_Atlas.npy', - '/home/mtllab/Documents/bicuculline/6/cells_transformed_to_Atlas.npy', - # '/home/mtllab/Documents/bicuculline/8/cells_transformed_to_Atlas.npy', - # '/home/mtllab/Documents/bicuculline/9/cells_transformed_to_Atlas.npy', - # '/home/mtllab/Documents/bicuculline/13/cells_transformed_to_Atlas.npy'] - '/home/mtllab/Documents/bicuculline/14/cells_transformed_to_Atlas.npy'] - -group1 = ['/home/mtllab/Documents/bicuculline/10/cells_transformed_to_Atlas.npy', - '/home/mtllab/Documents/bicuculline/11/cells_transformed_to_Atlas.npy', - '/home/mtllab/Documents/bicuculline/12/cells_transformed_to_Atlas.npy'] - - - -group1i = [fn.replace('cells_transformed_to_Atlas', 'intensities') for fn in group1]; -group2i = [fn.replace('cells_transformed_to_Atlas', 'intensities') for fn in group2]; - - - - -AnnotationFile = os.path.join(PathReg, 'annotation_25_right_fullWD.tif'); - - - -#ids, pc1, pc1i = stat.countPointsGroupInRegions(group1, intensityGroup = group1i, withIds = True, labeledImage = lbl.DefaultLabeledImageFile, withCounts = True); -#pc2, pc2i = stat.countPointsGroupInRegions(group2, intensityGroup = group2i, withIds = False, labeledImage = lbl.DefaultLabeledImageFile, withCounts = True); - -ids, pc1, pc1i = stat.countPointsGroupInRegions(group1, intensityGroup = group1i, returnIds = True, labeledImage = AnnotationFile, returnCounts = True, collapse=True); -pc2, pc2i = stat.countPointsGroupInRegions(group2, intensityGroup = group2i, returnIds = False, labeledImage = AnnotationFile, returnCounts = True, collapse=True); +ids, pc1, pc1i = stat.countPointsGroupInRegions(group1_roi, intensityGroup = group1i, returnIds = True, labeledImage = AnnotationFile, returnCounts = True, collapse=True); +pc2, pc2i = stat.countPointsGroupInRegions(group2_roi, intensityGroup = group2i, returnIds = False, labeledImage = AnnotationFile, returnCounts = True, collapse=True); pvals, psign = stat.tTestPointsInRegions(pc1, pc2, pcutoff = None, signed = True); pvalsi, psigni = stat.tTestPointsInRegions(pc1i, pc2i, pcutoff = None, signed = True, equal_var = True); -import ClearMap.Analysis.Tools.MultipleComparisonCorrection as FDR - - iid = pvalsi < 1; ids0 = ids[iid]; @@ -265,16 +206,16 @@ #make table - -dtypes = [('id','int64'),('mean1','f8'),('std1','f8'),('mean2','f8'),('std2','f8'),('pvalue', 'f8'),('qvalue', 'f8'),('psign', 'int64')]; -for i in range(len(group1)): +dtypes = [('id','int64'),('structureOrder', 'int64'),('mean1','f8'),('std1','f8'),('mean2','f8'),('std2','f8'),('pvalue', 'f8'),('qvalue', 'f8'),('psign', 'int64')]; +for i in range(len(group1_roi)): dtypes.append(('count1_%d' % i, 'f8')); -for i in range(len(group2)): +for i in range(len(group2_roi)): dtypes.append(('count2_%d' % i, 'f8')); dtypes.append(('name', 'a256')); table = numpy.zeros(ids0.shape, dtype = dtypes) table["id"] = ids0; +table["structureOrder"] = labelToStructureOrder(ids0); table["mean1"] = pc1i0.mean(axis = 1)/1000000; table["std1"] = pc1i0.std(axis = 1)/1000000; table["mean2"] = pc2i0.mean(axis = 1)/1000000; @@ -283,11 +224,11 @@ table["qvalue"] = qvalsi0; table["psign"] = psigni0; -for i in range(len(group1)): +for i in range(len(group1_roi)): table["count1_%d" % i] = pc10[:,i]; -for i in range(len(group2)): +for i in range(len(group2_roi)): table["count2_%d" % i] = pc20[:,i]; -table["name"] = lbl.labelToName(ids0); +table["name"] = labelToName(ids0); #sort by qvalue @@ -304,19 +245,18 @@ f.close(); ############################# - - #make table -dtypes = [('id','int64'),('mean1','f8'),('std1','f8'),('mean2','f8'),('std2','f8'),('pvalue', 'f8'),('qvalue', 'f8'),('psign', 'int64')]; -for i in range(len(group1)): +dtypes = [('id','int64'),('structureOrder', 'int64'),('mean1','f8'),('std1','f8'),('mean2','f8'),('std2','f8'),('pvalue', 'f8'),('qvalue', 'f8'),('psign', 'int64')]; +for i in range(len(group1_roi)): dtypes.append(('count1_%d' % i, 'f8')); -for i in range(len(group2)): +for i in range(len(group2_roi)): dtypes.append(('count2_%d' % i, 'f8')); dtypes.append(('name', 'a256')); table = numpy.zeros(ids0.shape, dtype = dtypes) table["id"] = ids0; +table["structureOrder"] = labelToStructureOrder(ids0); table["mean1"] = pc10.mean(axis = 1); table["std1"] = pc10.std(axis = 1); table["mean2"] = pc20.mean(axis = 1); @@ -325,11 +265,11 @@ table["qvalue"] = qvals0; table["psign"] = psigni0; -for i in range(len(group1)): +for i in range(len(group1_roi)): table["count1_%d" % i] = pc10[:,i]; -for i in range(len(group2)): +for i in range(len(group2_roi)): table["count2_%d" % i] = pc20[:,i]; -table["name"] = lbl.labelToName(ids0); +table["name"] = labelToName(ids0); #sort by qvalue @@ -337,10 +277,10 @@ tableSorted = table.copy(); tableSorted = tableSorted[ii]; -with open(os.path.join(baseDirectory, 'counts_table_slow.csv','w') as f: +with open(os.path.join(baseDirectory, 'counts_table_slow.csv'),'w') as f: f.write(', '.join([str(item) for item in table.dtype.names])); f.write('\n'); for sublist in tableSorted: f.write(', '.join([str(item) for item in sublist])); f.write('\n'); - f.close(); \ No newline at end of file + f.close(); diff --git a/ClearMap/Scripts/parameter_file_template.py b/ClearMap/Scripts/parameter_file_template.py index aa838a9..ad95219 100644 --- a/ClearMap/Scripts/parameter_file_template.py +++ b/ClearMap/Scripts/parameter_file_template.py @@ -18,7 +18,7 @@ from ClearMap.Analysis.Voxelization import voxelize from ClearMap.Analysis.Statistics import thresholdPoints from ClearMap.Utils.ParameterTools import joinParameter -from ClearMap.Analysis.Label import labelToName +from ClearMap.Analysis.Label import labelToName, labelToStructureOrder ######################### Data parameters diff --git a/ClearMap/Scripts/process_template.py b/ClearMap/Scripts/process_template.py index 385ea02..914a89f 100644 --- a/ClearMap/Scripts/process_template.py +++ b/ClearMap/Scripts/process_template.py @@ -91,16 +91,18 @@ ################## #With integrated weigths from the intensity file (here raw intensity): ids, counts = countPointsInRegions(points, labeledImage = AnnotationFile, intensities = intensities, intensityRow = 0); -table = numpy.zeros(ids.shape, dtype=[('id','int64'),('counts','f8'),('name', 'a256')]) +table = numpy.zeros(ids.shape, dtype=[('id','int64'), ('structureOrder', 'int64'), ('counts','f8'),('name', 'a256')]) table["id"] = ids; +table["structureOrder"] = labelToStructureOrder(ids); table["counts"] = counts; table["name"] = labelToName(ids); io.writeTable(os.path.join(BaseDirectory, 'Annotated_counts_intensities.csv'), table); #Without weigths (pure cell number): ids, counts = countPointsInRegions(points, labeledImage = AnnotationFile, intensities = None); -table = numpy.zeros(ids.shape, dtype=[('id','int64'),('counts','f8'),('name', 'a256')]) +table = numpy.zeros(ids.shape, dtype=[('id','int64'), ('structureOrder', 'int64'), ('counts','f8'),('name', 'a256')]) table["id"] = ids; +table["structureOrder"] = labelToStructureOrder(ids); table["counts"] = counts; table["name"] = labelToName(ids); io.writeTable(os.path.join(BaseDirectory, 'Annotated_counts.csv'), table);