[Eeglablist] Extracting ERP values after MPT clustering
Makoto Miyakoshi
mmiyakoshi at ucsd.edu
Thu Apr 21 10:30:14 PDT 2016
Dear Mohammed,
> Is there any alternative method to extract the ERP values after MPT
clustering or to extract the P values after run statistics by MPT GUI?
I have a utility tool for that. I gave it to Ozgur Balkan for which he
added me to his latest IEEE conf proc paper. If you are interested in
please let me know. I can send it to you privately. However, you MUST edit
it to adjust to your data (hopefully that is straightforward, but you have
to understand the data structure of MPT.)
> About the anatomical regions, do you know how can I do that?
I pasted below a part of code from my alpha version of groupSIFT (this name
is only tentative) plugin. Beware that there are more than 300 lines. No
claim no return, use it at your own risk-kind of thing.
Makoto
% --- Executes on button press in selectFilesButton.
function selectFilesButton_Callback(hObject, eventdata, handles)
% hObject handle to selectFilesButton (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Obtain multiple .set files
[allFiles, workingFolder] = uigetfile('*.set', 'MultiSelect', 'on');
if ~any(workingFolder)
disp('Cancelled.')
return
end
% Obtain user input file name
userInputFilename = get(handles.fileNameEdit, 'String');
% Move to the working folder
cd(workingFolder)
% Load empty dipolePairAndMeasure from eeglabroot/plugins/groupSIFT/+pr:
Thanks Clement!
dipolePairAndMeasureObj = pr.dipolePairAndMeasure;
% Load all .set data located under the working folder
ALLEEG = [];
for n = 1:length(allFiles)
loadName = allFiles{n};
EEG = pop_loadset('filename', loadName,'filepath',workingFolder,
'loadmode', 'info');
[ALLEEG, EEG, CURRENTSET] = eeg_store( ALLEEG, EEG, 0 );
end
% Compute how many IC combinations in all subjects
numICs = zeros(length(ALLEEG),1);
for n = 1:length(ALLEEG)
numICs(n) = size(ALLEEG(1,n).icaweights,1);
end
numIcCombination = sum(numICs.^2);
% Assign unique IDs to dipoles
if get(handles.rpdcButton,'value')==1
timeFreqSize = size(squeeze(ALLEEG(1).CAT.Conn.RPDC( 1,1,:,:)));
else
timeFreqSize = size(squeeze(ALLEEG(1).CAT.Conn.dDTF08(1,1,:,:)));
end
latencies = ALLEEG(1).CAT.Conn.erWinCenterTimes;
frequencies = ALLEEG(1).CAT.Conn.freqs;
dimensionLabels= ALLEEG(1).CAT.Conn.dims;
dipolePairAndMeasureObj.linearizedMeasure = zeros([numIcCombination
timeFreqSize(1)*timeFreqSize(2)]);
dipoleCounter = 0;
counter = 1;
for n = 1:length(ALLEEG)
numberOfSessiondipoles = length(ALLEEG(n).dipfit.model);
dipoleLocation = [];
dipoleResidualVariance = [];
for m = 1:numberOfSessiondipoles
dipoleLocation(m,:) = ALLEEG(n).dipfit.model(m).posxyz(1,:);
dipoleResidualVariance(m) = ALLEEG(n).dipfit.model(m).rv;
end;
dipoleId = (dipoleCounter+1):(dipoleCounter+length(dipoleLocation));
dipoleCounter = dipoleCounter + length(dipoleLocation);
for m = 1:length(dipoleLocation) % from location: confirm it with
ALLEEG(1,1).CAT.Conn: dims: {'var_to' 'var_from' 'freq' 'time'}
for k = 1:length(dipoleLocation) % to location
dipolePairAndMeasureObj.from.location = cat(1,
dipolePairAndMeasureObj.from.location,
ALLEEG(n).dipfit.model(m).posxyz(1,:));
dipolePairAndMeasureObj.from.residualVariance(end+1,1) =
ALLEEG(n).dipfit.model(m).rv;
dipolePairAndMeasureObj.to.location = cat(1,
dipolePairAndMeasureObj.to.location, ALLEEG(n).dipfit.model(k).posxyz(1,:));
if m == k
dipolePairAndMeasureObj.linearizedMeasure(counter,:) =
zeros(timeFreqSize(1)*timeFreqSize(2),1);
else
if get(handles.rpdcButton,'value')==1 % Again,
ALLEEG(1,1).CAT.Conn: dims: {'var_to' 'var_from' 'freq' 'time'}
tmpTimeFreqMatrix = squeeze(ALLEEG(n).CAT.Conn.RPDC(
k,m,:,:));
else
tmpTimeFreqMatrix =
squeeze(ALLEEG(n).CAT.Conn.dDTF08(k,m,:,:));
end
dipolePairAndMeasureObj.linearizedMeasure(counter,:) =
tmpTimeFreqMatrix(:);
end
dipolePairAndMeasureObj.sessionId(counter,1) = n;
dipolePairAndMeasureObj.fromDipoleId(counter,1) = dipoleId(m);
dipolePairAndMeasureObj.toDipoleId(counter,1) = dipoleId(k);
counter = counter + 1;
end
end
end
% Find unique dipoles
uniqueDipole = pr.dipole;
[~,fromId,fromIdReverse] = unique(dipolePairAndMeasureObj.fromDipoleId,
'stable');
% Extract coordinates + residual variance of the unique dipoles
uniqueDipole.location = dipolePairAndMeasureObj.from.location(fromId,:);
uniqueDipole.residualVariance =
dipolePairAndMeasureObj.from.residualVariance(fromId);
% Load headGrid cubes
headGrid = pr.headGrid;
% Set Gaussian smoothing kernel
% FWHM = 2.355*sigma See
https://en.wikipedia.org/wiki/Full_width_at_half_maximum
% 4.2 (FWHM==10mm, 12/23/2015)
% 8.5 (FWHM==20mm, 06/24/2015)
% 12.8 (FWHM==30mm, 09/05/2014)
% 17.1 (FWHM==40mm, 01/06/2015)
% Note that Gaussian is NOT truncated until reaching 150mm apart from
the center
userInputKernelSize = str2num(get(handles.gaussianSizeEdit, 'String'));
standardDeviationOfEstimatedDipoleLocation = userInputKernelSize/2.355; %
this calculates sigma in Gaussian equation
projectionParameter =
pr.projectionParameter(standardDeviationOfEstimatedDipoleLocation);
projectionParameter.numberOfStandardDeviationsToTruncatedGaussaian =
150/standardDeviationOfEstimatedDipoleLocation;
[~,~,gaussianWeightMatrix]=
pr.meanProjection.getProjectionMatrix(uniqueDipole, headGrid,
projectionParameter);
% Define valid anatomical ROIs as EEG sources agreed by Scott and Makoto in
Dec 2015; See below for the list
excludeRoiIdx = [21 22 37:42 71:78]; % 29 30 are insula... better to
include?
includeRoiIdx = setdiff(1:88, excludeRoiIdx);
roiLabels = pr.regionOfInterestFromAnatomy.getAllAnatomicalLabels;
includeRoiLabels = roiLabels(includeRoiIdx);
excludeRoiLabels = roiLabels(excludeRoiIdx);
roiLabels = includeRoiLabels;
% These regions are to be included
% 'Precentral_L'
% 'Precentral_R'
% 'Frontal_Sup_L'
% 'Frontal_Sup_R'
% 'Frontal_Sup_Orb_L'
% 'Frontal_Sup_Orb_R'
% 'Frontal_Mid_L'
% 'Frontal_Mid_R'
% 'Frontal_Mid_Orb_L'
% 'Frontal_Mid_Orb_R'
% 'Frontal_Inf_Oper_L'
% 'Frontal_Inf_Oper_R'
% 'Frontal_Inf_Tri_L'
% 'Frontal_Inf_Tri_R'
% 'Frontal_Inf_Orb_L'
% 'Frontal_Inf_Orb_R'
% 'Rolandic_Oper_L'
% 'Rolandic_Oper_R'
% 'Supp_Motor_Area_L'
% 'Supp_Motor_Area_R'
% 'Frontal_Sup_Medial_L'
% 'Frontal_Sup_Medial_R'
% 'Frontal_Med_Orb_L'
% 'Frontal_Med_Orb_R'
% 'Rectus_L'
% 'Rectus_R'
% 'Insula_L'
% 'Insula_R'
% 'Cingulum_Ant_L'
% 'Cingulum_Ant_R'
% 'Cingulum_Mid_L'
% 'Cingulum_Mid_R'
% 'Cingulum_Post_L'
% 'Cingulum_Post_R'
% 'Calcarine_L'
% 'Calcarine_R'
% 'Cuneus_L'
% 'Cuneus_R'
% 'Lingual_L'
% 'Lingual_R'
% 'Occipital_Sup_L'
% 'Occipital_Sup_R'
% 'Occipital_Mid_L'
% 'Occipital_Mid_R'
% 'Occipital_Inf_L'
% 'Occipital_Inf_R'
% 'Fusiform_L'
% 'Fusiform_R'
% 'Postcentral_L'
% 'Postcentral_R'
% 'Parietal_Sup_L'
% 'Parietal_Sup_R'
% 'Parietal_Inf_L'
% 'Parietal_Inf_R'
% 'SupraMarginal_L'
% 'SupraMarginal_R'
% 'Angular_L'
% 'Angular_R'
% 'Precuneus_L'
% 'Precuneus_R'
% 'Paracentral_Lobule_L'
% 'Paracentral_Lobule_R'
% 'Temporal_Sup_L'
% 'Temporal_Sup_R'
% 'Temporal_Pole_Sup_L'
% 'Temporal_Pole_Sup_R'
% 'Temporal_Mid_L'
% 'Temporal_Mid_R'
% 'Temporal_Pole_Mid_L'
% 'Temporal_Pole_Mid_R'
% 'Temporal_Inf_L'
% 'Temporal_Inf_R'
%
% These regions are to be excluded
% 'Olfactory_L'
% 'Olfactory_R'
% 'Hippocampus_L'
% 'Hippocampus_R'
% 'ParaHippocampal_L'
% 'ParaHippocampal_R'
% 'Amygdala_L'
% 'Amygdala_R'
% 'Caudate_L'
% 'Caudate_R'
% 'Putamen_L'
% 'Putamen_R'
% 'Pallidum_L'
% 'Pallidum_R'
% 'Thalamus_L'
% 'Thalamus_R'
%
% One can visualize these regions by running
visualizeAnatomicalRoiWithNHimasBlobs.m contained by the groupSIFT folder.
% Obtain dipole density and centroids in each ROI
numberOfRegionsOfInterest = length(roiLabels);
dipoleProbabilityInRegion = zeros(uniqueDipole.numberOfDipoles,
numberOfRegionsOfInterest);
roiCentroids = zeros(length(includeRoiLabels),3);
for i=1:numberOfRegionsOfInterest
regionOfInterest(i) = pr.regionOfInterestFromAnatomy(pr.headGrid,
roiLabels{i});
dipoleProbabilityInRegion(:,i) = gaussianWeightMatrix *
regionOfInterest(i).membershipProbabilityCube(headGrid.insideBrainCube);
% compute centroids of ROIs
xCube = regionOfInterest(i).headGrid.xCube;
yCube = regionOfInterest(i).headGrid.yCube;
zCube = regionOfInterest(i).headGrid.zCube;
membershipCube = regionOfInterest(i).membershipCube;
xCentroid = mean(xCube(membershipCube));
yCentroid = mean(yCube(membershipCube));
zCentroid = mean(zCube(membershipCube));
roiCentroids(i,:) = [xCentroid yCentroid zCentroid];
end
% Force the ROI centroids to be symmetric
leftHemisphereCentroids = roiCentroids(1:2:end,:);
xPositiveLeftHemisphereControids = bsxfun(@times, leftHemisphereCentroids,
[-1 1 1]);
meanAbsRoiCentroids =
round((xPositiveLeftHemisphereControids+roiCentroids(2:2:end,:))/2);
symmetricRoiCentroids = zeros(size(roiCentroids));
symmetricRoiCentroids(1:2:end,:) = bsxfun(@times, meanAbsRoiCentroids, [-1
1 1]);
symmetricRoiCentroids(2:2:end,:) = meanAbsRoiCentroids;
% generate subject list
fileNameList = {ALLEEG.filename}';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Compute dipole pair density (non-normalized) %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculate dipole densities in source regions and destination regions.
Reuse indices from 'from' dipole for 'to' dipole
fromDipoleProbabilityInRegion =
zeros(dipolePairAndMeasureObj.from.numberOfDipoles,
numberOfRegionsOfInterest);
toDipoleProbabilityInRegion =
zeros(dipolePairAndMeasureObj.to.numberOfDipoles,
numberOfRegionsOfInterest);
[~,~,toIdReverse] = unique(dipolePairAndMeasureObj.toDipoleId, 'stable');
for i = 1:numberOfRegionsOfInterest
fromDipoleProbabilityInRegion(:,i) =
dipoleProbabilityInRegion(fromIdReverse, i);
toDipoleProbabilityInRegion(:, i) =
dipoleProbabilityInRegion(toIdReverse, i);
end
% Prepare dipole pair density by multiplying from-density and to-density.
The first dimension is the sum of each subject's IC*IC, the second is
(reduced AAL regions)^2==5184.
dipolePairProbabilityOnFromRegionToRegionSubstrate =
zeros(dipolePairAndMeasureObj.from.numberOfDipoles, length(includeRoiIdx) *
length(includeRoiIdx));
for i=1:dipolePairAndMeasureObj.from.numberOfDipoles
dipolePairProbabilityOnFromRegionToRegionSubstrate(i, :) =
vec(fromDipoleProbabilityInRegion(i,:)' *
toDipoleProbabilityInRegion(i,:))';
end
% Calculate dipole pair density for each unique edge
numberOfsessions = max(dipolePairAndMeasureObj.sessionId); % session means
subjects
dipolePairDensity =
zeros(size(dipolePairProbabilityOnFromRegionToRegionSubstrate, 2),
numberOfsessions);
for i=1:size(dipolePairProbabilityOnFromRegionToRegionSubstrate, 2)
dipolePairProbabilityOnFromRegionToRegion =
dipolePairProbabilityOnFromRegionToRegionSubstrate(:,i);
sessionDensity = zeros(1,numberOfsessions);
for j=1:numberOfsessions
sessionDensity(j) =
sum(dipolePairProbabilityOnFromRegionToRegion(dipolePairAndMeasureObj.sessionId
== j));
end
% store dipole pair density
dipolePairDensity(i,:) = sessionDensity;
end
dipolePairDensity = reshape(dipolePairDensity, [length(includeRoiIdx)
length(includeRoiIdx) numberOfsessions]);
% calculate dipole-pair pairwise similarity (could be further optimized by
not calculating pair that do not share dipole-pair density in any ROIs)
correlationSimilaity =
1-squareform(pdist(dipolePairAndMeasureObj.linearizedMeasure',
'correlation'));
linearizedMeasureSimilarity =
estimate_mutual_information_from_correlation(correlationSimilaity);
% save dipolePairDensity and others that is needed for similarity test
save([workingFolder filesep userInputFilename '_dipolePairDensity'],
'dipolePairDensity', 'dipoleProbabilityInRegion', 'roiLabels',
'symmetricRoiCentroids', 'fileNameList')
% save convergenceData
linearizedMeasure = dipolePairAndMeasureObj.linearizedMeasure;
save([workingFolder filesep userInputFilename '_measureConvergence'],
'dipolePairProbabilityOnFromRegionToRegionSubstrate', 'linearizedMeasure',
'linearizedMeasureSimilarity');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Generate individual subject's connectivity map--This part is very slow
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prepare the map
numOfRegionsOfInterest = size(dipolePairDensity,1);
toRegionNumber = repmat((1:numOfRegionsOfInterest)', [1,
numOfRegionsOfInterest]); % first deimension (rows) contains toRegion
Numbers.
toRegionNumber = toRegionNumber(:);
fromRegionNumber = repmat( 1:numOfRegionsOfInterest,
[numOfRegionsOfInterest, 1]); % second deimension (columns) contains
fromRegion Numbers.
fromRegionNumber = fromRegionNumber(:);
% Double check the from and to!!
% ALLEEG(1,1).CAT.Conn: dims: {'var_to' 'var_from' 'freq' 'time'}
% figure; imagesc(toRegionNumber); title('To Indices')
% figure; imagesc(fromRegionNumber); title('From Indices')
if get(handles.rpdcButton,'value')==1
connectivityType = 'rPDC';
else
connectivityType = 'dDTF08';
end
for sessionIdx = unique(dipolePairAndMeasureObj.sessionId)'
effectiveConnectivityTimeFreq =
single(zeros(timeFreqSize(1)*timeFreqSize(2),
size(dipolePairDensity,1)*size(dipolePairDensity,2)));
tic
for fromIdx = 1:max(fromRegionNumber) % up to 72
for toIdx = 1:max(toRegionNumber) % up to 72
currentEdgeIdx = find((fromRegionNumber==fromIdx) &
(toRegionNumber==toIdx));
currentSessionIdx =
find(dipolePairAndMeasureObj.sessionId==sessionIdx);
dipolePairDensitiesInRoiPair =
dipolePairProbabilityOnFromRegionToRegionSubstrate(currentSessionIdx,
currentEdgeIdx);
% 02/10/2016 Makoto. Put it back after discussing with Nima and
Scott (during and after my presentation at lab meeting)
% 12/23/2015 Makoto. Disabled this normalization since this
will inflate near-zero values
normFactor = sum(dipolePairDensitiesInRoiPair);
dipolePairDensitiesInRoiPair = dipolePairDensitiesInRoiPair /
normFactor;
% This is the projected measure
projectedMeasure = bsxfun(@times,
dipolePairAndMeasureObj.linearizedMeasure(currentSessionIdx,:),
dipolePairDensitiesInRoiPair);
effectiveConnectivityTimeFreq(:,currentEdgeIdx) =
squeeze(sum(projectedMeasure));
end
end
effectiveConnectivityTimeFreq = single(effectiveConnectivityTimeFreq);
% use single to save time and disk space
effectiveConnectivityTimeFreq =
permute(reshape(effectiveConnectivityTimeFreq, timeFreqSize(1),
timeFreqSize(2), max(fromRegionNumber), max(fromRegionNumber)), [3 4 1 2]);
% Use most consistent dimensions for SIFT and most intuitive for users
save(sprintf([workingFolder filesep userInputFilename
'_%04.f_connectivity'], sessionIdx), 'effectiveConnectivityTimeFreq',
'connectivityType', 'latencies', 'frequencies', 'dimensionLabels',
'fileNameList');
timeLapse = toc;
disp(sprintf('%2.0d/%2.0d subjects done (%0.1d sec lapsed)',
sessionIdx, length(unique(dipolePairAndMeasureObj.sessionId)),
round(timeLapse)));
end
disp('Convert To Anatomical ROI done.')
On Thu, Apr 21, 2016 at 2:48 AM, Mohammed Jarjees <
m.jarjees.1 at research.gla.ac.uk> wrote:
> Dear Makoto,
>
>
>
> Thank you very much for your quick response.
>
>
>
> Sure, I am going to report this problem in Bugzilla.
>
>
>
> Is there any alternative method to extract the ERP values after MPT
> clustering or to extract the P values after run statistics by MPT GUI?
>
>
>
> About the anatomical regions, do you know how can I do that?
>
>
>
> Many Thanks
>
> Mohammed
>
>
>
> *From:* Makoto Miyakoshi [mailto:mmiyakoshi at ucsd.edu]
> *Sent:* 21 April 2016 03:38
> *To:* Mohammed Jarjees
> *Cc:* eeglablist at sccn.ucsd.edu
> *Subject:* Re: [Eeglablist] Extracting ERP values after MPT clustering
>
>
>
> Dear Mohammed,
>
>
>
> Sorry to hear the inconvenience. This is a good detailed question. I
> believe this should be posted for eeglab Bugzilla to be officially
> registered. Would you mind filing the report here? I appreciate your
> cooperation and patience.
>
> https://sccn.ucsd.edu/bugzilla/enter_bug.cgi
>
>
>
> That being said, I recently noticed that Measure Projection could have a
> potential problem of double dipping (because you cluster data using ERP,
> for example, and run statistics on that later). It's probably safer to use
> anatomical regions, not domains (but that may not be supported by GUI).
>
>
>
> Makoto
>
>
>
> On Mon, Apr 18, 2016 at 10:49 AM, Mohammed Jarjees <
> m.jarjees.1 at research.gla.ac.uk> wrote:
>
> Dear All,
>
>
>
> I am using Measure projection toolbox (MPT) to cluster my EEG (ERP) data.
> In this study I have 3 groups of subjects ( First group has 10 subjects,
> second group has 9 subjects and third group has 10 subjects) and I have 3
> conditions per group. I got 3 domains after MPT clustering. I extracted the
> ERP values from “sessionConditionCell” for each domains. The size of
> “sessionConditionCell” is 29*3 (subjects*conditions). I have used the
> following codes to extract the ERP values and to separate these values for
> each group and each condition.
>
>
>
> domainNumber = 1; % this should be 1,2 or 3%
>
> dipoleAndMeasure = STUDY.measureProjection.erp.object;
>
> domain =
> STUDY.measureProjection.erp.projection.domain(domainNumber);
>
> projection = STUDY.measureProjection.erp.projection;
>
> headGrid = STUDY.measureProjection.erp.headGrid;
>
> [linearProjectedMeasure sessionConditionCell groupId uniqeDatasetId
> dipoleDensity] =
> dipoleAndMeasure.getMeanProjectedMeasureForEachSession(headGrid,
> domain.membershipCube, projection.projectionParameter);
>
>
>
> for i =1 : 29 ;
>
> Condition_1 ( i , : )= sessionConditionCell { i, 1};
>
> Condition_2 ( i , : )= sessionConditionCell { i, 2};
>
> Condition_3 ( i , : )= sessionConditionCell { i, 3};
>
> end
>
>
>
> Group_1_Condition_1= Condition_1 ( 1:10, : ); Group_2_Condition_1=
> Condition_1 ( 11:19, : ); Group_3_Condition_1= Condition_1 ( 20:29, : );
>
> Group_1_Condition_2= Condition_2 ( 1:10, : ); Group_2_Condition_2=
> Condition_2 ( 11:19, : ); Group_3_Condition_2= Condition_2 ( 20:29, : );
>
> Group_1_Condition_3= Condition_3 ( 1:10, : ); Group_2_Condition_3=
> Condition_3 ( 11:19, : ); Group_3_Condition_3= Condition_3 ( 20:29, : );
>
>
>
> When I plotted the average ERP for each condition, each group and each
> domain, I got different in* some* results from the results that I got
> from MPT GUI. For example the results for group 1 at Domain 1 for 3
> conditions are same as the results from MPT GUI as well as the result for
> domain 2 and all three group and three condition are same as the MPT GUI.
> However, there are big different for group 2 & 3 at domain 1 as well as
> group 2 & 3 at domain 3 for all three conditions.
>
>
>
> Any Suggestion to solve this problem
>
>
>
> Best Regards
>
> Mohammed
>
>
> _______________________________________________
> Eeglablist page: http://sccn.ucsd.edu/eeglab/eeglabmail.html
> To unsubscribe, send an empty email to
> eeglablist-unsubscribe at sccn.ucsd.edu
> For digest mode, send an email with the subject "set digest mime" to
> eeglablist-request at sccn.ucsd.edu
>
>
>
>
>
> --
>
> Makoto Miyakoshi
> Swartz Center for Computational Neuroscience
> Institute for Neural Computation, University of California San Diego
>
--
Makoto Miyakoshi
Swartz Center for Computational Neuroscience
Institute for Neural Computation, University of California San Diego
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://sccn.ucsd.edu/pipermail/eeglablist/attachments/20160421/5c7d2a5d/attachment.html>
More information about the eeglablist
mailing list