source: plugins/fmri_plugins/network_analysis.m

Last change on this file was 163, checked in by tjniskan, 8 years ago
  • Minor changes to partial correlation analysis

M misclib/pcorr.m
M plugins/fmri_plugins/network_analysis.m
M aedes_revision.m

File size: 2.8 KB
Line 
1function network_analysis(DATA,ROI,AddInfo)
2%
3% This Aedes plugin calculates and displays partial correlations from ROIs
4%
5
6
7
8if AddInfo.isDataMixed
9  % Make data a 4D-matrix if it is loaded in aedes as a stack of 2D or 3D images
10  DATA2{1}=DATA{1};
11  DATA2{1}.FTDATA = zeros([size(DATA{1}.FTDATA,1),...
12    size(DATA{1}.FTDATA,2),size(DATA{1}.FTDATA,3),length(DATA)],...
13    class(DATA{1}.FTDATA));
14  for ii=1:length(DATA)
15    DATA2{1}.FTDATA(:,:,:,ii)=DATA{ii}.FTDATA;
16  end
17  DATA=DATA2;
18  CurrentVol=AddInfo.CurrentSlice;
19else
20  CurrentVol=AddInfo.CurrentVol;
21end
22
23CurrentVol=AddInfo.CurrentVol;
24
25% Check that there are ROIs defined at all
26if isempty(ROI)
27  % No ROIs defined, nothing to calculate
28  errordlg('No ROIs defined. At least three ROIs must be defined for calculating partial correlations.',...
29    'ROIs not defined.','modal');
30  return
31end
32
33% Check that at least 3 ROIs have been defined
34if length(ROI)<3
35        errordlg('At least three ROIs must be defined for calculating partial correlations.',...
36    'Too few ROIs','modal');
37  return
38end
39
40% Check that there are ROIs defined in this volume
41RoiInds = [];
42for ii=1:length(ROI)
43  if any(any(any(ROI(ii).voxels{1}(:,:,:,CurrentVol))))
44    RoiInds(end+1)=ii;
45  end
46end
47
48% Return if there were no ROIs in the current volume
49if isempty(RoiInds) || length(RoiInds)<3
50  errordlg('Too few ROIs defined in the current volume. At least three ROIs must be defined for calculating partial correlations.',...
51    'ROIs not defined in current volume.','modal');
52  return
53end
54
55X = zeros(size(DATA{1}.FTDATA,4),length(RoiInds));
56
57% Get time series from ROIs
58for kk=RoiInds
59        ind = repmat(ROI(kk).voxels{1}(:,:,:,CurrentVol),[1 1 1 size(DATA{1}.FTDATA,4)]);
60       
61        % Mean EPI time series
62        ts_data = reshape(double(DATA{1}.FTDATA(ind)),[],size(DATA{1}.FTDATA,4));
63        ts_data = ts_data.';
64        ts_data = mean(detrend(ts_data),2);%+repmat(mean(ts_data),size(ts_data,1),1),2);
65        %data_trend=aedes_trendest(double(ts_data),10000);
66        X(:,kk) = ts_data;
67end
68
69% Calculate partial correlations
70[PC,P,CC] = pcorr(X);
71
72% Linear indexes to unique correlations
73ind = find(tril(PC,-1)~=0).';
74
75% Display results
76fprintf('***********************************\n');
77fprintf('* Networks between ROIs\n')
78fprintf('***********************************\n');
79for ii=ind
80        [I,J] = ind2sub(size(PC),ii);
81        if P(ii)<=0.01
82                fprintf(2,'%s <-> %s (PCC=%.4f, CC=%.4f, PCC/CC=%.04f, P=%.4f)**\n',...
83                        ROI(J).label,ROI(I).label,PC(ii),CC(ii),PC(ii)/CC(ii),P(ii));
84        elseif P(ii) <= 0.05
85                fprintf(2,'%s <-> %s (PCC=%.4f, CC=%.4f, PCC/CC=%.04f, P=%.4f)*\n',...
86                        ROI(J).label,ROI(I).label,PC(ii),CC(ii),PC(ii)/CC(ii),P(ii));
87        else
88                fprintf(1,'%s <-> %s (PCC=%.4f, CC=%.4f, PCC/CC=%.04f, P=%.4f)\n',...
89                        ROI(J).label,ROI(I).label,PC(ii),CC(ii),PC(ii)/CC(ii),P(ii));
90        end
91end
92fprintf('***********************************\n\n');
93
94
95
96
Note: See TracBrowser for help on using the repository browser.

Powered by Trac 1.0.9.Copyright © Juha-Pekka Niskanen 2008