How to get coherence between stc files?

✉️ Mailing List Archive (read-only)
1

Dear mne experts,

In an example ?Compute coherence in source space using a MNE inverse
solution? mne.connectivity.spectral_connectivity is used to compute the
coherence between a seed in the left auditory cortex and the rest of the
brain based on single-trial MNE-dSPM inverse solutions. However, I would
like to use mne.connectivity.spectral_connectivity to compute the coherence
between two stc files measured in two different conditions.

I have tried to do this as follows:
stcs=np.concatenate((stc_1.data,stc_2.data),axis=0)
stcs=stcs.tolist()
indices=(np.arange(1,20484),np.arange(20485,40968)) (I would like to get
the coherence between each vertex in stc1 and the corresponding vertices in
stc2)
coh, freqs, times, n_epochs, n_tapers = spectral_connectivity(
    stcs, method='coh', mode='fourier', indices=indices,
    sfreq=sfreq, fmin=fmin, fmax=fmax, faverage=True, n_jobs=1)

However, I get:

Connectivity computation...

2

After this line:

stcs=np.concatenate((stc_1.data,stc_2.data),axis=0)

?
You will have stcs.shape == (40968, n_times).

However, for `spectral_connectivity`, the `data` argument needs to be of
shape (n_epochs, n_signals, n_times), i.e. you need multiple epochs/trials.
Here you appear to only have 1.

HTH,
Eric

3

Hi Eric,

Many thanks for your answer!

I managed to combine the data as follows:

data=np.concatenate((epochs_1.get_data(),epochs_2.get_data()),axis=0)

data.shape

Out[28]: (320, 308, 3721)

I am not sure how the parameter ?indices? should be set in my case. I would
like to calculate the coherence between each data point in the file 1 and
the corresponding data points in the file 2. Since I have 320 epochs after
concatenation, I tried:

indices=(np.arange(1,160),np.arange(161,320)) (i.e. 160 epochs in each
file).

However, this doesn?t seem to work:

?Connectivity computation...

    computing connectivity for 159 connections

    using t=0.000s..3.100s for estimation (3721 points)

    computing connectivity for the bands:

     band 1: 8.1Hz..12.9Hz (16 points)

     band 2: 13.2Hz..30.0Hz (53 points)

    connectivity scores will be averaged for each band

    using FFT with a Hanning window to estimate spectra

    the following metrics will be computed: Coherence

    computing connectivity for epoch 1

4

Hi Maria,

regarding the error, I couldn't replicate it myself and I don't see any
problems with your code (except that indices does not include 0 and 160).
My hunch is however, that the problem is that the allocation fails to
recognize the relevant channels. You could try dropping the irrelevant
channels from the epochs before concatenating the data to see if it works
(probably something like epochs.pick_types(meg=True)).

-Jaakko

5

data=np.concatenate((epochs_1.get_data(),epochs_2.get_data()),axis=0)

data.shape

Out[28]: (320, 308, 3721)

Following the nomenclature from the connectivity docs
<https://mne-tools.github.io/stable/generated/mne.connectivity.spectral_connectivity.html#mne.connectivity.spectral_connectivity&gt;,
your input data need to be formatted as (n_epochs, n_signals, n_times).
Keep in mind that the n_signals dimension is usually the spatial one, like
source space vertices, or in your case, channels. So the way you have this
set up, you have 308 spatial channels / signals, each with 320 epochs /
repeats.

I am not sure how the parameter ?indices? should be set in my case.

Indices should be into the (typically) spatial dimension, i.e. the
n_signals dimension.

I would like to calculate the coherence between each data point in the file

1 and the corresponding data points in the file 2.

"data point" is a bit ambiguous -- I assume you mean you want to estimate
the connectivity between *each channel* from the first dataset and *each
channel* in the second dataset....

Since I have 320 epochs after concatenation, I tried:

indices=(np.arange(1,160),np.arange(161,320)) (i.e. 160 epochs in each
file).

... which means this probably isn't set up to do what you want.

What you might want is this:

data=np.concatenate((epochs_1.get_data(),epochs_2.get_data()),axis=1)
data.shape
(160, 616, 3721)

?
This way, you have 616 signals of interest (308 signals / channels from one
condition, 308 from signals / channels from another), each with 160 epochs
/ repeats.

