angles

near_away_locus_analysis.m

+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Analyse motion near and away from locus
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+clear all
+global minTraj 
+minTraj  = 2;
+
+%% make sure you are in the parent folder of the results folder
+px = 110    % pixel size
+[fn pt] = uigetfile('*.mat', 'Select vbspt mat file');
+%%
+[fn2 pt2] = uigetfile('*.tif', 'Select locus image');
+im = imread(fullfile(pt2, fn2));
+%im = []
+%% read in vbspt result
+a = load(fullfile(pt, fn));
+
+a.options.dim = 2;
+X = VB3_readData(a.options);
+%%
+traj_state = a.Wbest.est2.sMaxP;
+start = 1;
+stop = length(X);
+longTrajs = [];
+% Loop over trajectories longer than minTraj
+trjColoring = '<d(t)>'
+lineWidth=1;
+
+
+q = zeros(stop-start+1,5);
+for i = start:stop
+    q(i,:) = [X{i}(1,1),X{i}(1,2),X{i}(end,1)-X{i}(1,1),X{i}(end,2)-X{i}(1,2),double(size(X{i},1))];
+   
+end
+
+t = multithresh(im-imgaussfilt(im,5),2)
+figure; imagesc(im-imgaussfilt(im,5)>t(2))
+hold on
+d = bwdist(im-imgaussfilt(im,5)>t(2));
+
+
+
+for k = 1:size(q,1)
+    startpos(k) = d(max(1,floor(q(k,2)/px)),max(1,floor(q(k,1)/px)));    
+end
+id = find(startpos<10);
+not_id = find(startpos>=10);
+
+%plot near locus tracks
+for k = 1:length(id)    
+    track_length(k) = length(X{id(k)}(:,1));
+   plot(X{id(k)}(:,1)/px,X{id(k)}(:,2)/px, '.'); hold on
+end
+hold off
+figure; imagesc(im-imgaussfilt(im,5)>t(2))
+hold on
+%plot away from locus tracks
+for k = 1:length(not_id)    
+    track_length_not_id(k) = length(X{not_id(k)}(:,1));
+   plot(X{not_id(k)}(:,1)/px,X{not_id(k)}(:,2)/px, '.'); hold on
+end
+hold off
+
+%find state occupancies for near locus
+for j = 1:length(id)
+    state_one(j) = sum(traj_state{id(j)}(:)==1);
+    state_two(j) = sum(traj_state{id(j)}(:)==2);
+    state_three(j) = sum(traj_state{id(j)}(:)==3);
+    state_four(j) = sum(traj_state{id(j)}(:)==4);  
+end
+total_ones = sum(state_one)
+total_twos = sum(state_two)
+total_threes = sum(state_three)
+total_fours = sum(state_four)
+
+total_sum = total_ones+total_twos+total_threes+total_fours;
+
+one_occupancy_near = total_ones/total_sum
+two_occupancy_near = total_twos/total_sum
+three_occupancy_near = total_threes/total_sum
+four_occupancy_near = total_fours/total_sum
+
+%find state occupancies away from locus
+
+for j = 1:length(not_id)
+    state_one_away(j) = sum(traj_state{not_id(j)}(:)==1);
+    state_two_away(j) = sum(traj_state{not_id(j)}(:)==2);
+    state_three_away(j) = sum(traj_state{not_id(j)}(:)==3);
+    state_four_away(j) = sum(traj_state{not_id(j)}(:)==4);  
+end
+total_ones_away = sum(state_one_away)
+total_twos_away = sum(state_two_away)
+total_threes_away = sum(state_three_away)
+total_fours_away = sum(state_four_away)
+
+total_sum_away = total_ones_away + total_twos_away + total_threes_away + total_fours_away;
+
+one_occupancy_away = total_ones_away/total_sum_away
+two_occupancy_away = total_twos_away/total_sum_away
+three_occupancy_away = total_threes_away/total_sum_away
+four_occupancy_away = total_fours_away/total_sum_away
+
+
+%investigate angles
+data = zeros(100000,3);
+for j = 1:length(id)
+    S = sum(track_length(1:j));
+    if track_length(j) >= 3
+        for k = 2 : track_length(j)-1
+        
+            s1 = traj_state{id(j)}(k-1);
+            s2 = traj_state{id(j)}(k);
+        
+            x1 = X{id(j)}(k-1,1);
+            x2 = X{id(j)}(k,1);
+            x3 = X{id(j)}(k+1,1);
+   
+            y1 = X{id(j)}(k-1,2);
+            y2 = X{id(j)}(k,2);
+            y3 = X{id(j)}(k+1,2);
+            
+            u = [x2-x1 y2-y1];
+            v = [x3-x2 y3-y2];
+        
+            %CosTheta = max(min(dot(u,v)/(norm(u)*norm(v)),1),-1);
+            %angle = real(acosd(CosTheta));
+            
+            CosTheta =(dot(u,v)/(norm(u)*norm(v)));
+            angle = real(acos(CosTheta) * 180/pi);
+            
+            data(k+S,1) = s1;
+            data(k+S,2) = s2;
+            data(k+S,3) = angle;            
+          
+        end
+    else
+            
+    end
+end    
+
+%delete zero rows  
+data(~any(data,2),:) = [];
+data = rmmissing(data);
+figure
+histogram(data(:,3))
+histogram(data(:,3),10)
+[y1,x1]=hist(data(:,3),10); hold on
+plot(x1,y1,'r')
+grid on;
+xlabel('Angle', 'FontSize', 10);
+ylabel('Count', 'FontSize', 10);
+%save data for statistical tests on R
+%colNames = {'Initial_State','Final_State','Angle'};
+%sTable = array2table(data,'VariableNames',colNames);
+%writetable(sTable,strcat('angles_',strtok(fn, '.'), '.csv')); 
+
+%plot histograms
+figure
+for k = 1 : 4
+    data_k = data((data(:,1)==k & data(:,2)==k),:);
+    
+    no_lines = length(data_k(:,3));
+    mu = mean(data_k(:,3));
+    sigma = std(data_k(:,3));
+    subplot(2,2,k);
+    %figure
+    histogram(data_k(:,3));
+    grid on;
+    
+    xlabel('Angle', 'FontSize', 10);
+    ylabel('Count', 'FontSize', 10);
+    title({'State',k,'Line no.=',no_lines,'Mean=',mu,'std =', sigma},...
+        'FontSize', 10);
+end
+
+%investigate angles away from locus
+data_away = zeros(100000,3);
+for j = 1:length(not_id)
+    S = sum(track_length_not_id(1:j));
+    if track_length_not_id(j) >= 3
+        for k = 2 : track_length_not_id(j)-1
+        
+            s1 = traj_state{not_id(j)}(k-1);
+            s2 = traj_state{not_id(j)}(k);
+        
+            x1 = X{not_id(j)}(k-1,1);
+            x2 = X{not_id(j)}(k,1);
+            x3 = X{not_id(j)}(k+1,1);
+   
+            y1 = X{not_id(j)}(k-1,2);
+            y2 = X{not_id(j)}(k,2);
+            y3 = X{not_id(j)}(k+1,2);
+            
+            u = [x2-x1 y2-y1];
+            v = [x3-x2 y3-y2];
+            
+            CosTheta =(dot(u,v)/(norm(u)*norm(v)));
+            angle = real(acos(CosTheta) * 180/pi);
+            
+            data_away(k+S,1) = s1;
+            data_away(k+S,2) = s2;
+            data_away(k+S,3) = angle;            
+          
+        end
+    else
+            
+    end
+end    
+
+%delete zero rows  
+data_away(~any(data_away,2),:) = [];
+data_away = rmmissing(data_away);
+figure
+histogram(data_away(:,3),50)
+[y,x]=hist(data_away(:,3),50); hold on
+plot(x,y,'r')
+figure
+histfit(data_away(:,3),20,'kernel')
+grid on;
+xlabel('Angle', 'FontSize', 10);
+ylabel('Count', 'FontSize', 10);
+
+
+
+%plot 1&2 states near locus
+Locus_image = imread(fn2);
+imshow(Locus_image, [])
+hold on
+cti = 1;
+%ctii = 1;
+for j = 1: length(id)
+    one_lines = find(traj_state{id(j)}(:)==1);
+    two_lines = find(traj_state{id(j)}(:)==2);
+    three_lines = find(traj_state{id(j)}(:)==3);
+      
+       if ~isempty(one_lines) && ~isempty(two_lines)%&& isempty(three_lines)
+       %plot(tracks_new{j}(:,1),tracks_new{j}(:,2),'y')
+       one_two_transitions{cti} = traj_state{id(j)};
+       cti = cti+1;
+           for k = 2:length(traj_state{id(j)})-1
+               if traj_state{id(j)}(k-1)==1 && traj_state{id(j)}(k)==1
+               %plot([X{id(j)}(k,1)/px; X{id(j)}(k+1,1)/px] ,[X{id(j)}(k,2)/px; X{id(j)}(k+1,2)/px],'g'); hold on
+               elseif traj_state{id(j)}(k-1)==1 && traj_state{id(j)}(k)==2
+               plot([X{id(j)}(k,1)/px;X{id(j)}(k+1,1)/px] ,[X{id(j)}(k,2)/px;X{id(j)}(k+1,2)/px],'Color','green','LineWidth',2); hold on
+               elseif traj_state{id(j)}(k-1)==2 && traj_state{id(j)}(k)==2
+               %plot([X{id(j)}(k,1)/px;X{id(j)}(k+1,1)/px] ,[X{id(j)}(k,2)/px;X{id(j)}(k+1,2)/px],'b'); hold on
+               elseif traj_state{id(j)}(k-1)==2 && traj_state{id(j)}(k)==1
+               plot([X{id(j)}(k,1)/px;X{id(j)}(k+1,1)/px] ,[X{id(j)}(k,2)/px;X{id(j)}(k+1,2)/px],'Color','red','LineWidth',2); hold on
+               end
+           end
+      
+ 
+    end
+       clear one_lines
+       clear two_lines
+
+end
+
+
+%count 1-2 and 2-1 transitions near locus
+ct_ii_i = 1;
+ct_i_ii = 1;
+for j = 1: length(id)
+    for k = 1:length(traj_state{id(j)})-1
+               if traj_state{id(j)}(k)==2 && traj_state{id(j)}(k+1)==1
+                  two_to_one_transition(ct_ii_i) = 1;
+                  ct_ii_i = ct_ii_i + 1;
+               elseif traj_state{id(j)}(k)== 1 && traj_state{id(j)}(k+1)== 2
+                  one_to_two_transition(ct_i_ii) = 1;
+                  ct_i_ii = ct_i_ii + 1;
+               else
+               end
+    end
+end
+total_two_to_one = sum(two_to_one_transition)
+total_one_to_two = sum(one_to_two_transition)
+% 
+% ct_ii_i = 1;
+% ct_i_ii = 1;
+% for j = 1: length(id)
+%     for k = 1:length(traj_state{id(j)})-1
+%                if traj_state{id(j)}(k)==2 && traj_state{id(j)}(k+1)==3
+%                   two_to_three_transition(ct_ii_i) = 1;
+%                   ct_ii_i = ct_ii_i + 1;
+%                elseif traj_state{id(j)}(k)== 3 && traj_state{id(j)}(k+1)== 2
+%                   three_to_two_transition(ct_i_ii) = 1;
+%                   ct_i_ii = ct_i_ii + 1;
+%                else
+%                end
+%     end
+% end
+% total_two_to_three = sum(two_to_three_transition)
+% total_three_to_two = sum(three_to_two_transition)
+% 
+% ct_ii_i = 1;
+% ct_i_ii = 1;
+% for j = 1: length(id)
+%     for k = 1:length(traj_state{id(j)})-1
+%                if traj_state{id(j)}(k)==1 && traj_state{id(j)}(k+1)==3
+%                   one_to_three_transition(ct_ii_i) = 1;
+%                   ct_ii_i = ct_ii_i + 1;
+%                elseif traj_state{id(j)}(k)== 3 && traj_state{id(j)}(k+1)== 1
+%                   three_to_one_transition(ct_i_ii) = 1;
+%                   ct_i_ii = ct_i_ii + 1;
+%                else
+%                end
+%     end
+% end
+% total_one_to_three = sum(one_to_three_transition)
+% total_three_to_one = sum(three_to_one_transition)
+
+
+for j = 1: length(id)
+    M(j) = length(X{id(j)})-1;
+end
+total_number_tracks = sum(M)
+
+fraction_one_to_two_near_locus = total_one_to_two/total_number_tracks
+
+fraction_two_to_one_near_locus = total_two_to_one/total_number_tracks
+
+                   
+%count 1-2 and 2-1 transitions away from locus
+ct_ii_i = 1;
+ct_i_ii = 1;
+for j = 1: length(not_id)
+    for k = 1:length(traj_state{not_id(j)})-1
+               if traj_state{not_id(j)}(k)==2 && traj_state{not_id(j)}(k+1)==1
+                  two_to_one_transition_away(ct_ii_i) = 1;
+                  ct_ii_i = ct_ii_i + 1;
+               elseif traj_state{not_id(j)}(k)== 1 && traj_state{not_id(j)}(k+1)== 2
+                  one_to_two_transition_away(ct_i_ii) = 1;
+                  ct_i_ii = ct_i_ii + 1;
+               else
+               end
+    end
+end
+total_two_to_one_away = sum(two_to_one_transition_away)
+total_one_to_two_away = sum(one_to_two_transition_away)                   
+                   
+                   
+                   
+  
+