%% chromosome similarity

%% Load theory
% addpath(genpath('/Users/moacarlsson/Desktop/Yeast reference files'))
% fastas=dir('/Users/moacarlsson/Desktop/Yeast reference files/*fasta');
% addpath(genpath('C:/Users/luisma/OneDrive - Chalmers/yeastdata/ReferenceScer/*fasta'))
% %fastas=dir('C:/Users/luisma/OneDrive - Chalmers/yeastdata/ReferenceScer/*fasta');
addpath(genpath('C:\Users\luisma\OneDrive - Chalmers\S.pombe genome proyect\Spombe_Genome\ref_fastas'))
fastas=dir('C:\Users\luisma\OneDrive - Chalmers\S.pombe genome proyect\Spombe_Genome\ref_fastas\*fasta');

for i=1:length(fastas)
fastaFile{1,i}=strjoin({fastas(i).folder,'/',fastas(i).name},"");
end

%quickly calculate theory for given single fasta file (fastaFile, nmbp, nmPerPx, psffac, islinear)
[theoryStructRev,theoryStruct,bT] = prep_thry_for_local_comp(fastaFile, 0.25, 1,254,1); 

%% settings

timestamp = datestr(clock(), 'yyyy-mm-dd_HH_MM_SS');

%% Data Prep
% make bitmask logical and add the chromosome name and theory check
for i=1:length(bT_newpsf)
bT_newpsf{i}.rawBitmask = logical(bT{i}.rawBitmask);
bT_newpsf{i}.name=theoryStruct_newpsf(i).name;
bT_newpsf{i}.theory=i;
end

for i=1:length(bT_newpsf_scer)
bT_newpsf_scer{i}.rawBitmask = logical(bT_newpsf_scer{i}.rawBitmask);
bT_newpsf_scer{i}.name=theoryStruct_newpsf_scer(i).name;
bT_newpsf_scer{i}.theory=i;
end
% Adds cushion 6px mask on the edges
bitmask_buff=0;
bT_newpsf=expand_mask(bT_newpsf,bitmask_buff,1);
bT_newpsf_scer=expand_mask(bT_newpsf_scer,bitmask_buff,1);
%
%% oS
sF=1
bT=bT_newpsf_pombe;
w=150;
[oSTheory{1},mpI1,mp1] = calc_overlap_mp(bT,sF, w,timestamp);
% In this one, lenB is static and lenA is stretched along it
% This means each column features each barcode/theory being static
% and everything else stretching to match it
% each row, however, represents the efforts of said barcode to match the
% others being theories

%% Exploratory
figure,tiledlayout(2,1)
nexttile
plot(mpI1{1})
nexttile
plot(mp1{1})

%% chr 1 mapping into 2
% Parameters
minRunLength = 10;   % minimum consecutive points
gapTolerance = 5;    % allowed missing points
extendLength = w;  % artificial extension length
tol = 2;             % slope tolerance
 
 
figure,tiledlayout(2,2)
nexttile

ccthreshold=0.75;
i=2;
x=1:5330;
y=mpI1{i}(x);
max_index=length(bT{i}.rawBarcode(bT{i}.rawBitmask))-w;
min_big=length(bT{i}.rawBarcode(bT{i}.rawBitmask))+1;
ceiling=min_big*2-w+1;
y(y>=min_big)=ceiling-y(y>=min_big);
%toplot
xplot=x(mp1{i}(x)>=ccthreshold);
yplot=y(xplot);
extend_chr_sim(xplot, yplot, 'minRunLength', minRunLength, 'gapTolerance', gapTolerance, 'extendLength', extendLength, 'tol', tol);
title('Chr1 vs Chr2');
xlabel('Chr1')
ylabel('Chr2')

%%
minRunLength = 10;   % minimum consecutive points
gapTolerance = 2;    % allowed missing points
extendLength = w;  % artificial extension length
tol = 1;             % slope tolerance
 ccthreshold=0.75;
