Dataset_Doron_Shin_et_al_PRhL1_Science_2020/codes/Histology analysis/hM4Di quantification/PlotAvgNormDreaddExpression.m

load('DreaddExpressionData.mat')
%%
AvgJS09norm = nanmean(JS09norm,2);

%%
rf = find(JS08(:,1) == 0);
AvgJS09norm = AvgJS08norm(rf-386:rf+386);

distance = JS08(rf-386:rf+386,1);

%%
TotalAvgNorm = [AvgR002norm, AvgR003norm, AvgR004norm, AvgY754norm, AvgY755norm, AvgJS09norm, ...,
    AvgJS12norm, AvgJS13norm, AvgJS14norm, AvgJS22norm; AvgJS08norm];
%%
plot(distance, TotalAvgNorm)
xlim([-1000 1000])
hold on
plot(distance, nanmean(TotalAvgNorm,2),'k', 'linew',2)
yval = get(gca,'ylim')
plot([-250 -250], yval, 'r--')
plot([750 750], yval, 'r--')


%% Quantification
interval = distance(2)-distance(1);
area = TotalAvgNorm.*interval;

TotalArea = nansum(area,1);

PRhBoundaries = [-250; 750];
PRhUpperIdx = dsearchn(distance, PRhBoundaries(1));
PRhLowerIdx = dsearchn(distance, PRhBoundaries(2));

PRhArea = nansum(area(PRhUpperIdx:PRhLowerIdx,:));
PRhAbove = nansum(area(1:PRhUpperIdx-1,:));
PRhBelow = nansum(area(PRhLowerIdx+1:end,:));

PRhPerc = 100*PRhArea./TotalArea
PRhAbovePerc = 100*PRhAbove./TotalArea
PRhBelowPerc = 100*PRhBelow./TotalArea

%Specificity = PRhPerc./NonPRhPerc


%%
group = [ones(11,1); ones(11,1)+1; ones(11,1)+2];
boxplot([PRhAbovePerc';PRhPerc';PRhBelowPerc'], group)
ylim([0 100])
ylabel('Expression (%)')
set(gca,'xticklabel', {'Above PRh', 'PRh', 'Below PRh'})