disp('Running MATLAB script file LobularGeneExpSigClust.m') ;
%
%    Explores SigClust Analysis of RNAseq Gene Expression 
%    Based on Lobular Freeze Breast Cancer Data
%    
%    Data in file TCGA_Lobular2015.mat
%    in Research/BioInf/GeneArray/JIVE/Qing-Lobular
%
%    These data are part of the TCGA project, and were first studied in:
%
%Ciriello, G., Gatza, M. L., Beck, A. H., Wilkerson, M. D., Rhie, S. K., Pastore, A., ... & Bowlby, R. (2015). Comprehensive molecular portraits of invasive lobular breast cancer. Cell, 163(2), 506-519.
%
%    and 
%
%HU, Y.-J., SUN, W., TZENG, J.-Y. & PEROU, C. M. (2015). Proper use of allele-specific expression improves
%statistical power for cis-eqtl mapping with rna-seq data. Journal of the American Statistical Association 110,
%962–974.
%
%The precise version here was used in:
%
%Feng, Q., Hannig, J., & Marron, J. S. (2015). Non-iterative Joint and Individual Variation Explained. arXiv %preprint arXiv:1512.04060.
%
%  This uses lines from Lobular1JIVE.m
%

ipart = 12 ;    %  1 - PCA scatterplot matrix of Full Data
                %  2 - Analyze LumAs alone
                %  3 - Analyze LumBs alone
                %  4 - Analyze Basals alone
                %  5 - Analyze Her2s alone
                %  6 - Analyze Basals vs rest split
                %  7 - Analyze Basals vs Her2 split
                %  8 - Analyze Her2 vs Luminals split
                %  9 - Analyze LumA vs LumB split
                %  10 - Analyze Her2 vs LumA split
                %  11 - Analyze Her2 vs LumB split
                %  12 - Analyze Basals vs rest (subsetted) split
                %  13 - Analyze Basals (subsetted) vs Her2 split
                %  14 - Analyze Her2 vs Luminals (subsetted) split
                %  15 - Analyze LumA (subsetted) vs LumB split
                %  16 - Analyze Her2 vs LumA (subsetted) split
                %  17 - Analyze Her2 vs LumB (subsetted) split


%  Read data, and set parameters
%
load ../TCGA_Lobular2015
    % Reads many things.  Will use:
    %        gemat
    %        subtype



%  Create Color and Symbol Matrices
%      will skip normal likes
%
mdata = [] ;
mcolor = [] ;
markerstr = '.' ;
stflag = [] ;
for i = 1:size(gemat,2) ;
  if strcmp(subtype{i},'LumA') ;
    stflag = [stflag 1] ;
    mdata = [mdata gemat(:,i)] ;
%    mcolor = [mcolor; [1 0 0]] ;
    mcolor = [mcolor; [0 0 1]] ;    %   Fixed 2020
    markerstr = strvcat(markerstr,'+') ;
  elseif strcmp(subtype{i},'LumB') ;
    stflag = [stflag 2] ;
    mdata = [mdata gemat(:,i)] ;
%    mcolor = [mcolor; [1 0 1]] ;
    mcolor = [mcolor; [0 1 1]] ;    %   Fixed 2020
    markerstr = strvcat(markerstr,'x') ;
  elseif strcmp(subtype{i},'Basal') ;
    stflag = [stflag 3] ;
    mdata = [mdata gemat(:,i)] ;
%    mcolor = [mcolor; [0 0 1]] ;
    mcolor = [mcolor; [1 0 0]] ;    %   Fixed 2020
    markerstr = strvcat(markerstr,'<') ;
  elseif strcmp(subtype{i},'Her2') ;
    stflag = [stflag 4] ;
    mdata = [mdata gemat(:,i)] ;
%    mcolor = [mcolor; [0 1 1]] ;
    mcolor = [mcolor; [1 0 1]] ;    %   Fixed 2020
    markerstr = strvcat(markerstr,'*') ;
  elseif strcmp(subtype{i},'Normal-like') ;
          %  Do nothing in these cases
