diff --git a/functions/03_Clustering/ConsensusClustering.m b/functions/03_Clustering/ConsensusClustering.m
index 4083034..79ad660 100755
--- a/functions/03_Clustering/ConsensusClustering.m
+++ b/functions/03_Clustering/ConsensusClustering.m
@@ -1,225 +1,225 @@
%% This function performs consensus clustering over a range of K values
% The goal is to provide a measure of how good each value of K is
%
% Inputs:
% - X is the data matrix (n_DP x n_DIM)
% - K_range is the range of K values to examine
% - Subsample_type defines how subsampling is done: across items (data
% points) if 'items', and across dimensions if 'dims', across subjects (all
% items of one subject) if 'subjects'
% - Subsample_fraction is the fraction of the original data points, or
% dimensions, to keep for a given fold
% - n_folds is the number of folds over which to run
% - DistType:
% - subject_labels is the assignment of every item in X to a subject, only
% needed for Subsample_type='subjects'
%
% Outputs:
% saves the results for every K in the consensus subfolder:
% Consensus_ordered - Consensus matrix for every fold (n_items x n_items)
%
%
% 14.6.2017 - Daniela:
% created based on Thomas' version, added subject subsampling,
% changed outputs to structs
% 29.05.2018 - Daniela:
% changed data saving, updated Build_Connectivity_Matrix, removed
% dimensions sampling
function ConsensusClustering(X,subject_labels,param)
K_range=param.K_vect;
n_folds=param.cons_n_folds;
n_rep=param.n_folds;
DistType=param.DistType;
Subsample_fraction=param.Subsample_fraction;
Subsample_type=param.Subsample_type;
if isfield(param,'MaxIter')
MaxIter=param.MaxIter;
else
MaxIter=100;
end
outDir_cons=param.outDir_cons;
% Number of data points
n_items = size(X,1);
disp(['n_items ',num2str(n_items),'...']);
% Number of dimensions
n_dims = size(X,2);
disp(['n_dims ',num2str(n_dims),'...']);
% Loop over all K values to assess
for k = 1:length(K_range)
disp(['Running consensus clustering for K = ',num2str(K_range(k)),'...']);
if exist(fullfile(outDir_cons,['consensusResults_' num2str(K_range(k)) '.mat']), 'file') && ...
~param.force_ConsensusClustering
disp('consensus clustering already done, skipping ...')
continue;
end
if exist(fullfile(outDir_cons,['Consensus_' num2str(K_range(k)) '.mat']),'file')
load(fullfile(outDir_cons,['Consensus_' num2str(K_range(k)) '.mat']));
else
% Connectivity matrix that will contain 0s or 1s depending on whether
% elements are clustered together or not
M_sum=zeros(n_items,n_items);
I_sum=zeros(n_items,n_items);
if strcmp(Subsample_type,'dims')
M_sum=zeros(n_dims,n_dims);
I_sum=zeros(n_dims,n_dims);
end
% Loops over all the folds to perform clustering for
for h = 1:n_folds
disp(['Fold ' num2str(h) ':'])
switch Subsample_type
case 'items'
% Number of items to subsample
n_items_ss = floor(Subsample_fraction*n_items);
% Does the subsampling
[~,tmp_ss] = datasample(X,n_items_ss,1,'Replace',false);
% Vector
I_vec = zeros(n_items,1);
I_vec(tmp_ss) = 1;
%Constructs the indicator matrix
I_sum=I_sum+I_vec*I_vec';
X_ss=X(I_vec>0,:);
tmp_ss=[];
% case 'dims'
% error('Subsampling according to dimensions is not implemented yet')
% this part of the routine is not testet, modifications
% may be needed before using it
% % Number of dimensions to subsample
% n_dims_ss = floor(Subsample_fraction*n_dims);
%
% % Does the subsampling
% [X_ss,tmp_ss] = datasample(X,n_dims_ss,2,'Replace',false);
%
% % Vector
% I_vec = zeros(n_dims,1);
% I_vec(tmp_ss) = 1;
%
% %Constructs the indicator matrix
% I_sum=I_sum+I_vec_s*I_vec_s';
case 'subjects'
if ~exist('subject_labels','var') || isempty(subject_labels)
error('subject labels needed for subjects subsampling!')
end
subject_list=unique(subject_labels);
n_subjects=length(subject_list);
% Number of items to subsample
n_subjects_ss = floor(Subsample_fraction*n_subjects);
% Does the subject subsampling
disp('Subject subsampling...'); tic
[subjects_ss,~] = datasample(subject_list,n_subjects_ss,1,'Replace',false);
% Vector
I_vec = zeros(n_items,1);
for iS=1:n_subjects_ss
I_vec(subject_labels==subjects_ss(iS)) = 1;
end
%Constructs the indicator matrix
I_sum=I_sum+I_vec*I_vec';
% subsampled data
disp('Data subsampling...');tic
X_ss=X(I_vec>0,:);
otherwise
errordlg('PROBLEM IN TYPE OF SUBSAMPLING');
end
% Does the clustering (for now, only with k-means), so that IDX
% contains the indices for each datapoint
disp('Clustering ...');tic
IDX = kmeans(X_ss,K_range(k),'Distance',DistType,'Replicates',n_rep,'MaxIter',MaxIter);
clear X_ss
% Builds the connectivity matrix
disp('Buiding connectivity matrix M ...');tic
M_sum=M_sum+Build_Connectivity_Matrix(IDX,find(I_vec>0),Subsample_type,n_items);
clear I_vec
clear X_ss
clear tmp_ss
clear IDX
end
disp('Computing Consensus Matrix ...');tic
% Constructs the consensus matrix for the considered K
Consensus = M_sum./I_sum;
if any(I_sum(:)==0)
warning([num2str(nnz(isnan(Consensus(:)))) ' (' ...
num2str(nnz(isnan(Consensus(:)))/length(Consensus(:))*100) ...
'%) items have not been selected during subsampling, you should increase the number of folds!']);
Consensus(isnan(Consensus))=0;
end
clear M_sum I_sum M I
disp('Saving consensus results (not ordered) ...');
save(fullfile(outDir_cons,['Consensus_' num2str(K_range(k))]),'Consensus','-v7.3');
end
disp('Ordering consensus matrix ...'); tic
tree = linkage(1-Consensus,'average');
toc
% Leaf ordering to create a nicely looking matrix
leafOrder = optimalleaforder(tree,1-Consensus);
toc
% Ordered consensus matrix
Consensus_ordered = Consensus(leafOrder,leafOrder);
toc
clear Consensus leafOrder
% computing CDF and AUC of consensus matrix
[CDF(k,:),AUC(k,:)] = ComputeClusteringQuality(Consensus_ordered,K_range(k));
toc
disp('Saving consensus results ...');
save(fullfile(outDir_cons,['Consensus_ordered_' num2str(K_range(k))]),'Consensus_ordered','-v7.3');
figure;imagesc(Consensus_ordered,[0 1]);colorbar;title(['k= ' num2str(K_range(k))]);
- savefig(fullfile(outDir_cons,['Consensus_ordered_' num2str(K_range(k))]));
+% savefig(fullfile(outDir_cons,['Consensus_ordered_' num2str(K_range(k))]));
print(fullfile(outDir_cons,['Consensus_ordered_' num2str(K_range(k))]), '-dpng','-painters');
close gcf
clearvars Consensus_ordered
end
save(fullfile(outDir_cons,'CDF'),'CDF');
save(fullfile(outDir_cons,'AUC'),'AUC');
figure;hold on;
for k = 1:length(K_range)
pl(k)=plot(0:0.01:1,CDF(k,:),'linewidth',2);
end
legend(pl,[num2str(K_range')],'location','southeast');
title('CDF');
print(fullfile(outDir_cons,'CDF'),'-depsc2','-painters');
figure;hold on;
plot(K_range,AUC,'-o','linewidth',2);
xlabel('K');
title('AUC');
print(fullfile(outDir_cons,'AUC'),'-depsc2','-painters');
end
\ No newline at end of file
diff --git a/functions/03_Clustering/getClusterConsensus.m b/functions/03_Clustering/getClusterConsensus.m
old mode 100755
new mode 100644
diff --git a/functions/03_Clustering/saveRegionTables.m b/functions/03_Clustering/saveRegionTables.m
old mode 100755
new mode 100644
diff --git a/functions/03_Clustering/saveSubjectMaps.m b/functions/03_Clustering/saveSubjectMaps.m
old mode 100755
new mode 100644
diff --git a/functions/04_Regression/computeTemporalCharacteristics.m b/functions/04_Regression/computeTemporalCharacteristics.m
index 36c4083..48012fe 100755
--- a/functions/04_Regression/computeTemporalCharacteristics.m
+++ b/functions/04_Regression/computeTemporalCharacteristics.m
@@ -1,323 +1,323 @@
% computes temporal characteristics of iCAPs
% occurrences of only one frame will not be counted
% Input: TC - nSubjects x 1 cell object
% cells: nClus x nTP time courses per iCAP
% clusteringResults - strcut with clustering results
% .AI_subject_labels
% .subject_labels
% .IDX
% [.scrub_labels] - only needed if param.excludeMotionFrames=1
% param - necessary field:
% .activityThres - threshold at which normalized time courses
% will be considered "active", default = 1
%
%
% Output: tempChar - structure containing fields with temporal
% characteristics (nClus x nSubs)
% * Characteristics of significant innovations
% .innov_counts_total - number of all significant innovations per
% subject
% .innov_counts_total_perc - percentage of significant
% innovations in all included time points
% .innov_counts - nClus x nSub matrix, number of significant
% innovations per iCAP per subject
% .innov_counts_percOfInnov - nClus x nSub matrix, percentage of
% innovations in this iCAP in all significant innovations
%
% * Thresholded time courses and overall characteristics
% (nSubjects x 1 cell objects)
% .TC_norm_thes - normalized and thresholded time courses, with
% occurrences of only one frame removed from the TC
% .TC_active - activity information, 1 if positive activity, -1
% if negative, 0 if none
% .coactiveiCAPs_total - 1 x nTP_subject with number of active
% iCAPs per TP
% .coactivation - nClus x nClus x nTP_subject with coactivation
% between every two iCAPs
%
% * Temporal characteristics of activity blocks in time courses
% (nClus x nSub matrices)
% .occurrences - number of activity blocks in thresholded time
% courses for each cluster in each subject
% .occurrences_pos - number of positive activity blocks in
% thresholded time courses
% .occurrences_neg - number of negative activity blocks in
% thresholded time courses
%
% .duration_total_counts - number of active time points per
% cluster per subject
% .duration_total_pos_counts - number of positively active time
% points per cluster per subject
% .duration_total_neg_counts - number of negatively active time
% points per cluster per subject
%
% .duration_total_perc - total duration of iCAP in percentage of
% the whole scan duration
% .duration_total_pos_perc - total positive duration of iCAP in
% percentage of the whole scan duration
% .duration_total_neg_perc - total negative duration of iCAP in
% percentage of the whole scan duration
%
% .duration_avg_counts - average duration of activity blocks
% (number of time points)
% .duration_avg_pos_counts - average duration of positive
% activity blocks (number of time points)
% .duration_avg_pos_counts - average duration of negative
% activity blocks (number of time points)
%
% * Co-activation characteristics of iCAPs time courses
% (nClus x nClus x nSub)
% .coupling_counts - coupling duration (number of time
% points)
% .coupling_jacc - coupling duration (percentage of
% total duration of both iCAPs)
%
% .coupling_sameSign - same-signed coupling (positive
% coupling) duration (number of time points)
% .coupling_diffSign - differently-signed coupling
% (negative coupling) duration (number of time points)
% .coupling_sameSign_jacc - same-signed coupling
% duration (percentage of total duration of both iCAPs)
% .coupling_diffSign_jacc - differently-signed
% coupling duration (percentage of tota duration of both
% iCAPs)
function tempChar = computeTemporalCharacteristics(TC,clusteringResults,param)
thres=param.activityThres;
AI_subject_labels=clusteringResults.AI_subject_labels;
subject_labels=clusteringResults.subject_labels;
IDX=clusteringResults.IDX;
% constants
nSub=length(TC);
nClus=size(TC{1},1);
for iS=1:nSub
nTP_sub(iS)=size(TC{iS},2);
end
nTP_all=length(AI_subject_labels);
if ~isfield(clusteringResults,'scrub_labels') || isempty(clusteringResults.scrub_labels) ...
- || ~isfield(param,'excludeMotionFrames') || ~param.excludeMotionFrames
+ || ~isfield(param,'excludeMotionFromTC') || ~param.excludeMotionFromTC
scrub_labels=ones(nTP_all,1);
else
scrub_labels=clusteringResults.scrub_labels;
end
% loop over all subjects
for iS = 1:nSub
% scrubbed time points
vols_iS=AI_subject_labels==iS;
tempChar.scrub_labels{iS}=scrub_labels(vols_iS)';
tempChar.nTP_sub(iS,1)=nnz(scrub_labels(vols_iS));
% innovation counts per subject
tempChar.innov_counts_total(1,iS)=nnz(subject_labels==iS); % number of all significant innovations per subject
tempChar.innov_counts_total_perc(1,iS)=tempChar.innov_counts_total(1,iS)/tempChar.nTP_sub(iS)*100; % percentage of significant innovations in all included time points
% normalize time courses
TC_norm{iS}=reshape(zscore(TC{iS}(:)),size(TC{iS},1),size(TC{iS},2));
TC_norm_thes=TC_norm;
TC_norm_thes{iS}(abs(TC_norm{iS})<thres)=0;
% remove occurrences of only one frame
for iC=1:size(TC_norm_thes{iS},1) % I am computing thresholded time courses for all iCAPs first because I will use them for the co-activation computation
activeComp=bwconncomp(TC_norm_thes{iS}(iC,:));
% splitting components with sign changes
activeComp=splitPosNegComps(activeComp,TC_norm_thes{iS}(iC,:));
% deleting occurrences of only one frame
TC_norm_thes{iS}(iC,:)=deleteOneFrameOcc(activeComp,TC_norm_thes{iS}(iC,:));
end
% number of co-active iCAPs per time point
tempChar.TC_norm_thes{iS}=TC_norm_thes{iS};
tempChar.TC_active{iS}=TC_norm_thes{iS};
tempChar.TC_active{iS}(tempChar.TC_active{iS}>0)=1;
tempChar.TC_active{iS}(tempChar.TC_active{iS}<0)=-1;
tempChar.coactiveiCAPs_total{iS}=sum(tempChar.TC_active{iS}~=0);
tempChar.coactiveiCAPs_total{iS}(tempChar.scrub_labels{iS}==0)=nan; % set not included time points to nan
% compute iCAPs-wise characteristics
for iC=1:size(TC_norm_thes{iS},1)
% compute the active components again
activeComp=bwconncomp(TC_norm_thes{iS}(iC,:));
activeComp=splitPosNegComps(activeComp,TC_norm_thes{iS}(iC,:));
% get signs of components
activeComp=getCompSign(activeComp,TC_norm_thes{iS}(iC,:));
for iA=1:activeComp.NumObjects
if length(activeComp.PixelIdxList{iA})<2
error('Something went wrong while suppressing occurrences of only one frame');
end
end
% compute occurrences and average duration
tempChar.occurrences(iC,iS)=activeComp.NumObjects;
tempChar.occurrences_pos(iC,iS)=nnz(activeComp.compSign>0);
tempChar.occurrences_neg(iC,iS)=nnz(activeComp.compSign<0);
tempChar.duration_total_counts(iC,iS)=nnz(TC_norm_thes{iS}(iC,:));
tempChar.duration_total_pos_counts(iC,iS)=nnz(TC_norm_thes{iS}(iC,:)>0);
tempChar.duration_total_neg_counts(iC,iS)=nnz(TC_norm_thes{iS}(iC,:)<0);
% compute duration in percentage of total scan time
tempChar.duration_total_perc(iC,iS)=tempChar.duration_total_counts(iC,iS)/tempChar.nTP_sub(iS)*100;
tempChar.duration_total_pos_perc(iC,iS)=tempChar.duration_total_pos_counts(iC,iS)/tempChar.nTP_sub(iS)*100;
tempChar.duration_total_neg_perc(iC,iS)=tempChar.duration_total_neg_counts(iC,iS)/tempChar.nTP_sub(iS)*100;
tempChar.duration_avg_counts(iC,iS)=tempChar.duration_total_counts(iC,iS)/activeComp.NumObjects;
tempChar.duration_avg_pos_counts(iC,iS)=tempChar.duration_total_pos_counts(iC,iS)/tempChar.occurrences_pos(iC,iS);
tempChar.duration_avg_neg_counts(iC,iS)=tempChar.duration_total_neg_counts(iC,iS)/tempChar.occurrences_neg(iC,iS);
% innovation counts
tempChar.innov_counts(iC,iS)=nnz(IDX==iC&subject_labels==iS);
tempChar.innov_counts_percOfInnov(iC,iS)=tempChar.innov_counts(iC,iS)/tempChar.innov_counts_total(1,iS)*100; % percentage of innovations in this iCAP in all significant innovations
% % innovation counts in time courses
% tempChar.innov_counts_TC(iC,iS)=nnz(diff(TC_norm{iS}(iC,:))~=0)-nnz(diff(tempChar.scrub_labels{iS}~=0)); % significant innovations in time courses, except those due to motion scrubbing
% tempChar.innov_counts_TC_perc(iC,iS)=tempChar.innov_counts_TC(iC,iS)/tempChar.nTP_sub(iS)*100; % percentage of innovations in this iCAP in included time points per subject
%
% % number of co-active iCAPs with given iCAP
% tempChar.coactiveiCAPs(iC,iS,:)=tempChar.coactiveiCAPs_total(iS,:);
% tempChar.coactiveiCAPs(iC,iS,tempChar.TC_active{iS}(iC,:)==0)=nan; % set to nan all time points where the given iCAP was not active
% tempChar.coactiveiCAPs_pos(iC,iS,:)=tempChar.coactiveiCAPs(iC,iS,:);
% tempChar.coactiveiCAPs_pos(iC,iS,tempChar.TC_active{iS}(iC,:)<=0)=nan; % set to nan all time points where the given iCAP was not positive
% tempChar.coactiveiCAPs_neg(iC,iS,:)=tempChar.coactiveiCAPs(iC,iS,:);
% tempChar.coactiveiCAPs_neg(iC,iS,tempChar.TC_active{iS}(iC,:)>=0)=nan; % set to nan all time points where the given iCAP was not negative
%
% compute signed co-activations with all other iCAPs
for iC2=1:nClus
% time points of co-activation of iCAP iC and iCAP iC2
tempChar.coupling{iS}(iC,iC2,:)=TC_norm_thes{iS}(iC,:) & ...
TC_norm_thes{iS}(iC2,:);
% percentage of co-activation with iCAP iC2, with respect to
% total activation of iCAP iC
tempChar.coupling_counts(iC,iC2,iS)=nnz(tempChar.coupling{iS}(iC,iC2,:));
tempChar.coupling_jacc(iC,iC2,iS)=nnz(tempChar.coupling{iS}(iC,iC2,:))/...
nnz(TC_norm_thes{iS}(iC,:)~=0 | TC_norm_thes{iS}(iC2,:)~=0);
% signed co-activation
tempChar.coupling_posPos_counts{iS}(iC,iC2,:)=TC_norm_thes{iS}(iC,:)>0 & ...
TC_norm_thes{iS}(iC2,:)>0;
tempChar.coupling_posNeg_counts{iS}(iC,iC2,:)=TC_norm_thes{iS}(iC,:)>0 & ...
TC_norm_thes{iS}(iC2,:)<0;
tempChar.coupling_negPos_counts{iS}(iC,iC2,:)=TC_norm_thes{iS}(iC,:)<0 & ...
TC_norm_thes{iS}(iC2,:)>0;
tempChar.coupling_negNeg_counts{iS}(iC,iC2,:)=TC_norm_thes{iS}(iC,:)<0 & ...
TC_norm_thes{iS}(iC2,:)<0;
tempChar.coupling_sameSign_counts(iC,iC2,iS)=(nnz(tempChar.coupling_posPos_counts{iS}(iC,iC2,:))+...
nnz(tempChar.coupling_negNeg_counts{iS}(iC,iC2,:)));
tempChar.coupling_diffSign_counts(iC,iC2,iS)=(nnz(tempChar.coupling_posNeg_counts{iS}(iC,iC2,:))+...
nnz(tempChar.coupling_negPos_counts{iS}(iC,iC2,:)));
% percentage of signed co-activation with iCAP iC2, with
% respect to total positive or negative activation of both
% iCAPs
tempChar.coupling_sameSign_jacc(iC,iC2,iS)=(nnz(tempChar.coupling_posPos_counts{iS}(iC,iC2,:))+...
nnz(tempChar.coupling_negNeg_counts{iS}(iC,iC2,:)))/...
nnz(TC_norm_thes{iS}(iC,:)~=0 | TC_norm_thes{iS}(iC2,:)~=0);
tempChar.coupling_diffSign_jacc(iC,iC2,iS)=(nnz(tempChar.coupling_posNeg_counts{iS}(iC,iC2,:))+...
nnz(tempChar.coupling_negPos_counts{iS}(iC,iC2,:)))/...
nnz(TC_norm_thes{iS}(iC,:)~=0 | TC_norm_thes{iS}(iC2,:)~=0);
% percentage of positive or negative couplings (if coupled)
tempChar.coupling_sameSign_perc(iC,iC2,iS)=tempChar.coupling_sameSign_counts(iC,iC2,iS)/tempChar.coupling_counts(iC,iC2,iS)*100;
tempChar.coupling_diffSign_perc(iC,iC2,iS)=tempChar.coupling_diffSign_counts(iC,iC2,iS)/tempChar.coupling_counts(iC,iC2,iS)*100;
if iC==iC2;
tempChar.coupling_counts(iC,iC2,iS)=nan;
tempChar.coupling_jacc(iC,iC2,iS)=nan;
tempChar.coupling_sameSign(iC,iC2,iS)=nan;
tempChar.coupling_diffSign(iC,iC2,iS)=nan;
tempChar.coupling_sameSign_jacc(iC,iC2,iS)=nan;
tempChar.coupling_diffSign_jacc(iC,iC2,iS)=nan;
tempChar.coupling_sameSign_perc(iC,iC2,iS)=nan;
tempChar.coupling_diffSign_perc(iC,iC2,iS)=nan;
end
end
end
% remove nans
tempChar.duration_avg_counts(isnan(tempChar.duration_avg_counts))=0;
tempChar.duration_avg_pos_counts(isnan(tempChar.duration_avg_pos_counts))=0;
tempChar.duration_avg_neg_counts(isnan(tempChar.duration_avg_neg_counts))=0;
end
% % compute measures in seconds
% tempChar.duration_total_sec=tempChar.duration_total_counts.*param.TR;
% tempChar.duration_avg_sec=tempChar.duration_avg_counts.*param.TR;
% tempChar.duration_total_pos_sec=tempChar.duration_total_pos_counts.*param.TR;
% tempChar.duration_avg_pos_sec=tempChar.duration_avg_pos_counts.*param.TR;
% tempChar.duration_total_neg_sec=tempChar.duration_total_neg_counts.*param.TR;
% tempChar.duration_avg_neg_sec=tempChar.duration_avg_neg_counts.*param.TR;
end
function activeComp=splitPosNegComps(activeComp,TC)
for iA=1:activeComp.NumObjects
% split connected components that include a sign change
sign_comp=sign(TC(activeComp.PixelIdxList{iA}));
sign_diff=find(diff(sign_comp));
if ~isempty(sign_diff)
activeComp.NumObjects=activeComp.NumObjects+length(sign_diff);
sign_diff=[sign_diff,length(activeComp.PixelIdxList{iA})];
% add additional components
for iN=1:length(sign_diff)-1
activeComp.PixelIdxList{end+1}=activeComp.PixelIdxList{iA}(sign_diff(iN)+1:sign_diff(iN+1));
end
% update existing component (only keep the first connected
% component)
activeComp.PixelIdxList{iA}=activeComp.PixelIdxList{iA}(1:sign_diff(1));
end
end
end
function TC=deleteOneFrameOcc(activeComp,TC)
% deleting occurrences of only one frame
for iA=1:activeComp.NumObjects
if length(activeComp.PixelIdxList{iA})<2
TC(activeComp.PixelIdxList{iA})=0;
end
end
end
function activeComp=getCompSign(activeComp,TC)
activeComp.compSign=[];
for iA=1:activeComp.NumObjects
% split connected components that include a sign change
sign_comp=sign(TC(activeComp.PixelIdxList{iA}));
sign_diff=find(diff(sign_comp));
if ~isempty(sign_diff)
activeComp=splitPosNegComps(activeComp,TC);
activeComp=getCompSign(activeComp,TC);
continue;
else
activeComp.compSign(iA)=sign_comp(1);
end
end
end
diff --git a/functions/04_Regression/getScrubInfo.m b/functions/04_Regression/getScrubInfo.m
deleted file mode 100755
index bc61734..0000000
--- a/functions/04_Regression/getScrubInfo.m
+++ /dev/null
@@ -1,29 +0,0 @@
-function [num_seg,startID_seg,endID_seg] = getScrubInfo(TM)
-
-
-% compute the number of segments that are together
-diff_TM=[0;diff(TM)]; % first TP is zero to fit to same length as TM
-
-if length(find(diff_TM==-1))==length(find(diff_TM==1))
- if ~TM(1)&&~TM(end) % if moved at beginning AND at the end
- num_seg=length(find(diff_TM==-1));
- startID_seg=find(diff_TM==1);
- endID_seg=find(diff_TM==-1)-1;
- else % if not moved neither at beginning nor at end
- num_seg=length(find(diff_TM==-1))+1;
- startID_seg=[1;find(diff_TM==1)];
- endID_seg=[find(diff_TM==-1)-1;length(TM)];
- end
-else % if the number of activations and deactivations is not the same, one block is at the end or the beginning ot the signal
- num_seg=max(length(find(diff_TM==-1)),length(find(diff_TM==1)));
-
- if ~TM(1) % moved at start
- startID_seg=find(diff_TM==1);
- endID_seg=[find(diff_TM==-1)-1;length(TM)];
- end
-
- if ~TM(end) % moved at end
- startID_seg=[1;find(diff_TM==1)];
- endID_seg=find(diff_TM==-1)-1;
- end
-end
\ No newline at end of file
diff --git a/functions/04_Regression/getTemporalCharacteristics.m b/functions/04_Regression/getTemporalCharacteristics.m
deleted file mode 100755
index 3a9c659..0000000
--- a/functions/04_Regression/getTemporalCharacteristics.m
+++ /dev/null
@@ -1,157 +0,0 @@
-%% This function computes the temporal characteristics (duration, occurences, etc.)
-% given a set of time courses
-%
-% Inputs:
-% - TC: cell array with iCAPs time courses for every subject
-% - clusteringResults: structure containing clustering results, fields:
-% .subject_labels
-% .time_labels
-% [.scrub_labels]
-% - param: structure with input parameters
-% .TR: TR
-%
-% Outputs:
-% - out: structure containing all computed characteristics:
-% .durations_total_n_subs: total number of frames an iCAP is active
-% in every subject
-% .durations_total_s_subs: total time (in seconds) an iCAP is active
-% in every subject
-%
-% .durations_avg(_n/_s): average duration of an iCAPs activity
-
-function [tempChar] = getTemporalCharacteristics(TC,clusteringResults,param)
-
- nSub=length(TC);
- for iS=1:nSub
- nTP_sub(iS)=nnz(clusteringResults.AI_subject_labels==iS);
- end
- subject_labels=clusteringResults.subject_labels;
- IDX=clusteringResults.IDX;
-
- % adding option of scrubbed time courses
- if isfield(param,'excludeMotionFrames') && param.excludeMotionFrames
- scrub_labels=iCAPsResults.scrub_labels;
- else
- scrub_labels=ones(nTP_sub*nSub,1);
- end
-
-
- tempChar = computeTemporalCharacteristics(iCAPsResults.Time_Courses,1,param,scrub_labels,subject_labels,IDX);
-% plotTCandInnov(iCAPsResults,param);
-
-
-end
-
-
-%%
-function plotTCandInnov(iCAPsResults,param)
- %% constants
- nSub=length(iCAPsResults.Time_Courses);
- nClus=size(iCAPsResults.Time_Courses{1},1);
- nTP_sub=size(iCAPsResults.Time_Courses{1},2);
-
- % adding option of scrubbed time courses
- if isfield(param,'excludeMotionFrames') && param.excludeMotionFrames
- scrub_labels=iCAPsResults.scrub_labels;
- else
- scrub_labels=ones(nTP_sub*nSub,1);
- end
-
- for iS=1:nSub
- vols_iS=(iS-1)*nTP_sub+1:iS*nTP_sub;
- nTP_sub_scrub(iS,1) = nnz(scrub_labels(vols_iS));
- TemporalMask{iS}=scrub_labels(vols_iS);
- end
-
-
- %% get data and iCAPs output directories
- if isfield(param,'data_title') % in case of new data saving structure (clustering in subfolders)
- % there is a main folder according to data+thresholding
- % and a separate one for clustering
- outDir_main=fullfile(param.PathData,'iCAPs_results',[param.data_title,'_',param.thresh_title]);
- outDir_iCAPs=fullfile(param.PathData,'iCAPs_results',[param.data_title,'_',param.thresh_title],param.iCAPs_title);
- else
- % before, there was a separate folder for
- % thresholding/data/clustering parameters
- outDir_main=fullfile(param.PathData,'iCAPs_results',[param.thresh_title,'_',param.iCAPs_title]);
- outDir_iCAPs=fullfile(param.PathData,'iCAPs_results',[param.thresh_title,'_',param.iCAPs_title]);
- end
-
-
- %% plotting time courses and innovation indices
- if ~exist(fullfile(outDir_iCAPs,'Time_Courses_and_Innovations'),'dir'); mkdir(fullfile(outDir_iCAPs,'Time_Courses_and_Innovations'));end;
- if ~exist(fullfile(outDir_iCAPs,'Time_Courses'),'dir'); mkdir(fullfile(outDir_iCAPs,'Time_Courses'));end;
-
- iCAPsResults.time_labels(iCAPsResults.time_labels>195)=-(iCAPsResults.time_labels(iCAPsResults.time_labels>195)-195);
-
- for iS = 1:nSub
- % normalize time courses
- iCAPsResults.Time_Courses_plot{iS,1}=reshape(zscore(iCAPsResults.Time_Courses{iS}(:)),size(iCAPsResults.Time_Courses{iS},1),size(iCAPsResults.Time_Courses{iS},2));
- end
-
-
- colors = cbrewer('div', 'Spectral', 10);
- cmap=cbrewer('div', 'RdYlBu',100);
- cmap=cmap(100:-1:1,:);
- % colors = [colors;colors([11:-1:1],:)];
- set(groot,'defaultAxesColorOrder',colors);
- set(groot,'defaultAxesFontSize',15);
-
- clims=[-3,3];
-
- for iS=1:10:nSub
- figure('position',[440 378 515 420]);
- imagesc(iCAPsResults.Time_Courses_plot{iS},clims);
- title(['subject ' num2str(iS)])
- colormap(cmap)
- c=colorbar;
- c.Label.String='amplitude';
- xlabel('time [frames]');
- % ylabel('iCAP');
- print(fullfile(outDir_iCAPs,'Time_Courses',['SUB' num2str(iS)]),'-dpng','-painters');
-
-
- for iC=1:nClus
- fig=figure('Position',[440 573 560 167]);
- hold on;
- grid on
- stemAmp=max(abs(iCAPsResults.Time_Courses_plot{iS,1}(iC,:)))/2;
- stemAmp=2;
- patch('Vertices',[0 -1; nTP_sub -1; nTP_sub 1; 0 1],'Faces',[1 2 3 4],'facecolor',0.5*[1 1 1],'edgeColor',0.5*[1 1 1],'FaceAlpha',0.3,'EdgeAlpha',0.3);
- pl(1)=stem(abs(iCAPsResults.time_labels(iCAPsResults.subject_labels==iS&iCAPsResults.IDX==iC)),...
- stemAmp*sign(iCAPsResults.time_labels(iCAPsResults.subject_labels==iS&iCAPsResults.IDX==iC)),'k')
- plot([0 200],[0 0],'k-')
- pl(2)=plot(iCAPsResults.Time_Courses_plot{iS}(iC,:),'linewidth',1.5,'color',colors(3,:));
- % title(['subject ' num2str(iS)])
- xlabel('time [frames]')
- % ylabel('amplitude')
- % set(gca,'ylim',2.5*[-stemAmp stemAmp]);
- set(gca,'ytick',[-20:2:20],'ylim',[min(-4,1.1*min(iCAPsResults.Time_Courses_plot{iS,1}(iC,:))) max(4,1.1*max(iCAPsResults.Time_Courses_plot{iS,1}(iC,:)))]);
- if ~exist(fullfile(outDir_iCAPs,'Time_Courses_and_Innovations'),'dir'); mkdir(fullfile(outDir,'Time_Courses_and_Innovations'));end;
- outFileName=fullfile(outDir_iCAPs,'Time_Courses_and_Innovations',['SUB' num2str(iS) ' - iCAP' num2str(iC)]);
- % figure(fig);
- % legend(pl,{'',''},'location','eastoutside');
- set(gca,'xlim',[1 195])
- print(outFileName,'-depsc2','-painters');
- close gcf
-
-
-
-% stemAmp=max(abs(iCAPsResults.Time_Courses_plot{iS}(iC,:)))/2;
-% fig=figure('position',[440 662 560 136]);
-% imagesc([1 nTP_sub]+1,[-5 5],TemporalMask{iS}','alphadata',0.1,[0,1]);colormap('gray');
-% hold on;
-% stem(abs(iCAPsResults.time_labels(iCAPsResults.subject_labels==iS&iCAPsResults.IDX==iC)),...
-% stemAmp*sign(iCAPsResults.time_labels(iCAPsResults.subject_labels==iS&iCAPsResults.IDX==iC)))
-% stairs(iCAPsResults.Time_Courses_plot{iS}(iC,:),'linewidth',1.5);
-% title(['SUB' num2str(iS) ' - iCAP' num2str(iC)])
-% xlabel('frame')
-% % ylabel('amplitude')
-% set(gca,'ylim',2.5*[-stemAmp stemAmp]);
-% outFileName=fullfile(outDir_iCAPs,'Time_Courses_and_Innovations',['SUB' num2str(iS) ' - iCAP' num2str(iC)]);
-% figure(fig);
-% print(outFileName,'-depsc2','-painters');
-% close gcf
- end
- end
-end