Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fixes #9 #10

Open
wants to merge 5 commits into
base: master
Choose a base branch
from
Open

Fixes #9 #10

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
71d17ff
Add StructureOrder to the tables
kannanuv Dec 9, 2016
ddd71ab
Merge remote-tracking branch 'refs/remotes/ChristophKirst/master' int…
kannanuv Dec 13, 2016
6cd85fb
Added structure order to the statistical analysis tables.
kannanuv Dec 13, 2016
40f86cf
Merge pull request #1 from kannanuv/structureOrderAdditionInCellCounts
kannanuv Dec 13, 2016
7ffaeca
Merge remote-tracking branch 'refs/remotes/ChristophKirst/master'
kannanuv Aug 25, 2017
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
13 changes: 11 additions & 2 deletions ClearMap/Analysis/Label.py
View file
Open in desktop
Expand Up @@ -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"""


Expand All @@ -77,6 +77,7 @@ class LabelInfo(object):
names = None;
acronyms = None;
colors = None;
structureOrder = None;
parents = None;
levels = None;
collapse = None;
Expand All @@ -92,14 +93,15 @@ 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();

slf.ids = [x.id for x in labels];
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};

Expand All @@ -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];

Expand Down Expand Up @@ -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];



Expand Down
184 changes: 62 additions & 122 deletions ClearMap/Scripts/analyze_template.py
View file
Open in desktop
@@ -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
##############################################
Expand All @@ -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
Expand All @@ -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];
Expand All @@ -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;
Expand All @@ -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
Expand All @@ -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);
Expand All @@ -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);
Expand All @@ -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];
Expand All @@ -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;
Expand All @@ -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
Expand All @@ -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);
Expand All @@ -325,22 +265,22 @@
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);
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();
f.close();