
%%% need to install matlab toolboxes:
%%%   FastICA http://research.ics.aalto.fi/ica/fastica
%%%   FSLNets (downloaded later than 26 August 2016) http://fsl.fmrib.ox.ac.uk/fsl/fslwiki/FSLNets
%%%   htmlTableToCell http://uk.mathworks.com/matlabcentral/fileexchange/34968-htmltabletocell
%%%   nearestSPD http://uk.mathworks.com/matlabcentral/fileexchange/42885-nearestspd
%%%   FDR https://www2.warwick.ac.uk/fac/sci/statistics/staff/academic-research/nichols/software/fdr
addpath ~/matlab ~/matlab/FSLNets ~/matlab/FastICA_25

%%% Import big text file with all imaging IDPs (except rfMRI IDPs), obtained from UK Biobank.
%%% This should be arranged with one column per IDP (starting with subject ID),
%%% using exactly the same IDP ordering as found in the IDPinfo.txt file found in
%%% http://www.fmrib.ox.ac.uk/ukbiobank
%%% The top line should list the IDP names: a horizontal version of the first column in IDPinfo.txt
%%% The second line should be the first subject's IDP values, etc.
IDPs=importdata('IDPs.txt');
ALL_IDPs=IDPs.data;
N=size(ALL_IDPs,1);

%%% load rfMRI partial correlation netmats and node-amplitudes (i.e. from separate files than the above)
%%% these files should come (from UK Biobank) with one line per subject, no header line, and
%%% for each line, start with the subject ID and then have the unwrapped netmat or node amplitudes.
grot=load('IDP_rfMRI_d25_partialcorr.txt'));  NET25=zeros(N,size(grot,2)-1)/0;
for i = 1:size(grot,1), NET25(find(ALL_IDPs(:,1)==grot(i,1)),:)=grot(i,2:end); end;
grot=load('IDP_rfMRI_d25_NodeAmplitudes.txt'));  NODEamps25=zeros(N,size(grot,2)-1)/0;
for i = 1:size(grot,1), NODEamps25(find(ALL_IDPs(:,1)==grot(i,1)),:)=grot(i,2:end); end;
grot=load('IDP_rfMRI_d100_partialcorr.txt'));  NET100=zeros(N,size(grot,2)-1)/0;
for i = 1:size(grot,1), NET100(find(ALL_IDPs(:,1)==grot(i,1)),:)=grot(i,2:end); end;
grot=load('IDP_rfMRI_d100_NodeAmplitudes.txt'));  NODEamps100=zeros(N,size(grot,2)-1)/0;
for i = 1:size(grot,1), NODEamps100(find(ALL_IDPs(:,1)==grot(i,1)),:)=grot(i,2:end); end;

%%% get list of IDPs and modality groups
! cat IDPinfo.txt | awk '{print $1}' > grot.txt
po=fopen('grot.txt','r');  IDPnames=textscan(po,'%s');  fclose(po);  IDPnames=IDPnames{1};
j=length(IDPnames)+1;
for i=1:size(NODEamps25,2),   IDPnames{j}=sprintf('rfMRI amplitudes (ICA25 node %d)',i);    j=j+1;  end;
for i=1:size(NODEamps100,2),  IDPnames{j}=sprintf('rfMRI amplitudes (ICA100 node %d)',i);   j=j+1;  end;
for i=1:size(NET25,2),        IDPnames{j}=sprintf('rfMRI connectivity (ICA25 edge %d)',i);  j=j+1;  end;
for i=1:size(NET100,2),       IDPnames{j}=sprintf('rfMRI connectivity (ICA100 edge %d)',i); j=j+1;  end;
IDP_modality_types(1:18)=0;      % QC
IDP_modality_types(19:43)=1;     % T1
%%%%IDP_modality_types(??????)=2;     % T2 - no IDPs yet.
IDP_modality_types(44:57)=3;     % SWI
IDP_modality_types(58:73)=4;     % tfMRI
IDP_modality_types(74:748)=5;    % dMRI
IDP_modality_types(749:824)=6;   % rfMRI amplitudes
IDP_modality_types(825:size(IDPnames,1))=7; % rfMRI netmats

