% Use the principal components of the filtered genes to select the initial
% motif pool. 
load filteredData.mat
[pc, score, latent, tsquare] = princomp(filteredY);
perc_exp = 100*latent/sum(latent);
for i=1:10
    subplot(5,2,i);
    hold off
    plot(pc(1:5,i));
    hold on
    plot(pc(6:10,i),':');
    label = sprintf('pc %d, explains %2.1f%% of var.', i, perc_exp(i));
    ylim([-1 1])
    xlabel(label);
end

% Note that the slopes for len 5,6,7 motifs were calculated USING ENTIRE GENE LIST
% It would be technically coherent to do the entire analysis again starting
% from slopes calculated using ONLY the filtered genes.
pcutoff=0.0001;
maxpc=5;
minpc=2;
pcmotifs = cell(1,maxpc-minpc+1);
motifpool = {};
for pcind = 2:3


    pcutoff=0.0001;
    load len5slopes.mat
    [pcScore,pcScoreVar,pcScoreP] = calcLinearScore(slopeLen5, bSigmaLen5, pc(:,pcind)');
    [pcScorePsorted, sortedInd]=sort(pcScoreP);
    n=length(find(pcScoreP<pcutoff));
    len5include = sortedInd(1:3:n-1);

    load len6slopes.mat
    [pcScore,pcScoreVar,pcScoreP] = calcLinearScore(slopeLen6, bSigma, pc(:,pcind)');
    [pcScorePsorted, sortedInd]=sort(pcScoreP);
    n=length(find(pcScoreP<pcutoff));
    len6include = sortedInd(1:3:n-1);

    load len7slopes.mat
    [pcScore,pcScoreVar,pcScoreP] = calcLinearScore(slopeLen7, bSigmaLen7, pc(:,pcind)');
    [pcScorePsorted, sortedInd]=sort(pcScoreP);
    n=length(find(pcScoreP<pcutoff));
    len7include = sortedInd(1:3:n-1);

    pcmotifs{pcind-minpc+1} =  [len5motifs(len5include) len6motifs(len6include) len7motifs(len7include)];
    motifpool = union(motifpool,pcmotifs{pcind-minpc+1})
end

save  len567pc23motifpool.mat motifpool