%    mcolor = [mcolor; [0 1 0]] ;
%    markerstr = strvcat(markerstr,'o') ;
  else ;
    disp('Error from Lobular1JIVE.m: unknown subtype') ;
    disp(subtype{i}) ;
    disp('Terminating Execution') ;  
  end ;
end ;
markerstr = markerstr(2:end) ;

n = size(mdata,2) ;
d = size(mdata,1) ;

disp(' ') ;
disp('Finished Data Read & Pre-processing') ;
disp(' ') ;


if ipart == 1 ;    %  PCA scatterplot matrix of Full Data

  figure(1) ;
  clf ;
  legendcellstr = {{'LumA' 'LumB' 'Basal' 'Her2'}} ;
  mlegendcolor = [[1 0 0]; ...
                  [1 0 1]; ...
                  [0 0 1]; ...
                  [0 1 1]] ;
  titlecellstr = {{'Overall PCA' 'Lobular Freeze GE Data' ...
                    ['n = ' num2str(n) ', d = ' num2str(d)]}} ;
  labelcellstr = {{'PC1 scores'; 'PC2 scores'; 'PC3 scores'; 'PC4 scores'}} ;
  savestr = ['LGESigCLustIp' num2str(ipart) 'AllPCA4x4'] ;
  paramstruct = struct('npcadiradd',4, ...
                       'icolor',mcolor, ...
                       'markerstr',markerstr, ...
                       'isubpopkde',1, ...
                       'datovlaymax',0.9, ...
                       'datovlaymin',0.1, ...
                       'legendcellstr',legendcellstr, ...
                       'mlegendcolor',mlegendcolor, ...
                       'titlecellstr',titlecellstr, ...
                       'labelcellstr',labelcellstr, ...
                       'savestr',savestr, ...
                       'iscreenwrite',1) ;
  scatplotSM(mdata,[],paramstruct) ;

  figure(2) ;
  clf ;
  savestr = ['LGESigCLustIp' num2str(ipart) 'AllPCA3x3'] ;
  paramstruct = struct('npcadiradd',3, ...
                       'icolor',mcolor, ...
                       'markerstr',markerstr, ...
                       'isubpopkde',1, ...
                       'datovlaymax',0.9, ...
                       'datovlaymin',0.1, ...
                       'legendcellstr',legendcellstr, ...
                       'mlegendcolor',mlegendcolor, ...
                       'titlecellstr',titlecellstr, ...
                       'labelcellstr',labelcellstr, ...
                       'savestr',savestr, ...
                       'iscreenwrite',1) ;
  scatplotSM(mdata,[],paramstruct) ;

  figure(3) ;
  clf ;
  savestr = ['LGESigCLustIp' num2str(ipart) 'AllPCA2x2'] ;
  paramstruct = struct('npcadiradd',2, ...
                       'icolor',mcolor, ...
                       'markerstr',markerstr, ...
                       'isubpopkde',1, ...
                       'datovlaymax',0.9, ...
                       'datovlaymin',0.1, ...
                       'legendcellstr',legendcellstr, ...
                       'mlegendcolor',mlegendcolor, ...
                       'titlecellstr',titlecellstr, ...
                       'labelcellstr',labelcellstr, ...
                       'savestr',savestr, ...
                       'iscreenwrite',1) ;
  scatplotSM(mdata,[],paramstruct) ;

  paramstruct = struct('npc',2, ...
                       'iscreenwrite',1, ...
                       'viout',[0 1]) ;
  outstruct = pcaSM(mdata,paramstruct) ;
  mdir = getfield(outstruct,'meigvec') ;
  savestr = ['LGESigCLustIp' num2str(ipart) 'AllPC12'] ;
  paramstruct = struct('icolor',mcolor, ...
                       'markerstr',markerstr, ...
                       'legendcellstr',legendcellstr, ...
                       'mlegendcolor',mlegendcolor, ...
                       'titlestr','Overall PCA, Lobular Freeze GE Data', ...
                       'xlabelstr','PC1 scores', ...
                       'ylabelstr','PC2 scores', ...
                       'ifigure',4, ...
                       'savestr',savestr, ...
                       'iscreenwrite',1) ;
  projplot2SM(mdata,mdir,paramstruct) ;


