disp('Running MATLAB script file Pancreatic2.m') ;
%
%    For analysis of Jen Jen Yeh's Pancreatic Cancer Data
%        2nd Generation (with paired normal data
%
%    Does some visualizations and tests
%
%    Note:  there was an inconsistency in the Class Labels,
%       addressed this by going with the labels in the 3rd sheet
%


ipart = 133 ;    %  0 - Read data, check inputs and save as:
%         Pancreatic2.mat
%          (everything else reads this)
%  1 - PCA Plot of full data set (class symbol & tissue color)
%  2 - DWD norm vs. rest, ortho PCs
%  3 - DWD Solid organ vs. rest, ortho PCs
%  4 - DWD liver normal vs. rest, ortho PCs, 100 gene list
%  5 - DWD liver normal vs. liver tumor, full data set, 100 gene list
%  11 - pancreas only, DWD tumor vs. rest, 100 gene list
%  12 - spleen only, DWD tumor vs. rest, 100 gene list
%  13 - lung only, DWD tumor vs. rest, 100 gene list
%  14 - liver only, DWD tumor vs. rest, 100 gene list
%  15 - lymph only, DWD tumor vs. rest, 100 gene list
%  16 - fat only, DWD tumor vs. rest, 100 gene list
%  17 - muscle only, DWD tumor vs. rest, 100 gene list
%  18 - bowel only, DWD tumor vs. rest, 100 gene list
%  19 - vessel only, DWD tumor vs. rest, 100 gene list
% 100 - Do DWD merge of resected with autopsy primary
% 101 - Before and after PCA plots of merged data
% 102 - DWD norm vs. rest, ortho PCs (merged data)
% 103 - DWD Solid organ vs. rest, ortho PCs (merged data)
% 104 - DWD liver normal vs. rest, ortho PCs, 100 gene list (merged data)
% 105 - DWD liver normal vs. liver tumor, full data set, 100 gene list (merged data)
% 111 - pancreas only, DWD tumor vs. rest, 100 gene list (merged data)
% 112 - spleen only, DWD tumor vs. rest, 100 gene list (merged data)
% 113 - lung only, DWD tumor vs. rest, 100 gene list (merged data)
% 114 - liver only, DWD tumor vs. rest, 100 gene list (merged data)
% 115 - lymph only, DWD tumor vs. rest, 100 gene list (merged data)
% 116 - fat only, DWD tumor vs. rest, 100 gene list (merged data)
% 117 - muscle only, DWD tumor vs. rest, 100 gene list (merged data)
% 118 - bowel only, DWD tumor vs. rest, 100 gene list (merged data)
% 119 - vessel only, DWD tumor vs. rest, 100 gene list (merged data)
% 121 - Study Liver mets vs. Primary, 500 gene list (merged data)
% 122 - Study Lung mets vs. Primary, 500 gene list (merged data)
% 123 - Study Lymph mets vs. Primary, 500 gene list (merged data)
% 124 - Study Fat mets vs. Primary, 500 gene list (merged data)
% 125 - Study Muscle mets vs. Primary, 500 gene list (merged data)
% 131 - Study Fat&Muscle mets vs. Primary, 500 gene list (merged data)
% 132 - Study Lung&Liver mets vs. Primary, 500 gene list (merged data)
% 133 - Study Lymph&Fat&Muscle&Lung&Liver mets vs. Primary, 500 gene list (merged data)



%  Check Matlab version
%
if version('-release') < 14 ;

    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    disp('!!!   Warning from Pancreatic2.m:      !!!') ;
    disp('!!!   A release of Matlab before 14    !!!') ;
    disp('!!!   is being used.  This can give    !!!') ;
    disp('!!!   unexpected results with data     !!!') ;
    disp('!!!   writing and reading              !!!') ;
    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
    pauseSM ;

end ;



if ipart == 0 ;    %  Do main data read, check inputs, and save as .mat

    %  Read in Data Values
    %
    infilestr = 'RAP.UNC.ALL.Marron.xls' ;

    [numeric,txt] = xlsread(infilestr,'RAP.UNC.all.pcl(Imputed)') ;
    %  read 1st sheet
    mdata = numeric(4:end,:) ;
    vResected = numeric(1,:) ;
    vPrimary = numeric(2,:) ;
    vNormal = numeric(3,:) ;
    SlideNames = txt(2,3:end) ;
    Samples = txt(3,3:end) ;
    Tissues = txt(4,3:end) ;
    AgilentIDs = txt(8:end,1) ;
    GeneNames = txt(8:end,2) ;
    numeric = [] ;
    txt = [] ;
    %  to save space

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


    disp(' ') ;
    disp('Check upper left data value = -0.2') ;
    mdata(1,1)

    disp(' ') ;
    disp('Check upper right data value = 1.79') ;
    mdata(1,end)

    disp(' ') ;
    disp('Check lower left data value = -0.28') ;
    mdata(end,1)

    disp(' ') ;
    disp('Check lower right data value = 0.11') ;
    mdata(end,end)


    disp(' ') ;
    disp('Check 1st vResected = 1') ;
    vResected(1)

    disp(' ') ;
    disp('Check last vResected = 0') ;
    vResected(end)


    disp(' ') ;
    disp('Check 1st vPrimary = 1') ;
    vPrimary(1)

    disp(' ') ;
    disp('Check last vPrimary = 0') ;
    vPrimary(end)


    disp(' ') ;
    disp('Check 1st vNormal = 0') ;
    vNormal(1)

    disp(' ') ;
    disp('Check last vNormal = 0') ;
    vNormal(end)


    disp(' ') ;
    disp('Check 1st SlideName = 251485019091_1_2') ;
    SlideNames{1}

    disp(' ') ;
    disp('Check last SlideName = 251485021225_1_1') ;
    SlideNames{end}


    disp(' ') ;
    disp('Check 1st Sample = 10R1750B Normal Liver (100%)') ;
    Samples{1}

    disp(' ') ;
    disp('Check last Sample = PDUCTTREErpt(2)') ;
    Samples{end}


    disp(' ') ;
    disp('Check 1st Tissue = LIVER') ;
    Tissues{1}

    disp(' ') ;
    disp('Check last Tissue = PANCREAS') ;
    Tissues{end}


    disp(' ') ;
    disp('Check 1st AgilentID = AGI_HUM1_OLIGO_A_23_P100001') ;
    AgilentIDs{1}

    disp(' ') ;
    disp('Check last AgilentID = AGI_HUM1_OLIGO_A_32_P99942') ;
    AgilentIDs{end}


    disp(' ') ;
    disp('Check 1st GeneName = Hs.27373 ^ NM_207446 ^ LOC400451 ^ hypothetical gene supported by AK075564; BC060873') ;
    GeneNames{1}

    disp(' ') ;
    disp('Check last GeneName = ^ THC2375058 ^  ^ ') ;
    GeneNames{end}


    %  Read in and check pair connections
    %
    [numeric,txt] = xlsread(infilestr,'Sheet1') ;
    %  read 1st sheet
    %  note this sheet actually contains ALL of the data
    %  even though Excel shows only the ones with pairs...
    PLSlideNames = txt(2:end,1) ;
    PLSamples = txt(2:end,2) ;
    PLTissues = txt(2:end,3) ;
    PLmatch = numeric(:,4) ;


    disp(' ') ;
    disp('Check 1st PLSlideName = 251485019091_1_2') ;
    %  taken from full Excel file
    PLSlideNames{1}

    disp(' ') ;
    disp('Check last PLSlideName = 251485021225_1_1') ;
    PLSlideNames{end}


    disp(' ') ;
    disp('Check 1st PLSample = 10R1750B Normal Liver (100%)') ;
    PLSamples{1}

    disp(' ') ;
    disp('Check last PLSample = PDUCTTREErpt(2)') ;
    PLSamples{end}


    disp(' ') ;
    disp('Check 1st PLTissue = LIVER') ;
    PLTissues{1}

    disp(' ') ;
    disp('Check last PLTissue = PANCREAS') ;
    PLTissues{end}


    disp(' ') ;
    disp('Check 1st match = NaN') ;
    PLmatch(1)

    disp(' ') ;
    disp('Check 2nd match = 1') ;
    PLmatch(2)

    disp(' ') ;
    disp('Check 8th match = 2') ;
    PLmatch(8)

    disp(' ') ;
    disp('Check last match = NaN') ;
    PLmatch(end)


    for i = 1:n ;
        if ~strcmp(PLSlideNames{i},SlideNames{i}) ;
            disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
            disp('!!!   Error from Pancreatic2.m    !!!') ;
            disp('!!!   mismatch of paired sample:  !!!') ;
            disp(['!!!   ' PLSlideNames{i}]) ;
            disp('!!!   Terminating Execution       !!!') ;
            disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
            return ;
        end ;
        if ~strcmp(PLSamples{i},Samples{i}) ;
            disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
            disp('!!!   Error from Pancreatic2.m    !!!') ;
            disp('!!!   mismatch of paired sample:  !!!') ;
            disp(['!!!   ' PLSamples{i}]) ;
            disp('!!!   Terminating Execution       !!!') ;
            disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
            return ;
        end ;
        if ~strcmp(PLTissues{i},Tissues{i}) ;
            disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
            disp('!!!   Error from Pancreatic2.m    !!!') ;
            disp('!!!   mismatch of paired sample:  !!!') ;
            disp(['!!!   ' PLTissues{i}]) ;
            disp('!!!   Terminating Execution       !!!') ;
            disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
            return ;
        end ;
    end ;

    flagNaN = isnan(PLmatch) ;
    PLmatch(flagNaN) = zeros(sum(flagNaN),1) ;
    %  reset NaNs to 0

    %  Construct idataconn, for showing connections between cases in scatplotSM
    %
    idataconn = [] ;
    vind = 1:n ;
    for ipair = 1:max(PLmatch) ;
        flagpair = (PLmatch == ipair) ;
        if sum(flagpair) == 2 ;
            idataconn = [idataconn; vind(flagpair)] ;
        else ;
            disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
            disp('!!!   Error from Pancreatic2.m      !!!') ;
            disp(['!!!   number of pair indices = ' num2str(sum(flagpair)) '    !!!']) ;
            disp('!!!   Terminating Execution         !!!') ;
            disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
            return ;
        end ;
    end ;


    disp(' ') ;
    disp('Check 1st idataconn = [2 7]') ;
    idataconn(1,:)

    disp(' ') ;
    disp('Check last idataconn = [164 165]') ;
    idataconn(end,:)



    %  Save results in .mat file
    %
    save pancreatic2.mat mdata vResected vPrimary vNormal SlideNames Samples ...
        Tissues AgilentIDs GeneNames idataconn  ;
    %  save these variables

    disp('Finished data save') ;




else ;    %  then do other parts of the program


    load pancreatic2 ;
    %  loads:
    %      mdata
    %      vResected
    %      vPrimary
    %      vNormal
    %      SlideNames
    %      Samples
    %      Tissues
    %      AgilentIDs
    %      GeneNames
    %      idataconn
    disp('  Finished Data Load') ;

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



    %  Set up generic legend
    %
    legendcellstr = {{'PANCREAS' 'SPLEEN' 'LUNG' 'LIVER' 'LYMPH NODE' ...
        'FAT' 'MUSCLE' 'BOWEL' 'VESSEL'}} ;
    mlegendcolor = [[0.6 0 0]; [1 0 0]; [1 0.8 0]; [0.8 1 0]; [0 1 0]; ...
        [0 1 0.8]; [0 0.8 1]; [0 0 1]; [0 0 0.6]] ;


    %  Create generic icolor & markerstr
    %
    icolor = [] ;
    markerstr = [] ;
    AutoMetsFlag = zeros(1,n) ;
    AutoPrimFlag = zeros(1,n) ;
    ResePrimFlag = zeros(1,n) ;
    NormalFlag = zeros(1,n) ;
    vTissueNum = zeros(1,n) ;

    for i = 1:n ;

        if vNormal(i) == 0 ;    %  then have a normal
            NormalFlag(i) = 1 ;
            AutoMetsFlag(i) = 0 ;
            AutoPrimFlag(i) = 0 ;
            ResePrimFlag(i) = 0 ;
            markerstr = strvcat(markerstr,'o') ;
        elseif vNormal(i) == 1 ;
            NormalFlag(i) = 0 ;
            if vResected(i) == 0 ;    %  then have a resected
                AutoMetsFlag(i) = 0 ;
                AutoPrimFlag(i) = 0 ;
                ResePrimFlag(i) = 1 ;
                markerstr = strvcat(markerstr,'+') ;
            elseif vResected(i) == 1 ;
                ResePrimFlag(i) = 0 ;
                if vPrimary(i) == 0 ;    %  then have a primary
                    AutoMetsFlag(i) = 0 ;
                    AutoPrimFlag(i) = 1 ;
                    markerstr = strvcat(markerstr,'*') ;
                elseif vPrimary(i) == 1 ;    %  then have a mets
                    AutoMetsFlag(i) = 1 ;
                    AutoPrimFlag(i) = 0 ;
                    markerstr = strvcat(markerstr,'x') ;
                else ;
                    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
                    disp('!!!   Error from Pancreatic2.m:   !!!') ;
                    disp('!!!   Bad value of vPrimary      !!!') ;
                    disp(['!!!   for i = ' num2str(i)]) ;
                    disp('!!!   Terminating Execution       !!!') ;
                    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
                    return ;
                end ;
            else ;
                disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
                disp('!!!   Error from Pancreatic2.m:   !!!') ;
                disp('!!!   Bad value of vResected      !!!') ;
                disp(['!!!   for i = ' num2str(i)]) ;
                disp('!!!   Terminating Execution       !!!') ;
                disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
                return ;
            end ;
        else ;
            disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
            disp('!!!   Error from Pancreatic2.m:   !!!') ;
            disp('!!!   Bad value of vNormal        !!!') ;
            disp(['!!!   for i = ' num2str(i)]) ;
            disp('!!!   Terminating Execution       !!!') ;
            disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
            return ;
        end ;


        if strcmp(Tissues{i},legendcellstr{1}{1}) ;
            icolor = [icolor; mlegendcolor(1,:)] ;
            vTissueNum(i) = 1 ;
        elseif strcmp(Tissues{i},legendcellstr{1}{2}) ;
            icolor = [icolor; mlegendcolor(2,:)] ;
            vTissueNum(i) = 2 ;
        elseif strcmp(Tissues{i},legendcellstr{1}{3}) ;
            icolor = [icolor; mlegendcolor(3,:)] ;
            vTissueNum(i) = 3 ;
        elseif strcmp(Tissues{i},legendcellstr{1}{4}) ;
            icolor = [icolor; mlegendcolor(4,:)] ;
            vTissueNum(i) = 4 ;
        elseif strcmp(Tissues{i},legendcellstr{1}{5}) ;
            icolor = [icolor; mlegendcolor(5,:)] ;
            vTissueNum(i) = 5 ;
        elseif strcmp(Tissues{i},legendcellstr{1}{6}) ;
            icolor = [icolor; mlegendcolor(6,:)] ;
            vTissueNum(i) = 6 ;
        elseif strcmp(Tissues{i},legendcellstr{1}{7}) ;
            icolor = [icolor; mlegendcolor(7,:)] ;
            vTissueNum(i) = 7 ;
        elseif strcmp(Tissues{i},legendcellstr{1}{8}) ;
            icolor = [icolor; mlegendcolor(8,:)] ;
            vTissueNum(i) = 8 ;
        elseif strcmp(Tissues{i},legendcellstr{1}{9}) ;
            icolor = [icolor; mlegendcolor(9,:)] ;
            vTissueNum(i) = 9 ;
        else ;
            disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
            disp('!!!   Error from Pancreatic2.m:   !!!') ;
            disp('!!!   Bad value of vTissue        !!!') ;
            disp(['!!!   for i = ' num2str(i)]) ;
            disp('!!!   Terminating Execution       !!!') ;
            disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
            return ;
        end ;

    end ;


    idataconncolor = icolor(idataconn(:,1),:) ;
    %  use color of left hand member of each pair

    AutoMetsFlag = logical(AutoMetsFlag) ;
    AutoPrimFlag = logical(AutoPrimFlag) ;
    ResePrimFlag = logical(ResePrimFlag) ;
    NormalFlag = logical(NormalFlag) ;


    LiverFlag = (vTissueNum == 4) ;
    PancreasFlag = (vTissueNum == 1) ;
    SolidFlag = (vTissueNum == 1) | ...
        (vTissueNum == 2) | ...
        (vTissueNum == 3) | ...
        (vTissueNum == 4) ;
    LymphFlag = (vTissueNum == 5) ;


    figure(1) ;
    clf ;



    if ipart == 1 ;    %  PCA Plot of full data set (class symbol & tissue colors)

        titlecellstr = {{['PCA, n = ' num2str(n) ', d = ' num2str(d)] 'Pancreatic data (all)' ...
            ['o - normal,   + - Resected'] ['* - AutoPrim,  x - AutoMets']}} ;

        labelcellstr = {{'PC 1' 'PC 2' 'PC 3' 'PC 4'}} ;

        savestr = ['plots/Pancreatic2ip00' num2str(ipart) 'PCAall'] ;

        paramstruct = struct('icolor',icolor, ...
            'npcadiradd',4, ...
            'markerstr',markerstr, ...
            'isubpopkde',1, ...
            'idataconn',idataconn, ...
            'idataconncolor',idataconncolor, ...
            'titlecellstr',titlecellstr, ...
            'legendcellstr',legendcellstr, ...
            'mlegendcolor',mlegendcolor, ...
            'labelcellstr',labelcellstr, ...
            'savestr',savestr, ...
            'iscreenwrite',1) ;

        scatplotSM(mdata,[],paramstruct) ;


    elseif ipart == 2 ;    %  DWD norm vs. rest, full data set (class symbol & tissue colors)

        mdir = DWD1SM(mdata(:,~NormalFlag),mdata(:,NormalFlag)) ;

        titlecellstr = {{'DWD Normal vs. Rest' 'Pancreatic data (all)' ...
            ['o - normal,   + - Resected'] ['* - AutoPrim,  x - AutoMets']}} ;

        labelcellstr = {{'DWD Direction' 'Ortho PC 1' 'Ortho PC 2' 'Ortho PC 3'}} ;

        savestr = ['plots/Pancreatic2ip00' num2str(ipart) 'DWDNormalVRestOrtho'] ;

        paramstruct = struct('icolor',icolor, ...
            'npcadiradd',-3, ...
            'markerstr',markerstr, ...
            'isubpopkde',1, ...
            'idataconn',idataconn, ...
            'idataconncolor',idataconncolor, ...
            'legendcellstr',legendcellstr, ...
            'mlegendcolor',mlegendcolor, ...
            'titlecellstr',titlecellstr, ...
            'labelcellstr',labelcellstr, ...
            'savestr',savestr, ...
            'iscreenwrite',1) ;

        scatplotSM(mdata,mdir,paramstruct) ;


    elseif ipart == 3 ;    %  DWD Solid organ vs. rest (class symbol & tissue colors)

        mdir = DWD1SM(mdata(:,~SolidFlag),mdata(:,SolidFlag)) ;

        titlecellstr = {{'DWD Solid Organ vs. Rest' 'Pancreatic data (all)' ...
            ['o - normal,   + - Resected'] ['* - AutoPrim,  x - AutoMets']}} ;

        labelcellstr = {{'DWD Direction' 'Ortho PC 1' 'Ortho PC 2' 'Ortho PC 3'}} ;

        savestr = ['plots/Pancreatic2ip00' num2str(ipart) 'DWDSolidOrgVRestOrtho'] ;

        paramstruct = struct('icolor',icolor, ...
            'npcadiradd',-3, ...
            'markerstr',markerstr, ...
            'isubpopkde',1, ...
            'idataconn',idataconn, ...
            'idataconncolor',idataconncolor, ...
            'legendcellstr',legendcellstr, ...
            'mlegendcolor',mlegendcolor, ...
            'titlecellstr',titlecellstr, ...
            'labelcellstr',labelcellstr, ...
            'savestr',savestr, ...
            'iscreenwrite',1) ;

        scatplotSM(mdata,mdir,paramstruct) ;


    elseif ipart == 4 ;    %  DWD liver normal vs. rest, ortho PCs

        mdir = DWD1SM(mdata(:,~(LiverFlag & NormalFlag)),mdata(:,(LiverFlag & NormalFlag))) ;

        titlestrsub = 'DWD Normal Liver vs. Rest' ;
        titlecellstr = {{titlestrsub 'Pancreatic data (all)' ...
            ['o - normal,   + - Resected'] ['* - AutoPrim,  x - AutoMets']}} ;

        labelcellstr = {{'DWD Direction' 'Ortho PC 1' 'Ortho PC 2' 'Ortho PC 3'}} ;

        savestrsub = 'DWDNorLivVRestOrtho' ;
        savestr = ['plots/Pancreatic2ip00' num2str(ipart) savestrsub] ;

        paramstruct = struct('icolor',icolor, ...
            'npcadiradd',-3, ...
            'markerstr',markerstr, ...
            'isubpopkde',1, ...
            'idataconn',idataconn, ...
            'idataconncolor',idataconncolor, ...
            'legendcellstr',legendcellstr, ...
            'mlegendcolor',mlegendcolor, ...
            'titlecellstr',titlecellstr, ...
            'labelcellstr',labelcellstr, ...
            'savestr',savestr, ...
            'iscreenwrite',1) ;

        scatplotSM(mdata,mdir,paramstruct) ;

        %  Make output text file
        %
        ngeneout = 100 ;
        [temp,visort] = sort(abs(mdir),1,'descend') ;
        fid = fopen(['plots/Pancreatic2ip00' num2str(ipart) 'GeneList' ...
            num2str(ngeneout) savestrsub '.txt'],'wt') ;
        %  'wt' is for "delete contents of this file and open
        %               for writing" (with 't' for "text").

        titlstr1 = 'Output from the MATLAB Script Pancreatic2.m,   ' ;
        titlstr1 = [titlstr1,date] ;
        cntbytes = fprintf(fid,'%1s\n\n\n',titlstr1) ;
        %  '%1s/n' says "print string, followed by a line feed"
        %       (Note, this only goes to next line)

        titlstr2 = ['     Top ' num2str(ngeneout) ' genes (in order), for ' titlestrsub ': '] ;
        cntbytes = fprintf(fid,'%1s\n\n',titlstr2) ;
        %  3 line feeds gives two blank lines

        forstr = '%s\t%s\t%s\t%s\t%s\n' ;
        cntbytes = fprintf(fid,forstr,'GeneRank','DWDwt','PvalT','AgilentID','GeneName') ;

        PvalTs = [] ;
        for i = 1:d ;
            iout = visort(i) ;
            DWDwt = abs(mdir(iout)) ;
            [H,P] = ttest2(mdata(iout,~(LiverFlag & NormalFlag)),mdata(iout,(LiverFlag & NormalFlag)),0.025,'both') ;
            %  do two sided 2 sample t-test,
            P = P / 2 ;
            %  but divide by 2, to get test on absolute value
            if i <= ngeneout ;    %  then print to file
                forstr = '%3d\t%6.4f\t%10.8f\t%s\t%s\n' ;
                %  Embedded text appears as written (including parens)
                cntbytes = fprintf(fid,forstr,i,DWDwt,P,AgilentIDs{iout},GeneNames{iout}) ;
                PvalTs = [PvalTs; P]  ;
            end ;
        end ;

        fclose('all') ;
        %  !!!! Careful !!!!  Need to do this to be able work with files after
        %  running matlab.  Alternative is:
        %  fclose(fid) ;



    elseif ipart == 5 ;    %  DWD liver normal vs. liver tumor, full data set

        mdir = DWD1SM(mdata(:,(LiverFlag & NormalFlag)),mdata(:,(LiverFlag & (~NormalFlag)))) ;

        titlestrsub = 'DWD Liver, Normal vs. Tumor' ;
        titlecellstr = {{titlestrsub 'Pancreatic data (all)' ...
            ['o - normal,   + - Resected'] ['* - AutoPrim,  x - AutoMets']}} ;

        labelcellstr = {{'DWD Direction' 'Ortho PC 1' 'Ortho PC 2' 'Ortho PC 3'}} ;

        savestrsub = 'DWDLivNorVTumOrtho' ;
        savestr = ['plots/Pancreatic2ip00' num2str(ipart) savestrsub] ;

        paramstruct = struct('icolor',icolor, ...
            'npcadiradd',-3, ...
            'markerstr',markerstr, ...
            'isubpopkde',1, ...
            'idataconn',idataconn, ...
            'idataconncolor',idataconncolor, ...
            'legendcellstr',legendcellstr, ...
            'mlegendcolor',mlegendcolor, ...
            'titlecellstr',titlecellstr, ...
            'labelcellstr',labelcellstr, ...
            'savestr',savestr, ...
            'iscreenwrite',1) ;

        scatplotSM(mdata,mdir,paramstruct) ;

        %  Make output text file
        %
        ngeneout = 100 ;
        [temp,visort] = sort(abs(mdir),1,'descend') ;
        fid = fopen(['plots/Pancreatic2ip00' num2str(ipart) 'GeneList' ...
            num2str(ngeneout) savestrsub '.txt'],'wt') ;
        %  'wt' is for "delete contents of this file and open
        %               for writing" (with 't' for "text").

        titlstr1 = 'Output from the MATLAB Script Pancreatic2.m,   ' ;
        titlstr1 = [titlstr1,date] ;
        cntbytes = fprintf(fid,'%1s\n\n\n',titlstr1) ;
        %  '%1s/n' says "print string, followed by a line feed"
        %       (Note, this only goes to next line)

        titlstr2 = ['     Top ' num2str(ngeneout) ' genes (in order), for ' titlestrsub ': '] ;
        cntbytes = fprintf(fid,'%1s\n\n',titlstr2) ;
        %  3 line feeds gives two blank lines

        forstr = '%s\t%s\t%s\t%s\t%s\n' ;
        cntbytes = fprintf(fid,forstr,'GeneRank','DWDwt','PvalT','AgilentID','GeneName') ;

        PvalTs = [] ;
        for i = 1:d ;
            iout = visort(i) ;
            DWDwt = abs(mdir(iout)) ;
            [H,P] = ttest2(mdata(iout,(LiverFlag & NormalFlag)),mdata(iout,(LiverFlag & (~NormalFlag))),0.025,'both') ;
            %  do two sided 2 sample t-test,
            P = P / 2 ;
            %  but divide by 2, to get test on absolute value
            if i <= ngeneout ;    %  then print to file
                forstr = '%3d\t%6.4f\t%10.8f\t%s\t%s\n' ;
                %  Embedded text appears as written (including parens)
                cntbytes = fprintf(fid,forstr,i,DWDwt,P,AgilentIDs{iout},GeneNames{iout}) ;
                PvalTs = [PvalTs; P]  ;
            end ;
        end ;

        fclose('all') ;
        %  !!!! Careful !!!!  Need to do this to be able work with files after
        %  running matlab.  Alternative is:
        %  fclose(fid) ;



    elseif  ipart == 11  | ...
            ipart == 12  | ...
            ipart == 13  | ...
            ipart == 14  | ...
            ipart == 15  | ...
            ipart == 16  | ...
            ipart == 17  | ...
            ipart == 18  | ...
            ipart == 19  ;    %  individual tissues only, DWD tumor vs. rest


        if ipart == 11 ;
            TissFlag = (vTissueNum == 1) ;
            TissTitStr = legendcellstr{1}{1} ;
            TissColor = mlegendcolor(1,:) ;
        elseif ipart == 12 ;
            TissFlag = (vTissueNum == 2) ;
            TissTitStr = legendcellstr{1}{2} ;
            TissColor = mlegendcolor(2,:) ;
        elseif ipart == 13 ;
            TissFlag = (vTissueNum == 3) ;
            TissTitStr = legendcellstr{1}{3} ;
            TissColor = mlegendcolor(3,:) ;
        elseif ipart == 14 ;
            TissFlag = (vTissueNum == 4) ;
            TissTitStr = legendcellstr{1}{4} ;
            TissColor = mlegendcolor(4,:) ;
        elseif ipart == 15 ;
            TissFlag = (vTissueNum == 5) ;
            TissTitStr = legendcellstr{1}{5} ;
            TissColor = mlegendcolor(5,:) ;
        elseif ipart == 16 ;
            TissFlag = (vTissueNum == 6) ;
            TissTitStr = legendcellstr{1}{6} ;
            TissColor = mlegendcolor(6,:) ;
        elseif ipart == 17 ;
            TissFlag = (vTissueNum == 7) ;
            TissTitStr = legendcellstr{1}{7} ;
            TissColor = mlegendcolor(7,:) ;
        elseif ipart == 18 ;
            TissFlag = (vTissueNum == 8) ;
            TissTitStr = legendcellstr{1}{8} ;
            TissColor = mlegendcolor(8,:) ;
        elseif ipart == 19 ;
            TissFlag = (vTissueNum == 9) ;
            TissTitStr = legendcellstr{1}{9} ;
            TissColor = mlegendcolor(9,:) ;
        end ;

        vind = 1:n ;
        %  all indices
        vindTiss = vind(TissFlag) ;
        %  indices for this tissue type
        idataconnTiss = [] ;
        idataconncolorTiss = [] ;
        for i = 1:size(idataconn,1) ;    %  Loop through data conns
            [minv1,ind1] = min(abs(idataconn(i,1) - vindTiss)) ;
            %  minimizer, and index (among this tissue)
            if minv1 == 0 ;    %  if 1st entry corresponds to this tissue type
                [minv2,ind2] = min(abs(idataconn(i,2) - vindTiss)) ;
                %  minimizer, and index (among this tissue)
                if minv2 == 0 ;    %  if 2nd entry corresponds to this tissue type
                    idataconnTiss = [idataconnTiss; [ind1 ind2]] ;
                    idataconncolorTiss = [idataconncolorTiss; TissColor] ;
                else ;
                    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
                    disp('!!!   Warning from Pancreatic2.m:              !!!') ;
                    disp('!!!   idataconn has mismatch of tissue types   !!!') ;
                    disp(['!!!   in row ' num2str(i)]) ;
                    disp(['!!!   case numbers ' num2str(idataconn(i,1)) num2str(idataconn(i,2))]) ;
                    disp('!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!') ;
                end ;
            end ;
        end ;

        mdir = DWD1SM(mdata(:,(TissFlag & NormalFlag)),mdata(:,(TissFlag & (~NormalFlag)))) ;

        titlestrsub = ['DWD ' TissTitStr ', Normal vs. Tumor'] ;
        titlecellstr = {{titlestrsub 'Indiv. Tissue Axes' ...
            ['o - normal,   + - Resected'] ['* - AutoPrim,  x - AutoMets']}} ;

        labelcellstr = {{'DWD Direction' 'Ortho PC 1' 'Ortho PC 2' 'Ortho PC 3'}} ;

        savestrsub = ['DWDNorVTum' TissTitStr] ;
        savestr = ['plots/Pancreatic2ip0' num2str(ipart) savestrsub] ;

        paramstruct = struct('icolor',icolor(TissFlag',:), ...
            'npcadiradd',-3, ...
            'markerstr',markerstr(TissFlag'), ...
            'isubpopkde',1, ...
            'idataconn',idataconnTiss, ...
            'idataconncolor',idataconncolorTiss, ...
            'titlecellstr',titlecellstr, ...
            'labelcellstr',labelcellstr, ...
            'savestr',savestr, ...
            'iscreenwrite',1) ;

        scatplotSM(mdata(:,TissFlag),mdir,paramstruct) ;

        %  Make output text file
        %
        ngeneout = 100 ;
        [temp,visort] = sort(abs(mdir),1,'descend') ;
        fid = fopen(['plots/Pancreatic2ip0' num2str(ipart) 'GeneList' ...
            num2str(ngeneout) savestrsub '.txt'],'wt') ;
        %  'wt' is for "delete contents of this file and open
        %               for writing" (with 't' for "text").

        titlstr1 = 'Output from the MATLAB Script Pancreatic2.m,   ' ;
        titlstr1 = [titlstr1,date] ;
        cntbytes = fprintf(fid,'%1s\n\n\n',titlstr1) ;
        %  '%1s/n' says "print string, followed by a line feed"
        %       (Note, this only goes to next line)

        titlstr2 = ['     Top ' num2str(ngeneout) ' genes (in order), for ' titlestrsub ': '] ;
        cntbytes = fprintf(fid,'%1s\n\n',titlstr2) ;
        %  3 line feeds gives two blank lines

        forstr = '%s\t%s\t%s\t%s\t%s\n' ;
        cntbytes = fprintf(fid,forstr,'GeneRank','DWDwt','PvalT','AgilentID','GeneName') ;

        PvalTs = [] ;
        for i = 1:d ;
            iout = visort(i) ;
            DWDwt = abs(mdir(iout)) ;
            [H,P] = ttest2(mdata(iout,(LiverFlag & NormalFlag)),mdata(iout,(LiverFlag & (~NormalFlag))),0.025,'both') ;
            %  do two sided 2 sample t-test,
            P = P / 2 ;
            %  but divide by 2, to get test on absolute value
            if i <= ngeneout ;    %  then print to file
                forstr = '%3d\t%6.4f\t%10.8f\t%s\t%s\n' ;
                %  Embedded text appears as written (including parens)
                cntbytes = fprintf(fid,forstr,i,DWDwt,P,AgilentIDs{iout},GeneNames{iout}) ;
                PvalTs = [PvalTs; P]  ;
            end ;
        end ;

        fclose('all') ;
        %  !!!! Careful !!!!  Need to do this to be able work with files after
        %  running matlab.  Alternative is:
        %  fclose(fid) ;


    elseif ipart == 100 ;    %  Do DWD adjustment, and save .mat file

        batchlabels = (vResected == 0) - ((vResected == 1) & (vPrimary == 0)) ;
        bflag = abs(batchlabels) > 0.5 ;
        savestr = ['plots/Pancreatic2ip' num2str(ipart) 'DWDBatchAdjust'] ;
        legcellstr = {{'Resected, Pancreas' 'Autopsy, Pancreas'}} ;
        paramstruct = struct('viplot',ones(4,1), ...
            'imeantype',-1, ...
            'savestr',savestr, ...
            'legcellstr',legcellstr, ...
            'iscreenwrite',1, ...
            'npc',4) ;
        adjdata = BatchAdjustSM(mdata(:,bflag),batchlabels(:,bflag),paramstruct) ;

        mdataadj = mdata ;
        mdataadj(:,bflag) = adjdata ;


        %  Save results in .mat file
        %
        save pancreatic2adj.mat mdataadj  ;
        %  save these variables

        %  Save results in Excel Spreadsheet
        %
        outfilestr = 'pancreatic2adj.xls' ;
        [Success,Message] = xlswrite(outfilestr,mdataadj,'MainAdjData')
        [Success,Message] = xlswrite(outfilestr,GeneNames,'GeneNames')
        [Success,Message] = xlswrite(outfilestr,SlideNames,'CaseNames')

        disp('Finished data save') ;



    elseif ipart >= 101 ;    %  load adjusted .mat file

        load pancreatic2adj ;
        %  loads:
        %      mdataadj
        disp('  Finished Adjusted Data Load') ;


        if ipart == 101 ;    %  PCA plots of merged data

            %  Compute common PCA Directions
            %
            paramstruct = struct('npc',4,...
                'viout',[0 1]) ;
            outstruct = pcaSM(mdata,paramstruct) ;
            mdir = getfield(outstruct,'meigvec') ;


            %  Compute DWD direction of shift
            %
            DWDdir = DWD1SM(mdata(:,vResected == 0),mdata(:,(vResected == 1) & (vPrimary == 0))) ;
            mdataproj = mdata - DWDdir * (DWDdir' * mdata) ;
            %  projection of data onto subspace ortho to DWDdir
            paramstruct = struct('npc',3,...
                'viout',[0 1]) ;
            outstruct = pcaSM(mdataproj,paramstruct) ;
            mdirortho = getfield(outstruct,'meigvec') ;
            mdirDWD = [DWDdir mdirortho] ;


            %  Before Adjustment
            %
            titlecellstr = {{['PCA, n = ' num2str(n) ', d = ' num2str(d)] 'Before Batch Adjustment' ...
                ['o - normal,   + - Resected'] ['* - AutoPrim,  x - AutoMets']}} ;

            labelcellstr = {{'PC 1' 'PC 2' 'PC 3' 'PC 4'}} ;

            savestr = ['plots/Pancreatic2ip' num2str(ipart) 'PCAallBefore'] ;

            paramstruct = struct('icolor',icolor, ...
                'npcadiradd',0, ...
                'markerstr',markerstr, ...
                'isubpopkde',1, ...
                'idataconn',idataconn, ...
                'idataconncolor',idataconncolor, ...
                'titlecellstr',titlecellstr, ...
                'legendcellstr',legendcellstr, ...
                'mlegendcolor',mlegendcolor, ...
                'labelcellstr',labelcellstr, ...
                'savestr',savestr, ...
                'iscreenwrite',1) ;

            scatplotSM(mdata,mdir,paramstruct) ;


            titlecellstr = {{['DWD & Ortho PCA, n = ' num2str(n) ', d = ' num2str(d)] 'Before Batch Adjustment' ...
                ['o - normal,   + - Resected'] ['* - AutoPrim,  x - AutoMets']}} ;

            labelcellstr = {{'DWD Resect. vs. Auto. Prim.' 'Ortho PC 1' 'Ortho PC 2' 'Ortho PC 3'}} ;

            savestr = ['plots/Pancreatic2ip' num2str(ipart) 'DWDnOrthoPCallBefore'] ;

            paramstruct = struct('icolor',icolor, ...
                'npcadiradd',0, ...
                'markerstr',markerstr, ...
                'isubpopkde',1, ...
                'idataconn',idataconn, ...
                'idataconncolor',idataconncolor, ...
                'titlecellstr',titlecellstr, ...
                'legendcellstr',legendcellstr, ...
                'mlegendcolor',mlegendcolor, ...
                'labelcellstr',labelcellstr, ...
                'savestr',savestr, ...
                'iscreenwrite',1) ;

            scatplotSM(mdata,mdirDWD,paramstruct) ;



            %  After Adjustment
            %
            titlecellstr = {{['PCA, n = ' num2str(n) ', d = ' num2str(d)] 'After Adjustment, new Axes' ...
                ['o - normal,   + - Resected'] ['* - AutoPrim,  x - AutoMets']}} ;

            labelcellstr = {{'PC 1' 'PC 2' 'PC 3' 'PC 4'}} ;

            savestr = ['plots/Pancreatic2ip' num2str(ipart) 'PCAallAfter'] ;

            paramstruct = struct('icolor',icolor, ...
                'npcadiradd',0, ...
                'markerstr',markerstr, ...
                'isubpopkde',1, ...
                'idataconn',idataconn, ...
                'idataconncolor',idataconncolor, ...
                'titlecellstr',titlecellstr, ...
                'legendcellstr',legendcellstr, ...
                'mlegendcolor',mlegendcolor, ...
                'labelcellstr',labelcellstr, ...
                'savestr',savestr, ...
                'iscreenwrite',1) ;

            scatplotSM(mdataadj,mdir,paramstruct) ;


            titlecellstr = {{['DWD & Ortho PCA, n = ' num2str(n) ', d = ' num2str(d)] 'After Batch Adjustment' ...
                ['o - normal,   + - Resected'] ['* - AutoPrim,  x - AutoMets']}} ;

            labelcellstr = {{'DWD Resect. vs. Auto. Prim.' 'Ortho PC 1' 'Ortho PC 2' 'Ortho PC 3'}} ;

            savestr = ['plots/Pancreatic2ip' num2str(ipart) 'DWDnOrthoPCallAfter'] ;

            paramstruct = struct('icolor',icolor, ...
                'npcadiradd',0, ...
                'markerstr',markerstr, ...
                'isubpopkde',1, ...
                'idataconn',idataconn, ...
                'idataconncolor',idataconncolor, ...
                'titlecellstr',titlecellstr, ...
                'legendcellstr',legendcellstr, ...
                'mlegendcolor',mlegendcolor, ...
                'labelcellstr',labelcellstr, ...
                'savestr',savestr, ...
                'iscreenwrite',1) ;

            scatplotSM(mdataadj,mdirDWD,paramstruct) ;



            %  After Adjustment, new PCA
            %
            titlecellstr = {{['PCA, n = ' num2str(n) ', d = ' num2str(d)] 'After Adjustment, new Axes' ...
                ['o - normal,   + - Resected'] ['* - AutoPrim,  x - AutoMets']}} ;

            labelcellstr = {{'PC 1' 'PC 2' 'PC 3' 'PC 4'}} ;

            savestr = ['plots/Pancreatic2ip' num2str(ipart) 'PCAallAfterNewAxes'] ;

            paramstruct = struct('icolor',icolor, ...
                'npcadiradd',4, ...
                'markerstr',markerstr, ...
                'isubpopkde',1, ...
                'idataconn',idataconn, ...
                'idataconncolor',idataconncolor, ...
                'titlecellstr',titlecellstr, ...
                'legendcellstr',legendcellstr, ...
                'mlegendcolor',mlegendcolor, ...
                'labelcellstr',labelcellstr, ...
                'savestr',savestr, ...
                'iscreenwrite',1) ;

            scatplotSM(mdataadj,[],paramstruct) ;




        elseif ipart == 102 ;


            %  Scott adds stuff here









        elseif  ipart == 121  | ...
                ipart == 122  | ...
                ipart == 123  | ...
                ipart == 124  | ...
                ipart == 125;

            if ipart == 121 ;    %  Study Liver mets vs. Primary, 500 gene list (merged data)
                TissueStr = 'Liver' ;
                TissueFlag = (vTissueNum == 4) ;
                %  liver is the 4th tissue
            elseif ipart == 122 ;    %  Study Lung mets vs. Primary, 500 gene list (merged data)
                TissueStr = 'Lung' ;
                TissueFlag = (vTissueNum == 3) ;
                %  lung is the 3rd tissue
            elseif ipart == 123 ;    %  Study Lymph mets vs. Primary, 500 gene list (merged data)
                TissueStr = 'Lymph' ;
                TissueFlag = (vTissueNum == 5) ;
                %  lymph is the 5th tissue
            elseif ipart == 124 ;    %  Study Fat mets vs. Primary, 500 gene list (merged data)
                TissueStr = 'Fat' ;
                TissueFlag = (vTissueNum == 6) ;
                %  fat is the 6th tissue
            elseif ipart == 125 ;    %  Study Muscle mets vs. Primary, 500 gene list (merged data)
                TissueStr = 'Muscle' ;
                TissueFlag = (vTissueNum == 7) ;
                %  muscle is the 7th tissue
            end ;

            redflag = PancreasFlag | TissueFlag ;
            mdatared = mdata(:,redflag) ;
            TissueNormalFlag = TissueFlag(redflag) & NormalFlag(redflag) ;
            PancreasNormalFlag = PancreasFlag(redflag) & NormalFlag(redflag) ;
            TissueTumorFlag = TissueFlag(redflag) & (~NormalFlag(redflag)) ;
            PancreasTumorFlag = PancreasFlag(redflag) & (~NormalFlag(redflag)) ;

            %  Lines from BatchAdjustSM.m
            %
            dirvec = DWD1SM(mdatared(:,PancreasNormalFlag),mdatared(:,TissueNormalFlag)) ;

            vprojp = mdatared(:,PancreasNormalFlag)' * dirvec ;
            vprojn = mdatared(:,TissueNormalFlag)' * dirvec ;

            meanprojp = mean(vprojp) ;
            meanprojn = mean(vprojn) ;

            adjdata = mdatared ;
            adjdata(:,PancreasFlag(redflag)) = mdatared(:,PancreasFlag(redflag)) ...
                - vec2matSM(meanprojp * dirvec,sum(PancreasFlag)) ...
                + vec2matSM(meanprojn * dirvec,sum(PancreasFlag)) ;
            adjdata(:,TissueFlag(redflag)) = mdatared(:,TissueFlag(redflag)) ;


            %  Make before and after PCA plots
            %
            %  First compute common axes
            %
            paramstruct = struct('npc',4,...
                'viout',[0 1]) ;
            outstruct = pcaSM(mdatared,paramstruct) ;
            mdir = getfield(outstruct,'meigvec') ;


            %  Before PCA plot
            %
            titlecellstr = {{[TissueStr ' & Pancreas PCA, n = ' num2str(n) ', d = ' num2str(d)] ...
                'Before Removing Tissue Effect' ...
                ['o - normal,   + - Resected'] ['* - AutoPrim,  x - AutoMets']}} ;
            labelcellstr = {{'PC 1' 'PC 2' 'PC 3' 'PC 4'}} ;
            savestr = ['plots/Pancreatic2ip' num2str(ipart) 'PCA' TissueStr '-PancreasBefore'] ;

            paramstruct = struct('icolor',icolor(redflag,:), ...
                'npcadiradd',0, ...
                'markerstr',markerstr(redflag), ...
                'isubpopkde',1, ...
                'titlecellstr',titlecellstr, ...
                'legendcellstr',legendcellstr, ...
                'mlegendcolor',mlegendcolor, ...
                'labelcellstr',labelcellstr, ...
                'savestr',savestr, ...
                'iscreenwrite',1) ;
            %  idataconn has deliberately been skipped
            %  may want to add back in later
            scatplotSM(mdatared,mdir,paramstruct) ;


            %  After PCA plot
            %
            titlecellstr = {{[TissueStr ' & Pancreas PCA, n = ' num2str(n) ', d = ' num2str(d)] ...
                'After Removing Tissue Effect' ...
                ['o - normal,   + - Resected'] ['* - AutoPrim,  x - AutoMets']}} ;
            labelcellstr = {{'PC 1' 'PC 2' 'PC 3' 'PC 4'}} ;
            savestr = ['plots/Pancreatic2ip' num2str(ipart) 'PCA' TissueStr '-PancreasAfter'] ;

            paramstruct = struct('icolor',icolor(redflag,:), ...
                'npcadiradd',0, ...
                'markerstr',markerstr(redflag), ...
                'isubpopkde',1, ...
                'titlecellstr',titlecellstr, ...
                'legendcellstr',legendcellstr, ...
                'mlegendcolor',mlegendcolor, ...
                'labelcellstr',labelcellstr, ...
                'savestr',savestr, ...
                'iscreenwrite',1) ;
            %  idataconn has deliberately been skipped
            %  may want to add back in later
            scatplotSM(adjdata,mdir,paramstruct) ;


            %  Show DWD direction Mets vs. Primary
            %
            DWDdir = DWD1SM(adjdata(:,PancreasTumorFlag),adjdata(:,TissueTumorFlag)) ;

            titlecellstr = {{['DWD Panc. vs. ' TissueStr ' Tumor, n = ' num2str(n) ', d = ' num2str(d)] ...
                'After Removing Tissue Effect' ...
                ['o - normal,   + - Resected'] ['* - AutoPrim,  x - AutoMets']}} ;
            labelcellstr = {{'DWD' 'OPC 1' 'OPC 2' 'OPC 3'}} ;
            savestr = ['plots/Pancreatic2ip' num2str(ipart) 'DWDPancv' TissueStr '-Tumors'] ;

            paramstruct = struct('icolor',icolor(redflag,:), ...
                'npcadiradd',-3, ...
                'markerstr',markerstr(redflag), ...
                'isubpopkde',1, ...
                'titlecellstr',titlecellstr, ...
                'legendcellstr',legendcellstr, ...
                'mlegendcolor',mlegendcolor, ...
                'labelcellstr',labelcellstr, ...
                'savestr',savestr, ...
                'iscreenwrite',1) ;
            %  idataconn has deliberately been skipped
            %  may want to add back in later
            scatplotSM(adjdata,DWDdir,paramstruct) ;


            %  Print gene list, i.e. make output text file
            %
            ngeneout = 500 ;
            [temp,visort] = sort(abs(DWDdir),1,'descend') ;
            fid = fopen(['plots/Pancreatic2ip' num2str(ipart) 'GeneList' ...
                num2str(ngeneout) 'PrimVs' TissueStr '.txt'],'wt') ;
            %  'wt' is for "delete contents of this file and open
            %               for writing" (with 't' for "text").

            titlstr1 = 'Output from the MATLAB Script Pancreatic2.m,   ' ;
            titlstr1 = [titlstr1,date] ;
            cntbytes = fprintf(fid,'%1s\n\n\n',titlstr1) ;
            %  '%1s/n' says "print string, followed by a line feed"
            %       (Note, this only goes to next line)

            titlstr2 = ['     Top ' num2str(ngeneout) ' genes (in absolute value order), for Primary vs. ' TissueStr ': '] ;
            cntbytes = fprintf(fid,'%1s\n\n',titlstr2) ;
            %  3 line feeds gives two blank lines

            forstr = '%s\t%s\t%s\t%s\t%s\n' ;
            cntbytes = fprintf(fid,forstr,'GeneRank','DWDwt','PvalT','AgilentID','GeneName') ;

            PvalTs = [] ;
            for i = 1:d ;
                iout = visort(i) ;
                DWDwt = abs(DWDdir(iout)) ;
                [H,P] = ttest2(adjdata(iout,PancreasTumorFlag),adjdata(iout,TissueTumorFlag),0.025,'both') ;
                %  do two sided 2 sample t-test,
                P = P / 2 ;
                %  but divide by 2, to get test on absolute value
                if i <= ngeneout ;    %  then print to file
                    forstr = '%3d\t%6.4f\t%10.8f\t%s\t%s\n' ;
                    %  Embedded text appears as written (including parens)
                    cntbytes = fprintf(fid,forstr,i,DWDwt,P,AgilentIDs{iout},GeneNames{iout}) ;
                    PvalTs = [PvalTs; P]  ;
                end ;
            end ;

            fclose('all') ;
            %  !!!! Careful !!!!  Need to do this to be able work with files after
            %  running matlab.  Alternative is:
            %  fclose(fid) ;





        elseif  ipart == 131  | ...
                ipart == 132  | ...
                ipart == 133  ;


            if ipart == 131 ;    %  Study Fat&Muscle mets vs. Primary, 500 gene list (merged data)
                TissueStr = 'Fat&Muscle' ;
                TissueFlag = (vTissueNum == 6) | (vTissueNum == 7) ;
                %  fat is the 6th tissue
                %  muscle is the 7th tissue
                vTissueIndex = [6 7] ;
                TissueCell = {'Fat' 'Muscle'} ;
            elseif ipart == 132 ;    %  Study Lung&Liver mets vs. Primary, 500 gene list (merged data)
                TissueStr = 'Lung&Liver' ;
                TissueFlag = (vTissueNum == 3) | (vTissueNum == 4) ;
                %  lung is the 3rd tissue
                %  liver is the 4th tissue
                vTissueIndex = [3 4] ;
                TissueCell = {'Lung' 'Liver'} ;
            elseif ipart == 133 ;    %  Study Lymph&Fat&Muscle&Lung&Liver mets vs. Primary, 500 gene list (merged data)
                TissueStr = 'Lymph&Fat&Muscle&Lung&Liver' ;
                TissueFlag = (vTissueNum == 5) | (vTissueNum == 6) | (vTissueNum == 7) | ...
                    (vTissueNum == 3) | (vTissueNum == 4) ;
                %  lymph is the 5th tissue
                %  fat is the 6th tissue
                %  muscle is the 7th tissue
                %  lung is the 3rd tissue
                %  liver is the 4th tissue
                vTissueIndex = [5 6 7 3 4] ;
                TissueCell = {'Lymph' 'Fat' 'Muscle' 'Lung' 'Liver'} ;
            end ;


            redflag = PancreasFlag | TissueFlag ;
            mdatared = mdata(:,redflag) ;
            TissueNormalFlag = TissueFlag(redflag) & NormalFlag(redflag) ;
            PancreasNormalFlag = PancreasFlag(redflag) & NormalFlag(redflag) ;
            TissueTumorFlag = TissueFlag(redflag) & (~NormalFlag(redflag)) ;
            PancreasTumorFlag = PancreasFlag(redflag) & (~NormalFlag(redflag)) ;


            %  Make before PCA plots
            %
            %  First compute common axes
            %
            paramstruct = struct('npc',4,...
                'viout',[0 1]) ;
            outstruct = pcaSM(mdatared,paramstruct) ;
            mdir = getfield(outstruct,'meigvec') ;


            %  Before PCA plot
            %
            titlecellstr = {{[TissueStr ' & Pancreas PCA'] ...
                'Before Removing Tissue Effect' ...
                ['o - normal,   + - Resected'] ['* - AutoPrim,  x - AutoMets']}} ;
            labelcellstr = {{'PC 1' 'PC 2' 'PC 3' 'PC 4'}} ;
            savestr = ['plots/Pancreatic2ip' num2str(ipart) 'PCA' TissueStr '-PancreasBefore'] ;

            paramstruct = struct('icolor',icolor(redflag,:), ...
                'npcadiradd',0, ...
                'markerstr',markerstr(redflag), ...
                'isubpopkde',1, ...
                'titlecellstr',titlecellstr, ...
                'legendcellstr',legendcellstr, ...
                'mlegendcolor',mlegendcolor, ...
                'labelcellstr',labelcellstr, ...
                'savestr',savestr, ...
                'iscreenwrite',1) ;
            %  idataconn has deliberately been skipped
            %  may want to add back in later
            scatplotSM(mdatared,mdir,paramstruct) ;


            %  Adjust each tissue normal to Pancreatic normal
            %
            i = 0 ;
            adjdata = mdatared ;
            % set tissue counter
            for iTI = vTissueIndex ;

                i = i + 1 ;
                %  increment tissue counter

                TissueFlag = (vTissueNum == iTI) ;
                IndivTissueNormalFlag = TissueFlag & NormalFlag ;
                IndivTissueNormalFlag = IndivTissueNormalFlag(redflag) ;


                %  Lines from BatchAdjustSM.m
                %
                dirvec = DWD1SM(mdatared(:,IndivTissueNormalFlag),mdatared(:,PancreasNormalFlag)) ;

                vprojp = mdatared(:,IndivTissueNormalFlag)' * dirvec ;
                vprojn = mdatared(:,PancreasNormalFlag)' * dirvec ;

                meanprojp = mean(vprojp) ;
                meanprojn = mean(vprojn) ;

                adjdata(:,TissueFlag(redflag)) = mdatared(:,TissueFlag(redflag)) ...
                    - vec2matSM(meanprojp * dirvec,sum(TissueFlag)) ...
                    + vec2matSM(meanprojn * dirvec,sum(TissueFlag)) ;
                adjdata(:,PancreasFlag(redflag)) = mdatared(:,PancreasFlag(redflag)) ;


                %  After PCA plot
                %
                titlecellstr = {{[TissueStr ' & Pancreas PCA'] ...
                    ['After Removing ' TissueCell{i} ' Effect'] ...
                    ['o - normal,   + - Resected'] ['* - AutoPrim,  x - AutoMets']}} ;
                labelcellstr = {{'PC 1' 'PC 2' 'PC 3' 'PC 4'}} ;
                savestr = ['plots/Pancreatic2ip' num2str(ipart) 'PCA' TissueCell{i} '-PancreasAfter'] ;

                paramstruct = struct('icolor',icolor(redflag,:), ...
                    'npcadiradd',0, ...
                    'markerstr',markerstr(redflag), ...
                    'isubpopkde',1, ...
                    'titlecellstr',titlecellstr, ...
                    'legendcellstr',legendcellstr, ...
                    'mlegendcolor',mlegendcolor, ...
                    'labelcellstr',labelcellstr, ...
                    'savestr',savestr, ...
                    'iscreenwrite',1) ;
                %  idataconn has deliberately been skipped
                %  may want to add back in later
                scatplotSM(adjdata,mdir,paramstruct) ;

            end ;


            %  Show DWD direction Mets vs. Primary
            %
            DWDdir = DWD1SM(adjdata(:,PancreasTumorFlag),adjdata(:,TissueTumorFlag)) ;

            titlecellstr = {{['DWD Panc. vs. ' TissueStr ' Tumor'] ...
                'After Removing Tissue Effect' ...
                ['o - normal,   + - Resected'] ['* - AutoPrim,  x - AutoMets']}} ;
            labelcellstr = {{'DWD' 'OPC 1' 'OPC 2' 'OPC 3'}} ;
            savestr = ['plots/Pancreatic2ip' num2str(ipart) 'DWDPancv' TissueStr '-Tumors'] ;

            paramstruct = struct('icolor',icolor(redflag,:), ...
                'npcadiradd',-3, ...
                'markerstr',markerstr(redflag), ...
                'isubpopkde',1, ...
                'titlecellstr',titlecellstr, ...
                'legendcellstr',legendcellstr, ...
                'mlegendcolor',mlegendcolor, ...
                'labelcellstr',labelcellstr, ...
                'savestr',savestr, ...
                'iscreenwrite',1) ;
            %  idataconn has deliberately been skipped
            %  may want to add back in later
            scatplotSM(adjdata,DWDdir,paramstruct) ;


            %  Print gene list, i.e. make output text file
            %
            ngeneout = 500 ;
            [temp,visort] = sort(abs(DWDdir),1,'descend') ;
            fid = fopen(['plots/Pancreatic2ip' num2str(ipart) 'GeneList' ...
                num2str(ngeneout) 'PrimVs' TissueStr '.txt'],'wt') ;
            %  'wt' is for "delete contents of this file and open
            %               for writing" (with 't' for "text").

            titlstr1 = 'Output from the MATLAB Script Pancreatic2.m,   ' ;
            titlstr1 = [titlstr1,date] ;
            cntbytes = fprintf(fid,'%1s\n\n\n',titlstr1) ;
            %  '%1s/n' says "print string, followed by a line feed"
            %       (Note, this only goes to next line)

            titlstr2 = ['     Top ' num2str(ngeneout) ' genes (in absolute value order), for Primary vs. ' TissueStr ': '] ;
            cntbytes = fprintf(fid,'%1s\n\n',titlstr2) ;
            %  3 line feeds gives two blank lines

            forstr = '%s\t%s\t%s\t%s\t%s\n' ;
            cntbytes = fprintf(fid,forstr,'GeneRank','DWDwt','PvalT','AgilentID','GeneName') ;

            PvalTs = [] ;
            for i = 1:d ;
                iout = visort(i) ;
                DWDwt = abs(DWDdir(iout)) ;
                [H,P] = ttest2(adjdata(iout,PancreasTumorFlag),adjdata(iout,TissueTumorFlag),0.025,'both') ;
                %  do two sided 2 sample t-test,
                P = P / 2 ;
                %  but divide by 2, to get test on absolute value
                if i <= ngeneout ;    %  then print to file
                    forstr = '%3d\t%6.4f\t%10.8f\t%s\t%s\n' ;
                    %  Embedded text appears as written (including parens)
                    cntbytes = fprintf(fid,forstr,i,DWDwt,P,AgilentIDs{iout},GeneNames{iout}) ;
                    PvalTs = [PvalTs; P]  ;
                end ;
            end ;

            fclose('all') ;
            %  !!!! Careful !!!!  Need to do this to be able work with files after
            %  running matlab.  Alternative is:
            %  fclose(fid) ;


            if ipart == 133 ;    %  then also print Excel spreadsheet

              outfilestr = 'pancreatic2adj5TissSubtr.xls' ;
              [Success,Message] = xlswrite(outfilestr,adjdata,'MainAdjData')
              [Success,Message] = xlswrite(outfilestr,GeneNames,'GeneNames')
              [Success,Message] = xlswrite(outfilestr,SlideNames(redflag),'CaseNames')


            end ;



        end ;



    end ;






end ;