i=2;
x=1:5330;
y=mpI1{i}(x);
max_index=length(bT{i}.rawBarcode(bT{i}.rawBitmask))-w;
min_big=length(bT{i}.rawBarcode(bT{i}.rawBitmask))+1;
ceiling=min_big*2-w+1;
y(y>=min_big)=ceiling-y(y>=min_big);
%toplot
xplot=x(mp1{i}(x)>=ccthreshold);
yplot=y(xplot);
 
chr_sim_marginals(bT_newpsf_pombe{1}, bT_newpsf_pombe{2}, xplot, yplot,minRunLength, gapTolerance, extendLength, tol)
%% chr 1 mapping into 3

nexttile


i=3;
x=1:5330;
y=mpI1{i}(x);
max_index=length(bT{i}.rawBarcode(bT{i}.rawBitmask))-150;
min_big=length(bT{i}.rawBarcode(bT{i}.rawBitmask))+1;
ceiling=min_big*2-150+1;
y(y>=min_big)=ceiling-y(y>=min_big);
%toplot
xplot=x(mp1{i}(x)>=ccthreshold);
yplot=y(xplot);
extend_chr_sim(xplot, yplot, 'minRunLength', minRunLength, 'gapTolerance', gapTolerance, 'extendLength', extendLength, 'tol', tol);
title('Chr1 vs Chr3');
xlabel('Chr1')
ylabel('Chr3')

%%
minRunLength = 10;   % minimum consecutive points
gapTolerance = 1;    % allowed missing points
extendLength = w;  % artificial extension length
tol = 1;             % slope tolerance
 ccthreshold=0.75;
i=3;
x=1:5330;
y=mpI1{i}(x);
max_index=length(bT{i}.rawBarcode(bT{i}.rawBitmask))-w;
min_big=length(bT{i}.rawBarcode(bT{i}.rawBitmask))+1;
ceiling=min_big*2-w+1;
y(y>=min_big)=ceiling-y(y>=min_big);
%toplot
xplot=x(mp1{i}(x)>=ccthreshold);
yplot=y(xplot);
 
chr_sim_marginals(bT_newpsf_pombe{1}, bT_newpsf_pombe{3}, xplot, yplot,minRunLength, gapTolerance, extendLength, tol)
%% chr 3 mapping into 2
nexttile

i=2;
x=5481:7733;
y=mpI1{i}(x);
x=x-5480;
max_index=length(bT{i}.rawBarcode(bT{i}.rawBitmask))-w;
min_big=length(bT{i}.rawBarcode(bT{i}.rawBitmask))+1;
ceiling=min_big*2-w+1;
y(y>=min_big)=ceiling-y(y>=min_big);
%toplot
xplot=x(mp1{i}(x)>=ccthreshold);
yplot=y(xplot);
extend_chr_sim(xplot, yplot, 'minRunLength', minRunLength, 'gapTolerance', gapTolerance, 'extendLength', extendLength, 'tol', tol);
title('Chr2 vs Chr3');
xlabel('Chr2')
ylabel('Chr3')

%%
minRunLength = 10;   % minimum consecutive points
gapTolerance = 1;    % allowed missing points
extendLength = w;  % artificial extension length
tol = 1;             % slope tolerance
 ccthreshold=0.75;
i=2;
x=5481:7733;
y=mpI1{i}(x);
x=x-5480;
max_index=length(bT{i}.rawBarcode(bT{i}.rawBitmask))-150;
min_big=length(bT{i}.rawBarcode(bT{i}.rawBitmask))+1;
ceiling=min_big*2-150+1;
y(y>=min_big)=ceiling-y(y>=min_big);
%toplot
xplot=x(mp1{i}(x)>=ccthreshold);
yplot=y(xplot);
 
chr_sim_marginals(bT_newpsf_pombe{3}, bT_newpsf_pombe{2}, xplot, yplot,minRunLength, gapTolerance, extendLength, tol)


