/
m3ha_compute_statistics.m
executable file
·470 lines (415 loc) · 20.1 KB
/
m3ha_compute_statistics.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
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
function statsTable = m3ha_compute_statistics (varargin)
%% Computes LTS and burst statistics for some indices (indOfInterest) in ltsOnsetTimeEachSwp, spikesPerLtsEachSwp, burstOnsetTimeEachSwp & spikesPerBurst
% Usage: statsTable = m3ha_compute_statistics (varargin)
%
% Explanation:
% TODO
%
% Example(s):
% statsTable = m3ha_compute_statistics
% statsTable = m3ha_compute_statistics('DataMode', 2)
% statsTable = m3ha_compute_statistics('PharmConditions', 1)
% statsTable = m3ha_compute_statistics('PharmConditions', 1:4)
% statsTable = m3ha_compute_statistics('PharmConditions', num2cell(1:4))
% statsTable = m3ha_compute_statistics('PharmConditions', num2cell(1:4), 'GIncrCondition', num2cell([100; 200; 400]))
%
% Outputs:
% statsTable - a table containing measures as row names
% and the following variables:
% measureTitle
% measureStr
% dataMode
% pharmCondition
% gIncrCondition
% vHoldCondition
% allValues
% nValues
% meanValue
% stdValue
% stderrValue
% errValue
% upper95Value
% lower95Value
% specified as a table
%
% Arguments:
% varargin - 'SwpInfo': a table of sweep info, with each row named by
% the matfile base containing the raw data
% must a 2D table with row names being file bases
% and with the fields:
% cellidrow - cell ID
% ltspeaktimes - LTS peak delays
% spikesperpeak - action potentials per LTS
% bursttimes - burst delays
% spikesperburst - action potentials per burst
% default == m3ha_load_sweep_info
% - 'DataMode': data mode
% must be one of:
% 0 - all data
% 1 - all of g incr = 100%, 200%, 400%
% 2 - all of g incr = 100%, 200%, 400%
% but exclude cell-pharm-g_incr sets
% containing problematic sweeps
% default == 0
% - 'PharmConditions': pharmacological condition(s)
% to restrict to
% must be empty or some of:
% 1 - control
% 2 - GAT1 blockade
% 3 - GAT3 blockade
% 4 - dual blockade
% or a cell array of them (will become 1st dimension)
% default == no restrictions
% - 'GIncrCondition': conductance amplitude condition(s) (%)
% to restrict to
% must be empty or some of: 25, 50, 100, 200, 400, 800
% or a cell array of them (will become 2nd dimension)
% default == no restrictions
% - 'VHoldConditions': holding potential condition(s) (mV)
% to restrict to
% must be empty or some of: -60, -65, -70
% or a cell array of them (will become 3rd dimension)
% default == no restrictions
%
% Requires:
% cd/argfun.m
% cd/array_fun.m
% cd/compute_stats.m
% cd/first_matching_field.m
% cd/force_column_cell.m
% cd/force_column_vector.m
% cd/match_row_count.m
% cd/m3ha_load_sweep_info.m
% cd/m3ha_select_sweeps.m
%
% Used by:
% cd/m3ha_compute_and_plot_statistics.m
% cd/m3ha_compute_and_plot_violin.m
% File History:
% 2016-08-19 Created
% 2016-08-29 Last Modified
% 2017-01-25 - Corrected errors to reflect t-confidence intervals
% (from the Gosset's t distribution)
% 2019-11-26 Improved code structure
% 2019-11-27 Now computes means of all measures from the data points
% for each cell (rather than for each sweep)
% 2019-11-27 Added 'PharmConditions', 'GIncrConditions' & 'VHoldConditions'
% as optional arguments
% 2019-12-04 Added other LTS features that might be of interest:
% 2019-12-04 Added other burst features that might be of interest
% 2020-02-16 Now makes sure the same cells are represented in all groups
% 2020-05-06 Now removes LTS amplitude values that are less than -100 mV
% TODO: Add 'MeasuresToCompute' as an optional argument
%% Hard-coded parameters
% Items to compute
measureTitle = {'LTS amplitude (mV)'; 'LTS maximum slope (V/s)'; ...
'LTS concavity (V^2/s^2)'; 'LTS prominence (mv)'; ...
'LTS width (ms)'; 'LTS onset time (ms)'; ...
'LTS time jitter (ms)'; ...
'LTS probability'; 'Spikes per LTS'; ...
'Maximum spike amplitude (mV)'; ...
'Minimum spike amplitude (mV)'; ...
'Spike frequency (Hz)'; 'Spike adaptation'; ...
'Burst onset time (ms)'; 'Burst time jitter (ms)'; ...
'Burst probability'; 'Spikes per burst'};
% Note: measureStr must match the variables defined in this file!
measureStr = {'ltsAmplitude'; 'ltsMaxSlope'; ...
'ltsConcavity'; 'ltsProminence'; ...
'ltsWidth'; 'ltsOnsetTime'; 'ltsTimeJitter'; ...
'ltsProbability'; 'spikesPerLts'; ...
'spikeMaxAmp'; 'spikeMinAmp'; ...
'spikeFrequency'; 'spikeAdaptation'
'burstOnsetTime'; 'burstTimeJitter'; ...
'burstProbability'; 'spikesPerBurst'};
% Columns defined in swpInfo
cellIdStr = 'cellidrow';
ltsAmplitudeStr = {'ltspeakval', 'ltsPeakValue'};
ltsMaxSlopeStr = {'maxslopeval', 'maxSlopeValue'};
ltsConcavityStr = {'peak2ndder', 'peak2ndDer'};
ltsProminenceStr = {'peakprom', 'peakProm'};
ltsWidthStr = {'peakwidth', 'peakWidth'};
ltsOnsetTimeStr = {'ltspeaktime', 'ltsPeakTime'};
spikesPerLtsStr = {'spikesperpeak', 'spikesPerPeak'};
spikeMaxAmpStr = {'maxspikeamp', 'maxSpikeAmp'};
spikeMinAmpStr = {'minspikeamp', 'minSpikeAmp'};
spikeFrequencyStr = {'spikefrequency', 'spikeFrequency'};
spikeAdaptationStr = {'spikeadaptation', 'spikeAdaptation'};
burstOnsetTimeStr = {'bursttime', 'burstTime'};
spikesPerBurstStr = {'spikesperburst', 'spikesPerBurst'};
%% Default values for optional arguments
swpInfoDefault = table.empty;
dataModeDefault = 0;
pharmConditionsDefault = [];
gIncrConditionsDefault = [];
vHoldConditionsDefault = [];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Deal with arguments
% Set up Input Parser Scheme
iP = inputParser;
iP.FunctionName = mfilename;
% Add parameter-value pairs to the Input Parser
addParameter(iP, 'SwpInfo', swpInfoDefault, ...
@(x) validateattributes(x, {'table'}, {'2d'}));
addParameter(iP, 'DataMode', dataModeDefault, ...
@(x) validateattributes(x, {'numeric'}, {'integer', 'scalar'}));
addParameter(iP, 'PharmConditions', pharmConditionsDefault, ...
@(x) validateattributes(x, {'numeric', 'cell'}, {'2d'}));
addParameter(iP, 'GIncrConditions', gIncrConditionsDefault, ...
@(x) validateattributes(x, {'numeric', 'cell'}, {'2d'}));
addParameter(iP, 'VHoldConditions', vHoldConditionsDefault, ...
@(x) validateattributes(x, {'numeric', 'cell'}, {'2d'}));
% Read from the Input Parser
parse(iP, varargin{:});
swpInfo = iP.Results.SwpInfo;
dataMode = iP.Results.DataMode;
pharmConditionsUser = iP.Results.PharmConditions;
gIncrConditionsUser = iP.Results.GIncrConditions;
vHoldConditionsUser = iP.Results.VHoldConditions;
% Count the items to compute
nMeasures = numel(measureStr);
%% Preparation
% Read the sweep info data
if isempty(swpInfo)
swpInfo = m3ha_load_sweep_info;
end
% Update the toUse column based on data mode
swpInfo = m3ha_select_sweeps('SwpInfo', swpInfo, 'Verbose', false, ...
'DataMode', dataMode);
%% Determine all possible cells
% Extract whether to use the sweep
toUse = swpInfo.toUse;
% Restrict swpInfo to those sweeps
swpInfoToUse = swpInfo(toUse, :);
% Extract the cell IDs
cellIdRow = swpInfoToUse.(cellIdStr);
% Find unique cell IDs
uniqueCellIds = unique(cellIdRow);
%% Create all conditions
if iscell(pharmConditionsUser) || iscell(gIncrConditionsUser) || ...
iscell(vHoldConditionsUser)
% Force as cell arrays of column vectors
[pharmConditionsUser, gIncrConditionsUser, vHoldConditionsUser] = ...
argfun(@(x) force_column_vector(x, 'ToLinearize', true, ...
'ForceCellOutput', true), ...
pharmConditionsUser, gIncrConditionsUser, vHoldConditionsUser);
% Force as column cell arrays
[pharmConditionsUser, gIncrConditionsUser, vHoldConditionsUser] = ...
argfun(@force_column_cell, ...
pharmConditionsUser, gIncrConditionsUser, vHoldConditionsUser);
% Count the number of values for each type of condition
nPharm = numel(pharmConditionsUser);
nGIncr = numel(gIncrConditionsUser);
nVHold = numel(vHoldConditionsUser);
% Mesh indices for all set of conditions
% Note: meshgrid returns a length(y) by length(x) by length(z) matrix
% so one needs to permute the first two dimensions
% to give a length(x) by length(y) by length(z) matrix
% TODO: meshgrid_custom.m
[indPharmCondition, indGIncrCondition, indVHoldCondition] = ...
meshgrid(1:nPharm, 1:nGIncr, 1:nVHold);
[indPharmCondition, indGIncrCondition, indVHoldCondition] = ...
argfun(@(x) permute(x, [2, 1]), ...
indPharmCondition, indGIncrCondition, indVHoldCondition);
% Mesh all set of conditions
pharmCondition = array_fun(@(x) pharmConditionsUser{x}, ...
indPharmCondition, 'UniformOutput', false);
gIncrCondition = array_fun(@(x) gIncrConditionsUser{x}, ...
indGIncrCondition, 'UniformOutput', false);
vHoldCondition = array_fun(@(x) vHoldConditionsUser{x}, ...
indVHoldCondition, 'UniformOutput', false);
% Set flag
manyConditionsFlag = true;
else
% There is only one set of conditions
pharmCondition = pharmConditionsUser;
gIncrCondition = gIncrConditionsUser;
vHoldCondition = vHoldConditionsUser;
% Set flag
manyConditionsFlag = false;
end
%% Compute statistics
if manyConditionsFlag
% Compute statistics for each set of conditions
% Note: This will return a nPharm x nGIncr x nVhold cell array
% where each element is either a nMeasures x 1 column cell vector
% or a nMeasures x 1 column numeric vector
[allValues, nValues, meanValue, stdValue, ...
stderrValue, errValue, upper95Value, lower95Value] = ...
array_fun(@(x, y, z) ...
m3ha_compute_statistics_helper(swpInfo, ...
x, y, z, uniqueCellIds, cellIdStr, measureStr, ...
ltsAmplitudeStr, ltsMaxSlopeStr, ltsConcavityStr, ...
ltsProminenceStr, ltsWidthStr, ltsOnsetTimeStr, ...
spikesPerLtsStr, spikeMaxAmpStr, spikeMinAmpStr, ...
spikeFrequencyStr, spikeAdaptationStr, ...
burstOnsetTimeStr, spikesPerBurstStr), ...
pharmCondition, gIncrCondition, vHoldCondition, ...
'UniformOutput', false);
% Reorganize cell arrays of cell arrays
% so that measures are grouped together
% Note: This will return a nMeasures x 1 cell array,
% where each element is a nPharm x nGIncr x nVhold cell array
allValues = extract_columns(allValues, transpose(1:nMeasures), ...
'TreatCnvAsColumns', true, 'OutputMode', 'single');
% Reorganize cell arrays of numeric vectors
% so that measures are grouped together
% Note: This will return a nMeasures x 1 cell array,
% where each element is a nPharm x nGIncr x nVhold numeric array
[nValues, meanValue, stdValue, ...
stderrValue, errValue, upper95Value, lower95Value] = ...
argfun(@(x) ...
array_fun(@(y) extract_elements(x, 'specific', 'Index', y), ...
transpose(1:nMeasures), 'UniformOutput', false), ...
nValues, meanValue, stdValue, ...
stderrValue, errValue, upper95Value, lower95Value);
else
% Compute statistics for this set of conditions
% Note: this will return either a nMeasures x 1 column cell vector
% or a nMeasures x 1 column numeric vector
[allValues, nValues, meanValue, stdValue, ...
stderrValue, errValue, upper95Value, lower95Value] = ...
m3ha_compute_statistics_helper(swpInfo, ...
pharmCondition, gIncrCondition, vHoldCondition, ...
uniqueCellIds, cellIdStr, measureStr, ...
ltsAmplitudeStr, ltsMaxSlopeStr, ltsConcavityStr, ...
ltsProminenceStr, ltsWidthStr, ltsOnsetTimeStr, ...
spikesPerLtsStr, spikeMaxAmpStr, spikeMinAmpStr, ...
spikeFrequencyStr, spikeAdaptationStr, ...
burstOnsetTimeStr, spikesPerBurstStr);
end
%% Output results
% Put in a cell array
[pharmCondition, gIncrCondition, vHoldCondition] = ...
argfun(@(x) {x}, pharmCondition, gIncrCondition, vHoldCondition);
% Match the number of rows
[dataMode, pharmCondition, gIncrCondition, vHoldCondition] = ...
argfun(@(x) match_row_count(x, nMeasures), ...
dataMode, pharmCondition, gIncrCondition, vHoldCondition);
% Create a statistics table
statsTable = table(measureTitle, measureStr, dataMode, ...
pharmCondition, gIncrCondition, vHoldCondition, ...
allValues, nValues, meanValue, stdValue, ...
stderrValue, errValue, upper95Value, lower95Value, ...
'RowNames', measureStr);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [allValues, nValues, meanValue, stdValue, ...
stderrValue, errValue, upper95Value, lower95Value] = ...
m3ha_compute_statistics_helper(swpInfo, ...
pharmConditions, gIncrConditions, vHoldConditions, ...
uniqueCellIds, cellIdStr, measureStr, ...
ltsAmplitudeStr, ltsMaxSlopeStr, ltsConcavityStr, ...
ltsProminenceStr, ltsWidthStr, ltsOnsetTimeStr, ...
spikesPerLtsStr, spikeMaxAmpStr, spikeMinAmpStr, ...
spikeFrequencyStr, spikeAdaptationStr, ...
burstOnsetTimeStr, spikesPerBurstStr)
%% Computes LTS and burst statistics for one condition
%% Select sweeps
% Restrict sweeps to use to the given conditions
swpInfo = m3ha_select_sweeps('SwpInfo', swpInfo, 'Verbose', false, ...
'PharmConditions', pharmConditions, ...
'GIncrConditions', gIncrConditions, ...
'VHoldConditions', vHoldConditions);
% Extract whether to use the sweep
toUse = swpInfo.toUse;
% Restrict swpInfo to those sweeps
swpInfoToUse = swpInfo(toUse, :);
% Extract the cell IDs
cellIdRow = swpInfoToUse.(cellIdStr);
%% Extract measures
% Extract the measures
[ltsAmplitudeEachSwp, ltsMaxSlopeEachSwp, ltsConcavityEachSwp, ...
ltsProminenceEachSwp, ltsWidthEachSwp, ltsOnsetTimeEachSwp, ...
spikesPerLtsEachSwp, spikeMaxAmpEachSwp, spikeMinAmpEachSwp, ...
spikeFrequencyEachSwp, spikeAdaptationEachSwp, ...
burstOnsetTimeEachSwp, spikesPerBurstEachSwp] = ...
argfun(@(x) first_matching_field(swpInfoToUse, x), ...
ltsAmplitudeStr, ltsMaxSlopeStr, ltsConcavityStr, ...
ltsProminenceStr, ltsWidthStr, ltsOnsetTimeStr, ...
spikesPerLtsStr, spikeMaxAmpStr, spikeMinAmpStr, ...
spikeFrequencyStr, spikeAdaptationStr, ...
burstOnsetTimeStr, spikesPerBurstStr);
% Determine whether each sweep has an LTS or has a burst
hasLts = ltsOnsetTimeEachSwp > 0;
hasBurst = burstOnsetTimeEachSwp > 0;
% Set features to be NaN if there is no LTS
% Note: This may not be necessary for most features but just in case
% This is necessary for concavity, prominence and width as of 20191204
% TODO: Fix and run m3ha_append_lts_properties.m
ltsAmplitudeEachSwp(~hasLts) = NaN;
ltsMaxSlopeEachSwp(~hasLts) = NaN;
ltsConcavityEachSwp(~hasLts) = NaN;
ltsProminenceEachSwp(~hasLts) = NaN;
ltsWidthEachSwp(~hasLts) = NaN;
ltsOnsetTimeEachSwp(~hasLts) = NaN;
spikesPerLtsEachSwp(~hasLts) = NaN;
spikeMaxAmpEachSwp(~hasLts) = NaN;
spikeMinAmpEachSwp(~hasLts) = NaN;
spikeFrequencyEachSwp(~hasLts) = NaN;
spikeAdaptationEachSwp(~hasLts) = NaN;
burstOnsetTimeEachSwp(~hasLts) = NaN;
spikesPerBurstEachSwp(~hasLts) = NaN;
% Remove LTS amplitude values that are less than -100 mV
ltsAmplitudeEachSwp(ltsAmplitudeEachSwp < -100) = NaN;
%% Compute LTS & burst measures for each cell
% Compute the LTS or burst probability for each cell
[ltsProbability, burstProbability] = ...
argfun(@(a) array_fun(@(b) compute_feature_prob(b, cellIdRow, a), ...
uniqueCellIds), ...
hasLts, hasBurst);
% Compute means of LTS and burst properties for each cell
[ltsAmplitude, ltsMaxSlope, ltsConcavity, ...
ltsProminence, ltsWidth, ltsOnsetTime, ...
spikesPerLts, spikeMaxAmp, spikeMinAmp, ...
spikeFrequency, spikeAdaptation, burstOnsetTime, spikesPerBurst] = ...
argfun(@(x) array_fun(@(y) nanmean(x(cellIdRow == y)), uniqueCellIds), ...
ltsAmplitudeEachSwp, ltsMaxSlopeEachSwp, ltsConcavityEachSwp, ...
ltsProminenceEachSwp, ltsWidthEachSwp, ltsOnsetTimeEachSwp, ...
spikesPerLtsEachSwp, spikeMaxAmpEachSwp, spikeMinAmpEachSwp, ...
spikeFrequencyEachSwp, spikeAdaptationEachSwp, ...
burstOnsetTimeEachSwp, spikesPerBurstEachSwp);
% Compute standard deviations of LTS and burst properties for each cell
[ltsTimeJitter, burstTimeJitter] = ...
argfun(@(x) array_fun(@(y) nanstd(x(cellIdRow == y)), uniqueCellIds), ...
ltsOnsetTimeEachSwp, burstOnsetTimeEachSwp);
%% Calculate overall LTS & burst statistics
% Collect all values for each measure
% Note: Each element of the cell array contains a nCells by 1 numeric vector
% that may contain NaNs
% Note: cellfun won't work here because of namespace differences
nMeasures = numel(measureStr);
allValues = cell(nMeasures, 1);
for iMeasure = 1:nMeasures
allValues{iMeasure} = eval(measureStr{iMeasure});
end
% Compute effective n's (number of values that are not NaNs)
nValues = array_fun(@(x) sum(~isnan(x)), allValues);
% Compute overall statistics
[meanValue, stdValue, stderrValue, errValue, upper95Value, lower95Value] = ...
argfun(@(x) array_fun(@(y) compute_stats(y, x, 'IgnoreNan', true), ...
allValues), ...
'mean', 'std', 'stderr', 'err', 'upper95', 'lower95');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function featProb = compute_feature_prob (cellId, cellIdVec, hasFeatVec)
%% Computes the probability of a feature occuring
% Determine whether each sweep is from this cell
isThisCell = cellIdVec == cellId;
% Count the number of sweeps for this cell
nSweepsThisCell = sum(isThisCell);
nSweepsHasFeatThisCell = sum(isThisCell & hasFeatVec);
% Compute the probability if there is a sweep
if nSweepsThisCell > 0
featProb = nSweepsHasFeatThisCell / nSweepsThisCell;
else
featProb = NaN;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%{
OLD CODE:
swpInfoStrs = {ltsAmplitudeStr; ltsMaxSlopeStr; ltsConcavityStr; ...
ltsProminenceStr; ltsWidthStr; ltsOnsetTimeStr; ...
spikesPerLtsStr; spikeMaxAmpStr; spikeMinAmpStr; ...
spikeFrequencyStr; spikeAdaptationStr; ...
burstOnsetTimeStr; spikesPerBurstStr};
swpInfoToUse.(x)
%}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%