-
Notifications
You must be signed in to change notification settings - Fork 0
/
30_hypermut.m
142 lines (104 loc) · 3.4 KB
/
30_hypermut.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
%
%parse in runs {{{
clear
F = [];
F.results = direc('LNP_posteriors/signature_subcohorts_mutrate_split/*/output/*.mat');
F = parsein(F, 'results', '.*/([A-Z]+)_(hi|lo)/output/(\d+)-(.*)\.mat$', {'sig' 'rate' 'ch96' 'context'});
F = makeapn(F);
F = sort_struct(F, {'sig' 'ch96' 'rate'});
[~, ui] = unique_combos(F.sig, F.ch96);
H = rmfields(reorder_struct(F, ui), {'results' 'rate'});
H.results_hi = F.results(ui);
H.results_lo = F.results(ui + 1);
%}}}
%
% Figure S5A {{{
cmap = distinguishable_colors(16);
S = [];
S.sig = H.sig;
S.ch96 = H.ch96;
S.CR95 = cell(slength(H), 1);
%which contexts to include?
%eyeballed from plots of posterior vs. prior (see below); really, we should do this properly
%and compute some sort of mutual information (KL div, difference of entropy, JS div?)
S.include = false(slength(S), 16);
I = [];
I.sig = {'APOBEC' 'POLE' 'POLE' 'UV' 'UV'}';
I.ch96 = {96 64 95 93:95 96}';
I.pent = {9 1:16 1:16 1:16 [14 16]}';
for i = 1:slength(I),
for j = find(strcmp(S.sig, I.sig{i}) & ismember(S.ch96, I.ch96{i}))',
S.include(j, I.pent{i}) = true;
end
end
%loop over contexts to include
for i = 1:slength(H),
if ~any(S.include(i, :)), continue; end
X_hi = load(H.results_hi{i});
X_lo = load(H.results_lo{i});
sig_tots_hi = sqrt(bsxfun(@plus, 1./X_hi.tau(1:16, 500:end), 1./X_hi.tau(end, 500:end)))';
sig_tots_lo = sqrt(bsxfun(@plus, 1./X_lo.tau(1:16, 500:end), 1./X_lo.tau(end, 500:end)))';
CL = cell(16, 1);
%contour code excised from 50_signature_analysis
for j = find(S.include(i, :)),
xx = sig_tots_lo(:, j); yy = sig_tots_hi(:, j);
[px py] = meshgrid(linspace(min(xx), max(xx), 100), linspace(min(yy), max(yy), 100));
[f, xi] = ksdensity([xx yy], [px(:) py(:)]);
C = contourc(px(1, 1:100), py(1:100, 1), reshape(f, 100, []), 40);
cidx = 1;
clidx = 1;
clines = NaN;
while true,
if cidx > size(C, 2), break; end
cidx0 = C(2, cidx);
clines = C(:, 1 + cidx:(cidx + cidx0));
mp = mean(inpolygon(xx, yy, clines(1, :), clines(2, :)));
if mp >= 0.5 && mp <= 0.8, break; end
cidx = 1 + cidx + cidx0;
clidx = clidx + 1;
end
CL{j} = clines;
end
S.CR95{i} = CL;
end
cols = [];
cols.sig = {'APOBEC' 'ESO' 'LUNG' 'MSI' 'POLE' 'POLE_MSI' 'UV' 'VANILLA'}';
cols.rgb = [1 0 0; 0.74 0.68 0.052; 0 1 1; 0 1 0; 0 0 1; 0 0.5 0.5; 1 0.65 0; 0 0 0];
S.colors = mapacross(S.sig, cols.sig, cols.rgb);
S.gfx = gobjects(slength(S), 1);
%scatterplot of joint sigma_hi vs. sigma_lo
figure(1); clf
hold on
for i = 1:slength(S),
if ~any(S.include(i, :)), continue; end
c = 0;
for j = find(S.include(i, :)),
x = S.CR95{i}{j}(1, :); y = S.CR95{i}{j}(2, :);
parea = polyarea(x, y);
%0.0043 is the minimum area of any joint
%alph = 0.5*0.0043/parea;
alph = max(0.1, min(0.25*0.0043/parea, 0.95));
pobj = patch(x, y, 0, ...
'FaceColor', S.colors(i, :), 'FaceAlpha', alph, ...
'EdgeColor', 'none');
if c == 0, S.gfx(i) = pobj; end
c = c + 1;
end
end
xlim([log(2) log(4.5)])
ylim([log(2) log(4.5)])
lt = 1:0.5:4.5;
ax = gca;
ax.XTick = log(lt);
ax.YTick = log(lt);
ax.XTickLabel = lt;
ax.YTickLabel = lt;
ax.Box = 'on';
line(xlim, xlim, 'LineStyle', '--', 'Color', 0.5*[1 1 1])
axis square
grid on
xlabel('exp(\sigma) non-hypermutants')
ylabel('exp(\sigma) hypermutants')
legend(S.gfx([4 8 13]), 'APOBEC', 'POLE', 'UV', 'Location', 'NorthWest')
print('figures/hypermut_sigmas.png', '-dpng', '-r600')
% }}}