%% Calculate nearest neighbor pairs for one condition
% For in vivo spine data analyzed using MapManager V_20171028 
% Author: Richard Roth (rhroth@stanford.edu)
% Date: Sept 2023

%     Requires the following functions:
%     getmyMapListValuesCondDist_Neighbour.m
%     segmentCall_withDist.m
%     removeNanRows.m

% (1) Retrive values
%     This requires myMapList and mySegmentTable loaded in workspace and

% (2) Plot Stat vs Stat of nearest neighbor
%     Generates simple correlation plot

% (3) Randomize spine pairing per segment then combining
%     Use these to generate 10,000 randomized spine pairings, by 
%     randmoizing only spines withing each segment.

% (4) Plot histograms
%     Generate histograms with distribution of seeing randomized spine
%     pairs with stat changes in same or opposite direction and plotting
%     actually observed value in original data

% (5) Obtain observed values per dendrite segment

%% (1) Retrive values
Condition = 'Gr1';
% Condition = 'Gr2';

% Make pool using same SessionConditions as used for other ways of ploting data
SessionConditions = {'day-3', 'day-2', 'day-1', 'day0', 'day1', 'day2', 'day3', 'day4'};

% Define sessions from which values should be averaged for neighbor analysis
NeighborConditions = {'day1', 'day2', 'day3', 'day4'};

IncThresh = 0.14; % Threshold value above which spines are considered "increased"
DecThresh = -0.14; % Threshold value below which spines are considered "decreased"

range=10; % range for histogram

StatName = 'nntb2utsmdtdmch1'; %Spine dsRed / Dendrite dsRed mean normalized
Channel = 1; %Only relevant for channel specific stats

SecStatName = 'pDist'; %Used for calculating position of spine (by default the average across all sessions is used for the position)

PersCond = 1; %Minimum number of sessions that spine needs to be present in
MinSpines = 1; %Minimum number of spines per dendrite for dendrite to be included

mynnb = getmyMapListValuesCondDist_Neighbour(myMapList, mySegmentTable, Condition, SessionConditions, NeighborConditions, StatName, SecStatName, Channel, PersCond, MinSpines);
clearvars PersCond MinSpines SecStatName Channel StatName NeighborConditions SessionConditions
%% (2) Plot Stat vs Stat of nearest neighbor
nnball=mynnb.nnball;
nnball_NaNremoved=nnball(~any(isnan(nnball),2),:);

figure
% scatter(nnball_NaNremoved(:,1),nnball_NaNremoved(:,2),'b','LineWidth',2)
scatter(nnball_NaNremoved(:,1),nnball_NaNremoved(:,2),'r','LineWidth',2)

P = polyfit(nnball_NaNremoved(:,1),nnball_NaNremoved(:,2),1);
yfit = P(1)*nnball_NaNremoved(:,1)+P(2);

hold on;
plot(nnball_NaNremoved(:,1),yfit,'m-','LineWidth',3);

xlim([-1 1])
ylim([-1 1])
xlabel('Spine Size Fold Change')
ylabel('Spine Size Fold Change in Nearest Neighbor')
title([Condition,'  Nearest Neighbor'])

% Plot appearance settings
ax = gca;
ax.Color = 'white';
ax.Box = 'off';
ax.TickDir = 'out';
ax.TickLength = [0.015 0.035];
ax.FontSize = 13;
ax.LineWidth = 1.75;
set(0,'defaultAxesColorOrder','remove')
set(gcf,'color','w');


%% (3) Randomize spine pairing per segment then combining
numofrand = 10000; %Number of randomizations

nnball3D = mynnb.nnball3D;
shuffle_inc=[];
shuffle_oposite=[];
shuffle_dec=[];
shuffle_same=[];
allspines=[];

for b=1:numofrand
nnball_perm=[];

for row=1:size(nnball3D,1)
currentnnball_perm=[];
currentnnball=nnball3D{row,1};
nnball_NaNremoved=currentnnball(~any(isnan(currentnnball),2),:);
idx = randperm(size(nnball_NaNremoved,1));
currentnnball_perm(idx,1)=nnball_NaNremoved(:,1);
currentnnball_perm(:,2)=nnball_NaNremoved(:,2);
currentnnball_perm(:,3)=nnball_NaNremoved(:,3);
nnball_perm = [nnball_perm ; currentnnball_perm];
end