elseif  ipart == 2  | ...
        ipart == 3  | ...
        ipart == 4  | ...
        ipart == 5  ;    %  Analyze a subtype alone
         
  if ipart == 2 ;    %  Analyze LumA alone
    mdatas = mdata(:,(stflag == 1)) ;
    icolor = 'r' ;
    markerstr = '+' ;
    subtstr = 'LumA' ;
    
  elseif ipart == 3 ;    %  Analyze LumB alone
    mdatas = mdata(:,(stflag == 2)) ;
    icolor = 'm' ;
    markerstr = 'x' ;
    subtstr = 'LumB' ;
         
  elseif ipart == 4 ;    %  Analyze Basal alone
    mdatas = mdata(:,(stflag == 4)) ;
    icolor = 'b' ;
    markerstr = '<' ;
    subtstr = 'Basal' ;
         
  elseif ipart == 5 ;    %  Analyze Her2 alone
    mdatas = mdata(:,(stflag == 3)) ;
    icolor = 'c' ;
    markerstr = '*' ;
    subtstr = 'Her2' ;

  end ;

  %  One Class Only PCA
  %
  figure(5) ;
  clf ;
  titlecellstr = {{[subtstr ' Only PCA'] 'Lobular Freeze GE Data' ...
                    ['n = ' num2str(size(mdatas,2)) ', d = ' num2str(d)]}} ;
  labelcellstr = {{'PC1 scores'; 'PC2 scores'; 'PC3 scores'}} ;
  savestr = ['LGESigCLustIp' num2str(ipart) subtstr 'OnlyPCA'] ;
  paramstruct = struct('npcadiradd',3, ...
                       'icolor',icolor, ...
                       'markerstr',markerstr, ...
                       'datovlaymax',0.9, ...
                       'datovlaymin',0.1, ...
                       'titlecellstr',titlecellstr, ...
                       'labelcellstr',labelcellstr, ...
                       'savestr',savestr, ...
                       'iscreenwrite',1) ;
  scatplotSM(mdatas,[],paramstruct) ;

  %  One Class Only SigClust
  %  
  figure(1) ;
  clf ;
  savestr = ['LGESigCLustIp' num2str(ipart) subtstr 'Only'] ;
  paramstruct = struct('vclass',0, ...
                       'nsim',1000, ...
                       'InitRandstate',74029378, ...
                       'InitRandnstate',23948755, ...
                       'SimRandstate',47093287, ...
                       'SimRandnstate',37028734, ...
                       'datastr',['Split ' subtstr 's Only'], ...
                       'BGSDsavestr',[savestr 'BGSD'], ...
                       'CovEsavestr',[savestr 'CovE'], ...
                       'pValsavestr',[savestr 'SigClust'], ...
                       'iscreenwrite',2) ;
  SigClustSM(mdatas,paramstruct) ;


