Skip to content

Commit 531e8b1

Browse files
author
Andy
committed
Change logic of fdr tests to return vector of significant net pairs instead of a single p threshold
1 parent d6adbf0 commit 531e8b1

8 files changed

Lines changed: 47 additions & 73 deletions

File tree

‎+nla/+net/+mcc/Base.m‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
classdef Base
22
methods (Abstract)
3-
p_max = correct(obj, net_atlas, input_struct, prob)
3+
[is_sig_vector, p_max] = correct(obj, net_atlas, input_struct, prob)
44
correction_label = createLabel(obj, net_atlas, input_struct, prob)
55
end
66
end

‎+nla/+net/+mcc/BenjaminiHochberg.m‎

Lines changed: 3 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -4,21 +4,11 @@
44
end
55

66
methods
7-
function p_max = correct(obj, net_atlas, input_struct, prob)
8-
[~, p_max] = nla.lib.fdr_bh(prob.v, input_struct.prob_max, 'pdep');
7+
function [is_sig_vector, p_max] = correct(obj, net_atlas, input_struct, prob)
8+
[is_sig_vector, p_max] = nla.lib.fdr_bh(prob.v, input_struct.prob_max, 'pdep');
99
end
1010
function correction_label = createLabel(obj, net_atlas, input_struct, prob)
11-
p_max = obj.correct(net_atlas, input_struct, prob);
12-
if p_max == 0
13-
correction_label = sprintf('FDR_{BH} produced no significant nets');
14-
else
15-
format_specs = "%g/%d tests";
16-
if isequal(input_struct.behavior_count, 1)
17-
format_specs = "%g/%d test";
18-
end
19-
correction_label = sprintf(strcat("FDR_{BH}(", format_specs, ")"), input_struct.prob_max * input_struct.behavior_count,...
20-
input_struct.behavior_count);
21-
end
11+
correction_label = sprintf("Benjamini-Hochberg (alpha = %g)",input_struct.prob_max);
2212
end
2313
end
2414
end

‎+nla/+net/+mcc/BenjaminiYekutieli.m‎

Lines changed: 3 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -4,21 +4,11 @@
44
end
55

66
methods
7-
function p_max = correct(obj, net_atlas, input_struct, prob)
8-
[~, p_max] = nla.lib.fdr_bh(prob.v, input_struct.prob_max, 'dep');
7+
function [is_sig_vector, p_max] = correct(obj, net_atlas, input_struct, prob)
8+
[is_sig_vector, p_max] = nla.lib.fdr_bh(prob.v, input_struct.prob_max, 'dep');
99
end
1010
function correction_label = createLabel(obj, net_atlas, input_struct, prob)
11-
p_max = obj.correct(net_atlas, input_struct, prob);
12-
if p_max == 0
13-
correction_label = sprintf('FDR_{BY} produced no significant nets');
14-
else
15-
format_specs = "%g/%d tests";
16-
if isequal(input_struct.behavior_count, 1)
17-
format_specs = "%g/%d test";
18-
end
19-
correction_label = sprintf(strcat("FDR_{BY}(", format_specs, ")"), input_struct.prob_max * input_struct.behavior_count,...
20-
input_struct.behavior_count);
21-
end
11+
correction_label = sprintf("Benjamini-Yekutieli (alpha = %g)",input_struct.prob_max);
2212
end
2313
end
2414
end

‎+nla/+net/+mcc/Bonferroni.m‎

Lines changed: 7 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -4,15 +4,18 @@
44
end
55