shuffle_inc=[shuffle_inc ; sum(nnball_perm(:,1) > IncThresh & nnball_perm(:,2) > IncThresh)];
shuffle_oposite=[shuffle_oposite ; sum(nnball_perm(:,1) < DecThresh & nnball_perm(:,2) > IncThresh) + sum(nnball_perm(:,1) > IncThresh & nnball_perm(:,2) < DecThresh)];
shuffle_dec=[shuffle_dec ; sum(nnball_perm(:,1) < DecThresh & nnball_perm(:,2) < DecThresh)];
shuffle_same=[shuffle_same ; sum(nnball_perm(:,1) > IncThresh & nnball_perm(:,2) > IncThresh) + sum(nnball_perm(:,1) < DecThresh & nnball_perm(:,2) < DecThresh)];
allspines = [allspines ; sum(~isnan(nnball_perm(:,1)) & ~isnan(nnball_perm(:,2)))]; %changed from [allspines ; sum(nnball_perm(:,1) >= 0)] on 20230831

end

nnball=mynnb.nnball;
nnball_NaNremoved=nnball(~any(isnan(nnball),2),:);

obs_inc=sum(nnball_NaNremoved(:,1) > IncThresh & nnball_NaNremoved(:,2) > IncThresh);
obs_oposite=sum(nnball_NaNremoved(:,1) < DecThresh & nnball_NaNremoved(:,2) > IncThresh)+sum(nnball_NaNremoved(:,1) > IncThresh & nnball_NaNremoved(:,2) < DecThresh);
obs_dec=sum(nnball_NaNremoved(:,1) < DecThresh & nnball_NaNremoved(:,2) < DecThresh);
obs_same=sum(nnball_NaNremoved(:,1) > IncThresh & nnball_NaNremoved(:,2) > IncThresh)+sum(nnball_NaNremoved(:,1) < DecThresh & nnball_NaNremoved(:,2) < DecThresh);
obs_all=sum(~isnan(nnball_NaNremoved(:,1)) & ~isnan(nnball_NaNremoved(:,2))); %changed from obs_all=sum(nnball_NaNremoved(:,1) > 0) on 20230831


%plot
figure
% scatter(nnball_perm(:,1),nnball_perm(:,2),'b','LineWidth',2)
scatter(nnball_perm(:,1),nnball_perm(:,2),'r','LineWidth',2)
nnball_perm_NaNremoved=nnball_perm(~any(isnan(nnball_perm),2),:);
P = polyfit(nnball_perm_NaNremoved(:,1),nnball_perm_NaNremoved(:,2),1);
yfit = P(1)*nnball_perm_NaNremoved(:,1)+P(2);
hold on;
plot(nnball_perm_NaNremoved(:,1),yfit,'m-','LineWidth',3);
ylim([-1 1])
xlim([-1 1])
xlabel('Spine Size Fold Change')
ylabel('Spine Size Fold Change in Nearest Neighbor')
title([Condition,'   Nearest Neighbour shuffled'])
% Plot appearance settings
ax = gca;
ax.Color = 'white';
ax.Box = 'off';
ax.TickDir = 'out';
ax.TickLength = [0.015 0.035];
ax.FontSize = 13;
ax.LineWidth = 1.75;
set(0,'defaultAxesColorOrder','remove')
set(gcf,'color','w');

clearvars P yfit nnball_perm_NaNremoved nnball_perm nnball_NaNremoved nnball allspines currentnnball currentnnball_perm idx row b numofrand nnball3D
%% (4) Plot histograms
% Generate histograms with distribution of seeing randomized spine
% pairs with stat changes in same or oposite direction and plotting
% actually observed value in original data
% Generates output variable xout and histout for exporting histogram to
% other software such as Prism

% Set histogram properties

numbins=51;

shuffeled_values = cell(4,1);
shuffeled_values{1} = shuffle_inc./obs_all*100;
shuffeled_values{2} = shuffle_dec./obs_all*100;
shuffeled_values{3} = shuffle_oposite./obs_all*100;
shuffeled_values{4} = shuffle_same./obs_all*100;

observered_values(1) = obs_inc./obs_all*100;
observered_values(2) = obs_dec./obs_all*100;
observered_values(3) = obs_oposite./obs_all*100;
observered_values(4) = obs_same./obs_all*100;

titles{1}='Both Increase';
titles{2}='Both Decrease';
titles{3}='Change in Opposite Direction';
titles{4}='Change in Same Direction';

figure
for j=1:4
val = shuffeled_values{j};
obs_plot = observered_values(j);
if min(shuffeled_values{j})-5 >0
minval=round(min(shuffeled_values{j})-5,0);
else
minval=0;
end
maximumval=minval+range;

binwidth = range ./ (numbins-1);