elseif  ipart == 6  | ...
        ipart == 7  | ...
        ipart == 8  | ...
        ipart == 9  | ...
        ipart == 10  | ...
        ipart == 11  ;    %  Analyze pairwise splits 


  if ipart == 6 ;    %  Analyze Basals vs rest split
    mdatas1 = mdata(:,(stflag == 3)) ;
    mdatas2 = mdata(:,(stflag == 1) | (stflag == 2) | (stflag == 4)) ;
    icolor1 = 'b' ;
    icolor1v = [0 0 1] ;
    icolor2 = 'k' ;
    icolor2v = [0 0 0] ;
    markerstr1 = '<' ;
    markerstr2 = 'o' ;
    subtstr = 'BasvRest' ;

  elseif ipart == 7 ;    %  Analyze Basals vs Her2 split
    mdatas1 = mdata(:,(stflag == 3)) ;
    mdatas2 = mdata(:,(stflag == 4)) ;
    icolor1 = 'b' ;
    icolor1v = [0 0 1] ;
    icolor2 = 'c' ;
    icolor2v = [0 1 1] ;
    markerstr1 = '<' ;
    markerstr2 = '*' ;
    subtstr = 'BasvHer2' ;

  elseif ipart == 8 ;    %  Analyze Her2 vs Luminals split
    mdatas1 = mdata(:,(stflag == 4)) ;
    mdatas2 = mdata(:,(stflag == 1) | (stflag == 2)) ;
    icolor1 = 'c' ;
    icolor1v = [0 1 1] ;
    icolor2 = 'k' ;
    icolor2v = [0 0 0] ;
    markerstr1 = '*' ;
    markerstr2 = 'o' ;
    subtstr = 'Her2vLums' ;

  elseif ipart == 9 ;    %  Analyze LumA vs LumB split
    mdatas1 = mdata(:,(stflag == 1)) ;
    mdatas2 = mdata(:,(stflag == 2)) ;
    icolor1 = 'r' ;
    icolor1v = [1 0 0] ;
    icolor2 = 'm' ;
    icolor2v = [1 0 1] ;
    markerstr1 = '+' ;
    markerstr2 = 'x' ;
    subtstr = 'LumAvLumB' ;

  elseif ipart == 10 ;    %  Analyze Her2 vs LumA
    mdatas1 = mdata(:,(stflag == 4)) ;
    mdatas2 = mdata(:,stflag == 1) ;
    icolor1 = 'c' ;
    icolor1v = [0 1 1] ;
    icolor2 = 'r' ;
    icolor2v = [1 0 0] ;
    markerstr1 = '*' ;
    markerstr2 = '+' ;
    subtstr = 'Her2vLumA' ;

  elseif ipart == 11 ;    %  Analyze Her2 vs LumB
    mdatas1 = mdata(:,(stflag == 4)) ;
    mdatas2 = mdata(:,stflag == 2) ;
    icolor1 = 'c' ;
    icolor1v = [0 1 1] ;
    icolor2 = 'm' ;
    icolor2v = [1 0 1] ;
    markerstr1 = '*' ;
    markerstr2 = 'x' ;
    subtstr = 'Her2vLumB' ;

  end ;

  mdatas = [mdatas1 mdatas2] ;
  vMD = mean(mdatas2,2) - mean(mdatas1,2) ;
  vMD = vMD / norm(vMD) ;
  n1 = size(mdatas1,2) ;
  n2 = size(mdatas2,2) ;
  vclass = [ones(1,n1) (2 * ones(1,n2))] ;
  icolor = [] ;
  markerstr = '.' ;
  for i = 1:n1 ;
    icolor = [icolor; icolor1v] ;
    markerstr = strvcat(markerstr,markerstr1) ;
  end ;
  for i = 1:n2 ;
    icolor = [icolor; icolor2v] ;
    markerstr = strvcat(markerstr,markerstr2) ;
  end ;
  markerstr = markerstr(2:end) ;

  %  Pairwise Classes PCA
  %
  figure(5) ;
  clf ;
  titlecellstr = {{[subtstr ' PCA'] 'Lobular Freeze GE Data' ...
                    ['n = ' num2str(n1 + n2) ', d = ' num2str(d)]}} ;
  labelcellstr = {{'PC1 scores'; 'PC2 scores'; 'PC3 scores'}} ;
  savestr = ['LGESigCLustIp' num2str(ipart) subtstr 'PCA'] ;
  paramstruct = struct('npcadiradd',3, ...
                       'icolor',icolor, ...
                       'markerstr',markerstr, ...
                       'isubpopkde',1, ...
                       'datovlaymax',0.9, ...
                       'datovlaymin',0.1, ...
                       'titlecellstr',titlecellstr, ...
                       'labelcellstr',labelcellstr, ...
                       'savestr',savestr, ...
                       'iscreenwrite',1) ;
  scatplotSM(mdatas,[],paramstruct) ;

  %  Pairwise Classes MD & PCA
  %
  figure(6) ;
  clf ;
  titlecellstr = {{[subtstr ' PCA'] 'Lobular Freeze GE Data' ...
                    ['n = ' num2str(n1 + n2) ', d = ' num2str(d)]}} ;
  labelcellstr = {{[subtstr ' MD']; 'OPC1 scores'}} ;
  savestr = ['LGESigCLustIp' num2str(ipart) subtstr 'MDnOPC'] ;
  paramstruct = struct('npcadiradd',-1, ...
                       'icolor',icolor, ...
                       'markerstr',markerstr, ...
                       'isubpopkde',1, ...
                       'datovlaymax',0.9, ...
                       'datovlaymin',0.1, ...
                       'titlecellstr',titlecellstr, ...
                       'labelcellstr',labelcellstr, ...
                       'savestr',savestr, ...
                       'iscreenwrite',1) ;
  scatplotSM(mdatas,vMD,paramstruct) ;

  %  Pairwise Classes SigClust
  %  
  figure(1) ;
  clf ;
  savestr = ['LGESigCLustIp' num2str(ipart) subtstr 'Test'] ;
  paramstruct = struct('vclass',vclass, ...
                       'nsim',1000, ...
                       'InitRandstate',74029378, ...
                       'InitRandnstate',23948755, ...
                       'SimRandstate',47093287, ...
                       'SimRandnstate',37028734, ...
                       'datastr',['Test ' subtstr], ...
                       'BGSDsavestr',[savestr 'BGSD'], ...
                       'CovEsavestr',[savestr 'CovE'], ...
                       'pValsavestr',[savestr 'SigClust'], ...
                       'iscreenwrite',2) ;
  SigClustSM(mdatas,paramstruct) ;