66
methods
7-
function p_max = correct(obj, net_atlas, input_struct, prob)
7+
function [is_sig_vector, p_max] = correct(obj, net_atlas, input_struct, prob)
8+
89
p_max = input_struct.prob_max / net_atlas.numNetPairs();
10+
is_sig_vector = prob.v < p_max;
911
end
1012
function correction_label = createLabel(obj, net_atlas, input_struct, prob)
11-
format_specs_tests = "%d tests";
13+
format_specs_tests = "%d tests)";
1214
if isequal(input_struct.behavior_count, 1)
13-
format_specs_tests = "%d test";
15+
format_specs_tests = "%d test)";
1416
end
15-
correction_label = sprintf(strcat("%g/%d net-pairs/", format_specs_tests), input_struct.prob_max * input_struct.behavior_count,...
17+
p_max = input_struct.prob_max / net_atlas.numNetPairs();
18+
correction_label = sprintf(strcat("P < %.2g (%.2g/%d net-pairs/", format_specs_tests), p_max, input_struct.prob_max * input_struct.behavior_count,...
1619
net_atlas.numNetPairs(), input_struct.behavior_count);
1720
end
1821
end

‎+nla/+net/+mcc/HolmBonferroni.m‎

Lines changed: 4 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -4,16 +4,10 @@
44
end
55

66
methods
7-
function p_max = correct(obj, net_atlas, input_struct, prob)
8-
[is_significant, adjusted_pvals, ~] = nla.lib.bonferroni_holm(prob.v, input_struct.prob_max);
9-
10-
%We need to generate a p_max threshold where all nets marked
11-
%'is_significant' are strictly below it. We can accomplish this
12-
%by finding the highest p value that passes and adding a small
13-
%delta that is unlikely to allow any additional net pairs
14-
%through the threshold
15-
DELTA = 1e-10;
16-
p_max = max(is_significant .* prob.v)+DELTA;
7+
function [is_sig_vector, p_max] = correct(obj, net_atlas, input_struct, prob)
8+
[is_sig_vector, adjusted_pvals, ~] = nla.lib.bonferroni_holm(prob.v, input_struct.prob_max);
9+
10+
p_max = max(is_sig_vector .* prob.v);
1711
end
1812
function correction_label = createLabel(obj, net_atlas, input_struct, prob)
1913
correction_label = sprintf("Holm-Bonferroni (alpha = %g)",input_struct.prob_max);

‎+nla/+net/+mcc/None.m‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -4,15 +4,16 @@
44
end
55

66
methods
7-
function p_max = correct(obj, net_atlas, input_struct, prob)
7+
function [is_sig_vector, p_max] = correct(obj, net_atlas, input_struct, prob)
88
p_max = input_struct.prob_max;
9+
is_sig_vector = prob.v < p_max;
910
end
1011
function correction_label = createLabel(obj, net_atlas, input_struct, prob)
11-
format_specs = "%g/%d tests";
12+
format_specs = "P < %.2g (%g/%d tests)";
1213
if isequal(input_struct.behavior_count, 1)
13-
format_specs = "%g/%d test";
14+
format_specs = "P < %.2g (%g/%d test)";
1415
end
15-
correction_label = sprintf(format_specs, input_struct.prob_max * input_struct.behavior_count,...
16+
correction_label = sprintf(format_specs, input_struct.prob_max, input_struct.prob_max * input_struct.behavior_count,...
1617
input_struct.behavior_count);
1718
end
1819
end

‎+nla/+net/+result/NetworkResultPlotParameter.m‎

Lines changed: 18 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,7 @@
2727
end
2828

2929
function result = plotProbabilityParameters(obj, edge_test_options, edge_test_result, test_method, plot_statistic,...
30-
plot_title, fdr_correction, significance_filter, ranking_method)
30+
plot_title, fdr_correction, effect_size_filter, ranking_method)
3131
% plot_title - this will be a string
3232
% plot_statistic - this is the stat that will be plotted
3333
% significance filter - this will be a boolean or some sort of object (like Cohen's D > D-value)
@@ -36,9 +36,9 @@
3636

3737
import nla.TriMatrix nla.TriMatrixDiag
3838
% We're going to use a default filter here
39-
if isequal(significance_filter, false)
40-
significance_filter = TriMatrix(obj.number_of_networks, "logical", TriMatrixDiag.KEEP_DIAGONAL);
41-
significance_filter.v = true(numel(significance_filter.v), 1);
39+
if isequal(effect_size_filter, false)
40+
effect_size_filter = TriMatrix(obj.number_of_networks, "logical", TriMatrixDiag.KEEP_DIAGONAL);
41+
effect_size_filter.v = true(numel(effect_size_filter.v), 1);
4242
end
4343

4444
% Adding on to the plot title if it's a -log10 plot
@@ -53,30 +53,26 @@
5353
if isstring(fdr_correction) || ischar(fdr_correction)
5454
fdr_correction = nla.net.mcc.(erase(fdr_correction, "-"))();
5555
end
56-
p_value_max = fdr_correction.correct(obj.network_atlas, obj.updated_test_options, statistic_input);
56+
57+
[is_sig_vector, p_value_max] = fdr_correction.correct(obj.network_atlas, obj.updated_test_options, statistic_input);
58+
5759
p_value_breakdown_label = fdr_correction.createLabel(obj.network_atlas, obj.updated_test_options,...
5860
statistic_input);
59-
60-
if isa(fdr_correction,'nla.net.mcc.HolmBonferroni')
61-
%Holm Bonferroni does not have a single static p threshold,
62-
%so use the custom label for this FDR that only includes
63-
%the alpha (original p_max from input struct), and not any
64-
%other p threshold
65-
name_label = sprintf("%s %s\n%s", obj.network_test_results.test_display_name, plot_title,...
61+
name_label = sprintf("%s %s\n%s", obj.network_test_results.test_display_name, plot_title,...
6662
p_value_breakdown_label);
67-
else
68-
69-
name_label = sprintf("%s %s\nP < %.2g (%s)", obj.network_test_results.test_display_name, plot_title,...
70-
p_value_max, p_value_breakdown_label);
71-
if p_value_max == 0
72-
name_label = sprintf("%s %s\nP = %.2g (%s)", obj.network_test_results.test_display_name, plot_title,...
73-
p_value_max, p_value_breakdown_label);
74-
end
75-
end
63+
% else
64+
%
65+
% name_label = sprintf("%s %s\nP < %.2g (%s)", obj.network_test_results.test_display_name, plot_title,...
66+
% p_value_max, p_value_breakdown_label);
67+
% if p_value_max == 0
68+
% name_label = sprintf("%s %s\nP = %.2g (%s)", obj.network_test_results.test_display_name, plot_title,...
69+
% p_value_max, p_value_breakdown_label);
70+
% end
71+
% end
7672

7773
% Filtering if there's a filter provided
7874
significance_plot = TriMatrix(obj.number_of_networks, "logical", TriMatrixDiag.KEEP_DIAGONAL);
79-
significance_plot.v = (statistic_input.v < p_value_max) & significance_filter.v;
75+
significance_plot.v = is_sig_vector & effect_size_filter.v;
8076

8177
% scale values very slightly for display so numbers just below
8278
% the threshold don't show up white but marked significant

‎+nla/+net/+result/NetworkTestResult.m‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -285,15 +285,15 @@ function createPValueTriMatrices(obj, number_of_networks, test_method)
285285
% I don't really know what these do and haven't really thought about it. Hence the bad naming.
286286
function [sig, name] = singleSigMat(obj, network_atlas, edge_test_options, p_value, mcc_method, title_prefix)
287287
mcc_method = nla.net.mcc.(mcc_method)();
288-
p_value_max = mcc_method.correct(network_atlas, edge_test_options, p_value);
288+
[is_sig_vector, p_value_max] = mcc_method.correct(network_atlas, edge_test_options, p_value);
289289
p_breakdown_labels = mcc_method.createLabel(network_atlas, edge_test_options, p_value);
290290

291291
sig = nla.TriMatrix(network_atlas.numNets(), 'double', nla.TriMatrixDiag.KEEP_DIAGONAL);
292-
sig.v = (p_value.v < p_value_max);
293-
name = sprintf("%s %s P < %.2g (%s)", title_prefix, obj.test_display_name, p_value_max, p_breakdown_labels);
294-
if p_value_max == 0
295-
name = sprintf("%s %s P = 0 (%s)", title_prefix, obj.test_display_name, p_breakdown_labels);
296-
end
292+
sig.v = is_sig_vector;
293+
name = sprintf("%s %s %s", title_prefix, obj.test_display_name, p_value_max, p_breakdown_labels);
294+
% if p_value_max == 0
295+
% name = sprintf("%s %s P = 0 (%s)", title_prefix, obj.test_display_name, p_breakdown_labels);
296+
% end
297297
end
298298

299299
function [number_of_tests, sig_count_mat, names] = appendSignificanceMatrix(...

0 commit comments

Comments
 (0)