session=1;
histo=[];
cumhist=[];
i=1;
histo(i)=(sum(val(:,session)<=((0.5*binwidth)+minval)))./sum(val(:,session)>0);
cumhist(i)=sum(val(:,session)<=((0.5*binwidth)+minval))./sum(val(:,session)>0);
for i = 2:(numbins-1)
histo(i)=(sum(val(:,session)<=(((i-0.5)*binwidth)+minval)) - sum(val(:,session)<=(((i-1.5)*binwidth)+minval)))./sum(val(:,session)>0);
cumhist(i)=sum(val(:,session)<=(((i-0.5)*binwidth)+minval)) ./ sum(val(:,session)>0);
end
i=numbins;
histo(i)=(sum(val(:,session)>(((i-1.5)*binwidth)+minval))) ./ sum(val(:,session)>0);
cumhist(i)=sum(val(:,session) > minval)./sum(val(:,session)>0);

%find approximate p-value by summing tails
if obs_plot > nanmean(val)
pval = 2*(sum(val>=obs_plot)+1)./10001;
else
pval = 2*(sum(val<=obs_plot)+1)./10001;
end

subplot(1,4,j)
p1=bar([minval:((range)./(numbins-1)):range+minval],histo(:)*100,1.2,'FaceColor',[1 0 0],'EdgeColor',[1 0 0]);
% p1=bar([minval:((range)./(numbins-1)):range+minval],histo(:)*100,1.2,'FaceColor',[0 0 1],'EdgeColor',[0 0 1]);
hold on
p2=plot([obs_plot obs_plot],[0 100],'--k','LineWidth',2);

ylim([0 30])
% xlim([minval maximumval])
xlim([0 10])
xlabel(['Fraction of Spines (%)'])
ylabel(['Count (%)'])
title([titles{j}])
legend([p2],['Observed Value=' num2str(round(obs_plot,1)) '%' newline 'p=' num2str(round(pval,3))])

% Plot appearance settings
ax = gca;
ax.LineWidth = 1.5;
ax.Color = 'white';
ax.Box = 'off';
ax.TickDir = 'out';
ax.TickLength = [0.015 0.035];
ax.FontSize = 13;
ax.LineWidth = 1.75;
set(0,'defaultAxesColorOrder','remove')
set(gcf,'color','w');



% Use this to copy histogram data to Prism etc. 
% Column 1: Both Increase, Column 2: Both Decrease, Column 3: Change in Oposite Direction, Column 4: Change in Same Direction
xout(:,j) = [minval:((range)./(numbins-1)):range+minval]';
histout(:,j) = histo'*100;
end

clearvars ax binwidth cumhist histo i j maximumval minval numbins p1 p2 pval session titles val


%% (5) Obtain observed values per dendrite segment

persegment_obs_inc=[];
persegment_obs_oposite=[];
persegment_obs_dec=[];
persegment_obs_same=[];
persegment_obs_all=[];

for i=1:size(mynnb.nnball3D,1)

currentnnb = mynnb.nnball3D{i,1};
    
current_obs_inc=sum(currentnnb(:,1) > IncThresh & currentnnb(:,2) > IncThresh);
current_obs_oposite=sum(currentnnb(:,1) < DecThresh & currentnnb(:,2) > IncThresh)+sum(currentnnb(:,1) > IncThresh & currentnnb(:,2) < DecThresh);
current_obs_dec=sum(currentnnb(:,1) < DecThresh & currentnnb(:,2) < DecThresh);
current_obs_same=sum(currentnnb(:,1) > IncThresh & currentnnb(:,2) > IncThresh)+sum(currentnnb(:,1) < DecThresh & currentnnb(:,2) < DecThresh);
current_obs_all=sum(~isnan(currentnnb(:,1)) & ~isnan(currentnnb(:,2)));


persegment_obs_inc = [persegment_obs_inc; current_obs_inc];
persegment_obs_oposite = [persegment_obs_oposite; current_obs_oposite];
persegment_obs_dec = [persegment_obs_dec; current_obs_dec];
persegment_obs_same = [persegment_obs_same; current_obs_same];
persegment_obs_all = [persegment_obs_all; current_obs_all];

end

persegment_obs_inc_percentage = persegment_obs_inc./persegment_obs_all*100;
persegment_obs_oposite_percentage = persegment_obs_oposite./persegment_obs_all*100;
persegment_obs_dec_percentage = persegment_obs_dec./persegment_obs_all*100;
persegment_obs_same_percentage = persegment_obs_same./persegment_obs_all*100;


clearvars currentnnb current_obs_inc current_obs_oposite current_obs_dec current_obs_same current_obs_all