elseif  ipart == 12  | ...
        ipart == 13  | ...
        ipart == 14  | ...
        ipart == 15  | ...
        ipart == 16  | ...
        ipart == 17  ;    %  Analyze pairwise splits, with subsetting
  


  if ipart == 12 ;    %  Analyze Basals vs rest split
    mdatas1 = mdata(:,(stflag == 3)) ;
    mdatas2 = mdata(:,(stflag == 1) | (stflag == 2) | (stflag == 4)) ;
    icolor1 = 'b' ;
    icolor1v = [0 0 1] ;
    icolor2 = 'k' ;
    icolor2v = [0 0 0] ;
    markerstr1 = '<' ;
    markerstr2 = 'o' ;
    subtstr = 'BasvRest' ;

  elseif ipart == 13 ;    %  Analyze Her2 vs Basals split
    mdatas1 = mdata(:,(stflag == 4)) ;
    mdatas2 = mdata(:,(stflag == 3)) ;
    icolor1 = 'c' ;
    icolor1v = [0 1 1] ;
    icolor2 = 'b' ;
    icolor2v = [0 0 1] ;
    markerstr2 = '*' ;
    markerstr1 = '<' ;
    subtstr = 'Her2vBas' ;

  elseif ipart == 14 ;    %  Analyze Her2 vs Luminals split
    mdatas1 = mdata(:,(stflag == 4)) ;
    mdatas2 = mdata(:,(stflag == 1) | (stflag == 2)) ;
    icolor1 = 'c' ;
    icolor1v = [0 1 1] ;
    icolor2 = 'k' ;
    icolor2v = [0 0 0] ;
    markerstr1 = '*' ;
    markerstr2 = 'o' ;
    subtstr = 'Her2vLums' ;

  elseif ipart == 15 ;    %  Analyze LumB vs LumA split
    mdatas1 = mdata(:,(stflag == 2)) ;
    mdatas2 = mdata(:,(stflag == 1)) ;
    icolor1 = 'm' ;
    icolor1v = [1 0 1] ;
    icolor2 = 'r' ;
    icolor2v = [1 0 0] ;
    markerstr1 = 'x' ;
    markerstr2 = '+' ;
    subtstr = 'LumBvLumA' ;

  elseif ipart == 16 ;    %  Analyze Her2 vs LumA
    mdatas1 = mdata(:,(stflag == 4)) ;
    mdatas2 = mdata(:,stflag == 1) ;
    icolor1 = 'c' ;
    icolor1v = [0 1 1] ;
    icolor2 = 'r' ;
    icolor2v = [1 0 0] ;
    markerstr1 = '*' ;
    markerstr2 = '+' ;
    subtstr = 'Her2vLumA' ;

  elseif ipart == 17 ;    %  Analyze Her2 vs LumB
    mdatas1 = mdata(:,(stflag == 4)) ;
    mdatas2 = mdata(:,stflag == 2) ;
    icolor1 = 'c' ;
    icolor1v = [0 1 1] ;
    icolor2 = 'm' ;
    icolor2v = [1 0 1] ;
    markerstr1 = '*' ;
    markerstr2 = 'x' ;
    subtstr = 'Her2vLumB' ;

  end ;


  %  Set Common Quantities
  %
  figure(1) ;    %  Figure for MD & OPC1
  clf ;

  figure(2) ;    %  Figure for Sig Clust
  clf ;

  n1 = size(mdatas1,2) ;
  n2 = size(mdatas2,2) ;
      %  note many things below assume n2 > n1
  vclass = [ones(1,n1) (2 * ones(1,n1))] ;
      %  setting up for subsetting
  icolor = [] ;
  markerstr = '.' ;
  for i = 1:n1 ;
    icolor = [icolor; icolor1v] ;
    markerstr = strvcat(markerstr,markerstr1) ;
  end ;
  for i = 1:n1 ;    %  setting up for subsetting
    icolor = [icolor; icolor2v] ;
    markerstr = strvcat(markerstr,markerstr2) ;
  end ;
  markerstr = markerstr(2:end) ;
  rng(68987656) ;


  for iax = 1:4 ;   %  Loop through subplots

    vflag = randperm(n2) ;
    vflag = vflag(1:n1) ; 
    mdata2ss = mdatas2(:,vflag) ;
    mdatas = [mdatas1 mdata2ss] ;


    %  Pairwise Classes SigClust
    %  
    paramstruct = struct('vclass',vclass, ...
                         'nsim',1000, ...
                         'iBGSDdiagplot',0, ...
                         'iCovEdiagplot',0, ...
                         'ipValplot',0, ...
                         'iscreenwrite',2) ;
    [pval, zscore] = SigClustSM(mdatas,paramstruct) ;
    disp(['  For run ' num2str(iax) ...
              ', p = ' num2str(pval) ...
              ', zscore = ' num2str(zscore)]) ;


    %  Make MD projections
    %
    figure(1) ;
    subplot(2,2,iax) ;
      vMD = mean(mdata2ss,2) - mean(mdatas1,2) ;
      vMD = vMD / norm(vMD) ;
      paramstruct = struct('icolor',mcolor, ...
                           'markerstr',markerstr, ...
                           'ifigure',0, ...
                           'xlabelstr','MD scores', ...
                           'iscreenwrite',1) ;
      projplot1SM(mdatas,mdir,paramstruct) ;
      hold on ;
        vax = axis ;
        text(vax(1) + 0.9 * (vax(2) - vax(1)), ...
             vax(3) + 0.1 * (vax(4) - vax(3)), ...
             ['p-value = ' num2str(pvalue)]) ;
        text(vax(1) + 0.9 * (vax(2) - vax(1)), ...
             vax(3) + 0.7 * (vax(4) - vax(3)), ...
             ['Z-score = ' num2str(zscore)]) ;
      hold off ;


    %  Make MD & OPC1 graphic
    %
    figure(2) ;
    subplot(2,2,iax) ;
      paramstruct = struct('npc',1, ...
                           'iscreenwrite',1, ...
                           'viout',[0 1]) ;
      mdatasc = mdatas - (vMD * (vMD' * mdatas)) ;
          %  subtract projection onto vMD
      outstruct = pcaSM(mdatasc,paramstruct) ;
      mdir = [vMD getfield(outstruct,'meigvec')] ;
      paramstruct = struct('icolor',mcolor, ...
                           'markerstr',markerstr, ...
                           'ifigure',0, ...
                           'xlabelstr','MD scores', ...
                           'ylabelstr','OPC1 scores', ...
                           'iscreenwrite',1) ;
      projplot2SM(mdatas,mdir,paramstruct) ;

  
  end ;


  figure(1) ;
    orient landscape ;
    savestr = ['LGESigCLustIp' num2str(ipart) subtstr 'MDnOPC'] ;
  print('-dpsc',[savestr '.ps']) ;

  figure(2) ;
    orient landscape ;
    savestr = ['LGESigCLustIp' num2str(ipart) subtstr 'SigClust'] ;
  print('-dpsc',[savestr '.ps']) ;


end ;