%%% load non-IDP data obtained from UK Biobank (including info from the supplied HTML doc files).
%%% The data text files should be arranged to have one line per subject (no header line),
%%% and with one column for each of the rows in the main initial table in the HTML.
%%% Hence the first column corresponds to variable number "0" in the HTML.
%%% The data text file should be preprocessed to only have numerical (and NaN) values, and be space-separated.
%%% Hence dates and/or times of importance should be converted to pure numerical form.
varsSUBS=[]; varsVARS=[];
%%% The kk list of numbers below corresponds to if you have multiple distinct sets of variables obtained
%%% from UK Biobank (ie multiple pairs of HTML and data text files).
kk=[ 5963 6034 6225 ];
for k = kk;
  varsORIG=load(sprintf('ukb%d.txt',k));
  varsSUBS=sort(unique([varsSUBS;varsORIG(:,1)])); % final list of subject IDs
  clear table;  table.idTableBy.plaintextPreceedingTable='Ukb2docs';  
  varsHTML=htmlTableToCell(sprintf('ukb%d.html',k),table);
  varsVARS=sort(unique([ varsVARS {varsHTML{3:end,2}}]));
end
%%% incrementally fill out full vars matrix and headers lists, overwriting earlier versions with later ones
vars=zeros(length(varsSUBS),length(varsVARS))/0;  clear varsHeader; clear varsSource;
for k = kk
  varsORIG=load(sprintf('ukb%d.txt',k));
  clear table;  table.idTableBy.plaintextPreceedingTable='Ukb2docs';
  varsHTML=htmlTableToCell(sprintf('ukb%d.html',k),table);
  [~,grotSa,grotSb]=intersect(varsSUBS,varsORIG(:,1));
  [~,grotVa,grotVb]=intersect(varsVARS,{varsHTML{3:end,2}});
  vars(grotSa,grotVa)=varsORIG(grotSb,1+grotVb);  % main variables data matrix
  varsHeader(grotVa')={varsHTML{2+grotVb,5}};     % variable names
  varsSource(grotVa')=k;                          % which of the original pairs of files was each variable taken from
end

%%% fill out empty variable names (occurs for repeat visits), and insert ". " before "Uses"
for i=size(varsHeader,2):-1:1
  if prod(size(varsHeader{i}))==0
    j=i-1;
    while prod(size(varsHeader{j}))==0
      j=j-1;
    end
    varsHeader{i}=[varsHeader{j} sprintf(' visit-%d',1+i-j)];
  end
  j=regexp(varsHeader{i},'Uses');
  if prod(size(j))>0
    grot=varsHeader{i};  varsHeader{i} = [ grot(1:j-1) '. ' grot(j:end) ];
  end
end

%%% discard non-imaged subjects
varsORIG=vars; vars=zeros(N,size(varsORIG,2))/0;
for i=1:N
  j=find(ALL_IDPs(i,1)==varsSUBS);
  if prod(size(j))==1
    vars(i,:)=varsORIG(j,:);
  end
end

%%% identify bad vars on the basis of too much missing data, etc
badvars=[];
for i=1:size(vars,2)
  Y=vars(:,i);
  if sum(abs(mod(Y(~isnan(Y)),1))) == 0  % is this an integer/categorical variable?
    Y(Y<0)=NaN;  % if so then, mostly, negative numbers should be treated as NaN, mostly
  end
  grotKEEP=~isnan(Y);  
  % do we have extreme outliers? This part currently ignored
  grot=(Y(grotKEEP)-median(Y(grotKEEP))).^2; grot=max(grot/mean(grot));  grot=0;
  if (sum(grotKEEP)/N>0.1) & (std(Y(grotKEEP))>0) & (max(sum(nets_class_vectomat(Y(grotKEEP))))/length(Y(grotKEEP))<0.95) & (grot<100)
    % the 3rd thing above is:  is the size of the largest equal-values-group too large?
    vars(:,i)=Y;
  else
    % remove this variable
    badvars=[badvars i];
  end
end
badvars=unique(badvars);  varskeep=setdiff([1:size(vars,2)],badvars);

%%% identify "bad" vars on the basis of redundancy
for i=varskeep
  for j=varskeep
    if j>i & corr(vars(:,i),vars(:,j),'rows','pairwise')>0.9999
      if sum(~isnan(vars(:,i))) > sum(~isnan(vars(:,j)))
        k=i;
      else
        k=j;
      end
      badvars=[badvars k];
    end
  end
end
badvars=unique(badvars);  varskeep=setdiff([1:size(vars,2)],badvars);

%%%%%%%%%%%%%%%%%%%% categories of non-IDP variables
% 1  age, sex
% 2  genetics
% 3  early life factors
% 10 lifestyle and environment - general
% 11 lifestyle and environment - exercise and work
% 12 lifestyle and environment - food and drink
% 13 lifestyle and environment - alcohol
% 14 lifestyle and environment - tobacco
% 20 physical measures - general
% 21 physical measures - bone density and sizes
% 22 physical measures - cardiac
% 30 blood assays
% 31 brain IDPs
% 32 cognitive phenotypes
% 50 health and medical history, health outcomes
% 99 misc, ignored

vt{1}=[31 34 22200];
vt{2}=[21000 22000:22125 22201:22325];
vt{3}=[129 130 1677 1687 1697 1737 1767 1777 1787 20022];
vt{10}=[189 2267 2277 4825 4836 6138 6142 20118 20119 20121 22501 22599 22606];
vt{11}=[1001 1011 767 777 796 806 816 826 845 864 874 884 894 904 914 924 943 971 981 991 1021 1050:10:1220 2624 2634 ...
        3426 3637 3647 6143 6162 6164 22604 22605 22607:22615 22620 22630 22631 22640:22655];
vt{12}=[1289:10:1389 1408:10:1548 2654 3680 6144 20084:20094 20098:20109 100010:10:100560 100760:10:104670];
vt{13}=[1558:10:1628 2664 3731 3859 4407 4418 4429 4440 4451 4462 5364 20095:20097 20117 100580:10:100740];
vt{14}=[1239:10:1279 2644 2867:10:2907 2926 2936 3436:10:3506 5959 6157 6158 6183 6194 20116 22506:22508];
vt{20}=[46:50 20015 1707 1717 1727 1747 1757 2306 3062 3063 3064 21001 21002 22400:22408 22427 23099:23130 23244:23289];
vt{21}=[77 78 4106 4125 4138 4143 23200:23243 23290:23320];
vt{22}=[93:95 102 4079 4080 4194 12673:12687 12697 12698 12702 21021 22420:22426];
vt{30}=[30000:10:30300];
vt{31}=[25000:25746 25761:25768];
vt{32}=[399 20016 20018 20023 20082 20129:20157 20159 20195 20196 20198 20200 20229 20230 20240 20244:20248];
vt{50}=[134:137 2178 2188 2207:10:2257 2296 2316 2335:10:2365 2415 2443:10:2473 2492 2844 2956:10:2986 3005 3079 3140 ...
        3393 3404 3414 3571 3606 3616 3627 3741 3751 3761 3773 3786 3799 3809 3894 3992 ...
        4012 4022 4041 4056 4067 4689 4700 4717 4728 4792 4803 4814 5408 5419 5430 5441 5452 5463 5474 5485 5496 5507 ...
        5518 5529 5540 5610 5832 5843 5855 5877 5890 5901 5912 5923 5934 5945 6119 6147 6148 6149 ...
        6150 6151 6152 6153 6154 6155 6159 6177 6179 6205 20001:20010 22126:22181 22502:22505 22616 22618 22619 40001:41253];
vt{99}=[21 52 53 54 120 132 757 1647 12139 12140 12141 12187 12188 12223 12224 12253 12254 12623 12624 12651 12652 ...
        12663 12664 12671 12695 12699 12700 12704 12706 12848 12851 12854 ...
        20012:20014 20024 20074 20075 20077 20078 20083 20115 20158 20201:20227 20249:20253 21003 22499 22500 ...
        22600:22603 22617 22660:22664 25747:25753 25780 40000 105010 105030];

vt_use=[ 3 10:14 20:22 30 32 ]; % by default use these in IDP-vars analyses
vt_not_use=[ 99 1 2 31 50 ];    % and don't use these (yet, at least)

%%% identify vars of special interest
dosI=nets_cellfind(varsVARS,'53-2.0'); scan_date=vars(:,dosI); scan_date(isnan(scan_date))=2015.7;
yobI=nets_cellfind(varsVARS,'34-0.0'); mobI=nets_cellfind(varsVARS,'52-0.0'); birth_date=vars(:,yobI)+(vars(:,mobI)-0.5)/12;
age=scan_date - birth_date;
sexI=nets_cellfind(varsVARS,'31-0.0');  sex=vars(:,sexI);
badvars=[ badvars dosI yobI mobI sexI ];

%%% discard some more vars
for i=[ vt{99} vt{1} vt{2} vt{31} vt{50} ]
  j=nets_cellfind(varsVARS,sprintf('%d-',i));
  if ~isempty(j)
    badvars=[badvars j];
  end
end
badvars=unique(badvars);  varskeep=setdiff([1:size(vars,2)],badvars);  varskeepOLD=varskeep;

%%% create clean nicely ordered list of kept-variables (varskeep)
varskeepNEW=[]; varskeepVT=[];
for ii=vt_use;
  for i = vt{ii}
    j=nets_cellfind(varsVARS,sprintf('%d-',i));
    for jj=j
      if ~isempty(find(varskeep==jj))
        varskeepNEW=[varskeepNEW jj];
        varskeepVT=[varskeepVT ii];
      end
    end
  end
end
varskeep=varskeepNEW;
for i=1:length(varskeep)
   disp( sprintf('%d %d %s',i,varskeepVT(i),varsHeader{varskeep(i)} ) )
end

%%% number of PCA components to input into CCA
Nkeep=200;

%%% create confounds
age2=nets_demean(age).^2;
conf=[ age age2 sex nets_demean(age).*nets_demean(sex) age2.*nets_demean(sex) ALL_IDPs(:,[12 13 18]) ];
%%% K will index the subset of subjects having full rfMRI data
K=sum(isnan(NET25),2)==0; KN=sum(K);  conf=nets_inormal(conf(K,:)); conf(isnan(conf))=0;  100*KN/N

%%% reduce netmats to top eigenvectors - actually don't, just keep everything!
NET=NET25(K,:);
  NET1=nets_demean(NET);     NET1=NET1/std(NET1(:)); 
  amNET=abs(mean(NET));  NET3=nets_demean(NET./repmat(amNET,size(NET,1),1));  NET3(:,amNET<0.1)=[]; NET3=NET3/std(NET3(:)); 
  %%%grot25=nets_svds([NET1 NET3],200); % or any other combination
  grot25=[NET1 NET3]; % or any other combination
NET=NET100(K,:);
  NET1=nets_demean(NET);     NET1=NET1/std(NET1(:)); 
  amNET=abs(mean(NET));  NET3=nets_demean(NET./repmat(amNET,size(NET,1),1));  NET3(:,amNET<0.1)=[]; NET3=NET3/std(NET3(:)); 
  %%%grot100=nets_svds([NET1 NET3],300); % or any other combination
  grot100=[NET1 NET3]; % or any other combination

%%% PCA-pre-reduce each modality-group in its own right first, then sort out missing data and do final PCA
grotUa=nets_inormal(ALL_IDPs(K,19:73));  % Gaussianise, but no pre-reduction for the smaller modalities
grot=nets_inormal(ALL_IDPs(K,74:end)); grotK=max(isnan(grot),[],2)==0; [grotU,~,~]=nets_svds(grot(grotK,:),80); % dMRI to dim=80
grotUb=zeros(sum(K),80)/0; grotUb(grotK,:)=grotU; 
[grotUc,~,~]= nets_svds(nets_inormal([NODEamps25(K,:) NODEamps100(K,:)]),20); % node amplitudes to dimensionality 20
[grotUd,~,~]= nets_svds([grot25 grot100],80); % netmats to dim=80
%%% following line originally had a bug in, with SVD outputs having very small stddev
%%% fortunately that makes little difference to the final ICA outputs
%NETd=[grotUa grotUb grotUc grotUd]; % original version
NETd=nets_normalise([grotUa grotUb grotUc grotUd]);
%%% deconfound
for i=1:size(NETd,2)
  grot=~isnan(NETd(:,i)); grotconf=nets_demean(conf(grot,:)); NETd(grot,i)=nets_demean(NETd(grot,i)-grotconf*(pinv(grotconf)*NETd(grot,i)));
end
%%% sort out missing data problem
NETdCOV=zeros(size(NETd,1));
for i=1:size(NETd,1)
  for j=i:size(NETd,1)
    grot=NETd([i j],:); grot=grot(:,sum(isnan(grot))==0); NETdCOV(i,j)=grot(1,:)*grot(2,:)'/size(grot,2); NETdCOV(j,i)=NETdCOV(i,j);
  end
end
NETdCOV2=nearestSPD(NETdCOV);  corr(NETdCOV(:),NETdCOV2(:))  % scatter(NETdCOV(:),NETdCOV2(:));
%%% compute IDP vertical eigenvectors for feeding into CCA
[uu1,dd]=eigs(NETdCOV2,Nkeep);    [~,grot]=eig(NETdCOV2); grot=flipud(diag(grot)); 100*sum(grot(1:Nkeep))/sum(grot)   % 98.5

%%% now regenerate NETd more completely (not reducing IDPs with any PCA)
NETdraw=nets_inormal([ALL_IDPs(K,19:end) NODEamps25(K,:) NODEamps100(K,:) NET25(K,:) NET100(K,:)]); NETd=NETdraw; NETdAGE=NETdraw;
for i=1:size(NETd,2)
  grot=~isnan(NETdraw(:,i));
  grotconf=nets_demean(conf(grot,:));          NETd(grot,i)=nets_demean(NETdraw(grot,i)-grotconf*(pinv(grotconf)*NETdraw(grot,i)));
  grotconf=nets_demean(conf(grot,[3 6 7 8]));  NETdAGE(grot,i)=nets_demean(NETdraw(grot,i)-grotconf*(pinv(grotconf)*NETdraw(grot,i)));
end

%%% create main working vars matrix and deconfound (and get main missing-data-adjusted subject-eigenvectors)
varsdraw=nets_inormal(vars(K,varskeep));  varsd=varsdraw; varsdAGE=varsdraw;
for i=1:size(varsd,2)
  grot=~isnan(varsd(:,i));
  grotconf=nets_demean(conf(grot,:));         varsd(grot,i)=nets_normalise(varsdraw(grot,i)-grotconf*(pinv(grotconf)*varsdraw(grot,i)));
  grotconf=nets_demean(conf(grot,[3 6 7 8])); varsdAGE(grot,i)=nets_normalise(varsdraw(grot,i)-grotconf*(pinv(grotconf)*varsdraw(grot,i)));
end
varsdCOV=zeros(size(varsd,1));
for i=1:size(varsd,1)
  for j=i:size(varsd,1)
    grot=varsd([i j],:); grot=grot(:,sum(isnan(grot))==0); varsdCOV(i,j)=grot(1,:)*grot(2,:)'/size(grot,2); varsdCOV(j,i)=varsdCOV(i,j);
  end
end
varsdCOV2=nearestSPD(varsdCOV);  corr(varsdCOV(:),varsdCOV2(:))  % scatter(varsdCOV(:),varsdCOV2(:));
[uu2,dd]=eigs(varsdCOV2,Nkeep);    [~,grot]=eig(varsdCOV2); grot=flipud(diag(grot)); 100*sum(grot(1:Nkeep))/sum(grot)   % 98.5

%%% combine across IDPs and non-IDP-vars
varsALL=[NETd varsd];
varsALLraw=[NETdraw varsdraw];

%%% do the CCA and permutation testing
Nperm=1000;
[ccaA,ccaB,ccaR,ccaU,ccaV,grotstats]=canoncorr(uu1,uu2); ccaR=[ccaR mean(ccaR)]
ccaRp=zeros(Nperm,size(ccaR,2));
for j=1:Nperm
  [ccaAr,ccaBr,ccaRp(j,1:end-1),ccaUr,ccaVr,grotstatsr]=canoncorr(uu1,uu2(randperm(size(uu2,1)),:)); ccaRp(j,end)=mean(ccaRp(j,1:end-1));
end
for i=1:size(ccaR,2);  % show corrected pvalues
  ccaRpval(i)=(1+sum(ccaRp(2:end,1)>=ccaR(i)))/Nperm;
end
ccaRpval
Ncca=sum(ccaRpval<0.05)  % number of significant CCA components
ccaRpval(1:Ncca)

%%% "spatial" ICA fed from top CCA modes
J=9;  ccaUV=ccaU(:,1:J)+ccaV(:,1:J);
ccaAAA=corr(ccaUV,varsALL,'rows','pairwise'); ccaAAA=0.5*log((1+ccaAAA)./(1-ccaAAA));  % OR equivalently:  ccaAAAa=ccaUV'*varsALL;
[icaSdb,icaAdb,icaWdb]=fastica(ccaAAA,'approach','symm','g','tanh','epsilon',1e-13,'maxNumIterations',3000,'lastEig',J);
icaREV=[-1 1 -1 -1 -1 -1 1 1 1];  icaSdb=diag(icaREV)*icaSdb;  % flip signs of some components for convenience of interpretation
%%% get ICA subject weights   Y = ccaU ccaAAA = icaU icaS  =>  icaU = ccaU * ccaAAA * icaS'
icaU=ccaUV(:,1:9) * ccaAAA * icaSdb'; % the right way

%%% univariate correlations between IDPs and SMs
[unicorrR,unicorrP]=corr(NETd,varsd,'rows','pairwise');               % Pearson
[unicorrRdAGE,unicorrPdAGE]=corr(NETdAGE,varsdAGE,'rows','pairwise'); % Pearson with no age-deconfounding
[unicorrRRaw,unicorrPRaw]=corr(NETdraw,varsdraw,'rows','pairwise');   % Pearson with no deconfounding
[FDR_pID,FDR_pN]=FDR(unicorrP(:),0.05)             % use more conservative threshold FDR_pN
[FDR_pIDRaw,FDR_pNRaw]=FDR(unicorrPRaw(:),0.05)    
BONF=0.05/prod(size(unicorrP))                     % even more conservative - Bonferroni

%%% ASPS-study replication
FHTweight=nets_cellfind(varsVARS,'21002-0');
FHTwaist=nets_cellfind(varsVARS,'48-0'); 
FHThip=nets_cellfind(varsVARS,'49-0'); 
FHTwhratio=vars(:,FHTwaist)./vars(:,FHThip);
FHTbmi=nets_cellfind(varsVARS,'21001-0'); 
FHTalc=nets_cellfind(varsVARS,'1558-2');       % alcohol intake frequency (scale needs inverting)
FHTsmok1=nets_cellfind(varsVARS,'2897-0');     % age started smoking 
FHTsmok2=nets_cellfind(varsVARS,'22507-0');    % age stopped smoking 
FHTsmok3=nets_cellfind(varsVARS,'20116-0');    % smoking status (never, previous, current)
%%% exclude FHTsmok1~FHTsmok and FHTweight~FHTbmi (r>.8)  %% and alternately try FHTsmok2 and 3
grotX=[age(K) sex(K) ALL_IDPs(K,[12 13 18]) Hypertension(K) BloodPressures(K,:) HighCholesterol(K) FHTwhratio(K,:) vars(K,[FHTbmi FHTsmok3]) -vars(K,FHTalc)];
[grotR,grotP]=corr(NETdraw(:,IDP_modality_types(19:end)==3),grotX,'rows','pairwise');
clear poop*; j=1; for i=find(IDP_modality_types(19:end)==3)
  [grot1,grot2]=robustfit(nets_normalise(grotX), nets_normalise(NETdraw(:,i)));
  poopT(j,:)=grot1(2:end);  poopP(j,:)=grot2.p(2:end); j=j+1;
end

