disp('Running MATLAB script file PanCan2.m') ;
%
%    For analysis of TCGA PanCan data set from Chris Fan
%    12 Cancer Types
%
%    This is all about finding an example where transformation gives a benefit
%    Focus is on only OV vs. UCEC
%    Keep things simple by:
%        - deleting genes with any missings
%        - keeping only 1000 most variable genes
%        - working with random subsets of cases 
%
%    The main input file is:    kviao.pancan12tu.txt
%        Note 1st three lines have text:
%            - Long names
%            - Cancer Type
%            - Platform
%        Remaining lines 4 - 12485 are Gene Expression
%


ipart = 6 ;    %  0 - Read data, reduce, store as .mat file
               %  1 - MargDistPlots (log2, raw, autotrans)
               %  2 - PCA  (log2, raw, autotrans)
               %  3 - DWD, w/ OPC, ROC  (log2, raw, autotrans)
               %  4 - DiProPerm  (log2, raw, autotrans)
               %  5 - Random subsets, n1 = n2 = 50
               %  6 - Random subsets, n1 = n2 = 30
               %  7 - 



savefilestr1 = 'PanCan1.mat' ;
    %  file saved from PanCan1.m, starting point here
savefilestr = 'PanCan2.mat' ;



if ipart == 0 ;    %  Read data, reduce, store as .mat file

  disp(' ') ;
  disp('    Start .mat file read') ;
  load(savefilestr1) ;
  %    Loads Variables:
  %        mdata:         Matrix of Gene Expression
  %        caGeneName:    Cell Array of Gene Names
  %        caCLID         Cell Array of Case identifiers
  %        caCanType      Cell Array of Cancer Types
  %        caPlatform     Cell Array of Platforms

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

  caGeneNamebase = {'BLCA' 'COAD' 'GBM' 'KIRC' 'LAML' 'LUAD' ...
                        'LUSC' 'OV' 'READ' 'HNSC' 'UCEC' 'BRCA'} ;

  %  Sort out Subtypes
  %
  vBLCAflag = zeros(1,n) ;
  vCOADflag = zeros(1,n) ;
  vGBMflag = zeros(1,n) ;
  vKIRCflag = zeros(1,n) ;
  vLAMLflag = zeros(1,n) ;
  vLUADflag = zeros(1,n) ;
  vLUSCflag = zeros(1,n) ;
  vOVflag = zeros(1,n) ;
  vREADflag = zeros(1,n) ;
  vHNSCflag = zeros(1,n) ;
  vUCECflag = zeros(1,n) ;
  vBRCAflag = zeros(1,n) ;
  mcolor = zeros(n,3) ;
  for i = 1:n ;
    if strcmp(caCanType{i},'BLCA') ;
      vBLCAflag(i) = 1 ;
    elseif strcmp(caCanType{i},'COAD') ;
      vCOADflag(i) = 1 ;
    elseif strcmp(caCanType{i},'GBM') ;
      vGBMflag(i) = 1 ;
    elseif strcmp(caCanType{i},'KIRC') ;
      vKIRCflag(i) = 1 ;
    elseif strcmp(caCanType{i},'LAML') ;
      vLAMLflag(i) = 1 ;
    elseif strcmp(caCanType{i},'LUAD') ;
      vLUADflag(i) = 1 ;
    elseif strcmp(caCanType{i},'LUSC') ;
      vLUSCflag(i) = 1 ;
    elseif strcmp(caCanType{i},'OV') ;
      vOVflag(i) = 1 ;
    elseif strcmp(caCanType{i},'READ') ;
      vREADflag(i) = 1 ;
    elseif strcmp(caCanType{i},'HNSC') ;
      vHNSCflag(i) = 1 ;
    elseif strcmp(caCanType{i},'UCEC') ;
      vUCECflag(i) = 1 ;
    elseif strcmp(caCanType{i},'BRCA') ;
      vBRCAflag(i) = 1 ;
    else ;
      disp('    Error from PanCan2.m') ;
      disp(['    Unknown Cancer Type:  ' caCanType{i}]) ;
      pauseSM ;
    end ;
  end ;
  vBLCAflag = logical(vBLCAflag) ;
  vCOADflag = logical(vCOADflag) ;
  vGBMflag = logical(vGBMflag) ;
  vKIRCflag = logical(vKIRCflag) ;
  vLAMLflag = logical(vLAMLflag) ;
  vLUADflag = logical(vLUADflag) ;
  vLUSCflag = logical(vLUSCflag) ;
  vOVflag = logical(vOVflag) ;
  vREADflag = logical(vREADflag) ;
  vHNSCflag = logical(vHNSCflag) ;
  vUCECflag = logical(vUCECflag) ;
  vBRCAflag = logical(vBRCAflag) ;

  disp(' ') ;
  disp(['    Number of BLCA cases is:  ' num2str(sum(vBLCAflag))]) ;
  disp(['    Number of COAD cases is:  ' num2str(sum(vCOADflag))]) ;
  disp(['    Number of GBM cases is:  ' num2str(sum(vGBMflag))]) ;
  disp(['    Number of KIRC cases is:  ' num2str(sum(vKIRCflag))]) ;
  disp(['    Number of LAML cases is:  ' num2str(sum(vLAMLflag))]) ;
  disp(['    Number of LUAD cases is:  ' num2str(sum(vLUADflag))]) ;
  disp(['    Number of LUSC cases is:  ' num2str(sum(vLUSCflag))]) ;
  disp(['    Number of OV cases is:  ' num2str(sum(vOVflag))]) ;
  disp(['    Number of READ cases is:  ' num2str(sum(vREADflag))]) ;
  disp(['    Number of HNSC cases is:  ' num2str(sum(vHNSCflag))]) ;
  disp(['    Number of UCEC cases is:  ' num2str(sum(vUCECflag))]) ;
  disp(['    Number of BRCA cases is:  ' num2str(sum(vBRCAflag))]) ;

  vkeepflag = vOVflag | vUCECflag ; 
      %  one for cases that are either OV or UCEC

  mdata = mdata(:,vkeepflag) ;
  vOVflag = vOVflag(vkeepflag) ;
      %  reduce to only OV and UCEC cases

  mflagmiss = isnan(mdata) ;
  vnmissgene = sum(mflagmiss,2) ;
  vflagmissgene = (vnmissgene > 0) ;
      %  one for genes with a missing
  mdata = mdata(~vflagmissgene,:) ;
  caGeneName = caGeneName(~vflagmissgene) ;
      %  reduce to genes with no missings

  vsd = std(mdata,0,2) ;
  [temp, vind] = sort(vsd,'descend') ;
  mdata = mdata(vind(1:1000),:) ;
  caGeneName = caGeneName(vind(1:1000)) ;
      %  reduce to top 1000 SD genes

  mdataR = 2.^mdata ;
      %  Inverse of log2 transform

  disp('Starting Shifted log transformations') ;
  mdataT = AutoTransQF(mdataR) ;
  disp('Transformations finished') ;
      %  auto log transformed version

  disp(' ') ;
  disp('    Start .mat file save') ;
  save(savefilestr,'mdata','caGeneName','vOVflag','mdataR','mdataT') ;
  disp('    Finished .mat file save') ;
  disp(' ') ;


