Calculate PAC using pactools on openneuro datasets

BIDS_PAC

""""""
Calculate PAC using pactools on openneuro datasets

""""""


"""
import os
import os.path as op
import openneuro
import mne_bids
from mne_bids import BIDSPath, read_raw_bids, print_dir_tree, make_report


# Open EEG dataset from openneuro
# `Parkinson's EEG dataset <https://openneuro.org/datasets/ds002778`_
dataset = 'ds002778'
subject = 'pd14'
# Download one subject's data from each dataset
target_dir = op.join('examples', dataset)
if not op.isdir(target_dir):
    os.makedirs(target_dir)
openneuro.download(dataset=dataset, target_dir=target_dir,
                   include=[f'sub-{subject}'])
bids_root = op.join('examples', dataset)
print_dir_tree(bids_root, max_depth=4)
print(make_report(bids_root))
datatype = 'eeg'
session = 'off'
bids_path = BIDSPath(root=bids_root, session=session, datatype=datatype)
print(bids_path.match())
task = 'rest'
suffix = 'eeg'
bids_path = BIDSPath(subject=subject, session=session, task=task,
                     suffix=suffix, datatype=datatype, root=bids_root)
print(bids_path)
bids_path
#display the details of the dataset
raw = read_raw_bids(bids_path=bids_path, verbose=False)
print(raw.info['subject_info'])
print(raw.info['line_freq'])
print(raw.info['sfreq'])
raw.plot()


#preprocessing using MNE tools to remove bad channels, average referencing and removing drifts

raw_drop=raw.drop_channels(['EXG1', 'EXG2', 'EXG3', 'EXG4','EXG5', 'EXG6', 'EXG7', 'EXG8',
                         'Status'])
raw_drop.plot()
raw_drop.set_eeg_reference(ref_channels='average')
raw_drop.filter(0.5, None, fir_design='firwin',phase='zero-double')


#picking c3 and c4 channel from the list of EEG channels to perform PAC

c3c4 = raw_drop.copy().pick_channels(['C3','C4'])
data, times = c3c4[['C3'], :]
data1,times1=c3c4[['C4'], :]


#calculate PAC for phase frequency of 13-30hz,bandwidth of 2hz and amplitude frequency of 50-100hz, bandwidth of 4hz.

import numpy as np
from pactools import Comodulogram, REFERENCES
import matplotlib.pyplot as plt

low_fq_width = 2
high_fq_width= 4
low_fq_range = np.linspace(13,30,100)
high_fq_range= np.linspace(50,100,100)
fig, ax = plt.subplots(1, 1)
fs=512

# Compute the comodulograms for c3, c4 and average c3 and c4 channel and plot them
#for ax, method in zip(axs, methods):
estimator_c3= Comodulogram (fs=fs, low_fq_range=low_fq_range,
                          low_fq_width=low_fq_width,
                            high_fq_range=high_fq_range,
                            high_fq_width=high_fq_width,
                            method='tort',progress_bar=True)
estimator_c4= Comodulogram (fs=fs, low_fq_range=low_fq_range,
                          low_fq_width=low_fq_width,
                            high_fq_range=high_fq_range,
                            high_fq_width=high_fq_width,
                            method='tort',progress_bar=True)
estimator_avg= Comodulogram (fs=fs, low_fq_range=low_fq_range,
                          low_fq_width=low_fq_width,
                            high_fq_range=high_fq_range,
                            high_fq_width=high_fq_width,
                            method='tort',progress_bar=True)
estimator_c3.fit(data[0])
estimator_c4.fit(data1[0])
estimator_avg.comod_=np.mean([estimator_c3.comod_,estimator_c4.comod_],axis=0)
fig,axs=plt.subplots(1,3,figsize=(15,5))
estimator_c3.plot(titles=['C3'], axs=[axs[0]])
estimator_c4.plot(titles=['C4'], axs=[axs[1]])
estimator_avg.plot(titles=['avg'], axs=[axs[2]])
 plt.show()
 print(f' c3 mean {np.mean([estimator_c3.comod_])}')
 print(f' c4 mean {np.mean([estimator_c4.comod_])}')
 print(f'c3_c4_avg mean {np.mean([estimator_avg.comod_])}')
 """