diff --git a/Analysis/Run_Clustering_Sim.m b/Analysis/Run_Clustering_Sim.m new file mode 100644 index 0000000..3b4be95 --- /dev/null +++ b/Analysis/Run_Clustering_Sim.m @@ -0,0 +1,117 @@ +function [CP2,Disp,Std_Clusters,idx,d,sfrac] = Run_Clustering_Sim(XONn,n_clusters,mask,brain_info,maskP,maskN,n_rep,idx_sep_seeds,SeedType) + + % Number of seeds + n_seeds = size(idx_sep_seeds,3); + + % Number of subjects + n_subjects = size(idx_sep_seeds,2); + + % Number of possible seed combinations + switch n_seeds + case 1 + n_combos = 1; + combvec = [1]; + case 2 + n_combos = 3; + combvec = [1 0; 0 1; 1 1]; + case 3 + n_combos = 7; + combvec = [1 0 0; 0 1 0; 0 0 1; 1 1 0; 0 1 1; 1 0 1; 1 1 1]; + otherwise + errordlg('PROBLEM AT SEED FRACTIONS'); + end + + % Will contain the fraction of frames linked to a given seed (if using + % the 'Intersection' method) + sfrac = zeros(n_subjects,n_clusters,n_combos); + + % Rows datapoints, columns variable (so here every 70700 activation is + % a datapoint) + % idx will contain 1462 elements (the index of the cluster to which the + % considered datapoint belongs + [idx,CP] = kmeans(XONn',n_clusters,'distance','correlation','replicates',n_rep,'empty','drop','maxiter',100); + + % idx2counts is of size K (number of clusters) and has the number of + % datapoints classified within a given cluster) + + % disp('idx2counts:'); + idx2counts = histc(idx, 1:max(idx)); + + % Output = Input(IX) + [~,IX] = sort(idx2counts,'descend'); + + % Size Kx70700 (location of each cluster); clusters are put with 'the + % most prominent one first' + CP = CP(IX,:); idx2 = idx; % order by occurrence + + % Changes the datapoint indexes so that they fit the new clusters order + for l=1:max(idx), idx2(idx==IX(l))=l; end + idx=idx2; + + CP2 = zeros(n_clusters,size(CP,2)); + Disp = zeros(1,n_clusters); + Std_Clusters = zeros(size(CP,2),n_clusters); + + % For each cluster index + for l=1:max(idx) + % Averages all data points belonging to one specific cluster, and + % stores the obtained pattern as a cell in CP2 + CP2(l,:) = mean(XONn(:,idx==l),2); %./ ( std(XON(:,idx==l),[],2) ) * sqrt(length(idx==l)); % Liu&Duyn + + % Measure of dispersion within the cluster considered + Disp(l) = mean(corr(CP2(l,:)',XONn(:,idx==l))); + + Std_Clusters(:,l) = std(XONn(:,idx==l),[],2); + end + + % d contains the correlation values of all frames to the CAPs + r = corr(CP2',XONn); + + d = zeros(1,length(idx)); + for k=1:max(idx) + d(idx==k) = r(k,idx==k); + end + + % Added part to compute the fraction of frames assigned to a given seed + % if using the Union option (in which a data point is retained + % as long as at least one seed region becomes significantly (de)active) + if strcmp(SeedType,'Union') + + % idx_all will contain the clustering indices as put on a whole + % temporal scale (time x subjects) + idx_all = zeros(size(idx_sep_seeds,1),n_subjects); + + % Index to properly fill in the matrix by adding up the number of + % frames per subject every time + tmp_loc = 1; + + + for s = 1:n_subjects + tmp = sum(squeeze(idx_sep_seeds(:,s,:)),2); + tmp(tmp >= 1) = 1; + tmp = logical(tmp); + idx_all(tmp,s) = idx(tmp_loc:(tmp_loc+sum(tmp)-1)); + tmp_loc = tmp_loc + sum(tmp); + end + + + + % I will compute my fractions for each seed combination of + % interest; for example, 3 seeds would yield 7 possible + % combinations + for s = 1:n_subjects + for t = 1:size(idx_sep_seeds,1) + + tmp = squeeze(idx_sep_seeds(t,s,:)); + + % If there is at least one seed active or deactive at + % the point of interest, we update the sfrac count at the + % appropriate CAP + if sum(tmp) > 0 + sfrac(s,idx_all(t,s),find(ismember(combvec,tmp','rows'))) = ... + sfrac(s,idx_all(t,s),find(ismember(combvec,tmp','rows'))) + 1; + end + end + end + end +end \ No newline at end of file diff --git a/CAP_PickState.m b/CAP_PickState.m new file mode 100644 index 0000000..70fc7aa --- /dev/null +++ b/CAP_PickState.m @@ -0,0 +1,23 @@ +function [State2] = CAP_PickState(State,TM) + + % Picks a random number between 0 and 1 + T = randi([0,100000])/100000; + + % Our probabilities of interest (knowing the start state) + POI = TM(State,:); + CP = cumsum(POI); + + if T <= CP(1) + State2 = 1; + elseif T <= CP(2) + State2 = 2; + elseif T <= CP(3) + State2 = 3; + elseif T <= CP(4) + State2 = 4; + elseif T <= CP(5) + State2 = 5; + else + errordlg('WTH'); + end +end \ No newline at end of file diff --git a/DefaultData/2019_03_03_BCT/degrees_dir.m b/DefaultData/2019_03_03_BCT/degrees_dir.m index aa77875..087f3f4 100755 --- a/DefaultData/2019_03_03_BCT/degrees_dir.m +++ b/DefaultData/2019_03_03_BCT/degrees_dir.m @@ -1,31 +1,31 @@ function [id,od,deg] = degrees_dir(CIJ) %DEGREES_DIR Indegree and outdegree % % [id,od,deg] = degrees_dir(CIJ); % % Node degree is the number of links connected to the node. The indegree % is the number of inward links and the outdegree is the number of % outward links. % % Input: CIJ, directed (binary/weighted) connection matrix % % Output: id, node indegree % od, node outdegree % deg, node degree (indegree + outdegree) % % Notes: Inputs are assumed to be on the columns of the CIJ matrix. % Weight information is discarded. % % % Olaf Sporns, Indiana University, 2002/2006/2008 % ensure CIJ is binary... -CIJ = double(CIJ~=0); +% CIJ = double(CIJ~=0); % compute degrees id = sum(CIJ,1); % indegree = column sum of CIJ od = sum(CIJ,2)'; % outdegree = row sum of CIJ deg = id+od; % degree = indegree+outdegree diff --git a/Script_ArtificialData.m b/Script_ArtificialData.m new file mode 100644 index 0000000..efb5c40 --- /dev/null +++ b/Script_ArtificialData.m @@ -0,0 +1,55 @@ +%% This script generates artificial data to test the outputs of the toolbox + +% Number of subjects +S = 10; + +% Number of voxels +V = 1000; + +% Number of time points +T = 200; + +% Number of CAPs +K = 4; + +N_seed_voxels = 50; + +% CAP patterns +CAP = [[ones(250,1);zeros(750,1)], [zeros(250,1);ones(250,1);zeros(500,1)],... + [zeros(500,1);ones(250,1);zeros(250,1)], [zeros(750,1);ones(250,1)]]; + +% The first 50 voxels are the seed voxels +CAP = [ones(N_seed_voxels,4);CAP]; +CAP = [zeros(V+N_seed_voxels,1),CAP]; + +% Transition probabilities betwen CAPs and baseline (row/column 1 = +% baseline, the other ones = CAPs +TM = [0.8,0.2,0,0,0; ... + 0,0,1,0,0;... + 0,0,0,0.5,0.5;... + 0.4,0,0,0.6,0;... + 0,0,0,0.2,0.8]; + +FD = zeros(T,S); +mask = logical(ones(V+N_seed_voxels,1)); +brain_info = []; +Seed = logical([ones(N_seed_voxels,1);zeros(V,1)]); + + +% Generates the state in which we are +for s = 1:S + s + State(s,1) = 1; + + for t = 2:T + State(s,t) = CAP_PickState(State(s,t-1),TM); + end +end + +% Generates the associated time courses +for s = 1:S + for t = 1:T + TC{s}(t,:) = CAP(:,State(s,t))'; + end +end + diff --git a/Script_onepop.m b/Script_onepop.m index 884b4a0..0611b5c 100644 --- a/Script_onepop.m +++ b/Script_onepop.m @@ -1,109 +1,109 @@ %% This is an example script to run the CAPs analyses without the GUI % In this script, we assume one population of subjects only - +addpath(genpath(pwd)); %% 1. Loading the data files % Data: cell array, each cell of size n_TP x n_masked_voxels TC = % Mask: n_voxels x 1 logical vector mask = % Header: the header (obtained by spm_vol) of one NIFTI file with proper % data dimension and .mat information brain_info = % Framewise displacement: a n_TP x n_subj matrix with framewise % displacement information FD = % Seed: a n_masked_voxels x n_seed logical vector with seed information Seed = %% 2. Specifying the main parameters % Threshold above which to select frames -T = 1; +T = 0.8; % Selection mode ('Threshold' or 'Percentage') SelMode = 'Threshold'; % Threshold of FD above which to scrub out the frame and also the t-1 and % t+1 frames (if you want another scrubbing setting, directly edit the % code) Tmot = 0.5; % Type of used seed information: select between 'Average','Union' or % 'Intersection' SeedType = 'Average'; % Contains the information, for each seed (each row), about whether to % retain activation (1 0) or deactivation (0 1) time points SignMatrix = [1,0]; % Percentage of positive-valued voxels to retain for clustering Pp = 100; % Percentage of negative-valued voxels to retain for clustering Pn = 100; % Number of repetitions of the K-means clustering algorithm n_rep = 50; % Percentage of frames to use in each fold of consensus clustering Pcc = 80; % Number of folds we run consensus clustering for N = 50; %% 3. Selecting the frames to analyse % Xon will contain the retained frames, and Indices will tag the time % points associated to these frames, for each subject (it contains a % subfield for retained frames and a subfield for scrubbed frames) -[Xon,~,Indices] = CAP_find_activity(TC,seed,T,FD,Tmot,SelMode,SeedType,SignMatrix); +[Xon,~,Indices] = CAP_find_activity(TC,Seed,T,FD,Tmot,SelMode,SeedType,SignMatrix); %% 4. Consensus clustering (if wished to determine the optimum K) % This specifies the range of values over which to perform consensus % clustering: if you want to run parallel consensus clustering processes, % you should feed in different ranges to each call of the function -K_range = 2:20; +K_range = 2:6; % Have each of these run in a separate process on the server =) [Consensus] = CAP_ConsensusClustering(Xon,K_range,'items',Pcc/100,N,'correlation'); % Calculates the quality metrics -[~,Qual] = ComputeClusteringQuality(Consensus,[]); +[~,PAC] = ComputeClusteringQuality(Consensus,[]); % Qual should be inspected to determine the best cluster number(s) % You should fill this with the actual value -K_opt = +K_opt = 4; %% 5. Clustering into CAPs -[CAP,~,~,idx] = Run_Clustering(cell2mat(Xon),... +[CAP,~,~,idx] = Run_Clustering_Sim(cell2mat(Xon),... K_opt,mask,brain_info,Pp,Pn,n_rep,[],SeedType); %% 6. Computing metrics % The TR of your data in seconds -TR = +TR = 2; [ExpressionMap,Counts,Entries,Avg_Duration,Duration,TransitionProbabilities,... From_Baseline,To_Baseline,Baseline_resilience,Resilience,Betweenness,... InDegree,OutDegree,SubjectEntries] = Compute_Metrics_simpler(idx,... Indices.kept.active,Indices.scrubbedandactive,K_opt,TR); \ No newline at end of file