1515
1616
1717# Base Directory, usually where your experiment is saved:
18- baseDirectory = '/d2/studies/ClearMap/IA_iDISCO/analysisFiles '
19-
18+ baseDirectory = '/d2/studies/ClearMap/ROC_iDISCO/batchDirectory '
19+ sampleName = 'ROC_iDISCO'
2020
2121
2222# Voxel-based statistics:
2323#########################
2424
2525#Load the data (heat maps generated previously )
26- group1 = [os .path .join (baseDirectory , 'cells_heatmap_vox15_IA1_LB.tif' ),
27- os .path .join (baseDirectory , 'cells_heatmap_vox15_IA2_NP.tif' ),
28- os .path .join (baseDirectory , 'cells_heatmap_vox15_IA2_LT.tif' ),
29- os .path .join (baseDirectory , 'cells_heatmap_vox15_IA2_LB.tif' )
26+ group1 = [os .path .join (baseDirectory , 'cells_heatmap_vox15_ROC9.tif' ),
27+ os .path .join (baseDirectory , 'cells_heatmap_vox15_ROC10.tif' ),
28+ os .path .join (baseDirectory , 'cells_heatmap_vox15_ROC11.tif' ),
29+ os .path .join (baseDirectory , 'cells_heatmap_vox15_ROC12.tif' ),
30+ os .path .join (baseDirectory , 'cells_heatmap_vox15_ROC13.tif' ),
31+ os .path .join (baseDirectory , 'cells_heatmap_vox15_ROC14.tif' ),
32+ os .path .join (baseDirectory , 'cells_heatmap_vox15_ROC15.tif' ),
33+ os .path .join (baseDirectory , 'cells_heatmap_vox15_ROC16.tif' )
3034 ]
3135
32- group2 = [os .path .join (baseDirectory , 'cells_heatmap_vox15_IA1_RT.tif' ),
33- os .path .join (baseDirectory , 'cells_heatmap_vox15_IA1_RB.tif' ),
34- os .path .join (baseDirectory , 'cells_heatmap_vox15_IA1_LT.tif' ),
35- os .path .join (baseDirectory , 'cells_heatmap_vox15_IA2_RT.tif' ),
36- os .path .join (baseDirectory , 'cells_heatmap_vox15_IA2_RB.tif' )
36+ group2 = [os .path .join (baseDirectory , 'cells_heatmap_vox15_ROC1.tif' ),
37+ os .path .join (baseDirectory , 'cells_heatmap_vox15_ROC2.tif' ),
38+ os .path .join (baseDirectory , 'cells_heatmap_vox15_ROC4.tif' ),
39+ os .path .join (baseDirectory , 'cells_heatmap_vox15_ROC6.tif' ),
40+ os .path .join (baseDirectory , 'cells_heatmap_vox15_ROC7.tif' ),
41+ os .path .join (baseDirectory , 'cells_heatmap_vox15_ROC8.tif' ),
42+ os .path .join (baseDirectory , 'cells_heatmap_vox15_ROC17.tif' ),
43+ os .path .join (baseDirectory , 'cells_heatmap_vox15_ROC18.tif' ),
44+
3745 ]
3846
3947
5058g2a = numpy .mean (g2 ,axis = 0 );
5159g2s = numpy .std (g2 ,axis = 0 );
5260
53- io .writeData (os .path .join (baseDirectory , 'Yoked_full_mean_horizontal_vox15.mhd' ), g1a )# rsp.sagittalToCoronalData(g1a));
54- io .writeData (os .path .join (baseDirectory , 'Yoked_full_std_horizontal_vox15.mhd' ), g1s )# rsp.sagittalToCoronalData(g1s));
55-
56- io .writeData (os .path .join (baseDirectory , 'IA30_full_mean_horizontal_vox15.mhd' ), g2a )# rsp.sagittalToCoronalData(g2a));
57- io .writeData (os .path .join (baseDirectory , 'IA30_full_std_horizontal_vox15.mhd' ), g2s )# rsp.sagittalToCoronalData(g2s));
58-
59-
61+ io .writeData (os .path .join (baseDirectory , 'Group1_full_mean_horizontal.mhd' ), g1a )# rsp.sagittalToCoronalData(g1a));
62+ io .writeData (os .path .join (baseDirectory , 'Group1_full_std_horizontal.mhd' ), g1s )# rsp.sagittalToCoronalData(g1s));
6063
64+ io .writeData (os .path .join (baseDirectory , 'Group2_full_mean_horizontal.mhd' ), g2a )# rsp.sagittalToCoronalData(g2a));
65+ io .writeData (os .path .join (baseDirectory , 'Group2_full_std_horizontal.mhd' ), g2s )# rsp.sagittalToCoronalData(g2s));
6166
6267
6368#Generate the p-values map
6671pvals , psign = stat .tTestVoxelization (g1 .astype ('float' ), g2 .astype ('float' ), signed = True , pcutoff = 0.05 );
6772
6873#color the p-values according to their sign (defined by the sign of the difference of the means between the 2 groups)
69- pvalsc = stat .colorPValues (pvals , psign , positive = [0 , 1 ], negative = [1 , 0 ]);
70- io .writeData (os .path .join (baseDirectory , 'pvalues_full_vox15 .tif' ), pvalsc .astype ('float32' ))#, rsp.sagittalToCoronalData(pvalsc.astype('float32')));
74+ pvalsc = stat .colorPValues (pvals , psign , positive = [1 , 0 ], negative = [0 , 1 ]);
75+ io .writeData (os .path .join (baseDirectory , sampleName + '_pvalues .tif' ), pvalsc .astype ('float32' ))#, rsp.sagittalToCoronalData(pvalsc.astype('float32')));
7176
7277
7378
7479#############################################################################
75-
7680# Regions-based statistics:
7781###########################
7882
100104PathReg = '/d1/studies/ClearMap/ClearMap_resources/Parameter_files/' ;
101105AnnotationFile = '/d2/studies/ClearMap/IA_iDISCO/Annotations_Horizontal_8-302_Full.tif' ;
102106
103-
104-
105- #ids, pc1, pc1i = stat.countPointsGroupInRegions(group1, intensityGroup = group1i, withIds = True, labeledImage = lbl.DefaultLabeledImageFile, withCounts = True);
106- #pc2, pc2i = stat.countPointsGroupInRegions(group2, intensityGroup = group2i, withIds = False, labeledImage = lbl.DefaultLabeledImageFile, withCounts = True);
107-
107+ #Calculate stats for regions. Set collapse=True to use collapsed regions (e.g. collapse cortical layers into single region)
108108ids , pc1 , pc1i = stat .countPointsGroupInRegions (group1 , intensityGroup = group1i , returnIds = True , labeledImage = AnnotationFile , returnCounts = True , collapse = None );
109109pc2 , pc2i = stat .countPointsGroupInRegions (group2 , intensityGroup = group2i , returnIds = False , labeledImage = AnnotationFile , returnCounts = True , collapse = None );
110110
156156table ["name" ] = lbl .labelToName (ids0 );
157157
158158
159- #sort by qvalue
159+ #sort by qvalue & write intensity table
160160ii = numpy .argsort (pvalsi0 );
161161tableSorted = table .copy ();
162162tableSorted = tableSorted [ii ];
163163
164- with open (os .path .join (baseDirectory , 'counts-intensity_table_6 .csv' ),'w' ) as f :
164+ with open (os .path .join (baseDirectory , 'counts-intensity_table .csv' ),'w' ) as f :
165165 f .write (', ' .join ([str (item ) for item in table .dtype .names ]));
166166 f .write ('\n ' );
167167 for sublist in tableSorted :
198198table ["name" ] = lbl .labelToName (ids0 );
199199
200200
201- #sort by qvalue
201+ #sort by qvalue & write counts table
202202ii = numpy .argsort (pvals0 );
203203tableSorted = table .copy ();
204204tableSorted = tableSorted [ii ];
205205
206- with open (os .path .join (baseDirectory , 'counts_table_6 .csv' ),'w' ) as f :
206+ with open (os .path .join (baseDirectory , 'counts_table .csv' ),'w' ) as f :
207207 f .write (', ' .join ([str (item ) for item in table .dtype .names ]));
208208 f .write ('\n ' );
209209 for sublist in tableSorted :
210210 f .write (', ' .join ([str (item ) for item in sublist ]));
211211 f .write ('\n ' );
212212 f .close ();
213213
214-
215- #############################################################################
216-
217-
218- baseDirectory = '/home/mtllab/Documents/bicuculline'
219-
220-
221- group2 = [#'/home/mtllab/Documents/bicuculline/2/cells_heatmap.tif',
222- # '/home/mtllab/Documents/bicuculline/3/cells_transformed_to_Atlas.npy',
223- '/home/mtllab/Documents/bicuculline/7/cells_transformed_to_Atlas.npy' ,
224- '/home/mtllab/Documents/bicuculline/6/cells_transformed_to_Atlas.npy' ,
225- # '/home/mtllab/Documents/bicuculline/8/cells_transformed_to_Atlas.npy',
226- # '/home/mtllab/Documents/bicuculline/9/cells_transformed_to_Atlas.npy',
227- # '/home/mtllab/Documents/bicuculline/13/cells_transformed_to_Atlas.npy']
228- '/home/mtllab/Documents/bicuculline/14/cells_transformed_to_Atlas.npy' ]
229-
230- group1 = ['/home/mtllab/Documents/bicuculline/10/cells_transformed_to_Atlas.npy' ,
231- '/home/mtllab/Documents/bicuculline/11/cells_transformed_to_Atlas.npy' ,
232- '/home/mtllab/Documents/bicuculline/12/cells_transformed_to_Atlas.npy' ]
233-
234-
235-
236- group1i = [fn .replace ('cells_transformed_to_Atlas' , 'intensities' ) for fn in group1 ];
237- group2i = [fn .replace ('cells_transformed_to_Atlas' , 'intensities' ) for fn in group2 ];
238-
239-
240-
241-
242- AnnotationFile = os .path .join (PathReg , 'annotation_25_right_fullWD.tif' );
243-
244-
245-
246- #ids, pc1, pc1i = stat.countPointsGroupInRegions(group1, intensityGroup = group1i, withIds = True, labeledImage = lbl.DefaultLabeledImageFile, withCounts = True);
247- #pc2, pc2i = stat.countPointsGroupInRegions(group2, intensityGroup = group2i, withIds = False, labeledImage = lbl.DefaultLabeledImageFile, withCounts = True);
248-
249- ids , pc1 , pc1i = stat .countPointsGroupInRegions (group1 , intensityGroup = group1i , returnIds = True , labeledImage = AnnotationFile , returnCounts = True , collapse = True );
250- pc2 , pc2i = stat .countPointsGroupInRegions (group2 , intensityGroup = group2i , returnIds = False , labeledImage = AnnotationFile , returnCounts = True , collapse = True );
251-
252-
253- pvals , psign = stat .tTestPointsInRegions (pc1 , pc2 , pcutoff = None , signed = True );
254- pvalsi , psigni = stat .tTestPointsInRegions (pc1i , pc2i , pcutoff = None , signed = True , equal_var = True );
255-
256- import ClearMap .Analysis .Tools .MultipleComparisonCorrection as FDR
257-
258-
259- iid = pvalsi < 1 ;
260-
261- ids0 = ids [iid ];
262- pc1i0 = pc1i [iid ];
263- pc2i0 = pc2i [iid ];
264- pc10 = pc1 [iid ];
265- pc20 = pc2 [iid ];
266- psigni0 = psigni [iid ];
267- pvalsi0 = pvalsi [iid ];
268- qvalsi0 = FDR .estimateQValues (pvalsi0 );
269- psign0 = psign [iid ];
270- pvals0 = pvals [iid ];
271- qvals0 = FDR .estimateQValues (pvals0 );
272-
273-
274- #make table
275-
276- dtypes = [('id' ,'int64' ),('mean1' ,'f8' ),('std1' ,'f8' ),('mean2' ,'f8' ),('std2' ,'f8' ),('pvalue' , 'f8' ),('qvalue' , 'f8' ),('psign' , 'int64' )];
277- for i in range (len (group1 )):
278- dtypes .append (('count1_%d' % i , 'f8' ));
279- for i in range (len (group2 )):
280- dtypes .append (('count2_%d' % i , 'f8' ));
281- dtypes .append (('name' , 'a256' ));
282-
283- table = numpy .zeros (ids0 .shape , dtype = dtypes )
284- table ["id" ] = ids0 ;
285- table ["mean1" ] = pc1i0 .mean (axis = 1 )/ 1000000 ;
286- table ["std1" ] = pc1i0 .std (axis = 1 )/ 1000000 ;
287- table ["mean2" ] = pc2i0 .mean (axis = 1 )/ 1000000 ;
288- table ["std2" ] = pc2i0 .std (axis = 1 )/ 1000000 ;
289- table ["pvalue" ] = pvalsi0 ;
290- table ["qvalue" ] = qvalsi0 ;
291-
292- table ["psign" ] = psigni0 ;
293- for i in range (len (group1 )):
294- table ["count1_%d" % i ] = pc10 [:,i ];
295- for i in range (len (group2 )):
296- table ["count2_%d" % i ] = pc20 [:,i ];
297- table ["name" ] = lbl .labelToName (ids0 );
298-
299-
300- #sort by qvalue
301- ii = numpy .argsort (pvalsi0 );
302- tableSorted = table .copy ();
303- tableSorted = tableSorted [ii ];
304-
305- with open (os .path .join (baseDirectory , 'counts-intensity_table_slow.csv' ),'w' ) as f :
306- f .write (', ' .join ([str (item ) for item in table .dtype .names ]));
307- f .write ('\n ' );
308- for sublist in tableSorted :
309- f .write (', ' .join ([str (item ) for item in sublist ]));
310- f .write ('\n ' );
311- f .close ();
312-
313- #############################
314-
315-
316- #make table
317-
318- dtypes = [('id' ,'int64' ),('mean1' ,'f8' ),('std1' ,'f8' ),('mean2' ,'f8' ),('std2' ,'f8' ),('pvalue' , 'f8' ),('qvalue' , 'f8' ),('psign' , 'int64' )];
319- for i in range (len (group1 )):
320- dtypes .append (('count1_%d' % i , 'f8' ));
321- for i in range (len (group2 )):
322- dtypes .append (('count2_%d' % i , 'f8' ));
323- dtypes .append (('name' , 'a256' ));
324-
325- table = numpy .zeros (ids0 .shape , dtype = dtypes )
326- table ["id" ] = ids0 ;
327- table ["mean1" ] = pc10 .mean (axis = 1 );
328- table ["std1" ] = pc10 .std (axis = 1 );
329- table ["mean2" ] = pc20 .mean (axis = 1 );
330- table ["std2" ] = pc20 .std (axis = 1 );
331- table ["pvalue" ] = pvals0 ;
332- table ["qvalue" ] = qvals0 ;
333-
334- table ["psign" ] = psigni0 ;
335- for i in range (len (group1 )):
336- table ["count1_%d" % i ] = pc10 [:,i ];
337- for i in range (len (group2 )):
338- table ["count2_%d" % i ] = pc20 [:,i ];
339- table ["name" ] = lbl .labelToName (ids0 );
340-
341-
342- #sort by qvalue
343- ii = numpy .argsort (pvals0 );
344- tableSorted = table .copy ();
345- tableSorted = tableSorted [ii ];
346-
347- with open (os .path .join (baseDirectory , 'counts_table_slow.csv' ,'w' ) as f :
348- f .write (', ' .join ([str (item ) for item in table .dtype .names ]));
349- f .write ('\n ' );
350- for sublist in tableSorted :
351- f .write (', ' .join ([str (item ) for item in sublist ]));
352- f .write ('\n ' );
353- f .close ();
0 commit comments