업데이트 'README.md'

README.md

@@ -2,7 +2,7 @@
 
 
 
-version 1.0.2
+version 2.0.1 (2020-02-11. 15:01:33)
 
 ##### https://gin.g-node.org/hiobeen/Mouse_EEG_ChronicSleepRestriction_Kim_et_al
 
@@ -27,7 +27,7 @@ Kim, B., Kocsis, B., Hwang, E., Kim, Y., Strecker, R. E., McCarley, R. W., & Cho
 
 Kim, B., Hwang, E., Strecker, R. E., Choi, J. H., & Kim, Y. (2020). Differential modulation of NREM sleep regulation and EEG topography by chronic sleep restriction in mice. Scientific reports, 10(1), 1-13 [https://www.nature.com/articles/s41598-019-54790-y].
 
-See also: Hwang, E., Han, H. B., Kim, J. Y., & Choi, J. H. (2020). High-density EEG of auditory steady-state responses during stimulation of basal forebrain parvalbumin neurons. Scientific Data, 7(1), 1-9 [https://www.nature.com/articles/s41598-019-54790-y].
+See also: Hwang, E., Han, H. B., Kim, J. Y., & Choi, J. H. (2020). High-density EEG of auditory steady-state responses during stimulation of basal forebrain parvalbumin neurons. Scientific Data, 7(1), 1-9 [https://www.nature.com/articles/s41597-020-00621-z].
 * Dataset URL: https://gin.g-node.org/hiobeen/Mouse_hdEEG_ASSR_Hwang_et_al
 
 Visit our website [https://www.jeelab.net/] for full puiblication list.
@@ -86,9 +86,14 @@ from scipy.io import loadmat, savemat
 from scipy.signal import butter, sosfiltfilt, find_peaks
 from scipy.stats import pearsonr
 from scipy.interpolate import interp1d
+#!pip install scipy==1.1.0
 from scipy.misc import imresize as resize
 def fft_half(x, fs=500): return np.fft.fft(x.reshape(-1))[:int(len(x)/2)]/len(x), np.linspace(0,fs/2,int(len(x)/2))
 def find(x): return np.argwhere( x.flatten() ).flatten()
+def find_idx(target_range, continuous_vector):
+    start_point = np.where(continuous_vector<target_range[0])[0][-1]+1
+    end_point = np.where(continuous_vector<target_range[-1])[0][-1]
+    return (start_point, end_point)
 def smooth(x, span): return np.convolve(x,np.ones(int(span))/int(span), mode='same')
 def interp_nan(x):
     placeholder = np.arange(x.size)
@@ -99,8 +104,11 @@ def bandpower(x, band, fs=500):
 def find_seg(ind): # Handy function
   ind_diff = np.diff(np.concatenate((np.ones(1)*-1,ind),axis=0))
   discontinuous = np.where(ind_diff!=1)[0]
-  data_seg = np.zeros((len(discontinuous)-1,2),dtype='int')
-  for idx in range(len(discontinuous)-1): data_seg[idx,:]=[ind[discontinuous[idx]] ,ind[ discontinuous[idx+1]-1]]
+  if len(discontinuous)>1:
+    data_seg = np.zeros((len(discontinuous)-1,2),dtype='int')
+    for idx in range(len(discontinuous)-1): data_seg[idx,:]=[ind[discontinuous[idx]],ind[discontinuous[idx+1]-1]]
+  else:
+    data_seg = np.array([])
   return data_seg
 def bandpass_filt(x, band=(8,18), fs=500, order=10):
   #if x.shape[0]>x.shape[1]: x=x.transpose()
@@ -118,12 +126,12 @@ After downloading dataset, you can import EEG dataset using the function load_da
 
 ```python
 # File loading demo
-
 """ Important notice: Change these directories to fit your local environment """
-bids_path = YOUR_OWN_PATH_1 # Ex: 'D:/eeg_files/data_BIDS/'
-save_directory = YOUR_OWN_PATH_2 # Ex: 'D:/eeg_files/output/'
+bids_path = 'D:/SciData_Bowon_Sleep/data_BIDS_0812/'
+save_directory = 'D:/output_ipynb/'; 
 if not os.path.isdir( save_directory ): os.mkdir( save_directory )
-figure_directory =  '%s/figures/'%save_directory; if not os.path.isdir( figure_directory ): os.mkdir( figure_directory )
+figure_directory =  '%s/figures/'%save_directory; 
+if not os.path.isdir( figure_directory ): os.mkdir( figure_directory )
 figure_dpi = 300
 
 step_names = ['BL', 'SR1', 'SR2', 'SR3', 'SR4', 'SR5', 'R1', 'R2', 'R3']
@@ -134,7 +142,18 @@ def load_data(sbjIdx=0, sessIdx=0, verbose = True):
     
     # Load screw EEG data
     print( '** Data loading ... Subject %d, Session %d (Step: %s) **'%(sbjIdx+1, sessIdx+1,step_names[sessIdx]))
-    EEG_screw = load( fname_screw+'.set', preload = True, verbose = False)
+    print( fname_screw )
+    EEG_screw = load( fname_screw+'.set', preload = False, verbose = False)
+    # Author's note by Hio-Been Han
+    #
+    # If you have a problem with the line above (in case of mne==0.24.1)
+    # Please modify 'line 355-356' at (eeglab.py), by eliminating following two lines
+    ## (Line 355 at eeglab.py) # self.set_annotations(annot) 
+    ## (Line 356 at eeglab.py) # _check_boundary(annot, None)
+    #
+    EEG_screw.times = np.linspace(EEG_screw.first_time,
+                            EEG_screw.n_times/EEG_screw.info['sfreq'], 
+                            EEG_screw.n_times)
     if ~hasattr(EEG_screw, 'data'): EEG_screw.data = EEG_screw.get_data()
     EEG_screw.data = EEG_screw.data/CAL[0] # Voltage value calibration
     if verbose: print( 'Screw data shape: [%d channels x %d times]'%EEG_screw.data.shape )
@@ -146,7 +165,10 @@ def load_data(sbjIdx=0, sessIdx=0, verbose = True):
     
     # Load hdEEG data if exist
     if os.path.exists(fname_hdeeg+'.set'):
-        EEG_hdeeg = load( fname_hdeeg+'.set', preload=True, verbose = False)
+        EEG_hdeeg = load( fname_hdeeg+'.set', preload=False, verbose = False)
+        EEG_hdeeg.times = np.linspace(EEG_hdeeg.first_time,
+                                    EEG_hdeeg.n_times/EEG_hdeeg.info['sfreq'], 
+                                    EEG_hdeeg.n_times)
         if ~hasattr(EEG_hdeeg, 'data'): EEG_hdeeg.data = EEG_hdeeg.get_data()
         EEG_hdeeg.data = EEG_hdeeg.data/CAL[1] # Voltage value calibration
         if verbose: print( 'hdEEG data shape: [%d channels x %d times]'%EEG_hdeeg.data.shape )
@@ -167,7 +189,7 @@ def load_data(sbjIdx=0, sessIdx=0, verbose = True):
 ```
 
 Now we're all set for EEG signal processing.
-   
+
 ### 3-2) Sleep wave analysis (SWA, Spindle)
 
 Several analysis can be done on this data such as spectral power, slow wave activity, sleep spindles, traveling wave, interaction among different oscillatory activities. Here, we provide examples of the sleep wave analyses including SWA (slow wave activity) and spindle. For that, example EEG data from the dataset (subject 1, session 1) is loaded. 
@@ -179,6 +201,7 @@ EEG_screw, EEG_hdeeg, sleep_score = load_data(sbjIdx=0, sessIdx=0, verbose = Tru
 ```
 
     ** Data loading ... Subject 1, Session 1 (Step: BL) **
+    D:/SciData_Bowon_Sleep/data_BIDS_0812/sub-01/ses-01/eeg/sub-01-BL_screw
     Screw data shape: [4 channels x 43200000 times]
     Screw bad channels: []
     hdEEG data shape: [36 channels x 3600000 times]
@@ -298,11 +321,11 @@ Similarly, detection and visualization of spindle events can be done by followin
 ```python
 """ Spindle Detection Function """
 from scipy import signal
-def detect_spindle( data, fs=500 ):
-    spindlef = bandpass_filt( data, band=(8,18), fs=fs )
-    thr_updown = ( np.std(spindlef)*5, np.std(spindlef)*2.5 )
-    thr_duration = 0.30 # sec
-    thr_interval = 0.10 # sec
+def detect_spindle( data, fs=500, param_band = (8,18), param_thr=(5, 2.5, .3, .1) ):
+    spindlef = bandpass_filt( data, band=param_band, fs=fs )
+    thr_updown = ( np.std(spindlef)*param_thr[0], np.std(spindlef)*param_thr[1] )
+    thr_duration = param_thr[2] # sec
+    thr_interval = param_thr[3] # sec
     hilb_amptd = np.abs( signal.hilbert( spindlef ))
     amptd_envelop = smooth( hilb_amptd, fs/5 )
     big_enough = np.where(amptd_envelop>thr_updown[1])[0]
@@ -374,6 +397,176 @@ plt.show();
 ![png](figures/output_14_0.png)
 
 
+#### 3. Multi-channel detection of Spindle and SWA activities
+
+The detection functions above were designed for single-channel data only, however, we have 36-channel hdEEG data. To detect sleep-related EEG activities in multi-channel, following two functions detect_swa_multichan() and detect_spindles_multichan() are introduced by slightly modifying the original functions. FYI, threshold parameters for detection are adjusted to cover channel-variant voltage range.
+
+
+
+```python
+"""  Define functions: Detection in multiple channels """
+# (1) SWA
+def detect_swa_multichan(data, fs=500, t=EEG_hdeeg.times, shift_band = (.05, 6),
+               thr_duration = (0.12, 0.70), thr_std = .7, thr_xcorr = 0.7):
+    
+    # Data filtering & downshifting
+    data_filt_all, data_shift_all = data.copy(), data.copy()
+    for chanIdx in range(data.shape[0]):
+        data_filt_all[chanIdx,:] = bandpass_filt( data[chanIdx,:], band=shift_band, fs=fs )
+        data_shift_all[chanIdx,:] = data_filt_all[chanIdx,:].reshape(-1) - get_upper_envelope(data_filt_all[chanIdx,:],t)
+    std_data = np.std( data_shift_all.flatten() )
+    
+    idx_swas = []
+    for chanIdx in range(data.shape[0]):
+        data_filt = data_filt_all[chanIdx,:]
+        data_shift = data_shift_all[chanIdx,:]
+    
+        # Get upper/lower indeces near zero
+        cross_negative = np.insert( find( data_shift<0 ), 0,0)
+        cross_positive = find( data_shift>=0 )
+        ups = cross_negative[ find(np.diff(cross_negative)>1) ]+1
+        downs = cross_positive[ find(np.diff(cross_positive)>1) ]
+        if downs[0] < ups[0]:
+            min_length = min(downs.shape[0], ups.shape[0]);
+            updowns=np.array( [list(downs),list(ups) ])
+        else:
+            min_length = min(downs.shape[0], ups.shape[0]-1);
+            updowns=np.array( [downs[0:min_length-1],ups[1:min_length] ])
+    
+        # Main inspection loop
+        result = []
+        for j in range(min_length-1):
+            data_piece = data_shift[updowns[0][j]:updowns[1][j]]
+            min_where = np.argmin(data_piece)
+            n =len(data_piece)
+      
+            # Rejection criteria
+            flag_long = n < (thr_duration[0]*fs)
+            flag_short = n >= (thr_duration[1]*fs)
+            flag_amptd = np.max(np.abs(data_piece)) < (std_data*thr_std)
+            shift_squared = np.diff(data_filt[updowns[0][j]:updowns[1][j]-1])*np.diff(data_filt[updowns[0][j]+1:updowns[1][j]])
+            n_inflection = find(shift_squared<0).shape[0]
+            flag_inflection = n_inflection>2
+            reject = flag_long | flag_short | flag_amptd | flag_inflection
+
+            # Calculate correlation coefficient
+            if not reject:
+                min_where = np.argmin(data_piece)
+                templateL = np.arange(0,min_where) / (min_where-1)
+                templateR = np.arange(0,n-min_where)[::-1] / (n-min_where) 
+                template = np.concatenate((templateL,templateR),axis=0)*np.min(data_piece)
+                template[find(np.isnan(template))]=0
+                r, _ = pearsonr( template, data_piece )
+                if r > thr_xcorr: 
+                    result.append([updowns[0][j],updowns[1][j]])
+        idx_swas.append(np.array(result, dtype=int))
+    return idx_swas
+
+# (2) Spindles
+def detect_spindle_multichan( data, fs=500, param_band = (8,18), param_thr=(14, 1.6, .5, .1) ):
+    spindlef, amptd_envelop = data.copy(), data.copy()
+    for chanIdx in range(data.shape[0]):
+        spindlef[chanIdx,:] = bandpass_filt( data[chanIdx,:], band=param_band, fs=fs )
+        amptd_envelop[chanIdx,:] = smooth(np.abs(signal.hilbert(spindlef[chanIdx,:])), fs/5 )
+        
+    std_all = np.std( spindlef.flatten() )
+    thr_updown = ( std_all*param_thr[0], std_all*param_thr[1] )
+    thr_duration = param_thr[2] # sec
+    thr_interval = param_thr[3] # sec
+    big_enough_all = amptd_envelop>thr_updown[1]
+    idx_spindles = []
+    for chanIdx in range(data.shape[0]):
+        idx_merged = np.array([])
+        big_enough = np.where(big_enough_all[chanIdx,:])[0]
+        if len(big_enough)>0:
+            # Duration Thresholding
+            start_end = find_seg( big_enough )
+            if start_end.shape[0]:
+                long_enough = start_end[np.where( (start_end[:,1]-start_end[:,0])>(thr_duration*fs))[0],:]
+                # Amplitude Thresholding
+                include_flag = []
+                for i in range(long_enough.shape[0]):
+                  data_to_inspect = data[chanIdx,range(long_enough[i,0],long_enough[i,1])]
+                  if len( np.where(data_to_inspect>thr_updown[0])[0])==0:
+                    include_flag.append(i)
+                strong_enough = long_enough[include_flag,:]
+                # Merging short interval events
+                idx_merged = strong_enough.copy()
+                interval = strong_enough[1:,0] - strong_enough[0:-1,1]
+                id_to_merge = np.where(interval < thr_interval*fs)[0]
+                for j in range(len(id_to_merge))[::-1]:
+                  idx_merged[id_to_merge[j],1]= idx_merged[id_to_merge[j]+1,1]
+                  idx_merged = np.delete(idx_merged, id_to_merge[j]+1,axis=0)
+        idx_spindles.append( idx_merged )
+    return idx_spindles
+```
+
+Now we have multi-channel detection functions for both SWA and Spindles, let's try it on example dataset.
+
+
+```python
+# Sample data selection
+_, EEG_hdeeg, sleep_score = load_data(2, 6, verbose=False) # selection of arbitrary data
+data_filtered = EEG_hdeeg.data.copy() # for visualization
+for chanIdx in range(data_filtered.shape[0]): data_filtered[chanIdx,:]=bandpass_filt(data_filtered[chanIdx,:],band=(1,50))            
+
+# Multi-channel detection of sleep-related EEG activities 
+print('Calculating multi-channel spindle and swas... it takes time (3-5 minutes)')
+idx_spindles = detect_spindle_multichan( EEG_hdeeg.data )
+idx_swas = detect_swa_multichan( EEG_hdeeg.data )
+```
+
+    ** Data loading ... Subject 3, Session 7 (Step: R1) **
+    D:/SciData_Bowon_Sleep/data_BIDS_0812/sub-03/ses-07/eeg/sub-03-R1_screw
+    Calculating multi-channel spindle and swas... it takes time (3-5 minutes)
+    
+
+After running detection scripts, visualization of multi-channel data can be performed as follow. In case of bad channel, we're gonna skip drawing it.
+
+
+```python
+# Visualization 
+plt.figure(figsize=(12,12)); plt.subplot(111)
+spacing=5100
+y_center = np.linspace( -spacing, spacing, int(data_filtered.shape[0]) )
+xlims = (953.5,961.5); t_idx = find_idx( xlims, EEG_hdeeg.times ) # selection of arbitrary time window
+for chanIdx in range(data_filtered.shape[0]):
+    bad_channel = bool( list( EEG_hdeeg.bad_channels ).count(chanIdx) )
+    if not bad_channel:
+        # Plot raw data
+        x_epoch = EEG_hdeeg.times[t_idx[0]:t_idx[1]]
+        y_epoch = data_filtered[chanIdx,t_idx[0]:t_idx[1]]
+        line1 = plot( x_epoch, y_epoch-y_center[chanIdx], color=[.2,.2,.2] )
+        # Mark Spindles
+        spindle_idx = idx_spindles[chanIdx]
+        for evIdx in range(spindle_idx.shape[0]):
+            if (spindle_idx[evIdx,0]>=t_idx[0])&(spindle_idx[evIdx,1]<=t_idx[1]):
+                x = EEG_hdeeg.times[ spindle_idx[evIdx,0]:spindle_idx[evIdx,1] ]
+                y = data_filtered[chanIdx,spindle_idx[evIdx,0]:spindle_idx[evIdx,1]]
+                line2 = plot(x,y-y_center[chanIdx],'m')
+        # Mark SWAs
+        swa_idx = idx_swas[chanIdx]
+        for evIdx in range(swa_idx.shape[0]):
+            if (swa_idx[evIdx,0]>=t_idx[0])&(swa_idx[evIdx,1]<=t_idx[1]):
+                x = EEG_hdeeg.times[ swa_idx[evIdx,0]:swa_idx[evIdx,1] ]
+                y = data_filtered[chanIdx,swa_idx[evIdx,0]:swa_idx[evIdx,1]]
+                line3 = plot(x,y-y_center[chanIdx],'c')
+
+        plt.xticks(color=[0,0,0]); plt.yticks(color=[0,0,0])
+        plt.yticks(y_center, EEG_hdeeg.ch_names[::-1]);
+        plt.ylim((-1.10*spacing,1.10*spacing))
+    
+plt.title('Example multi-channel SWA & Spindle detection')
+plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
+plt.gca().set_facecolor((1,1,1))
+xlabel('Time (sec)'); ylabel('Voltage (uV)')
+plt.savefig( '%sFig_QQ_ExampleMultiDetection.jpg'%(figure_directory), dpi=figure_dpi );
+```
+
+
+![png](figures/output_20_0.png)
+
+
 ### 3-3) Time-frequency transformation (Spectrogram)
 
 In many cases, EEG signals are transformed into time-frequency representation. In our demonstration, we'll calculate delta-band (1-4 Hz) and theta-band (4-13 Hz) powers to compare NREM vs. REM sleep EEGs. Function get_spectrogram() which takes multi-channel time domain EEG data (channel * time) to produce time-frequency domain EEG data (time * frequency * channel) is provided below. 
@@ -454,7 +647,7 @@ for chanIdx in range(1):
     
 
 
-![png](figures/output_16_1.png)
+![png](figures/output_22_1.png)
 
 
 ### 3-4) Topographical reconstruction of HD-EEG 
@@ -497,7 +690,7 @@ plt.savefig( '%sFig_QQ_TutorialLabeling.jpg'%(figure_directory), dpi=figure_dpi
 ```
 
 
-![png](figures/output_18_0.png)
+![png](figures/output_24_0.png)
 
 
 Next, we need to extract band power information of each channel to draw topographical map. To easily extract specific band power with given state, we made a handy function bandpower_state() as below. It slices in temporal and frequency domain for corresponding point of sleep state and frequency bin, respectively. 
@@ -632,7 +825,7 @@ class bi_interp2:
 """ (2) Function for Topography plot """
 from matplotlib import cm
 from pandas import read_csv
-montage_file = bids_path+'montage.csv'
+montage_file = 'D:/SciData_Bowon_Sleep/montage.csv'
 from mpl_toolkits.mplot3d import Axes3D
 
 def get_boundary():
@@ -722,7 +915,7 @@ plt.savefig( '%sFig_QQ_TutorialTopography.jpg'%(figure_directory), dpi=figure_dp
 ```
 
 
-![png](figures/output_24_0.png)
+![png](figures/output_30_0.png)
 
 
 ### 3-5) Example trace figures
@@ -731,6 +924,13 @@ plt.savefig( '%sFig_QQ_TutorialTopography.jpg'%(figure_directory), dpi=figure_dp
 Below figures depict some example cases of EEG analysis (e.g., drawing multichannel raw data, rich visualization of SWA and spindle) contained in this dataset. For the purpose of demonstration, few pieces of data were arbitrarily selected. 
 
 
+```python
+example_info_trace = { 'subject':1, 'day':6 }
+example_info_swa = { 'subject':0, 'day':0, 'event':249 }
+example_info_spindle = { 'subject':2, 'day':6, 'event':3 }
+```
+
+
 ```python
 # Figure: Example signal visualization of screw/hdeeg data
 
@@ -787,7 +987,7 @@ for chanIdx in range(y.shape[0]):
 
 state_indicator_ypos = spacing*1.23
 plot( x[wake_toi], np.ones((wake_toi.shape[0],1))*state_indicator_ypos, '.', color=[.1, 1, .1])
-plot( x[sleep_toi], np.ones((sleep_toi.shape[0],1))*state_indicator_ypos,2 '.', color=[.5,.5,.5])
+plot( x[sleep_toi], np.ones((sleep_toi.shape[0],1))*state_indicator_ypos, '.', color=[.5,.5,.5])
 plt.gca().set_facecolor((1,1,1))
 plt.savefig( '%sFig_QQ_ExampleTrace.jpg'%(figure_directory), dpi=figure_dpi )
 plt.show();
@@ -803,7 +1003,7 @@ plt.show();
     
 
 
-![png](figures/output_26_1.png)
+![png](figures/output_33_1.png)
 
 
 
@@ -860,7 +1060,7 @@ plt.show();
     
 
 
-![png](figures/output_27_1.png)
+![png](figures/output_34_1.png)
 
 
 
@@ -888,8 +1088,8 @@ for chanIdx in range(y_epoch.shape[0]):
     plt.plot(x_epoch-t[spindle_idx[evIdx,0]], y_epoch[chanIdx,:]-shift, color=color_space[chanIdx,], linewidth=1);
     plt.xlabel('Time (sec)', color=[0,0,0])
     plt.xticks(color=[0,0,0]); plt.yticks(color=[0,0,0])
-    plt.yticks(y_center, EEG_hdeeg.ch_names[::-1]);
-    plt.ylim((-1.05*spacing,1.20*spacing))
+plt.yticks(y_center, EEG_hdeeg.ch_names[::-1]);
+plt.ylim((-1.05*spacing,1.20*spacing))
 plot( [0,0],plt.ylim(), color=(1,0,0), linestyle = '--' )
 plot( [0,0]+t[spindle_idx[evIdx,1]]-t[spindle_idx[evIdx,0]],plt.ylim(), color=(1,0,0), linestyle = '--' )
 plt.text(0, 1.12*spacing, 'Onset > ', horizontalalignment = 'right', color=(1,0,0))
@@ -903,8 +1103,15 @@ plt.show();
 
 ```
 
+    ** Data loading ... Subject 3, Session 7 (Step: R1) **
+    Screw data shape: [4 channels x 43200000 times]
+    Screw bad channels: []
+    hdEEG data shape: [36 channels x 3600000 times]
+    hdEEG bad channels: [9, 24, 27, 30, 34, 35]
+    
+
 
-![png](figures/output_28_0.png)
+![png](figures/output_35_1.png)
 
 
 ### 3-6) Obtaining grand-averaged results 
@@ -1137,21 +1344,19 @@ plt.show()
     
 
 
-![png](figures/output_33_1.png)
+![png](figures/output_40_1.png)
 
 
     REM theta
     
 
 
-![png](figures/output_33_3.png)
+![png](figures/output_40_3.png)
 
 
 Now, do it on your own!
 
 
 ```python
-# Code block
-   
 
-```
+```