Hopefully from there you can come up with how to get connectivity you want.
For example, if you want each channel only with its corresponding channel
from the other condition (308 estimates), this would be something like:

(np.arange(308), np.arange(308) + 308)

?
If you wanted all channels from condition 1 with all channels from
condition 2 (94864 estimates), this would be something like:

(np.repeat(np.arange(308), 308), np.tile(np.arange(308) + 308, 308))

?
And 308 is a somewhat rare number of data channels -- you might want to
check to make sure you have only data channels picked as Jaakko suggests.

HTH,
Eric
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.nmr.mgh.harvard.edu/pipermail/mne_analysis/attachments/20170616/e93812c2/attachment.html

6

Hi,

Many thanks for your answers! This clarified the matter for me.

However, instead of calculating coherence in sensor space, I would like to
calculate coherence in source space (i.e. I would need to get coherence
between each vertex from the first dataset and their counterparts from the
second dataset). Therefore, I have applied inverse operator for each epoch
of condition 1 and condition 2:

stcs1 = apply_inverse_epochs(epochs1, inverse_operator1, lambda2, method,
                            pick_ori="normal")
stcs2 = apply_inverse_epochs(epochs2, inverse_operator2, lambda2, method,
                            pick_ori="normal")

This gives me a list of 160 stcs that are as follows:
stcs1[0]
Out[25]: <SourceEstimate | 20482 vertices, subject : ah, tmin :
-99.9936018905 (ms), tmax : 2899.81445482 (ms), tstep : 0.833280015754
(ms), data size : 20482 x 3601>

Then, I could combine these two list of stcs and give the resulting list as
an argument to mne.connectivity.spectral_connectivity:
coh, freqs, times, n_epochs, n_tapers = spectral_connectivity(
    stcs, method='coh', mode='fourier', indices=indices,
    sfreq=sfreq, fmin=fmin, fmax=fmax, faverage=True, n_jobs=1)

However, I would need to estimate coherence between each vertex from the
first dataset and each vertex in the second dataset. As far as I
understand, ?indices? can only determine the vertices between which to
compute coherence within a dataset but not between two datasets.

The documentation says:

"By default, the connectivity between all signals is computed (only
connections corresponding to the lower-triangular part of the connectivity
matrix). If one is only interested in the connectivity between some
signals, the ?indices? parameter can be used. For example, to compute the
connectivity between the signal with index 0 and signals ?2, 3, 4? (a total
of 3 connections) one can use the following:

indices = (np.array([0, 0, 0]), # row indices
           np.array([2, 3, 4])) # col indices

con_flat = spectral_connectivity(data, method='coh',
                                 indices=indices, ...)"

Could you please let me know if there is any way in mne to calculate
coherence between two datasets in source space?

Best,
Maria

2017-06-16 18:32 GMT+03:00 Eric Larson <larson.eric.d at gmail.com>:

data=np.concatenate((epochs_1.get_data(),epochs_2.get_data()),axis=0)

data.shape

Out[28]: (320, 308, 3721)

Following the nomenclature from the connectivity docs
<https://mne-tools.github.io/stable/generated/mne.connectivity.spectral_connectivity.html#mne.connectivity.spectral_connectivity&gt;,
your input data need to be formatted as (n_epochs, n_signals, n_times).
Keep in mind that the n_signals dimension is usually the spatial one, like
source space vertices, or in your case, channels. So the way you have this
set up, you have 308 spatial channels / signals, each with 320 epochs /
repeats.

I am not sure how the parameter ?indices? should be set in my case.

Indices should be into the (typically) spatial dimension, i.e. the
n_signals dimension.

I would like to calculate the coherence between each data point in the

file 1 and the corresponding data points in the file 2.

"data point" is a bit ambiguous -- I assume you mean you want to estimate
the connectivity between *each channel* from the first dataset and *each
channel* in the second dataset....

Since I have 320 epochs after concatenation, I tried:

indices=(np.arange(1,160),np.arange(161,320)) (i.e. 160 epochs in each
file).

... which means this probably isn't set up to do what you want.

What you might want is this:

data=np.concatenate((epochs_1.get_data(),epochs_2.get_data()),axis=1)
data.shape
(160, 616, 3721)

?
This way, you have 616 signals of interest (308 signals / channels from
one condition, 308 from signals / channels from another), each with 160
epochs / repeats.