else ;    %  Work with .mat version of data

  disp(' ') ;
  disp('    Start .mat file read') ;
  load(savefilestr) ;
  %    Loads Variables:
  %        mdata:         Matrix of Gene Expression
  %        caGeneName     cell array of Gene Names 
  %        vOVflag       Indicator of OV cases (others are UCEC)
  %        mdataR         Matrix of Gene Expression (Raw coubts scale)
  %        mdataT         Matrix of Gene Expression (Auto transformed)

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

  %  Create icolor & markerstr
  %
  icolor = [] ;
  markerstr = [] ;
  for i = 1:n ;
    if vOVflag(i) ;
      icolor = [icolor; [1 0 1]] ;
      markerstr = strvcat(markerstr,'o') ;
    else ;
      icolor = [icolor; [0 0.7 0]] ;
      markerstr = strvcat(markerstr,'+') ;
    end ;
  end ;

  legendcellstr = {{'OV' 'UCEC'}} ;
  mlegendcolor = [[1 0 1]; ...
                  [0 0.7 0]] ;


  if ipart == 1 ;    %  MargDistPlots (log2, raw, autotrans)

    figure(1) ;
    clf ;
    varnamecellstr = {caGeneName} ;
    titlecellstr = {{'Summary Stats (means)' ...
                     'Raw PanCan OV & UCEC' ...
                     'Top 1000 Non-Missing Genes' ...
                     ['n = ' num2str(n) ',  d = ' num2str(d)]}} ;
    savestr = ['PanCan2ip' num2str(ipart) 'MargDistMeanRaw'] ;
    paramstruct = struct('istat',1, ...
                         'varnamecellstr',varnamecellstr, ...
                         'icolor',icolor, ...
                         'markerstr',markerstr, ...
                         'isubpopkde',1, ...
                         'datovlaymax',0.9, ...
                         'datovlaymin',0.4, ...
                         'titlecellstr',titlecellstr, ...
                         'savestr',savestr) ;
    MargDistPlotSM(mdataR,paramstruct) ;

    figure(2) ;
    clf ;
    varnamecellstr = {caGeneName} ;
    titlecellstr = {{'Summary Stats (means)' ...
                     'Log2 PanCan OV & UCEC' ...
                     'Top 1000 Non-Missing Genes' ...
                     ['n = ' num2str(n) ',  d = ' num2str(d)]}} ;
    savestr = ['PanCan2ip' num2str(ipart) 'MargDistMeanLog2'] ;
    paramstruct = struct('istat',1, ...
                         'varnamecellstr',varnamecellstr, ...
                         'icolor',icolor, ...
                         'markerstr',markerstr, ...
                         'isubpopkde',1, ...
                         'datovlaymax',0.9, ...
                         'datovlaymin',0.4, ...
                         'titlecellstr',titlecellstr, ...
                         'savestr',savestr) ;
    MargDistPlotSM(mdata,paramstruct) ;


    figure(3) ;
    clf ;
    varnamecellstr = {caGeneName} ;
    titlecellstr = {{'Summary Stats (means)' ...
                     'AutoTran PanCan OV & UCEC' ...
                     'Top 1000 Non-Missing Genes' ...
                     ['n = ' num2str(n) ',  d = ' num2str(d)]}} ;
    savestr = ['PanCan2ip' num2str(ipart) 'MargDistMeanAutoTran'] ;
    paramstruct = struct('istat',1, ...
                         'varnamecellstr',varnamecellstr, ...
                         'icolor',icolor, ...
                         'markerstr',markerstr, ...
                         'isubpopkde',1, ...
                         'datovlaymax',0.9, ...
                         'datovlaymin',0.4, ...
                         'titlecellstr',titlecellstr, ...
                         'savestr',savestr) ;
    MargDistPlotSM(mdataT,paramstruct) ;

    figure(4) ;
    clf ;
    varnamecellstr = {caGeneName} ;
    titlecellstr = {{'Summary Stats (Skewness)' ...
                     'Raw PanCan OV & UCEC' ...
                     'Top 1000 Non-Missing Genes' ...
                     ['n = ' num2str(n) ',  d = ' num2str(d)]}} ;
    savestr = ['PanCan2ip' num2str(ipart) 'MargDistSkewRaw'] ;
    paramstruct = struct('istat',1, ...
                         'varnamecellstr',varnamecellstr, ...
                         'viplot',[(1:7), ((d - 7):d)], ...
                         'icolor',icolor, ...
                         'markerstr',markerstr, ...
                         'isubpopkde',1, ...
                         'datovlaymax',0.9, ...
                         'datovlaymin',0.4, ...
                         'titlecellstr',titlecellstr, ...
                         'savestr',savestr) ;
    MargDistPlotSM(mdataR,paramstruct) ;


    figure(5) ;
    clf ;
    varnamecellstr = {caGeneName} ;
    titlecellstr = {{'Summary Stats (Skewness)' ...
                     'Log2 PanCan OV & UCEC' ...
                     'Top 1000 Non-Missing Genes' ...
                     ['n = ' num2str(n) ',  d = ' num2str(d)]}} ;
    savestr = ['PanCan2ip' num2str(ipart) 'MargDistSkewLog2'] ;
    paramstruct = struct('istat',1, ...
                         'varnamecellstr',varnamecellstr, ...
                         'viplot',[(1:7), ((d - 7):d)], ...
                         'icolor',icolor, ...
                         'markerstr',markerstr, ...
                         'isubpopkde',1, ...
                         'datovlaymax',0.9, ...
                         'datovlaymin',0.4, ...
                         'titlecellstr',titlecellstr, ...
                         'savestr',savestr) ;
    MargDistPlotSM(mdata,paramstruct) ;


    figure(6) ;
    clf ;
    varnamecellstr = {caGeneName} ;
    titlecellstr = {{'Summary Stats (Skewness)' ...
                     'AutoTran PanCan OV & UCEC' ...
                     'Top 1000 Non-Missing Genes' ...
                     ['n = ' num2str(n) ',  d = ' num2str(d)]}} ;
    savestr = ['PanCan2ip' num2str(ipart) 'MargDistSkewAutoTran'] ;
    paramstruct = struct('istat',1, ...
                         'varnamecellstr',varnamecellstr, ...
                         'viplot',[(1:7), ((d - 7):d)], ...
                         'icolor',icolor, ...
                         'markerstr',markerstr, ...
                         'isubpopkde',1, ...
                         'datovlaymax',0.9, ...
                         'datovlaymin',0.4, ...
                         'titlecellstr',titlecellstr, ...
                         'savestr',savestr) ;
    MargDistPlotSM(mdataT,paramstruct) ;


  elseif ipart == 2 ;    %  PCA  (log2, raw, autotrans)

    figure(1) ;
    clf ;
    titlecellstr = {{'PCA Scatterplot' ...
                     'Raw PanCan OV & UCEC' ...
                     'Top 1000 Non-Missing Genes' ...
                     ['n = ' num2str(n) ',  d = ' num2str(d)]}} ;
    labelcellstr = {{'PC 1 Scores' 'PC 2 Scores' 'PC 3 Scores' 'PC 4 Scores'}} ;
    savestr = ['PanCan2ip' num2str(ipart) 'PCA-Raw'] ;
    paramstruct = struct('npcadiradd',4, ...
                         'icolor',icolor, ...
                         'markerstr',markerstr, ...
                         'legendcellstr',legendcellstr, ...
                         'mlegendcolor',mlegendcolor, ...
                         'isubpopkde',1, ...
                         'titlecellstr',titlecellstr, ...
                         'labelcellstr',labelcellstr, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    scatplotSM(mdataR,[],paramstruct) ;

    figure(2) ;
    clf ;
    titlecellstr = {{'PCA Scatterplot' ...
                     'log2 PanCan OV & UCEC' ...
                     'Top 1000 Non-Missing Genes' ...
                     ['n = ' num2str(n) ',  d = ' num2str(d)]}} ;
    labelcellstr = {{'PC 1 Scores' 'PC 2 Scores' 'PC 3 Scores' 'PC 4 Scores'}} ;
    savestr = ['PanCan2ip' num2str(ipart) 'PCA-log2'] ;
    paramstruct = struct('npcadiradd',4, ...
                         'icolor',icolor, ...
                         'markerstr',markerstr, ...
                         'legendcellstr',legendcellstr, ...
                         'mlegendcolor',mlegendcolor, ...
                         'isubpopkde',1, ...
                         'titlecellstr',titlecellstr, ...
                         'labelcellstr',labelcellstr, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    scatplotSM(mdata,[],paramstruct) ;

    figure(3) ;
    clf ;
    titlecellstr = {{'PCA Scatterplot' ...
                     'AutoTran PanCan OV & UCEC' ...
                     'Top 1000 Non-Missing Genes' ...
                     ['n = ' num2str(n) ',  d = ' num2str(d)]}} ;
    labelcellstr = {{'PC 1 Scores' 'PC 2 Scores' 'PC 3 Scores' 'PC 4 Scores'}} ;
    savestr = ['PanCan2ip' num2str(ipart) 'PCA-AutoTran'] ;
    paramstruct = struct('npcadiradd',4, ...
                         'icolor',icolor, ...
                         'markerstr',markerstr, ...
                         'legendcellstr',legendcellstr, ...
                         'mlegendcolor',mlegendcolor, ...
                         'isubpopkde',1, ...
                         'titlecellstr',titlecellstr, ...
                         'labelcellstr',labelcellstr, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    scatplotSM(mdataT,[],paramstruct) ;


  elseif ipart == 3 ;    %  DWD, w/ OPC, ROC  (log2, raw, autotrans)

    %  Compute DWD directions
    %
    DWDlog2 = DWD2XQ(mdata(:,vOVflag),mdata(:,~vOVflag)) ;
    DWDraw = DWD2XQ(mdataR(:,vOVflag),mdataR(:,~vOVflag)) ;
    DWDtran = DWD2XQ(mdataT(:,vOVflag),mdataT(:,~vOVflag)) ;

    %  Make DWD scatterplot
    %
    figure(1) ;
    clf ;
    titlecellstr = {{'DWD & Ortho PCs' ...
                     'Raw PanCan OV & UCEC' ...
                     'Top 1000 Non-Missing Genes' ...
                     ['n = ' num2str(n) ',  d = ' num2str(d)]}} ;
    labelcellstr = {{'DWD OV v UCEC' 'OPC 1 Scores' 'OPC 2 Scores' 'OPC 3 Scores'}} ;
    savestr = ['PanCan2ip' num2str(ipart) 'DWDnOrthoPCraw'] ;
    paramstruct = struct('npcadiradd',-3, ...
                         'icolor',icolor, ...
                         'markerstr',markerstr, ...
                         'isubpopkde',1, ...
                         'legendcellstr',legendcellstr, ...
                         'mlegendcolor',mlegendcolor, ...
                         'titlecellstr',titlecellstr, ...
                         'labelcellstr',labelcellstr, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    scatplotSM(mdataR,DWDraw,paramstruct) ;

    figure(2) ;
    clf ;
    titlestr = ['ROC curve, Raw PanCan OV & UCEC'] ;
    savestr = ['PanCan2ip' num2str(ipart) 'ROC-AUCraw'] ;
    paramstruct = struct('titlestr',titlestr, ...
                         'xlabelstr','P{OV <= testpoint}', ...
                         'ylabelstr','P{UCEC <= testpoint}', ...
                         'savestr',savestr) ;
    ROCcurveSM(mdataR(:,vOVflag)' * DWDraw, ...
               mdataR(:,~vOVflag)' * DWDraw,paramstruct) ;

    figure(3) ;
    clf ;
    titlecellstr = {{'DWD & Ortho PCs' ...
                     'Log2 PanCan OV & UCEC' ...
                     'Top 1000 Non-Missing Genes' ...
                     ['n = ' num2str(n) ',  d = ' num2str(d)]}} ;
    labelcellstr = {{'DWD OV v UCEC' 'OPC 1 Scores' 'OPC 2 Scores' 'OPC 3 Scores'}} ;
    savestr = ['PanCan2ip' num2str(ipart) 'DWDnOrthoPClog2'] ;
    paramstruct = struct('npcadiradd',-3, ...
                         'icolor',icolor, ...
                         'markerstr',markerstr, ...
                         'isubpopkde',1, ...
                         'legendcellstr',legendcellstr, ...
                         'mlegendcolor',mlegendcolor, ...
                         'titlecellstr',titlecellstr, ...
                         'labelcellstr',labelcellstr, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    scatplotSM(mdata,DWDlog2,paramstruct) ;

    figure(4) ;
    clf ;
    titlestr = ['ROC curve, Log2 PanCan OV & UCEC'] ;
    savestr = ['PanCan2ip' num2str(ipart) 'ROC-AUClog2'] ;
    paramstruct = struct('titlestr',titlestr, ...
                         'xlabelstr','P{OV <= testpoint}', ...
                         'ylabelstr','P{UCEC <= testpoint}', ...
                         'savestr',savestr) ;
    ROCcurveSM(mdata(:,vOVflag)' * DWDlog2, ...
               mdata(:,~vOVflag)' * DWDlog2,paramstruct) ;

    figure(5) ;
    clf ;
    titlecellstr = {{'DWD & Ortho PCs' ...
                     'Autotrans PanCan OV & UCEC' ...
                     'Top 1000 Non-Missing Genes' ...
                     ['n = ' num2str(n) ',  d = ' num2str(d)]}} ;
    labelcellstr = {{'DWD OV v UCEC' 'OPC 1 Scores' 'OPC 2 Scores' 'OPC 3 Scores'}} ;
    savestr = ['PanCan2ip' num2str(ipart) 'DWDnOrthoPCtran'] ;
    paramstruct = struct('npcadiradd',-3, ...
                         'icolor',icolor, ...
                         'markerstr',markerstr, ...
                         'isubpopkde',1, ...
                         'legendcellstr',legendcellstr, ...
                         'mlegendcolor',mlegendcolor, ...
                         'titlecellstr',titlecellstr, ...
                         'labelcellstr',labelcellstr, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    scatplotSM(mdataT,DWDtran,paramstruct) ;

    figure(6) ;
    clf ;
    titlestr = ['ROC curve, AutoTrans PanCan OV & UCEC'] ;
    savestr = ['PanCan2ip' num2str(ipart) 'ROC-AUCtran'] ;
    paramstruct = struct('titlestr',titlestr, ...
                         'xlabelstr','P{OV <= testpoint}', ...
                         'ylabelstr','P{UCEC <= testpoint}', ...
                         'savestr',savestr) ;
    ROCcurveSM(mdataT(:,vOVflag)' * DWDtran, ...
               mdataT(:,~vOVflag)' * DWDtran,paramstruct) ;


  elseif ipart == 4 ;    %  DiProPerm  (log2, raw, autotrans)

    seed = 97825167 ;

    figure(1) ;
    clf ;
    savestr = ['PanCan2ip' num2str(ipart) 'DiProPermDWDraw'] ;
    title2str = ['Raw PanCan, n = ' num2str(n) ',  d = ' num2str(d)] ;
    paramstruct = struct('nsim',100, ...
                         'nreport',5, ...
                         'seed',seed, ...
                         'icolor',mlegendcolor, ...
                         'markerstr',['o'; '+'], ...
                         'title1str','DWD DiProPerm OV & UCEC', ...
                         'title2str',title2str, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    DiProPermSM(mdataR(:,vOVflag),mdataR(:,~vOVflag),paramstruct) ;

    figure(2) ;
    clf ;
    savestr = ['PanCan2ip' num2str(ipart) 'DiProPermMDraw'] ;
    title2str = ['Raw PanCan, n = ' num2str(n) ',  d = ' num2str(d)] ;
    paramstruct = struct('idir',2, ...
                         'nsim',100, ...
                         'nreport',5, ...
                         'seed',seed, ...
                         'icolor',mlegendcolor, ...
                         'markerstr',['o'; '+'], ...
                         'title1str','Mean Diff DiProPerm OV & UCEC', ...
                         'title2str',title2str, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    DiProPermSM(mdataR(:,vOVflag),mdataR(:,~vOVflag),paramstruct) ;

    figure(3) ;
    clf ;
    savestr = ['PanCan2ip' num2str(ipart) 'DiProPermDWDlog2'] ;
    title2str = ['Log2 PanCan, n = ' num2str(n) ',  d = ' num2str(d)] ;
    paramstruct = struct('nsim',100, ...
                         'nreport',5, ...
                         'seed',seed, ...
                         'icolor',mlegendcolor, ...
                         'markerstr',['o'; '+'], ...
                         'title1str','DWD DiProPerm OV & UCEC', ...
                         'title2str',title2str, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    DiProPermSM(mdata(:,vOVflag),mdata(:,~vOVflag),paramstruct) ;


    figure(4) ;
    clf ;
    savestr = ['PanCan2ip' num2str(ipart) 'DiProPermMDlog2'] ;
    title2str = ['Log2 PanCan, n = ' num2str(n) ',  d = ' num2str(d)] ;
    paramstruct = struct('idir',2, ...
                         'nsim',100, ...
                         'nreport',5, ...
                         'seed',seed, ...
                         'icolor',mlegendcolor, ...
                         'markerstr',['o'; '+'], ...
                         'title1str','Mean Diff DiProPerm OV & UCEC', ...
                         'title2str',title2str, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    DiProPermSM(mdata(:,vOVflag),mdata(:,~vOVflag),paramstruct) ;

    figure(5) ;
    clf ;
    savestr = ['PanCan2ip' num2str(ipart) 'DiProPermDWDtran'] ;
    title2str = ['Auto Trans PanCan, n = ' num2str(n) ',  d = ' num2str(d)] ;
    paramstruct = struct('nsim',100, ...
                         'nreport',5, ...
                         'seed',seed, ...
                         'icolor',mlegendcolor, ...
                         'markerstr',['o'; '+'], ...
                         'title1str','DWD DiProPerm OV & UCEC', ...
                         'title2str',title2str, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    DiProPermSM(mdataT(:,vOVflag),mdataT(:,~vOVflag),paramstruct) ;


    figure(6) ;
    clf ;
    savestr = ['PanCan2ip' num2str(ipart) 'DiProPermMDtran'] ;
    title2str = ['Auto Trans PanCan, n = ' num2str(n) ',  d = ' num2str(d)] ;
    paramstruct = struct('idir',2, ...
                         'nsim',100, ...
                         'nreport',5, ...
                         'seed',seed, ...
                         'icolor',mlegendcolor, ...
                         'markerstr',['o'; '+'], ...
                         'title1str','Mean Diff DiProPerm OV & UCEC', ...
                         'title2str',title2str, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    DiProPermSM(mdataT(:,vOVflag),mdataT(:,~vOVflag),paramstruct) ;


  elseif  ipart == 5  | ...
          ipart == 6  ;    %  Random subsets 

    if ipart == 5 ;
      nsubset = 50 ;
    elseif ipart == 6 ;
      nsubset = 30 ;
    end ;

    seed = 51096755 ;
%  Note:  forgot to set this with:    rng(seed)

    vind = 1:n ;
    vindOV = vind(vOVflag) ;
    vindUCEC = vind(~vOVflag) ;
        %  overall indices in each type
    nOV = sum(vOVflag) ;
    nUCEC = sum(~vOVflag) ;
    vindOVs = vindOV(randperm(nOV,nsubset)) ;
    vindUCECs = vindUCEC(randperm(nUCEC,nsubset)) ;
        %  random subsets of overall indices

    mdataOVs = mdata(:,vindOVs) ;
    mdataROVs = mdataR(:,vindOVs) ;
    mdataTOVs = mdataT(:,vindOVs) ;
    mdataUCECs = mdata(:,vindUCECs) ;
    mdataRUCECs = mdataR(:,vindUCECs) ;
    mdataTUCECs = mdataT(:,vindUCECs) ;

    icolors = [(ones(nsubset,1) * mlegendcolor(1,:)); ...
              (ones(nsubset,1) * mlegendcolor(2,:))] ;

    markerstrs = [] ;
    markerstrs2 = [] ;
    for i = 1:nsubset ;
      markerstrs = strvcat(markerstrs,'o') ;
      markerstrs2 = strvcat(markerstrs2,'+') ;
    end ;
    markerstrs = strvcat(markerstrs,markerstrs2) ;


    %  Make PCA scatterplots
    %
    figure(1) ;
    clf ;
    titlecellstr = {{'PCA Scatterplot' ...
                     'Raw PanCan OV & UCEC' ...
                     'Top 1000 Non-Missing Genes' ...
                     ['n1 = n2 = ' num2str(nsubset) ',  d = ' num2str(d)]}} ;
    labelcellstr = {{'PC 1 Scores' 'PC 2 Scores' 'PC 3 Scores' 'PC 4 Scores'}} ;
    savestr = ['PanCan2ip' num2str(ipart) 'PCA-Raw-ns' num2str(nsubset)] ;
    paramstruct = struct('npcadiradd',4, ...
                         'icolor',icolors, ...
                         'markerstr',markerstrs, ...
                         'legendcellstr',legendcellstr, ...
                         'mlegendcolor',mlegendcolor, ...
                         'isubpopkde',1, ...
                         'titlecellstr',titlecellstr, ...
                         'labelcellstr',labelcellstr, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    scatplotSM([mdataROVs mdataRUCECs],[],paramstruct) ;

    figure(2) ;
    clf ;
    titlecellstr = {{'PCA Scatterplot' ...
                     'Log2 PanCan OV & UCEC' ...
                     'Top 1000 Non-Missing Genes' ...
                     ['n1 = n2 = ' num2str(nsubset) ',  d = ' num2str(d)]}} ;
    labelcellstr = {{'PC 1 Scores' 'PC 2 Scores' 'PC 3 Scores' 'PC 4 Scores'}} ;
    savestr = ['PanCan2ip' num2str(ipart) 'PCA-Log2-ns' num2str(nsubset)] ;
    paramstruct = struct('npcadiradd',4, ...
                         'icolor',icolors, ...
                         'markerstr',markerstrs, ...
                         'legendcellstr',legendcellstr, ...
                         'mlegendcolor',mlegendcolor, ...
                         'isubpopkde',1, ...
                         'titlecellstr',titlecellstr, ...
                         'labelcellstr',labelcellstr, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    scatplotSM([mdataOVs mdataUCECs],[],paramstruct) ;

    figure(3) ;
    clf ;
    titlecellstr = {{'PCA Scatterplot' ...
                     'Auto Trans PanCan OV & UCEC' ...
                     'Top 1000 Non-Missing Genes' ...
                     ['n1 = n2 = ' num2str(nsubset) ',  d = ' num2str(d)]}} ;
    labelcellstr = {{'PC 1 Scores' 'PC 2 Scores' 'PC 3 Scores' 'PC 4 Scores'}} ;
    savestr = ['PanCan2ip' num2str(ipart) 'PCA-Trans-ns' num2str(nsubset)] ;
    paramstruct = struct('npcadiradd',4, ...
                         'icolor',icolors, ...
                         'markerstr',markerstrs, ...
                         'legendcellstr',legendcellstr, ...
                         'mlegendcolor',mlegendcolor, ...
                         'isubpopkde',1, ...
                         'titlecellstr',titlecellstr, ...
                         'labelcellstr',labelcellstr, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    scatplotSM([mdataTOVs mdataTUCECs],[],paramstruct) ;


    %  Do DWD & Ottho PCA
    %
    %  Compute DWD directions
    %
    DWDlog2 = DWD2XQ(mdataOVs,mdataUCECs) ;
    DWDraw = DWD2XQ(mdataROVs,mdataRUCECs) ;
    DWDtran = DWD2XQ(mdataTOVs,mdataTUCECs) ;

    %  Make DWD scatterplot
    %
    figure(4) ;
    clf ;
    titlecellstr = {{'DWD & Ortho PCs' ...
                     'Raw PanCan OV & UCEC' ...
                     'Top 1000 Non-Missing Genes' ...
                     ['n1 = n2 = ' num2str(nsubset) ',  d = ' num2str(d)]}} ;
    labelcellstr = {{'DWD OV v UCEC' 'OPC 1 Scores' 'OPC 2 Scores' 'OPC 3 Scores'}} ;
    savestr = ['PanCan2ip' num2str(ipart) 'DWDnOrthoPCraw-ns' num2str(nsubset)] ;
    paramstruct = struct('npcadiradd',-3, ...
                         'icolor',icolors, ...
                         'markerstr',markerstrs, ...
                         'isubpopkde',1, ...
                         'legendcellstr',legendcellstr, ...
                         'mlegendcolor',mlegendcolor, ...
                         'titlecellstr',titlecellstr, ...
                         'labelcellstr',labelcellstr, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    scatplotSM([mdataROVs mdataRUCECs],DWDraw,paramstruct) ;

    figure(5) ;
    clf ;
    titlecellstr = {{'DWD & Ortho PCs' ...
                     'Log2 PanCan OV & UCEC' ...
                     'Top 1000 Non-Missing Genes' ...
                     ['n1 = n2 = ' num2str(nsubset) ',  d = ' num2str(d)]}} ;
    labelcellstr = {{'DWD OV v UCEC' 'OPC 1 Scores' 'OPC 2 Scores' 'OPC 3 Scores'}} ;
    savestr = ['PanCan2ip' num2str(ipart) 'DWDnOrthoPClog2-ns' num2str(nsubset)] ;
    paramstruct = struct('npcadiradd',-3, ...
                         'icolor',icolors, ...
                         'markerstr',markerstrs, ...
                         'isubpopkde',1, ...
                         'legendcellstr',legendcellstr, ...
                         'mlegendcolor',mlegendcolor, ...
                         'titlecellstr',titlecellstr, ...
                         'labelcellstr',labelcellstr, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    scatplotSM([mdataOVs mdataUCECs],DWDlog2,paramstruct) ;

    figure(6) ;
    clf ;
    titlecellstr = {{'DWD & Ortho PCs' ...
                     'Auto Trans PanCan OV & UCEC' ...
                     'Top 1000 Non-Missing Genes' ...
                     ['n1 = n2 = ' num2str(nsubset) ',  d = ' num2str(d)]}} ;
    labelcellstr = {{'DWD OV v UCEC' 'OPC 1 Scores' 'OPC 2 Scores' 'OPC 3 Scores'}} ;
    savestr = ['PanCan2ip' num2str(ipart) 'DWDnOrthoPCtran-ns' num2str(nsubset)] ;
    paramstruct = struct('npcadiradd',-3, ...
                         'icolor',icolors, ...
                         'markerstr',markerstrs, ...
                         'isubpopkde',1, ...
                         'legendcellstr',legendcellstr, ...
                         'mlegendcolor',mlegendcolor, ...
                         'titlecellstr',titlecellstr, ...
                         'labelcellstr',labelcellstr, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    scatplotSM([mdataTOVs mdataTUCECs],DWDtran,paramstruct) ;


    %  Do ROC analysis
    %
    figure(7) ;
    clf ;
    titlestr = ['ROC curve, Raw PanCan OV & UCEC'] ;
    savestr = ['PanCan2ip' num2str(ipart) 'ROC-AUCraw-ns' num2str(nsubset)] ;
    paramstruct = struct('titlestr',titlestr, ...
                         'xlabelstr','P{OV <= testpoint}', ...
                         'ylabelstr','P{UCEC <= testpoint}', ...
                         'savestr',savestr) ;
    ROCcurveSM(mdataROVs' * DWDraw, mdataRUCECs' * DWDraw, paramstruct) ;

    figure(8) ;
    clf ;
    titlestr = ['ROC curve, Log2 PanCan OV & UCEC'] ;
    savestr = ['PanCan2ip' num2str(ipart) 'ROC-AUClog2-ns' num2str(nsubset)] ;
    paramstruct = struct('titlestr',titlestr, ...
                         'xlabelstr','P{OV <= testpoint}', ...
                         'ylabelstr','P{UCEC <= testpoint}', ...
                         'savestr',savestr) ;
    ROCcurveSM(mdataOVs' * DWDlog2, mdataUCECs' * DWDlog2, paramstruct) ;

    figure(9) ;
    clf ;
    titlestr = ['ROC curve, Auto Trans PanCan OV & UCEC'] ;
    savestr = ['PanCan2ip' num2str(ipart) 'ROC-AUCtrans-ns' num2str(nsubset)] ;
    paramstruct = struct('titlestr',titlestr, ...
                         'xlabelstr','P{OV <= testpoint}', ...
                         'ylabelstr','P{UCEC <= testpoint}', ...
                         'savestr',savestr) ;
    ROCcurveSM(mdataTOVs' * DWDtran, mdataTUCECs' * DWDtran, paramstruct) ;


    %  Do DiProPerms
    %
    figure(10) ;
    clf ;
    savestr = ['PanCan2ip' num2str(ipart) 'DiProPermDWDraw-ns' num2str(nsubset)] ;
    title2str = ['Raw PanCan, n1 = n2 = ' num2str(nsubset)  ',  d = ' num2str(d)] ;
    paramstruct = struct('nsim',100, ...
                         'nreport',5, ...
                         'seed',seed, ...
                         'icolor',mlegendcolor, ...
                         'markerstr',['o'; '+'], ...
                         'title1str','DWD DiProPerm OV & UCEC', ...
                         'title2str',title2str, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    DiProPermSM(mdataROVs,mdataRUCECs,paramstruct) ;

    figure(11) ;
    clf ;
    savestr = ['PanCan2ip' num2str(ipart) 'DiProPermMDraw-ns' num2str(nsubset)] ;
    title2str = ['Raw PanCan, n1 = n2 = ' num2str(nsubset)  ',  d = ' num2str(d)] ;
    paramstruct = struct('idir',2, ...
                         'nsim',100, ...
                         'nreport',5, ...
                         'seed',seed, ...
                         'icolor',mlegendcolor, ...
                         'markerstr',['o'; '+'], ...
                         'title1str','Mean Diff DiProPerm OV & UCEC', ...
                         'title2str',title2str, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    DiProPermSM(mdataROVs,mdataRUCECs,paramstruct) ;

    figure(12) ;
    clf ;
    savestr = ['PanCan2ip' num2str(ipart) 'DiProPermDWDlog2-ns' num2str(nsubset)] ;
    title2str = ['Log2 PanCan, n1 = n2 = ' num2str(nsubset)  ',  d = ' num2str(d)] ;
    paramstruct = struct('nsim',100, ...
                         'nreport',5, ...
                         'seed',seed, ...
                         'icolor',mlegendcolor, ...
                         'markerstr',['o'; '+'], ...
                         'title1str','DWD DiProPerm OV & UCEC', ...
                         'title2str',title2str, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    DiProPermSM(mdataOVs,mdataUCECs,paramstruct) ;

    figure(13) ;
    clf ;
    savestr = ['PanCan2ip' num2str(ipart) 'DiProPermMDlog2-ns' num2str(nsubset)] ;
    title2str = ['Log2 PanCan, n1 = n2 = ' num2str(nsubset)  ',  d = ' num2str(d)] ;
    paramstruct = struct('idir',2, ...
                         'nsim',100, ...
                         'nreport',5, ...
                         'seed',seed, ...
                         'icolor',mlegendcolor, ...
                         'markerstr',['o'; '+'], ...
                         'title1str','Mean Diff DiProPerm OV & UCEC', ...
                         'title2str',title2str, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    DiProPermSM(mdataOVs,mdataUCECs,paramstruct) ;

    figure(14) ;
    clf ;
    savestr = ['PanCan2ip' num2str(ipart) 'DiProPermDWDtran-ns' num2str(nsubset)] ;
    title2str = ['Auto Trans PanCan, n1 = n2 = ' num2str(nsubset)  ',  d = ' num2str(d)] ;
    paramstruct = struct('nsim',100, ...
                         'nreport',5, ...
                         'seed',seed, ...
                         'icolor',mlegendcolor, ...
                         'markerstr',['o'; '+'], ...
                         'title1str','DWD DiProPerm OV & UCEC', ...
                         'title2str',title2str, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    DiProPermSM(mdataTOVs,mdataTUCECs,paramstruct) ;

    figure(15) ;
    clf ;
    savestr = ['PanCan2ip' num2str(ipart) 'DiProPermMDtran-ns' num2str(nsubset)] ;
    title2str = ['AutoTrans PanCan, n1 = n2 = ' num2str(nsubset)  ',  d = ' num2str(d)] ;
    paramstruct = struct('idir',2, ...
                         'nsim',100, ...
                         'nreport',5, ...
                         'seed',seed, ...
                         'icolor',mlegendcolor, ...
                         'markerstr',['o'; '+'], ...
                         'title1str','Mean Diff DiProPerm OV & UCEC', ...
                         'title2str',title2str, ...
                         'savestr',savestr, ...
                         'iscreenwrite',1) ;
    DiProPermSM(mdataTOVs,mdataTUCECs,paramstruct) ;

    %  Do Mean Difference ROC
    %
    MDlog2 = mean(mdataOVs,2) - mean(mdataUCECs,2) ;
    MDraw = mean(mdataROVs,2) - mean(mdataRUCECs,2) ;
    MDtran = mean(mdataTOVs,2) - mean(mdataTUCECs,2) ;
    MDlog2 = MDlog2 / norm(MDlog2) ;
    MDraw = MDraw / norm(MDraw) ;
    MDtran = MDtran / norm(MDtran) ;

    figure(16) ;
    clf ;
    titlestr = ['MD ROC curve, Raw PanCan OV & UCEC'] ;
    savestr = ['PanCan2ip' num2str(ipart) 'MDROC-AUCraw-ns' num2str(nsubset)] ;
    paramstruct = struct('titlestr',titlestr, ...
                         'xlabelstr','P{OV <= testpoint}', ...
                         'ylabelstr','P{UCEC <= testpoint}', ...
                         'savestr',savestr) ;
    ROCcurveSM(mdataROVs' * MDraw, mdataRUCECs' * MDraw, paramstruct) ;

    figure(17) ;
    clf ;
    titlestr = ['MD ROC curve, Log2 PanCan OV & UCEC'] ;
    savestr = ['PanCan2ip' num2str(ipart) 'MDROC-AUClog2-ns' num2str(nsubset)] ;
    paramstruct = struct('titlestr',titlestr, ...
                         'xlabelstr','P{OV <= testpoint}', ...
                         'ylabelstr','P{UCEC <= testpoint}', ...
                         'savestr',savestr) ;
    ROCcurveSM(mdataOVs' * MDlog2, mdataUCECs' * MDlog2, paramstruct) ;

    figure(18) ;
    clf ;
    titlestr = ['MD ROC curve, Autotrans PanCan OV & UCEC'] ;
    savestr = ['PanCan2ip' num2str(ipart) 'MDROC-AUCtran-ns' num2str(nsubset)] ;
    paramstruct = struct('titlestr',titlestr, ...
                         'xlabelstr','P{OV <= testpoint}', ...
                         'ylabelstr','P{UCEC <= testpoint}', ...
                         'savestr',savestr) ;
    ROCcurveSM(mdataTOVs' * MDtran, mdataTUCECs' * MDtran, paramstruct) ;


  elseif ipart == 6 ;    %  


  end ;    %  of inner ipart loop


end ;    %  of outer ipart loop



