-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathwriteClusterInformationIntoHDF5.m
executable file
·176 lines (170 loc) · 8.88 KB
/
writeClusterInformationIntoHDF5.m
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
% calculate distances to nearest neighbr and number of close neighbrs and
% write these into the hdf5 files
clear
close all
pixelsize = 100/19.5; % 100 microns are 19.5 pixels
strains = {'N2','npr1'};
wormnums = {'40'};
intensityThresholds_g = containers.Map({'40','HD','1W'},{60, 40, 100});
maxBlobSize_r = 2.5e5;
maxBlobSize_g = 1e4;
minSkelLength = 850;
maxSkelLength = 1500;
clusterCutOff = 500;
for strainCtr = 1:length(strains)
for wormnum = wormnums
%% load data
filenames_r = importdata(['datalists/' strains{strainCtr} '_' wormnum{1} '_r_list.txt']);
filenames_g = importdata(['datalists/' strains{strainCtr} '_' wormnum{1} '_g_list.txt']);
numFiles = length(filenames_g);
for fileCtr = 1:numFiles % can be parfor
filename_r = filenames_r{fileCtr};
trajData_r = h5read(filename_r,'/trajectories_data');
blobFeats_r = h5read(filename_r,'/blob_features');
skelData_r = h5read(filename_r,'/skeleton');
filename_g = filenames_g{fileCtr};
trajData_g = h5read(filename_g,'/trajectories_data');
% filter green channel data by blob size and intensity
blobFeats_g = h5read(filename_g,'/blob_features');
trajData_g.filtered = filterIntensityAndSize(blobFeats_g,pixelsize,...
intensityThresholds_g(wormnum{1}),maxBlobSize_g);
% filter red channel data by blob size and intensity
if contains(filename_r,'55')||contains(filename_r,'54')
intensityThreshold_r = 80;
else
intensityThreshold_r = 40;
end
trajData_r.filtered = filterIntensityAndSize(blobFeats_r,pixelsize,...
intensityThreshold_r,maxBlobSize_r);
% filter by skeleton length
trajData_r.filtered = trajData_r.filtered&logical(trajData_r.is_good_skel)&...
filterSkelLength(skelData_r,pixelsize,minSkelLength,maxSkelLength);
%% calculate stats - red channel files
try
neighbr_distances_skeleton = h5read(filename_r,'/neighbr_distances_skeleton');
catch
disp(['Calculating skeleton distances from ' filename_r ])
neighbr_distances_skeleton ...
= calculateNeighbrDistanceFromSkeleton(skelData_r,trajData_r,trajData_g,pixelsize);
assert(size(neighbr_distances_skeleton,1)==length(trajData_r.frame_number))
assert(size(neighbr_distances_skeleton,2)==size(skelData_r,2));
assert(size(neighbr_distances_skeleton,3)==10)
% write stats to hdf5-file
h5create(filename_r,'/neighbr_distances_skeleton',size(neighbr_distances_skeleton),...
'Datatype','single')
h5write(filename_r,'/neighbr_distances_skeleton',single(neighbr_distances_skeleton))
end
try
neighbr_distances = h5read(filename_r,'/neighbr_distances');
min_neighbr_dist = h5read(filename_r,'/min_neighbr_dist');
num_close_neighbrs = h5read(filename_r,'/num_close_neighbrs');
catch
disp(['Calculating cluster status from ' filename_r ])
[ neighbr_distances, min_neighbr_dist, num_close_neighbrs ] ...
= calculateNeighbrDistance(trajData_r,trajData_g,pixelsize,clusterCutOff);
% check lengths
assert(size(neighbr_distances,1)==length(trajData_r.frame_number))
assert(size(neighbr_distances,2)==10)
assert(length(min_neighbr_dist)==length(trajData_r.frame_number))
assert(length(num_close_neighbrs)==length(trajData_r.frame_number))
% write stats to hdf5-file
try
h5create(filename_r,'/neighbr_distances',size(neighbr_distances),...
'Datatype','single')
h5write(filename_r,'/neighbr_distances',single(neighbr_distances))
catch
disp('/neighbr_distances already exists')
end
try
h5create(filename_r,'/min_neighbr_dist',...
size(min_neighbr_dist),'Datatype','single')
h5write(filename_r,'/min_neighbr_dist',...
single(min_neighbr_dist))
catch
disp('/min_neighbr_dist already exists')
end
try
h5create(filename_r,'/num_close_neighbrs',...
size(num_close_neighbrs),'Datatype','uint16')
h5write(filename_r,'/num_close_neighbrs',...
uint16(num_close_neighbrs))
catch
disp('/num_close_neighbrs already exists')
end
end
%% calculate stats - green channel files
try
neighbr_distances = h5read(filename_g,'/neighbr_distances');
min_neighbr_dist = h5read(filename_g,'/min_neighbr_dist');
num_close_neighbrs = h5read(filename_g,'/num_close_neighbrs');
catch
disp(['Calculating cluster status from ' filename_g ])
[ neighbr_distances, min_neighbr_dist, num_close_neighbrs ] ...
= calculateNeighbrDistance(trajData_g,trajData_r,pixelsize,clusterCutOff);
% check lengths
assert(size(neighbr_distances,1)==length(trajData_g.frame_number))
assert(size(neighbr_distances,2)==10)
assert(length(min_neighbr_dist)==length(trajData_g.frame_number))
assert(length(num_close_neighbrs)==length(trajData_g.frame_number))
% write stats to hdf5-file
h5create(filename_g,'/neighbr_distances',size(neighbr_distances),...
'Datatype','single')
h5write(filename_g,'/neighbr_distances',single(neighbr_distances))
h5create(filename_g,'/min_neighbr_dist',...
size(min_neighbr_dist),'Datatype','single')
h5write(filename_g,'/min_neighbr_dist',...
single(min_neighbr_dist))
h5create(filename_g,'/num_close_neighbrs',...
size(num_close_neighbrs),'Datatype','uint16')
h5write(filename_g,'/num_close_neighbrs',...
uint16(num_close_neighbrs))
end
end
end
end
%% first data set - green channel files only, one more strain, and slightly different filters
strains = {'npr1','HA','N2'};
wormnums = {'40','HD'};
intensityThresholds_g('40') = 50;
for strainCtr = 1:length(strains)
for wormnum = wormnums
%% load data
filenames_g = importdata(['datalists/' strains{strainCtr} '_' wormnum{1} '_list.txt']);
numFiles = length(filenames_g);
for fileCtr = 1:numFiles % can be parfor
filename_g = filenames_g{fileCtr};
trajData_g = h5read(filename_g,'/trajectories_data');
% filter green channel data by blob size and intensity
blobFeats_g = h5read(filename_g,'/blob_features');
trajData_g.filtered = filterIntensityAndSize(blobFeats_g,pixelsize,...
intensityThresholds_g(wormnum{1}),maxBlobSize_g);
%% calculate stats - green channel files
try
neighbr_distances = h5read(filename_g,'/neighbr_distances');
min_neighbr_dist = h5read(filename_g,'/min_neighbr_dist');
num_close_neighbrs = h5read(filename_g,'/num_close_neighbrs');
catch
disp(['Calculating cluster status from ' filename_g ])
[ neighbr_distances, min_neighbr_dist, num_close_neighbrs ] ...
= calculateNeighbrDistance(trajData_g,[],pixelsize,clusterCutOff);
% check lengths
assert(size(neighbr_distances,1)==length(trajData_g.frame_number))
assert(size(neighbr_distances,2)==10)
assert(length(min_neighbr_dist)==length(trajData_g.frame_number))
assert(length(num_close_neighbrs)==length(trajData_g.frame_number))
% write stats to hdf5-file
h5create(filename_g,'/neighbr_distances',size(neighbr_distances),...
'Datatype','single')
h5write(filename_g,'/neighbr_distances',single(neighbr_distances))
h5create(filename_g,'/min_neighbr_dist',...
size(min_neighbr_dist),'Datatype','single')
h5write(filename_g,'/min_neighbr_dist',...
single(min_neighbr_dist))
h5create(filename_g,'/num_close_neighbrs',...
size(num_close_neighbrs),'Datatype','uint16')
h5write(filename_g,'/num_close_neighbrs',...
uint16(num_close_neighbrs))
end
end
end
end