Hopefully from there you can come up with how to get connectivity you
want. For example, if you want each channel only with its corresponding
channel from the other condition (308 estimates), this would be something
like:

(np.arange(308), np.arange(308) + 308)

?
If you wanted all channels from condition 1 with all channels from
condition 2 (94864 estimates), this would be something like:

(np.repeat(np.arange(308), 308), np.tile(np.arange(308) + 308, 308))

?
And 308 is a somewhat rare number of data channels -- you might want to
check to make sure you have only data channels picked as Jaakko suggests.

HTH,
Eric

_______________________________________________
Mne_analysis mailing list
Mne_analysis at nmr.mgh.harvard.edu
Mne_analysis Info Page

The information in this e-mail is intended only for the person to whom it
is
addressed. If you believe this e-mail was sent to you in error and the
e-mail
contains patient information, please contact the Partners Compliance
HelpLine at
MyComplianceReport.com: Compliance and Ethics Reporting . If the e-mail was sent to you in
error
but does not contain patient information, please contact the sender and
properly
dispose of the e-mail.

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.nmr.mgh.harvard.edu/pipermail/mne_analysis/attachments/20170617/b89efa94/attachment.html

7

Concatenate stcs1 and stcs2 along axis=1 (spatial dimension), similar to
before. The logic for vertices here is the same as it was for channels in
your previous example. The connectivity code has no notion of datasets,
only single signals (axis=1) -- so whatever you want to compute
connectivity between (channels, sources, label means, whatever) needs to be
defined along that axis (with the epochs along axis=0 being meaningfully
paired across signals).

Eric

Hi,

Many thanks for your answers! This clarified the matter for me.

However, instead of calculating coherence in sensor space, I would like to
calculate coherence in source space (i.e. I would need to get coherence
between each vertex from the first dataset and their counterparts from the
second dataset). Therefore, I have applied inverse operator for each epoch
of condition 1 and condition 2:

stcs1 = apply_inverse_epochs(epochs1, inverse_operator1, lambda2, method,
                            pick_ori="normal")
stcs2 = apply_inverse_epochs(epochs2, inverse_operator2, lambda2, method,
                            pick_ori="normal")

This gives me a list of 160 stcs that are as follows:
stcs1[0]
Out[25]: <SourceEstimate | 20482 vertices, subject : ah, tmin :
-99.9936018905 (ms), tmax : 2899.81445482 (ms), tstep : 0.833280015754
(ms), data size : 20482 x 3601>

Then, I could combine these two list of stcs and give the resulting list as
an argument to mne.connectivity.spectral_connectivity:
coh, freqs, times, n_epochs, n_tapers = spectral_connectivity(
    stcs, method='coh', mode='fourier', indices=indices,
    sfreq=sfreq, fmin=fmin, fmax=fmax, faverage=True, n_jobs=1)

However, I would need to estimate coherence between each vertex from the
first dataset and each vertex in the second dataset. As far as I
understand, ?indices? can only determine the vertices between which to
compute coherence within a dataset but not between two datasets.

The documentation says:

"By default, the connectivity between all signals is computed (only
connections corresponding to the lower-triangular part of the connectivity
matrix). If one is only interested in the connectivity between some
signals, the ?indices? parameter can be used. For example, to compute the
connectivity between the signal with index 0 and signals ?2, 3, 4? (a total
of 3 connections) one can use the following:

indices = (np.array([0, 0, 0]), # row indices
           np.array([2, 3, 4])) # col indices

con_flat = spectral_connectivity(data, method='coh',
                                 indices=indices, ...)"

Could you please let me know if there is any way in mne to calculate
coherence between two datasets in source space?

Best,
Maria

2017-06-16 18:32 GMT+03:00 Eric Larson <larson.eric.d at gmail.com>:

8

Hi,

Unfortunately, I think that I can?t do that so easily since
apply_inverse_epochs gives a list of stcs and
len(stcs) is 160, so it has only one dimension. If I try
data=np.concatenate((stcs_1,stcs_2),axis=1)

I get:
IndexError: axis 1 out of bounds [0, 1)

Only this works:
data=np.concatenate((stcs_1,stcs_2),axis=0)

In [45]: data.shape
Out[45]: (320,)

data[0]
<SourceEstimate | 20482 vertices, subject : ah, tmin : -99.9936018905
(ms), tmax : 2899.81445482 (ms), tstep : 0.833280015754 (ms), data size :
20482 x 3601>

