Upload New File

demo_code.m

+% ConsensusOPLS Script example
+%
+% Goal: Discriminant Analysis of NCI-60 cancer cell lines from two origins (Colon vs Ovary)
+% Data: 3 data blocks evaluated on the same observations
+%       Metabolomics: 14 x 150 variables
+%       Microarray: 14 x 200 variables
+%       Proteomics: 14 x 100 variables
+% 
+% Please cite: 
+% J. Boccard, D.N. Rutledge.
+% A consensus OPLS-DA strategy for multiblock Omics data fusion
+% Analytica Chimica Acta, 769, 30-39 (2013).
+%
+%
+%% Data preparation
+
+load('demo_data.mat');
+
+% Unit Variance scaling of each data matrix
+collection(1)=matrix2saisir(zscore(MetaboData));
+collection(2)=matrix2saisir(zscore(MicroData));
+collection(3)=matrix2saisir(zscore(ProteoData));
+
+BlockNames={'Metabolomics', 'Microarray', 'Proteomics'};
+
+collection(1).i=MetaboVarNames;
+collection(2).i=MicroVarNames;
+collection(3).i=ProteoVarNames;
+
+
+%% Compute a consensusOPLS model with leave-one-out cross-validation
+
+% Definition of parameters
+LVsPred=1; % Number of predictive component(s)
+LVsOrtho=1; % Maximum number of orthogonal components
+CVfolds=14; % Number of cross-validation folds
+
+% This is the main function to compute the consensusOPLS model
+RVConsensusOPLSModelCV=RVConsensusOPLS(collection,Y,LVsPred,LVsOrtho,CVfolds,'nfold','da',0);
+
+disp('RVConsensusOPLS Results ');
+disp(['R2: ' num2str(RVConsensusOPLSModelCV.koplsModel.R2Yhat(RVConsensusOPLSModelCV.cv.OrthoLVsOptimalNum+1))]);
+disp(['Q2: ' num2str(RVConsensusOPLSModelCV.cv.Q2Yhat(RVConsensusOPLSModelCV.cv.OrthoLVsOptimalNum+1))]);
+disp(['DQ2: ' num2str(RVConsensusOPLSModelCV.cv.DQ2Yhat(RVConsensusOPLSModelCV.cv.OrthoLVsOptimalNum+1))]);
+disp(' ');
+disp('Confusion Matrix:');
+disp(RVConsensusOPLSModelCV.da.confusionMatrix);
+
+%% Plot the results
+
+% Consensus Score plot
+figure; plot(RVConsensusOPLSModelCV.koplsModel.T(:,1),RVConsensusOPLSModelCV.koplsModel.To(:,1),'k.');
+title('ConsensusOPLS Score plot');
+axis([-0.5 0.5 -1 1]);
+xlabel('Predictive component')
+ylabel('Orthogonal component')
+text(RVConsensusOPLSModelCV.koplsModel.T(:,1), RVConsensusOPLSModelCV.koplsModel.To(:,1), ObsNames(:,1), 'HorizontalAlignment','left','VerticalAlignment','bottom','FontSize',8)
+
+% Block contributions to the predictive component
+figure; bar(RVConsensusOPLSModelCV.koplsModel.lambda(:,1));
+set(gca,'xticklabel',BlockNames)
+ylabel('Weight');
+title('Block contributions to the predictive component');
+
+% Block contributions to the first orthogonal component
+figure; bar(RVConsensusOPLSModelCV.koplsModel.lambda(:,2));
+set(gca,'xticklabel',BlockNames)
+ylabel('Weight');
+title('Block contributions to the first orthogonal component');
+
+% Block contributions predictive vs. orthogonal
+figure; scatter(RVConsensusOPLSModelCV.koplsModel.lambda(:,1),RVConsensusOPLSModelCV.koplsModel.lambda(:,2));
+axis([0 0.5 0 0.5]);
+xlabel('Predictive component')
+ylabel('Orthogonal component')
+text(RVConsensusOPLSModelCV.koplsModel.lambda(:,1), RVConsensusOPLSModelCV.koplsModel.lambda(:,2), BlockNames, 'HorizontalAlignment','left','VerticalAlignment','bottom','FontSize',8)
+title('Block contributions');
+
+% Loading plots (one for each table)
+figure;
+title('ConsensusOPLS Loading plot');
+xlabel('Predictive component');
+ylabel('Orthogonal component');
+hold all;
+for i=1:length(collection)
+    scatter(RVConsensusOPLSModelCV.koplsModel.loadings{i,1},RVConsensusOPLSModelCV.koplsModel.loadings{i,2},'.','DisplayName',['Block ' num2str(i)]);
+end
+legend(gca,'show')
+
+
+%% Permutations
+[PermRes]=RVConsensusOPLSPerm(collection,Y,1000,LVsPred,LVsOrtho);
+