MATLAB script for correlation maps: Difference between revisions
Jump to navigation
Jump to search
imported>Psych 204 Created page with '% want to make a parameter map on subjectB's brain that shows the % correlation in each voxel for a particular voxel chosen from subjectA's % brain % load some data: clear all f…' |
imported>Psych 204 No edit summary |
||
(One intermediate revision by the same user not shown) | |||
Line 1: | Line 1: | ||
% want to make a parameter map on subjectB's brain that shows the | % want to make a parameter map on subjectB's brain that shows the | ||
% correlation in each voxel for a particular voxel chosen from subjectA's | % correlation in each voxel for a particular voxel chosen from subjectA's | ||
% brain | % brain | ||
% load some data: | % load some data: | ||
clear all | clear all | ||
format long g | format long g | ||
%set directory for the two subjects to be compared | %set directory for the two subjects to be compared | ||
% make sure you run this file from cjimmy | % make sure you run this file from cjimmy | ||
if ispc | if ispc | ||
subjADir='Y:\New_Localizers\dr021108-1p5mm-3mm'; | subjADir='Y:\New_Localizers\dr021108-1p5mm-3mm'; | ||
subjBDir='Y:\New_Localizers\kw011708-1p5mm-3mm'; | subjBDir='Y:\New_Localizers\kw011708-1p5mm-3mm'; | ||
else | else | ||
% kevin subjects | % kevin subjects | ||
% subjADir='/biac2/kgs/projects/New_Localizers/dr021108-1p5mm-3mm'; | % subjADir='/biac2/kgs/projects/New_Localizers/dr021108-1p5mm-3mm'; | ||
Line 23: | Line 23: | ||
% nathan's prosopagnosic subjects | % nathan's prosopagnosic subjects | ||
% subjADir='/biac2/kgs/projects/Prosopagnosia/fMRIData/112909whFMRI'; | % subjADir='/biac2/kgs/projects/Prosopagnosia/fMRIData/112909whFMRI'; | ||
% subjADir='/biac2/kgs/projects/Prosopagnosia/fMRIData/020209_mr_fmri'; | % subjADir='/biac2/kgs/projects/Prosopagnosia/fMRIData/020209_mr_fmri'; | ||
% subjBDir='/biac2/kgs/projects/Prosopagnosia/fMRIData/52309AF'; | % subjBDir='/biac2/kgs/projects/Prosopagnosia/fMRIData/52309AF'; | ||
% subjADir = | % subjADir = | ||
% '/biac2/kgs/projects/Kids/fmri/localizer/adult_jc_27yo_052408'; | % '/biac2/kgs/projects/Kids/fmri/localizer/adult_jc_27yo_052408'; | ||
subjADir = '/biac2/kgs/projects/Kids/fmri/localizer/adult_al_22yo_051108/'; | |||
% subjBDir = | % subjBDir = | ||
% '/biac2/kgs/projects/Kids/fmri/localizer/adult_kw_25yo_090308'; | % '/biac2/kgs/projects/Kids/fmri/localizer/adult_kw_25yo_090308'; | ||
subjBDir = '/biac2/kgs/projects/Kids/fmri/localizer/adult_dy_25yo_041908'; | subjBDir = '/biac2/kgs/projects/Kids/fmri/localizer/adult_dy_25yo_041908'; | ||
end | end | ||
% %choose an ROI in subject A and B | % %choose an ROI in subject A and B | ||
% roiA='lPPA'; | |||
% roiA='lPPA'; | % roiA='rFFA'; | ||
% roiA='rFFA'; | % alana | ||
% alana | % roiB='AllRois_Combined'; | ||
% roiB='AllRois_Combined'; | |||
% roiA='lPPA_IOvAC_p3'; | |||
% roiA='lPPA_IOvAC_p3'; | roiA='lFFA_MBvAC_p3'; | ||
roiA='lFFA_MBvAC_p3 | |||
roiB = 'allROIS'; | |||
%want to get the motion corrected data type which is usually 3 but might be | |||
%anything depending on the experimenter and experiments. | |||
dt = 3; %'MotionComp_RefScan5'; | |||
%then pick the scans we will use. | |||
% mtloc,eccbias1,eccbias2,4xretinotopy,objloc1,objloc2 | |||
% so to get our scans in order we have | |||
% scan = [8 9 1 2 3]; | |||
scan = [1 2 3 4]; | |||
% .. for subject A brain | |||
%go to subjects data directory | |||
cd(subjADir); | |||
%initialize our hidden inplane set to the right data type,scan(s), and roi | |||
% hiA = initHiddenInplane(dt, scan, roisA); | |||
hiA = initHiddenInplane(dt, scan); | |||
% | % get the roi data from subject 1 | ||
%then we | %then we grab the multivoxel data from the ROI | ||
% | %this should give us mvA which is a struct with lots of useful stuff | ||
% | %in particular | ||
% | %mvA.trials has information about conditions and order of stimuli etc | ||
%mvA. | %mvA.coord has the coordinates for each voxel in the roi should be 3 x | ||
% | %numVoxels | ||
%for | %mvA.tSeries has the time series for each voxel in % signal change | ||
% | %mvA.roi had roi info | ||
mvA = mv_init(hiA, roiA, scan, dt,1); | |||
%then we apply our glm to each of the voxels in the roi | |||
%we will probably want to set these parameters ourselves so that we know | |||
%what happened | |||
%this gives us | |||
%mvA.glm has the betas, the design matrix and so on | |||
%for a block design there will be as many betas as there are conditions | |||
%for an event related design there will be timepoints x conditions x voxels | |||
%betas | |||
% mvA{i} = mv_applyGlm(mvA{i}); | |||
% | % now do the same for subject B | ||
% | % set scan numbers for B | ||
% | % scan=[1 2 3 4 5]; | ||
scan=[1 2 3 4]; | |||
% .. for all ROIs in subj B | |||
% instead invoke hidden inplane with multiple rois. then loop through the | |||
%rois when comparing to subjectA | |||
%go to subject b directory | |||
cd(subjBDir); | |||
% get dataTpe for subjectB | |||
dt=3; | |||
%load data with all rois | |||
hiB = initHiddenInplane(dt,scan,roiB); | |||
mvB = mv_init(hiB, roiB, scan, dt, 1); | |||
% mvB{i} = mv_applyGlm(mvB{i}); | |||
% can get the size of the brain map in the view for a particular scan by | |||
% | % doing | ||
sz=dataSize(hiB,1); | |||
% | % which should be x,y, slice | ||
% | % when we use mv_init on our rois we end up with coordinates in the field | ||
% mvA.roi.coords that have a particular numbering convention. we want to | |||
% transform that into the space that the view uses which we hope is done | |||
% with the function | |||
% | % so we pick an arbitrary voxel from subject A's roi | ||
voxelInA = 10; | |||
% then we compute the correlation with each voxel in subject B's roi | |||
% vector for the correlations | |||
corrs=zeros(size(mvA.tSeries,2),size(mvB.tSeries,2)); | % vector for the correlations | ||
coordinates=zeros(3,size(mvB.roi.coords,1)); | corrs=zeros(size(mvA.tSeries,2),size(mvB.tSeries,2)); | ||
for j = 1:size(mvA.tSeries, 2) | coordinates=zeros(3,size(mvB.roi.coords,1)); | ||
for j = 1:size(mvA.tSeries, 2) | |||
for k = 1:size(mvB.tSeries, 2); %number of voxels in the ROI in B | for k = 1:size(mvB.tSeries, 2); %number of voxels in the ROI in B | ||
%get correlation between the time series of the two voxels | %get correlation between the time series of the two voxels | ||
Line 146: | Line 144: | ||
end | end | ||
% | % | ||
end | end | ||
% get an average correlation | % get an average correlation | ||
meancorr = mean(corrs,1); | meancorr = mean(corrs,1); | ||
% function [roiInd, coords] = roiIndices(view, coords, preserveCoords); | % function [roiInd, coords] = roiIndices(view, coords, preserveCoords); | ||
% in this case we think it should be | % in this case we think it should be | ||
[roiInd, coords] = roiIndices(hiB, mvB.roi.coords, 1); | [roiInd, coords] = roiIndices(hiB, mvB.roi.coords, 1); | ||
% [oldcoords viewcoords] = roiIndices(view, ourroi, 0) | % [oldcoords viewcoords] = roiIndices(view, ourroi, 0) | ||
% with the coords variable containing the mapping to the inplane. | % with the coords variable containing the mapping to the inplane. | ||
% then we propose to make an array of the same size as the view | % then we propose to make an array of the same size as the view | ||
ourmap = zeros(sz); | ourmap = zeros(sz); | ||
% then we will assign our correlations for each voxel to its location in | % then we will assign our correlations for each voxel to its location in | ||
% the matrix | % the matrix | ||
for i = 1:size(coords,2) | for i = 1:size(coords,2) | ||
ourmap(coords(1,i),coords(2,i),coords(3,i))=meancorr(i); | |||
end | end | ||
map = {ourmap,[]}; | map = {ourmap,[]}; | ||
mapName = 'cor_meanallFFA_KWall'; | mapName = 'cor_meanallFFA_KWall'; | ||
mapUnits = 'r'; | mapUnits = 'r'; | ||
co ={[],[]}; | co ={[],[]}; | ||
pth = ['Inplane/GLMs/' mapName]; | pth = ['Inplane/GLMs/' mapName]; | ||
save(pth, 'map', 'mapName', 'mapUnits'); | |||
% save our parameter map | % save our parameter map | ||
% -> load parameter map in mrvista and see what comes out. | % -> load parameter map in mrvista and see what comes out. |
Latest revision as of 09:24, 9 December 2009
% want to make a parameter map on subjectB's brain that shows the % correlation in each voxel for a particular voxel chosen from subjectA's % brain
% load some data: clear all format long g
%set directory for the two subjects to be compared % make sure you run this file from cjimmy if ispc subjADir='Y:\New_Localizers\dr021108-1p5mm-3mm'; subjBDir='Y:\New_Localizers\kw011708-1p5mm-3mm'; else % kevin subjects % subjADir='/biac2/kgs/projects/New_Localizers/dr021108-1p5mm-3mm'; % % subjBDir='/biac2/kgs/projects/New_Localizers/kw011708-1p5mm-3mm'; % subjBDir='/biac2/kgs/projects/New_Localizers/kgs020408-1p5mm-3mm';
% location of code home = '/biac2/kgs/projects/SameDifferent/'
% nathan's prosopagnosic subjects % subjADir='/biac2/kgs/projects/Prosopagnosia/fMRIData/112909whFMRI'; % subjADir='/biac2/kgs/projects/Prosopagnosia/fMRIData/020209_mr_fmri'; % subjBDir='/biac2/kgs/projects/Prosopagnosia/fMRIData/52309AF'; % subjADir = % '/biac2/kgs/projects/Kids/fmri/localizer/adult_jc_27yo_052408'; subjADir = '/biac2/kgs/projects/Kids/fmri/localizer/adult_al_22yo_051108/'; % subjBDir = % '/biac2/kgs/projects/Kids/fmri/localizer/adult_kw_25yo_090308'; subjBDir = '/biac2/kgs/projects/Kids/fmri/localizer/adult_dy_25yo_041908'; end
% %choose an ROI in subject A and B % roiA='lPPA'; % roiA='rFFA'; % alana % roiB='AllRois_Combined';
% roiA='lPPA_IOvAC_p3'; roiA='lFFA_MBvAC_p3';
roiB = 'allROIS';
%want to get the motion corrected data type which is usually 3 but might be %anything depending on the experimenter and experiments. dt = 3; %'MotionComp_RefScan5'; %then pick the scans we will use. % mtloc,eccbias1,eccbias2,4xretinotopy,objloc1,objloc2 % so to get our scans in order we have % scan = [8 9 1 2 3]; scan = [1 2 3 4];
% .. for subject A brain %go to subjects data directory cd(subjADir); %initialize our hidden inplane set to the right data type,scan(s), and roi % hiA = initHiddenInplane(dt, scan, roisA); hiA = initHiddenInplane(dt, scan);
% get the roi data from subject 1
%then we grab the multivoxel data from the ROI %this should give us mvA which is a struct with lots of useful stuff %in particular %mvA.trials has information about conditions and order of stimuli etc %mvA.coord has the coordinates for each voxel in the roi should be 3 x %numVoxels %mvA.tSeries has the time series for each voxel in % signal change %mvA.roi had roi info mvA = mv_init(hiA, roiA, scan, dt,1);
%then we apply our glm to each of the voxels in the roi %we will probably want to set these parameters ourselves so that we know %what happened %this gives us %mvA.glm has the betas, the design matrix and so on %for a block design there will be as many betas as there are conditions %for an event related design there will be timepoints x conditions x voxels %betas % mvA{i} = mv_applyGlm(mvA{i});
% now do the same for subject B % set scan numbers for B % scan=[1 2 3 4 5]; scan=[1 2 3 4];
% .. for all ROIs in subj B % instead invoke hidden inplane with multiple rois. then loop through the %rois when comparing to subjectA %go to subject b directory cd(subjBDir); % get dataTpe for subjectB dt=3; %load data with all rois hiB = initHiddenInplane(dt,scan,roiB); mvB = mv_init(hiB, roiB, scan, dt, 1); % mvB{i} = mv_applyGlm(mvB{i});
% can get the size of the brain map in the view for a particular scan by % doing
sz=dataSize(hiB,1); % which should be x,y, slice
% when we use mv_init on our rois we end up with coordinates in the field % mvA.roi.coords that have a particular numbering convention. we want to % transform that into the space that the view uses which we hope is done % with the function
% so we pick an arbitrary voxel from subject A's roi voxelInA = 10;
% then we compute the correlation with each voxel in subject B's roi
% vector for the correlations corrs=zeros(size(mvA.tSeries,2),size(mvB.tSeries,2)); coordinates=zeros(3,size(mvB.roi.coords,1)); for j = 1:size(mvA.tSeries, 2) for k = 1:size(mvB.tSeries, 2); %number of voxels in the ROI in B %get correlation between the time series of the two voxels r=corrcoef(mvA.tSeries(:,j), mvB.tSeries(:,k)); %records the correlations in a % record correlations corrs(j,k) = r(2); % in order to get the same number of timeseries as anatomical % coordinates must set preserveCoord variable in mv_init to be 1 % some time series will be duplicated but we don't care % coordinates(:,k) = mvB.roi.coords(:,k)
end % end
% get an average correlation
meancorr = mean(corrs,1);
% function [roiInd, coords] = roiIndices(view, coords, preserveCoords);
% in this case we think it should be [roiInd, coords] = roiIndices(hiB, mvB.roi.coords, 1);
% [oldcoords viewcoords] = roiIndices(view, ourroi, 0) % with the coords variable containing the mapping to the inplane.
% then we propose to make an array of the same size as the view ourmap = zeros(sz);
% then we will assign our correlations for each voxel to its location in % the matrix for i = 1:size(coords,2) ourmap(coords(1,i),coords(2,i),coords(3,i))=meancorr(i); end
map = {ourmap,[]}; mapName = 'cor_meanallFFA_KWall'; mapUnits = 'r'; co ={[],[]};
pth = ['Inplane/GLMs/' mapName]; save(pth, 'map', 'mapName', 'mapUnits'); % save our parameter map
% -> load parameter map in mrvista and see what comes out.