However, I can concatenate data of stc files from two lists of stcs:
data=np.concatenate((stcs_1[0].data,stcs_2[0].data),axis=0)
data.shape
Out[39]: (40964, 3601) (vertices=40964 and time=3601)

Thereafter, I tried to create a new stc as follows:
stc_3=stc_1
stc_3.data = data

but got:
AttributeError: can't set attribute

I wonder if there is any way to replace ?data? attribute of stc file? I
haven?t managed to find any example about this.

If I could create new stc files of the concatenated data, I could probably
create a list of them and give it as an argument to
mne.connectivity.spectral_connectivity. In this case, the indices would be
indices = (np.array([1, 2, ? ,20482]), np.array([20483, ? 40964]))

Best,
Maria

2017-06-18 0:20 GMT+03:00 Eric Larson <larson.eric.d at gmail.com>:

Concatenate stcs1 and stcs2 along axis=1 (spatial dimension), similar to
before. The logic for vertices here is the same as it was for channels in
your previous example. The connectivity code has no notion of datasets,
only single signals (axis=1) -- so whatever you want to compute
connectivity between (channels, sources, label means, whatever) needs to be
defined along that axis (with the epochs along axis=0 being meaningfully
paired across signals).

Eric

Hi,

Many thanks for your answers! This clarified the matter for me.

However, instead of calculating coherence in sensor space, I would like to
calculate coherence in source space (i.e. I would need to get coherence
between each vertex from the first dataset and their counterparts from the
second dataset). Therefore, I have applied inverse operator for each epoch
of condition 1 and condition 2:

stcs1 = apply_inverse_epochs(epochs1, inverse_operator1, lambda2, method,
                            pick_ori="normal")
stcs2 = apply_inverse_epochs(epochs2, inverse_operator2, lambda2, method,
                            pick_ori="normal")

This gives me a list of 160 stcs that are as follows:
stcs1[0]
Out[25]: <SourceEstimate | 20482 vertices, subject : ah, tmin :
-99.9936018905 (ms), tmax : 2899.81445482 (ms), tstep : 0.833280015754
(ms), data size : 20482 x 3601>

Then, I could combine these two list of stcs and give the resulting list
as an argument to mne.connectivity.spectral_connectivity:
coh, freqs, times, n_epochs, n_tapers = spectral_connectivity(
    stcs, method='coh', mode='fourier', indices=indices,
    sfreq=sfreq, fmin=fmin, fmax=fmax, faverage=True, n_jobs=1)

However, I would need to estimate coherence between each vertex from the
first dataset and each vertex in the second dataset. As far as I
understand, ?indices? can only determine the vertices between which to
compute coherence within a dataset but not between two datasets.

The documentation says:

"By default, the connectivity between all signals is computed (only
connections corresponding to the lower-triangular part of the connectivity
matrix). If one is only interested in the connectivity between some
signals, the ?indices? parameter can be used. For example, to compute the
connectivity between the signal with index 0 and signals ?2, 3, 4? (a total
of 3 connections) one can use the following:

indices = (np.array([0, 0, 0]), # row indices
           np.array([2, 3, 4])) # col indices

con_flat = spectral_connectivity(data, method='coh',
                                 indices=indices, ...)"

Could you please let me know if there is any way in mne to calculate
coherence between two datasets in source space?

Best,
Maria

2017-06-16 18:32 GMT+03:00 Eric Larson <larson.eric.d at gmail.com>:

data=np.concatenate((epochs_1.get_data(),epochs_2.get_data()),axis=0)

data.shape

Out[28]: (320, 308, 3721)

Following the nomenclature from the connectivity docs
<https://mne-tools.github.io/stable/generated/mne.connectivity.spectral_connectivity.html#mne.connectivity.spectral_connectivity&gt;,
your input data need to be formatted as (n_epochs, n_signals, n_times).
Keep in mind that the n_signals dimension is usually the spatial one, like
source space vertices, or in your case, channels. So the way you have this
set up, you have 308 spatial channels / signals, each with 320 epochs /
repeats.

I am not sure how the parameter ?indices? should be set in my case.

Indices should be into the (typically) spatial dimension, i.e. the
n_signals dimension.

I would like to calculate the coherence between each data point in the

file 1 and the corresponding data points in the file 2.

