forked from rejniaklab/LinG3D
-
Notifications
You must be signed in to change notification settings - Fork 0
/
linG3DAll.m
executable file
·212 lines (177 loc) · 8.74 KB
/
linG3DAll.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
function linG3DAll(pathData,numClones,IsGradient,xmin,xmax,...
ymin,ymax,tmin,tmax,fileStep,toPrint)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% This is a companion code for the paper "LinG3D: Visualizing the %%
%% Spatio-Temporal Dynamics of Clonal Evolution" by A. Hu, A.M.E. %%
%% Ojwang', K.D. Olumoyin, and K.A. Rejniak %%
%% This code generates the 3D lineage tree of all clones %%
%% %%
%% The following parameters need to be specified: %%
%% pathData -- directory with input data %%
%% numClones -- total number of clones in the data %%
%% IsGradient -- 1 to draw drug in the background, 0 not to draw %%
%% xmin,xmax,ymin,ymax -- dimensions of the spacial domain %%
%% tmin, tmax -- dimensions of the temporal domain %%
%% fileStep -- frequency of the sampled data %%
%% toPrint -- 1 to save the generated figure, 0 not to save %%
%% It requires the following data in the pathData/data/ directory: %%
%% cell_history.txt -- file with info about each cell %%
%% cellID_##.txt -- cell IDs in a file with index number ## %%
%% cellXY_##.txt -- cell coordinates in a file with index ## %%
%% drug.txt -- concentration of a drug for background %%
%% for the examples discussed in the paper use: %%
%% example 1: pathData='exampleB05'; numClones=9; %%
%% example 2: pathData='exampleB005'; numClones=147; %%
%% example 3: pathData='exampleExp'; numClones=10; %%
%% %%
%% January 10, 2024 %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
timeStep=(tmax-tmin)/(2.5*(xmax-xmin)); %frequency of sampled temporal data
% draw 3D figure
figure
grid
axis([xmin,xmax,0,tmax/timeStep,ymin,ymax])
axis equal
view(3)
hold on
col=DefineColorPalette; Ncol=size(col,1); % define color palette
dataDirectory='/data/'; % directory with cell and drug data
% draw background with drug gradient
if IsGradient==1
drug=load([pathData,dataDirectory,'drug.txt']);
DrawBackground(drug,tmax,timeStep,xmin,xmax,ymin,ymax)
end
% load cell history file
hist=load([pathData,dataDirectory,'cell_history.txt']);
% [cell ID, clone ID, mother ID, birth iter, div/death iter]
% inspect all cells for a given clone
for cloneNum=0:numClones
disp(['clone=',num2str(cloneNum),' of ',num2str(numClones)]);
% define matrix of line segments (3D branches) to draw
matrix_to_draw=zeros(1,6); % [x1,k1,y1,x2,k2,y2]
Nmatrix=0;
indLast=find(hist(:,2)==cloneNum);
for ii=1:length(indLast) % for every cell with index in indLast
if mod(ii,100)==0 disp('... calculating'); end
cellNum=hist(indLast(ii),1); % cell ID
mothNum=hist(indLast(ii),3); % mother ID
strtNum=hist(indLast(ii),4); % cell birth
endNum =hist(indLast(ii),5);
endNum=max(tmax,min(endNum,tmax)); % cell div/death/tmax
% find all appearances of the cellNum
kkStart=fileStep*floor(strtNum/fileStep); % initial file number
kkEnd =fileStep*floor(endNum/fileStep); % final file number
for kk=kkEnd:-fileStep:kkStart+fileStep % inspect all files
% cell ID and cell XY from the first file
fileMeID =load([pathData,dataDirectory,'cellID_',num2str(kk),'.txt']);
fileMeXY =load([pathData,dataDirectory,'cellXY_',num2str(kk),'.txt']);
indMe =find(fileMeID==cellNum); % find current indices of cellID
% cell ID and cell XY from the second file
fileMe2ID=load([pathData,dataDirectory,'cellID_',num2str(kk-fileStep),'.txt']);
fileMe2XY=load([pathData,dataDirectory,'cellXY_',num2str(kk-fileStep),'.txt']);
indMe2=find(fileMe2ID==cellNum); % find current indices of cellID
if isempty(indMe)
elseif isempty(indMe2)
while kkStart<hist(mothNum,4) % find file with the grand-mother cell
mothNum=hist(mothNum,3);
end
fileMe2ID=load([pathData,dataDirectory,'cellID_',num2str(kkStart),'.txt']);
fileMe2XY=load([pathData,dataDirectory,'cellXY_',num2str(kkStart),'.txt']);
indMe2=find(fileMe2ID==mothNum); % find current indices of mother cellID
if isempty(indMe2)
else
Nmatrix=Nmatrix+1; % save branch to draw [x1,t1,y1,x2,t2,y2]
matrix_to_draw(Nmatrix,1:6)=[fileMeXY(indMe,1),kkStart+fileStep,...
fileMeXY(indMe,2),fileMe2XY(indMe2,1),kkStart,fileMe2XY(indMe2,2)];
end
else
Nmatrix=Nmatrix+1; % save branch to draw [x1,t1,y1,x2,t2,y2]
matrix_to_draw(Nmatrix,1:6)=[fileMe2XY(indMe2,1),kk-fileStep,...
fileMe2XY(indMe2,2),fileMeXY(indMe,1),kk,fileMeXY(indMe,2)];
end
end
end
% drawing clones
NumCol=mod(cloneNum,Ncol)+1; % clone color
title(['clone=',num2str(cloneNum)])
axis equal
axis([xmin,xmax,tmin,tmax/timeStep,ymin,ymax])
view(3)
ylabel(['iterations/time x',num2str(timeStep)])
hYLabel = get(gca,'YLabel');
set(hYLabel,'rotation',-35) %,'VerticalAlignment','middle')
hold on
for ii=1:Nmatrix
plot3([matrix_to_draw(ii,1),matrix_to_draw(ii,4)],...
[matrix_to_draw(ii,2),matrix_to_draw(ii,5)]/timeStep,...
[matrix_to_draw(ii,3),matrix_to_draw(ii,6)],'color',...
col(NumCol,1:3),'linewidth',1.5)
if mod(ii,100)==0 pause(0.1); end
end
view(3)
axis equal
axis([xmin,xmax,tmin,tmax/timeStep,ymin,ymax])
pause(0.01)
end%for cloneNum
title('final 3D traces of all clones')
view(3)
axis equal
axis([xmin,xmax,tmin,tmax/timeStep,ymin,ymax])
if toPrint==1
print('-djpeg100','-r300',[pathData,'/tree_clones_combined.jpg'])
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function DrawBackground(drug,tmax,timeStep,xmin,xmax,ymin,ymax)
drugmin=min(min(drug)); drugmax=max(max(drug)); drugstep=(drugmax-drugmin)/4;
kk=tmax/timeStep;
[Nx,Ny]=size(drug); hgx=(xmax-xmin)/Nx; hgy=(ymax-ymin)/Ny;
for ii=1:Nx
for jj=1:Ny
if (drug(ii,jj)>=drugmin)&&(drug(ii,jj)<drugmin+drugstep)
patch([xmin+(ii-1)*hgx,xmin+ii*hgx,xmin+ii*hgx,xmin+(ii-1)*hgx,...
xmin+(ii-1)*hgx],[kk,kk,kk,kk,kk],[ymin+(jj-1)*hgy,ymin+(jj-1)*hgy,...
ymin+jj*hgy,ymin+jj*hgy,ymin+(jj-1)*hgy],'b','edgecolor','none')
hold on
elseif (drug(ii,jj)>=drugmin+drugstep)&&(drug(ii,jj)<drugmin+2*drugstep)
patch([xmin+(ii-1)*hgx,xmin+ii*hgx,xmin+ii*hgx,xmin+(ii-1)*hgx,...
xmin+(ii-1)*hgx],[kk,kk,kk,kk,kk],[ymin+(jj-1)*hgy,ymin+(jj-1)*hgy,...
ymin+jj*hgy,ymin+jj*hgy,ymin+(jj-1)*hgy],'c','edgecolor','none')
hold on
elseif (drug(ii,jj)>=drugmin+2*drugstep)&&(drug(ii,jj)<drugmin+3*drugstep)
patch([xmin+(ii-1)*hgx,xmin+ii*hgx,xmin+ii*hgx,xmin+(ii-1)*hgx,...
xmin+(ii-1)*hgx],[kk,kk,kk,kk,kk],[ymin+(jj-1)*hgy,ymin+(jj-1)*hgy,...
ymin+jj*hgy,ymin+jj*hgy,ymin+(jj-1)*hgy],'y','edgecolor','none')
hold on
else
patch([xmin+(ii-1)*hgx,xmin+ii*hgx,xmin+ii*hgx,xmin+(ii-1)*hgx,...
xmin+(ii-1)*hgx],[kk,kk,kk,kk,kk],[ymin+(jj-1)*hgy,ymin+(jj-1)*hgy,...
ymin+jj*hgy,ymin+jj*hgy,ymin+(jj-1)*hgy],'r','edgecolor','none')
hold on
end
end
end
alpha(0.25)
hold on
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function col=DefineColorPalette
col=[255,0,255;255,0,0;0,255,255;0,0,255;0,255,0;0,0,0;255,191,0;...
255,255,0;191,255,0;128,128,0;255,182,193;0,191,255;0,128,255;...
250,235,215;128,0,255;154,205,50;255,0,128;102,0,0;102,77,0;...
0,102,102;204,204,255;255,204,255;153,204,255;255,153,153;0,153,0;...
0,153,153;153,0,77;255,228,225;128,0,0;102,102,153;153,255,204;...
218,112,214;255,128,0;192,192,192;128,128,128;75,0,130;165,42,42;...
216,191,216;220,20,60;245,222,179;255,99,71;255,127,80;205,92,92;...
240,128,128;233,150,122;250,128,114;255,160,122;255,69,0;...
255,140,0;255,165,0;255,215,0;184,134,11;218,165,32;0,100,0;...
255,240,245;188,143,143;255,248,220;50,205,50;144,238,144;...
152,251,152;143,188,143;0,250,154;0,255,127;46,139,87;102,205,170;...
60,179,113;32,178,170;47,79,79;0,128,128;0,139,139;240,230,140;...
245,245,220;224,255,255;0,206,209;255,228,181;255,20,147;...
175,238,238;127,255,212;176,224,230;95,158,160;70,130,180;...
100,149,237;222,184,135;30,144,255;238,232,170;189,183,107;...
107,142,35;124,252,0;127,255,0;173,255,47;...
178,34,34;221,160,221;255,235,205;]/255;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%