% Metabolomics_toolbox is required. (Available at
% https://github.com/artedison/Edison_Lab_Shared_Metabolomics_UGA)
% PLS-toolbox is required for OPLS-DA (version 881, Eigenvector Research,
% Inc., Manson, WA).
%% %%%%%%%%%%%%%%% Spectra processing %%%%%%%%%%%%%%%
%% Load Runorder file
load('RunOrders.mat');
%% Load spectra into workspace
spectra_cL=Load1D('Sheep_heart_cpmg_SZ_L_03172020','bruker');%cpmg,L
spectra_cL([1:4,6,8,9,16,20:22,31:33,40,45,46])=[]; %remove buffer
spectra_cL=spectra_cL(order_cL(:,1));
  spectra_cR=Load1D('Sheep_heart_cpmg_SZ_R_06222020','bruker');%cpmg R
  spectra_cR([1:13,15,17,27,36,37,44:46,48,51])=[];   %R; remove bad ones and buffer
  spectra_cR=spectra_cR(order_cR);
 spectra_cS=Load1D('Sheep_heart_cpmg_SZ_S_06302020','bruker');%cpmg S
 spectra_cS([1:3,5,6,11,20,30,38])=[];   %S; remove bad ones and buffer
 spectra_cS=spectra_cS(order_cS);
%% Combine 3 cpmgs
spectra_cAll=[spectra_cL,spectra_cR,spectra_cS];

%% Y vector
[GroupCode,~,y]=unique(RunOrder.group);

y=repmat(y,3,1);
y_run=repmat([1,2,3],29,1);
y_run=y_run(:);
%% Set up your X matrix and visualize the spectra to make sure they look correct

[X_cAll,ppm_cAll,XTitles_cAll]=Setup1D(spectra_cAll);

figure, plotr(ppm_cAll,X_cAll)
set(gca,'XDir','rev')
%% Remove water region
     
    %XR_cL=remove_region(X_cL,ppm_cL,4.817,5.211); %cpmg,L
    %XR_cS=remove_region(X_cS,ppm_cS,4.787,5.149);%cpmg,S
    %XR_cR=remove_region(X_cR,ppm_cR,4.8,5.15);%cpmg,R
    XR_cAll=remove_region(X_cAll,ppm_cAll,4.786,5.211);% cpmg,all
    Uplim=max(find(ppm_cAll<-0.5));Downlim=min(find(ppm_cAll>10));
    XR_cAll(:,[1:Uplim,Downlim:end])=[];
    ppmR_cAll=ppm_cAll;
    ppmR_cAll(:,[1:Uplim,Downlim:end])=[];
%% baseline correction-visualization
% [A,S]=showBaseline(X_cAll, ppm_cAll, 200, 3);%1000,300 for XR_All
% %% baseline correction
% XB_All=CorrectBl(X_All, A, S);
% %check with all spectra
% figure;plotr(ppm_All,XB_All)
% hold on
% plot(ppm_All,zeros(1,size(XB_All,2)),'r')
% hold off
% %%
% XNew_All=remove_region(XB_All,ppm_All,4.79,5.15);% noesy,all
%     Uplim=max(find(ppm_All<-0.5));Downlim=min(find(ppm_All>10));
%     XNew_All(:,[1:Uplim,Downlim:end])=[];
%     ppmR_All=ppm_All;
%     ppmR_All(:,[1:Uplim,Downlim:end])=[];
%% Alignment of peaks
tic
XAL1_1=star_align1D(XR_cAll,ppmR_cAll,'mean','RAFFT');
XAL2_1=star_align1D(XR_cAll,ppmR_cAll,'median','RAFFT');
XAL3_1=star_align1D(XR_cAll,ppmR_cAll,'mean','PAFFT');% for XNew_All
XAL4_1=star_align1D(XR_cAll,ppmR_cAll,'median','PAFFT');%cpmgR
XAL5_1=star_align1D(XR_cAll,ppmR_cAll,'mean','CCOW'); % for noesyL,cpmgL, noesyAll_rerun
XAL6_1=star_align1D(XR_cAll,ppmR_cAll,'median','CCOW');%for noesyR, cpmgS, cpmgAll, noesyAll
toc
figure;
ax1=subplot(3,2,1);
plotr(ppmR_cAll,XAL1_1);
ax2=subplot(3,2,2);
plotr(ppmR_cAll,XAL2_1);
ax3=subplot(3,2,3);
plotr(ppmR_cAll,XAL3_1);
ax4=subplot(3,2,4);
plotr(ppmR_cAll,XAL4_1);
ax5=subplot(3,2,5);
plotr(ppmR_cAll,XAL5_1);
ax6=subplot(3,2,6);
plotr(ppmR_cAll,XAL6_1);
linkaxes([ax1,ax2,ax3,ax4,ax5,ax6]);
%% Check normalization 
% XALrm_cAll=remove_region(XAL5_1,ppmR_All,-0.05,0.05); % remove DSS peak
% normcheck(XALrm_cAll)
%reference:
% Craig, A., et al., Scaling and normalization effects in NMR spectroscopic metabonomic 
%data sets. Anal Chem, 2006. 78(7): p. 2262-2267.
%% Normalization by PQN
%XALN_w_cS=XAL6_1./repmat(w_S,1,length(ppmR_cS));
XALN_cAll=normalize(XAL6_1,ppmR_cAll,'PQN');
normcheck(XALN_cAll) % good on noesyAll,cpmgS

figure, plotr(ppmR_cAll,XALN_cAll)
 %% Remove outlier
XALN_cAllo=XALN_cAll;
XALN_cAllo(47,:)=[];
y_o=y([1:46,48:end],:);
%% Spectra colored by group, along with the median spectrum for each group
figure
hold on
p1=plot(ppmR_cAll,mean(XALN_cAll(y==1,:)),'b','LineWidth',1.5);
plot(ppmR_cAll,XALN_cAll(y==1,:),'b');
p2=plot(ppmR_cAll,mean(XALN_cAll(y==2,:)),'r','LineWidth',1.5);
plot(ppmR_cAll,XALN_cAll(y==2,:),'r');
p3=plot(ppmR_cAll,mean(XALN_cAll(y==3,:)),'g','LineWidth',1.5);
plot(ppmR_cAll,XALN_cAll(y==3,:),'g');
p4=plot(ppmR_cAll,mean(XALN_cAll(y==4,:)),'m','LineWidth',1.5);
plot(ppmR_cAll,XALN_cAll(y==4,:),'m');
hold off
legend([p1 p2 p3 p4],{'CORT','DCA control','DCA cort','control'})
set(gca,'XDir','rev')
%% %%%%%%%%%%%%%%% Binning %%%%%%%%%%%%%%%
%% Peak-picking
 %% Optimize Peak Picking (threshold; for representative spectrum)
        matrix = XALN_cAll;
        ppm = ppmR_cAll;
             optimize_Peakpick1D(matrix,ppm,'max',0.67:0.01:0.7,'Complex');   
             doToAllOpenFigs('set(gca,''xlim'',[0.7,1.7])')    
             
             %% Do the actual Peak Picking (for representative spectrum)
             peaks = struct();       
             [peaks.ints, peaks.shifts]= Peakpick1D(matrix,ppm,'max',0.67,'Complex');
 
               %%
     figure,hold    
     plotr(ppm,matrix, 'k','LineWidth',1.5)
     plotr(peaks(1).shifts,peaks(1).ints, 'co')   
      clearvars matrix ppm 
%% Automatic Binning (Bucketing)
%% Generate buckets using a range of both params
        matrix = XALN_cAll;
        ppm = ppmR_cAll;
        
    sb = 0.004:0.001:0.006;
    sl = 0.3:0.1:0.5;
  
    [optOB_out] = optimize_optBucket(matrix,ppm,sb,sl);
     
    %% Filter out the bins with no peaks
       [optOB_out] = filterBuckets_Peaks_opt(ppm,optOB_out, peaks);        
    %% Plot the results of optOB
         [optOB_out] =  plotOptBucket_optResult(matrix,ppm,optOB_out,[3.6584    4.0], [min(matrix(:)) max( matrix(:,(ppm>3.6584  & ppm<4.0)), [],'all' ) ]);
         
           %% First peak picked every spectra independently, so that each spectrum contributes to the peak shape between the boundaries
        peakthresh = 0.15;
        mode = 'Complex';
              [itpeaks] = Peakpick1D_per_spectra(matrix,ppm,peakthresh,mode);  
              %%  Expand buckets to each of the bins boundaries
    [optOB_out] = expandBucketBounds (optOB_out, matrix, ppm, 'plotResult');   

%%  Manual Refinement of the boundaries            
 
    %% [optOB_out] = refineBuckets(matrix,ppm,optOB_out,5);   
        
     [optOB_out] = refineBuckets_GG2(matrix,ppm,optOB_out,itpeaks,peaks, 'expandedBuckets');    

     %%
     %save('post_buckets_3Aug2020.mat')
%%
%% This calculates the number of peaks in each bin in the original Peakpick1D_per_spectra data
        % Calculates the max peak within each bin and its chemical shift 
        % gap fills both ppm and intensities of peaks that were no present
        % in the bin or not picked (not detected/below baseline)

 buckets = optOB_out.refinedBuckets.refinedBuckets;
 ppeak_struct = itpeaks;
        
[itpeaks] = peaks_per_bin (matrix, ppm,buckets, ppeak_struct);
%% Align metric
 buckets = itpeaks.max_iteration.sorted_bucket_list;
 ppeak_struct = itpeaks;
 metadata = [];
 
itpeaks = align_metric ( matrix, ppm, ppeak_struct, buckets, metadata);
%% %%%%%%%%%%%%%%% Statistical analysis %%%%%%%%%%%%%%%
%% On full-resolution spectra
%% Check the scaling
varcheck(XALN_cAllo)
%% scaling
%XALNrm_All=remove_region(XALN_All,ppmR_All,-0.05,0.05); % remove DSS peak
XALSN1_1=scale(XALN_cAllo,'pareto'); %cpmg_All_rerun
varcheck(XALSN1_1)
 XALSN2_1=scale(XALN_cAllo,'log',0.05);
 varcheck(XALSN2_1)
XALSN3_1=scale(XALN_cAllo,'logoff'); 
varcheck(XALSN3_1)
%% Principle Component Analysis- lets look at 5 components

PCA_1=nipalsPCA(XALSN1_1,5);
% Score plot of the PCA
%figure,VisScoresLabel(XALSN1_1,PCA_1,[1 2],y);
VisScores(XALN_cAllo,PCA_1,[1 2],'Y',y_o,'conf_ellipse',1,'showlabels',1);
%% Visualize the loadings plot associated with the first 2 componants of the PCA

VisLoadings1D(XALN_cAllo,PCA_1.loadings(1,:),ppmR_cAll);
title('1st');
VisLoadings1D(XALN_cAllo,PCA_1.loadings(2,:),ppmR_cAll);
title('2nd');

%% Partial Least Squares analysis

PLS_1=plsCV(y_o(y_o==4|y_o==3,:),XALSN1_1(y_o==4|y_o==3,:),7,'da',10);
% Score plot of the PLS
VisScores(XALSN1_1(y_o==4|y_o==3,:),PLS_1,[1 2],'Y',y_o(y_o==4|y_o==3,:),'conf_ellipse',1,'showlabels',true);
%% Loadings plot of the first componant of the PLS
VisLoadings1D(XALN_cAllo,PLS_1.loadings(1,:),ppmR_cAll);
title('1st loading')
VisLoadings1D(XALN_cAllo,PLS_1.loadings(2,:),ppmR_cAll);
title('2nd loading')
%% Combine CORT and control groups for PLS
y_comb=y_o;
y_comb(y_comb==4)=1;
y_comb(y_comb==3)=2;
%%
PLS_1=plsCV(y_comb,XALSN1_1,7,'da',10);
%(y_o==1|y_o==3,:)
% Score plot of the PLS
%figure,VisScoresLabel(XALSN1_1([1:27,29:46,48:end],:),PLS_1,[1 2],y([1:27,29:46,48:end],:));
VisScores(XALSN1_1,PLS_1,[1 2],'Y',y_comb,'conf_ellipse',1,'showlabels',true);
%% Loadings plot of the first componant of the PLS
VisLoadings1D(XALN_cAllo,PLS_1.loadings(1,:),ppmR_cAll);
title('1st loading')
VisLoadings1D(XALN_cAllo,PLS_1.loadings(2,:),ppmR_cAll);
title('2nd loading')
%% OPLS-DA
X_sub12=XALSN1_1(y_o==2|y_o==1,:);
Y_sub12=y_o(y_o==2|y_o==1,:);
%%
%figure,VisLoadings1D(XALN_cAllo,plsmodel_24_4.loads{2,1}(:,1),ppmR_cAll);
VisLoadings1D(XALN_cAllo,plsmodel_12_1.wts(:,1),ppmR_cAll);
%% STOCSY
STOCSY_RMB_threshold(XALN_cAll,ppmR_cAll,0);
 %% T-test on the entire spectrum
 [H,p,sigPvals] = ttestColorSpectra(XALN_cAll,ppmR_cAll,y,[1 3],'',0.05);
%%
ttestshow(ppmR_cAll,XALN_cAll(y==2,:),XALN_cAll(y==4,:),0.05)

%% Make sub-matrix from bucketed matrix
XALSN_cbin=scale(itpeaks.max_iteration.gap_f_ints,'range');
Xbinc_All=itpeaks.max_iteration.gap_f_ints;
Xbinc_24=itpeaks.max_iteration.gap_f_ints(y==4|y==2,:);
Xbinc_12=itpeaks.max_iteration.gap_f_ints(y==1|y==2,:);
Xbinc_34=itpeaks.max_iteration.gap_f_ints(y==4|y==3,:);
Xbinc_13=itpeaks.max_iteration.gap_f_ints(y==1|y==3,:);
Xbinc_14=itpeaks.max_iteration.gap_f_ints(y==1|y==4,:);
Xbinc_23=itpeaks.max_iteration.gap_f_ints(y==2|y==3,:);
%% Combine CORT and control groups for PLS
y_comb=y;
y_comb(y_comb==4)=1;
y_comb(y_comb==3)=2;
Ysub_24=y(y==4|y==2,:);
Ysub_13=y(y==1|y==3,:);
Ysub_12=y(y==1|y==2,:);
Ysub_34=y(y==4|y==3,:);
Ysub_14=y(y==4|y==1,:);
Ysub_23=y(y==2|y==3,:);
%% OSC-PLS-DA by PLS-toolbox (version 881, Eigenvector Research,Inc., Manson, WA). %%%

%% Check VIP score
%vip_scores= vip(plsmodel_c13bin); %extract vip scores from model, need pls_toolbox
vip_scores_1=vip_scores_c12bin(:,1); %Choose component that samples separate on (Component 1)
[VIP_scores_sorted1,vip_index]=sort(vip_scores_1,'descend'); %Sort VIP scores
% %Use chemical shifts for Y-axis

%Titles_Metabolites_VIP=num2str(itpeaks.max_iteration.gap_f_mean_ppm(vip_index),4);
Titles_metabolites_ppm=string(num2str(itpeaks.max_iteration.gap_f_mean_ppm,4));
Titles_metabolites_ano=strings(length(vip_scores_1),1);
Titles_metabolites_ano(matchPPMs(Table_Annotation.ano_ppm,itpeaks.max_iteration.gap_f_mean_ppm),:)=Table_Annotation.annotation;
Titles_metabolites=strcat(Titles_metabolites_ano,Titles_metabolites_ppm);
Titles_Metabolites_VIP=Titles_metabolites(vip_index);
Yaxis=[1:30]; %Make Yaxis
figure,hold;
plot(VIP_scores_sorted1(1:30),Yaxis,'o');
yticks(Yaxis)
yticklabels(Titles_Metabolites_VIP(1:30,:))
set( gca, 'YGrid', 'on' );
set(gca,'YDir','rev')

% %Use chemical names for Y-axis

% Titles_metabolites=string(num2str(itpeaks.max_iteration.gap_f_mean_ppm,4));
% Titles_metabolites(matchPPMs(Table_Annotation.ano_ppm,itpeaks.max_iteration.gap_f_mean_ppm),:)=Table_Annotation.annotation;
% Titles_Metabolites_VIP=Titles_metabolites(vip_index); %Sort titles in same order as VIP scores
% [uTitles_Metabolites_VIP,ind_u,~]=unique(Titles_Metabolites_VIP,'stable');
% % Remember to replace with least overlapped buckets!
% % How many metabolites on the plot?
% Yaxis=[1:15]; %Make Yaxis
% figure,hold;
% plot(VIP_scores_sorted1(ind_u(1:15)),Yaxis,'o');
% yticks(Yaxis)
% yticklabels(uTitles_Metabolites_VIP(1:15,:))
% set( gca, 'YGrid', 'on' );
% set(gca,'YDir','rev')
%% PLot VIP scores for DCA control vs. control
vip_scores_1=vip_scores_c24bin(:,1); %Choose component that samples separate on (Component 1)
[VIP_scores_sorted1,vip_index]=sort(vip_scores_1,'descend');
Titles_Metabolites_VIP=Titles_metabolites(vip_index);
Ind_c24bin=[1:3,5,6,8,10,12:15,17:19,22]';
Yaxis=[1:15]; %Make Yaxis
Color_c24=mean(Xbinc_24(Ysub_24==2,:))-mean(Xbinc_24(Ysub_24==4,:)); %DCA control - control
Color_c24=Color_c24(vip_index);
Color_c24=Color_c24(Ind_c24bin);

VIP_c24_sorted=VIP_scores_sorted1(Ind_c24bin);
Titles_Metabolites_VIP_c24bin=Titles_Metabolites_VIP(Ind_c24bin,:);
figure,hold on
plot(VIP_c24_sorted(Color_c24>0),Yaxis(Color_c24>0),'r.','MarkerSize',16);
plot(VIP_c24_sorted(Color_c24<0),Yaxis(Color_c24<0),'g.','MarkerSize',16);
yticks(Yaxis)
yticklabels(Titles_Metabolites_VIP_c24bin)
set( gca, 'YGrid', 'on','YDir','rev' );
%% Write fold-change table for DCA control vs. control
Tidx_c24bin=vip_index(Ind_c24bin);
FC_c24bin=mean(itpeaks.max_iteration.gap_f_ints(y==2,Tidx_c24bin))./mean(itpeaks.max_iteration.gap_f_ints(y==4,Tidx_c24bin));
[~,P_c24bin]=ttest2(itpeaks.max_iteration.gap_f_ints(y==2,:),itpeaks.max_iteration.gap_f_ints(y==4,:));
PFDR_c24bin=mafdr(P_c24bin,'BHFDR',true);
DCActrl_cTable=table(Titles_metabolites_ano(Tidx_c24bin),FC_c24bin',P_c24bin(Tidx_c24bin)',PFDR_c24bin(Tidx_c24bin)');

%% PLot VIP scores for DCA CORT vs. CORT
vip_scores_1=vip_scores_c13bin(:,1); %Choose component that samples separate on (Component 1)
[VIP_scores_sorted1,vip_index]=sort(vip_scores_1,'descend');
Titles_Metabolites_VIP=Titles_metabolites(vip_index);
Ind_c13bin=[1:6,8,9,11,13:17,19]';
Titles_Metabolites_VIP_c13bin=Titles_Metabolites_VIP(Ind_c13bin,:);
Color_c13=mean(Xbinc_13(Ysub_13==3,:))-mean(Xbinc_13(Ysub_13==1,:)); %DCA CORT - CORT
Color_c13=Color_c13(vip_index);
Color_c13=Color_c13(Ind_c13bin);
VIP_c13_sorted=VIP_scores_sorted1(Ind_c13bin);
Yaxis=[1:15];
figure,hold;
plot(VIP_c13_sorted(Color_c13>0),Yaxis(Color_c13>0),'r.','MarkerSize',16);
plot(VIP_c13_sorted(Color_c13<0),Yaxis(Color_c13<0),'g.','MarkerSize',16);
yticks(Yaxis)
yticklabels(Titles_Metabolites_VIP_c13bin)
set( gca, 'YGrid', 'on' ,'YDir','rev');
hold off
%% Write fold-change table for DCA CORT vs. CORT
Tidx_c13bin=vip_index(Ind_c13bin);
FC_c13bin=mean(itpeaks.max_iteration.gap_f_ints(y==3,Tidx_c13bin))./mean(itpeaks.max_iteration.gap_f_ints(y==1,Tidx_c13bin));
[~,P_c13bin]=ttest2(itpeaks.max_iteration.gap_f_ints(y==3,:),itpeaks.max_iteration.gap_f_ints(y==1,:));
PFDR_c13bin=mafdr(P_c13bin,'BHFDR',true);
DCACORT_cTable=table(Titles_metabolites_ano(Tidx_c13bin),FC_c13bin',P_c13bin(Tidx_c13bin)',PFDR_c13bin(Tidx_c13bin)');

%% PLot VIP scores for DCA  vs. non-DCA
vip_scores_1=vip_scores_cCombBin(:,1); %Choose component that samples separate on (Component 1)
[VIP_scores_sorted1,vip_index]=sort(vip_scores_1,'descend');
Titles_Metabolites_VIP=Titles_metabolites(vip_index);
Ind_cCombBin=[1:3,5:9,11:14,16:18]';
Color_cCombBin=mean(Xbinc_All(y_comb==2,:))-mean(Xbinc_All(y_comb==1,:)); %DCA  - non-DCA
Color_cCombBin=Color_cCombBin(vip_index);
Color_cCombBin=Color_cCombBin(Ind_cCombBin);
VIP_cCombBin_sorted=VIP_scores_sorted1(Ind_cCombBin);
Titles_Metabolites_VIP_cCombBin=Titles_Metabolites_VIP(Ind_cCombBin,:);
Yaxis=[1:15];
figure,hold;
plot(VIP_cCombBin_sorted(Color_cCombBin>0),Yaxis(Color_cCombBin>0),'r.','MarkerSize',16);
plot(VIP_cCombBin_sorted(Color_cCombBin<0),Yaxis(Color_cCombBin<0),'g.','MarkerSize',16);
yticks(Yaxis)
yticklabels(Titles_Metabolites_VIP_cCombBin)
set( gca, 'YGrid', 'on' ,'YDir','rev');
%% Write fold-change table for DCA vs. non-DCA
Tidx_cCombBin=vip_index(Ind_cCombBin);
FC_cCombBin=mean(itpeaks.max_iteration.gap_f_ints(y_comb==2,Tidx_cCombBin))./mean(itpeaks.max_iteration.gap_f_ints(y_comb==1,Tidx_cCombBin));
[~,P_cCombBin]=ttest2(itpeaks.max_iteration.gap_f_ints(y_comb==2,:),itpeaks.max_iteration.gap_f_ints(y_comb==1,:));
PFDR_cCombBin=mafdr(P_cCombBin,'BHFDR',true);
DCA_cTable=table(Titles_metabolites_ano(Tidx_cCombBin),FC_cCombBin',PFDR_cCombBin(Tidx_cCombBin)');
%% PLot VIP scores for DCA Control vs. CORT
vip_scores_1=vip_scores_c12bin(:,1); %Choose component that samples separate on (Component 1)
[VIP_scores_sorted1,vip_index]=sort(vip_scores_1,'descend');
Titles_Metabolites_VIP=Titles_metabolites(vip_index);
Ind_c12bin=[2,3,5,7:15,18,21,23]';
Titles_Metabolites_VIP_c12bin=Titles_Metabolites_VIP(Ind_c12bin,:);
Color_c12=mean(Xbinc_12(Ysub_12==2,:))-mean(Xbinc_12(Ysub_12==1,:)); %DCA control - CORT
Color_c12=Color_c12(vip_index);
Color_c12=Color_c12(Ind_c12bin);
VIP_c12_sorted=VIP_scores_sorted1(Ind_c12bin);
Yaxis=[1:15];
figure,hold;
plot(VIP_c12_sorted(Color_c12>0),Yaxis(Color_c12>0),'r.','MarkerSize',16);
plot(VIP_c12_sorted(Color_c12<0),Yaxis(Color_c12<0),'g.','MarkerSize',16);
yticks(Yaxis)
yticklabels(Titles_Metabolites_VIP_c12bin)
set( gca, 'YGrid', 'on' ,'YDir','rev');
%% PLot VIP scores for DCA CORT vs. control
vip_scores_1=vip_scores_c34bin(:,1); %Choose component that samples separate on (Component 1)
[VIP_scores_sorted1,vip_index]=sort(vip_scores_1,'descend');
Titles_Metabolites_VIP=Titles_metabolites(vip_index);
Ind_c34bin=[1,3,5,6,8,10:12,14:17,19:21]';
Titles_Metabolites_VIP_c34bin=Titles_Metabolites_VIP(Ind_c34bin,:);
Color_c34=mean(Xbinc_34(Ysub_34==3,:))-mean(Xbinc_34(Ysub_34==4,:)); %DCA cort - control
Color_c34=Color_c34(vip_index);
Color_c34=Color_c34(Ind_c34bin);
VIP_c34_sorted=VIP_scores_sorted1(Ind_c34bin);
Yaxis=[1:15];
figure,hold;
plot(VIP_c34_sorted(Color_c34>0),Yaxis(Color_c34>0),'r.','MarkerSize',16);
plot(VIP_c34_sorted(Color_c34<0),Yaxis(Color_c34<0),'g.','MarkerSize',16);
yticks(Yaxis)
yticklabels(Titles_Metabolites_VIP_c34bin)
set( gca, 'YGrid', 'on' ,'YDir','rev');
%% STOCSY
STOCSY_RMB_threshold(XALN_cAll,ppmR_cAll,0);
%% STOCSY for Physiological parameters
% MAP
PInfo_MAP_g=[];
  for ii=1:size(PInfo.MAP,1) 
        Add_peak_g(ii,:)= create_pseudo_peak(0,2,PInfo.MAP(ii)); 
  end
    PInfo_MAP_g = [PInfo_MAP_g Add_peak_g(:,1:end-1)];
PInfo_MAP_g=repmat(PInfo_MAP_g,3,1);
Joint_MAP=[XALN_cAll,PInfo_MAP_g./20];
Joint_MAP=Joint_MAP(repmat(logical(PInfo.MAP),3,1),:);
%figure,plotr([ppmR_cAll,10+(1:size(PInfo_MAP_g,2))./1000],Joint_MAP)
STOCSY_RMB_threshold(Joint_MAP,[ppmR_cAll,10+(1:size(PInfo_MAP_g,2))./1000],0.9);
% Pearson correlation with bins
corr(repmat(PInfo.MAP(PInfo.MAP~=0),3,1),itpeaks.max_iteration.gap_f_ints(repmat(PInfo.MAP~=0,3,1),:));
% Spearman correlation
STOCSY_spearman(10.1,Joint_MAP,[ppmR_cAll,10+(1:size(PInfo_MAP_g,2))./1000]);
%%
MAP=repmat(PInfo.MAP,3,1);
p=polyfit(MAP(MAP~=0),itpeaks.max_iteration.gap_f_ints(MAP~=0,38),1);
y=polyval(p,MAP(MAP~=0));
figure,hold on
plot(MAP(MAP~=0),y)
plot(MAP(MAP~=0),itpeaks.max_iteration.gap_f_ints(MAP~=0,38),'*')
text(45,0.08,['R=',num2str(corr(MAP(MAP~=0),itpeaks.max_iteration.gap_f_ints(MAP~=0,38)),4)],'FontSize',12)
xlabel('MAP')
ylabel('Alanine')
%%
% HR
PInfo_HR_g=[];
  for ii=1:size(PInfo.HR,1) 
        Add_peak_g(ii,:)= create_pseudo_peak(0,2,PInfo.HR(ii)); 
  end
    PInfo_HR_g = [PInfo_HR_g Add_peak_g(:,1:end-1)];
PInfo_HR_g=repmat(PInfo_HR_g,3,1);
Joint_HR=[XALN_cAll,PInfo_HR_g./50];
Joint_HR=Joint_HR(repmat(logical(PInfo.HR),3,1),:);
figure,plotr([ppmR_cAll,10+(1:size(PInfo_HR_g,2))./1000],Joint_HR)
%STOCSY_RMB_threshold(Joint_HR,[ppmR_cAll,10+(1:size(PInfo_HR_g,2))./1000],0.9);

% QRSAvg46
PInfo_QRS1_g=[];
  for ii=1:size(PInfo.QRSAvg46,1) 
        Add_peak_g(ii,:)= create_pseudo_peak(0,2,PInfo.QRSAvg46(ii)); 
  end
    PInfo_QRS1_g = [PInfo_QRS1_g Add_peak_g(:,1:end-1)];
PInfo_QRS1_g=repmat(PInfo_QRS1_g,3,1);
Joint_QRS1=[XALN_cAll,PInfo_QRS1_g./50];
Joint_QRS1=Joint_QRS1(repmat(logical(PInfo.QRSAvg46),3,1),:);
figure,plotr([ppmR_cAll,10+(1:size(PInfo_QRS1_g,2))./1000],Joint_QRS1)
%STOCSY_RMB_threshold(Joint_QRS1,[ppmR_cAll,10+(1:size(PInfo_QRS1_g,2))./1000],0.9);

% QRSavglast3
PInfo_QRS2_g=[];
  for ii=1:size(PInfo.QRSavglast3,1) 
        Add_peak_g(ii,:)= create_pseudo_peak(0,2,PInfo.QRSavglast3(ii)); 
  end
    PInfo_QRS2_g = [PInfo_QRS2_g Add_peak_g(:,1:end-1)];
PInfo_QRS2_g=repmat(PInfo_QRS2_g,3,1);
Joint_QRS2=[XALN_cAll,PInfo_QRS2_g./50];
Joint_QRS2=Joint_QRS2(repmat(logical(PInfo.QRSavglast3),3,1),:);
figure,plotr([ppmR_cAll,10+(1:size(PInfo_QRS2_g,2))./1000],Joint_QRS2)
%STOCSY_RMB_threshold(Joint_QRS2,[ppmR_cAll,10+(1:size(PInfo_QRS2_g,2))./1000],0);
%% %%%%%%%%%%%%%% Annotation %%%%%%%%
% Set up
spectra_eDCtrl=Load1D('SheepHeart_before2D_SZ_08202020','bruker');
spectra_eDCtrl=spectra_eDCtrl(3); %1D (before 2D) spectrum of DCA_control extraction
[X_eDCtrl,ppm_eDCtrl,XTitles_eDCtrl]=Setup1D(spectra_eDCtrl);
figure, plotr(ppm_eDCtrl,X_eDCtrl)
set(gca,'XDir','rev')
Uplim=max(find(ppm_eDCtrl<-0.5));Downlim=min(find(ppm_eDCtrl>10));
XR_eDCtrl=X_eDCtrl;
XR_eDCtrl(:,[1:Uplim,Downlim:end])=[];
ppmR_eDCtrl=ppm_eDCtrl;
ppmR_eDCtrl(:,[1:Uplim,Downlim:end])=[];
%% plot figure with one representative HRMAS spectrum
% Align extraction spectrum with all HRMAS is bad (distort spectra)
figure,plotr(ppmR_eDCtrl,XR_eDCtrl*1e-8)
hold on
plotr(ppmR_cAll,XALN_cAll([13,44],:))
legend('Extraction','HRMAS\_cpmg\_1401L','HRMAS\_cpmg\_3022R')