"data point" is a bit ambiguous -- I assume you mean you want to estimate
the connectivity between *each channel* from the first dataset and *each
channel* in the second dataset....

Since I have 320 epochs after concatenation, I tried:

indices=(np.arange(1,160),np.arange(161,320)) (i.e. 160 epochs in each
file).

... which means this probably isn't set up to do what you want.

What you might want is this:

data=np.concatenate((epochs_1.get_data(),epochs_2.get_data()),axis=1)
data.shape
(160, 616, 3721)

?
This way, you have 616 signals of interest (308 signals / channels from
one condition, 308 from signals / channels from another), each with 160
epochs / repeats.

Hopefully from there you can come up with how to get connectivity you
want. For example, if you want each channel only with its corresponding
channel from the other condition (308 estimates), this would be something
like:

(np.arange(308), np.arange(308) + 308)

?
If you wanted all channels from condition 1 with all channels from
condition 2 (94864 estimates), this would be something like:

(np.repeat(np.arange(308), 308), np.tile(np.arange(308) + 308, 308))

?
And 308 is a somewhat rare number of data channels -- you might want to
check to make sure you have only data channels picked as Jaakko suggests.

HTH,
Eric

_______________________________________________
Mne_analysis mailing list
Mne_analysis at nmr.mgh.harvard.edu
Mne_analysis Info Page

The information in this e-mail is intended only for the person to whom it
is
addressed. If you believe this e-mail was sent to you in error and the
e-mail
contains patient information, please contact the Partners Compliance
HelpLine at
MyComplianceReport.com: Compliance and Ethics Reporting . If the e-mail was sent to you
in error
but does not contain patient information, please contact the sender and
properly
dispose of the e-mail.

_______________________________________________
Mne_analysis mailing list
Mne_analysis at nmr.mgh.harvard.edu
Mne_analysis Info Page

The information in this e-mail is intended only for the person to whom it
is
addressed. If you believe this e-mail was sent to you in error and the
e-mail
contains patient information, please contact the Partners Compliance
HelpLine at
MyComplianceReport.com: Compliance and Ethics Reporting . If the e-mail was sent to you in
error
but does not contain patient information, please contact the sender and
properly
dispose of the e-mail.

_______________________________________________
Mne_analysis mailing list
Mne_analysis at nmr.mgh.harvard.edu
Mne_analysis Info Page

The information in this e-mail is intended only for the person to whom it
is
addressed. If you believe this e-mail was sent to you in error and the
e-mail
contains patient information, please contact the Partners Compliance
HelpLine at
MyComplianceReport.com: Compliance and Ethics Reporting . If the e-mail was sent to you in
error
but does not contain patient information, please contact the sender and
properly
dispose of the e-mail.

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://mail.nmr.mgh.harvard.edu/pipermail/mne_analysis/attachments/20170618/0a862153/attachment-0001.html

9

Hi,

from having discussed with Maria, I gathered she doesn?t want to compute connectivity at all. She just wants to compute the coherence (across time, so one value for each vertex) between two experimental conditions.

Going by Erics recommendations, it may be that the spectral_connectivity function may be able to be coaxed into doing that:

# Apply the inverse operator on each epoch
stcs_condition1 = apply_inverse_epochs(epochs_condition1, inverse_operator, ?)
stcs_condition2 = apply_inverse_epochs(epochs_condition1, inverse_operator, ?)

# Concatenate the epochs into a big 3D array
data_condition1 = np.data([s.data for s in secs_condition1])
data_condition2 = np.data([s.data for s in secs_condition2])

# The number of vertices is the size of one of the data arrays along the 2nd dimension
n_vertices = data_condition1.shape[1]

# Concatenate the two conditions along the spatial axis (what Eric suggested)
data = np.concatenate([data_condition1, data_condition2], axis=1)

# We are going to tell the spectral_connectivity function to compute coherence for each vertex between both conditions.
# So vertex 0 should be compared to vertex 0 + n_vertices, and vertex 1 to vertex 1 + n_vertices, etc.
indices = [np.arange(n_vertices), np.arange(n_vertices) + n_vertices]

# Finally, call the spectral_connectivity function
con, freqs, times, n_epochs, n_tapers = spectral_connectivity(data, indices=indices, ...)

I?ve no idea whether that will work, but it is my best guess :slight_smile: Maybe someone with more experience in computing coherence can give better guidance.

Marijn.