上传文件至 'ECoG/code'

ECoG/code/EphysSparseTWFbatchAnalysis.m

@@ -0,0 +1,132 @@
+%% this batch file is meant to be run in conjunction with TWF_ECOG_sparse_analysis.py
+% to clean TWF data using Alain's DSS algorithm
+
+% first, define some datapaths and constants, and load the list of data files.
+
+
+dirRoot='/ephys/2';
+if ~isdir(dirRoot)
+    dirRoot='/home/colliculus/ephys/2';
+end;
+addpath([dirRoot '/Cecilia_TWF/Analysis/NoiseTools'])
+%addpath('/home/jan/jan/matlab')
+addpath('/home/colliculus/behaviourBoxes/software/ratCageProgramsV2')
+dataroot=[dirRoot '/Cecilia_TWF/']
+%
+blankChans=[8,57,64]; % these channels are not connected on teh Viventi ecog electrodes and should therefore be excluded from analysis
+goodChans=1:64; goodChans(blankChans)=[];
+
+%condstr='trained_ITD'
+condstr='ITD';
+
+analysisPath=[dirRoot '/Cecilia_TWF/Analysis/Analysis_AC_ECoG_TWF_' condstr '_SparseData']
+%savePath = '/ephys/2/Cecilia_TWF/Analysis/Analysis_AC_ECoG_TWF_'+cond+'_SparseData/swp_stim_allData'
+filepaths=importdata([analysisPath, '/AC_ECoG_TWF_' condstr '_SparseData_filesToAnalyse.csv'])
+
+filepath=filepaths{1};
+%% 
+for ff=1:length(filepaths) 
+    filepath=filepaths{ff};
+    disp(['denoising ' filepath ' ...'])
+    % read in the LFP filtered ephys data
+    [fpath, fname,ext]=fileparts(filepath);
+    nextDataFile=[dataroot fpath '/LFP_' fname '.ephys'];
+    data=readEphysFiles(nextDataFile);
+%     if size(data(1).signal,1) == nsamples
+%         disp('... already cleaned.')
+%         continue % file already done
+%     end
+    % read in the stimulus grid
+    %stims=importdata([dataroot filepath]); % load itd codes
+    stim=readtable([dataroot filepath]);
+    %hasPresilence=sum(strcmp('preSilence_s_',stim.Properties.VariableNames));
+    firstITDpar=find(strcmp([condstr(end-2:end) '0'],stim.Properties.VariableNames));
+    itds=stim{:,firstITDpar:firstITDpar+3};
+    % sometimes the last sweep is too short and has been dropped in the
+    % python preprocessing. If that is the case, the last ITD needs to be
+    % dropped too.
+    while size(itds,1) > length(data)
+        itds(end,:)=[];
+    end
+    cond=1+(itds>0)*[8 4 2 1]'; % 1 to 16
+
+    %% nt_dss1 wants a time x channel x trials matrix. Let's build one
+%     if hasPresilence
+%         % we cut out the preSilence, if there is one
+%         presamples=round(stim.preSilence_s_(1)*data(1).sampleRate);
+%     else
+%         presamples=0;
+%     end
+    nsamples=1300;
+    for trial=1:length(data)
+        if size(data(trial).signal,1) < nsamples
+            nsamples=size(data(trial).signal,1);
+        end
+    end;
+
+    x=zeros(nsamples,length(goodChans), length(data));
+    for trial=1:length(data)
+        x(:,:,trial)=data(trial).signal(1:nsamples,goodChans);
+    end
+    %% now try denoising
+
+    NKEEP=11;
+    % DSS to enhance repeatability
+    c0=nt_cov(x)/size(x,3);
+    c1=zeros(size(x,2));
+    for iCond=1:16
+        c1=c1+nt_cov(mean(x(:,:,find(cond==iCond)),3));
+    end
+    c1=c1/16;
+    [todss,pwr0,pwr1]=nt_dss0(c0,c1);
+    fromdss=pinv(todss);
+    x_clean=nt_mmat(x,todss(:,1:NKEEP)*fromdss(1:NKEEP,:)); % 'cleaned' data
+    
+%     % DSS to enhance between-condition differences
+%     z=nt_mmat(x,todss(:,1:NKEEP));
+%     zz=[];
+%     for iCond=1:16
+%         zz(:,:,iCond)=mean(z(:,:,find(iCond==cond)),3);
+%     end
+%     [todss2,pwr0,pwr1]=nt_dss1(zz);
+%     fromdss=pinv(todss2);
+%     x_clean=nt_mmat(x,todss2(:,1:NKEEP)*fromdss(1:NKEEP,:)); % 'cleaned' data
+
+%     figure(1); clf
+%     plot(pwr1./pwr0, '.-'); xlabel('component'); ylabel('score'); title('DSS repeatibility')
+
+
+
+    % 
+    % %% what do the results look like?
+    % for chan=1:size(x,2)
+    %     figure(2); clf;
+    %     raw=mean(squeeze(x(:,chan,:)),2);
+    %     filtered=mean(squeeze(x_clean(:,chan,:)),2);
+    %     chan=chan+1;
+    %     subplot(1,2,1);
+    %     plot(raw); title(sprintf('channel %d',chan))
+    %     subplot(1,2,2);
+    %     plot(filtered);
+    %     pause
+    % end;
+    % 
+    % %%
+    % figure(2); clf;
+    % 
+    % subplot 121;
+    % plot(mean(x,3)); title('data');
+    % subplot 122;
+    % nt_bsplot(z(:,20,:)); title('recovered'); 
+
+    %% save the results back to the LFP file 
+    for trial=1:length(data)
+        data(trial).signal=x_clean(:,:,trial);
+    end
+    disp('saving cleaned data.');
+    nextDataFile=[dataroot fpath '/CLEAN_' fname '.ephys'];    
+    writeEphysFiles(nextDataFile,data);
+end;
+
+disp('All done!');
+

ECoG/code/UniveratiateECOGanalysis.py

@@ -0,0 +1,605 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Dec 10 15:02:37 2019
+
+@author: Cecilia and Jan
+
+
+The overall workflow is as follows:
+    1) ITD TWF data files are LFP filtered and written back to disk
+    2) In Matlab, LFP filtered data files are cleaned using Alain's DSS algorithm
+    3) Cleaned datafiles are read back in and LFP response amplitudes are computed
+    4) LFP response amplitudes are read back in and analysed with linear regression.
+
+"""
+
+#%% Step 0: import libraries and set paths
+import os
+from glob import glob
+import csv
+
+def zscore(X):
+    return (X-X.mean())/X.std()
+
+
+blankChans=[8,57,64] # these channels aren't connected ont eh Viventi electodes
+## create a datalog of all recorded filenames
+#PATH = r'Z:\ephys\2\Cecilia_TWF'
+PATH = r'Z:\home\colliculus\ephys\2\Cecilia_TWF'
+if not os.path.exists(PATH):
+    PATH='/ephys/2/Cecilia_TWF'
+if not os.path.exists(PATH):
+    PATH='/twinkle/ephys/2/Cecilia_TWF'
+if not os.path.exists(PATH):
+    PATH='/home/colliculus/ephys/2/Cecilia_TWF'    
+EXT = "*.csv"
+
+
+import statsmodels.api as sm
+import os   
+import sys
+        
+#    codeDir='//twinkle.bms.cityu.edu.hk/colliculus/behaviourBoxes/software/ratCageProgramsV2' #directory at ephys station
+codeDir='/home/colliculus/behaviourBoxes/software/ratCageProgramsV2' #directory at ephys station
+
+if not os.path.exists(codeDir):
+    #codeDir='d:/jan/behavbox/ratCageProgramsV2' #directory at ephys station
+    codeDir = 'Z:/behaviourBoxes/software/ratCageProgramsV2'  ## directory at CC's desktop
+    
+if not os.path.exists(codeDir):
+    codeDir = 'Z:/home/colliculus/behaviourBoxes/software/ratCageProgramsV2' ## directory at CC's laptop
+    
+if not os.path.exists(codeDir):
+    raise BaseException('No code directory')    
+
+sys.path.append(codeDir)
+
+import numpy as np
+import pandas as pd
+from matplotlib import pyplot as plt
+import RZ2ephys as ep
+import ntpath 
+import copy   
+#% Pre-step 1: read in tables listing relevant datafiles provided by Cecilia. Put into variable "filepaths"
+filepaths=[];
+#analysisPath=PATH+'/Analysis/Analysis_AC_ECoG_TWF_ITD_SparseData'
+# cond='ITD'
+#cond='ILD'
+cond='trained_ITD'
+# cond='trained_ILD'
+analysisPath=PATH+'/Analysis/Analysis_AC_ECoG_TWF_'+cond+'_SparseData'
+savePath = r'Z:\ephys\2\Cecilia_TWF\Analysis\Analysis_AC_ECoG_TWF_'+cond+'_SparseData\swp_stim_allData'
+if not os.path.exists(savePath):
+    savePath = '/ephys/2/Cecilia_TWF/Analysis/Analysis_AC_ECoG_TWF_'+cond+'_SparseData/swp_stim_allData'
+if not os.path.exists(savePath):
+    savePath = '/home/colliculus/ephys/2/Cecilia_TWF/Analysis/Analysis_AC_ECoG_TWF_'+cond+'_SparseData/swp_stim_allData'
+#with open(analysisPath+'/AC_ECoG_TWF_ITD_SparseData_filesToAnalyse.csv', newline='') as myFile:
+with open(analysisPath+'/AC_ECoG_TWF_'+cond+'_SparseData_filesToAnalyse.csv', newline='') as myFile:
+#with open('Z:/ephys/2/Cecilia_TWF/Analysis/Analysis_AC_ECoG_TWF_trained_ITD_SparseData/AC_ECoG_TWF_trained_ITD_SparseData_datalog.csv', newline='') as myFile:
+    reader = csv.reader(myFile)
+    for row in reader:                  
+        filepaths.append(row[0])
+        
+#filePath=filepaths[2]
+
+#%% Step 1: Loop through data files, calculate LFP and write back to disk    
+
+for filePath in filepaths: 
+#%               
+    filename = os.path.splitext(filePath)[0]
+    fname = ntpath.basename(filename)
+    LFPfilePath, LFPfileName = os.path.split(filePath)
+    LFPfileName='LFP_'+LFPfileName[:-4]
+    
+    print("===   ===   ===")        
+    print('Processing data from ' + fname)
+#    swp, stim = ep.readEphysFile((PATH+'/'+filePath)[:-4])
+    swp, stim = ep.readEphysFile(filePath[:-4])     
+    print('Reading in data done!')
+    
+    print('Start to match swp and stim.')
+    #stim=stim.loc[:, ~stim.columns.str.contains('timeStamp')] 
+    # get rid of timeStamps. We won't need them and they are all unique
+    
+    #% there is likely to be a pre-sweep. For analysis purposes we throw that out
+    if len(swp)-stim.shape[0] == 1:
+        swp=swp[1:]
+        print('Matching done!')
+    else:        
+        print('Warning: unexpected that number of sweeps does not equal 1-row of stim')
+    
+    #% due to a bug in early versions of the recording program
+    # the last sweep in some files may be too short.
+    # If so, drop it. We use as criterion that the last sweep should at least have 
+    # of 0.8 s
+    
+    print('Taking care of the last swp.')
+    if swp[-1].signal.shape[0]/swp[-1].sampleRate  < 0.8:
+        swp=swp[:-1]
+        stim=stim[:-1]
+    print('Trimming done!')
+    
+#    print('Applying AMUA filter.')            
+#    AMUAFilters=swp[0].AMUAFilterCoeffs()
+#    for ii in range(len(swp)):
+#        swp[ii]=swp[ii].calcAMUA(AMUAFilters)
+    
+    print('Applying LFP filter.')
+    lpfCoefs=swp[1].LFPfilterCoeffs()
+    for ss in range(1,len(swp)):
+         swp[ss]=swp[ss].calcLFP(lpfCoefs,downsample=4)       
+         
+#% save LFP data 
+    
+    #ep.writeEphysFile(PATH+'/'+LFPfilePath+'/'+LFPfileName,swp)
+    ep.writeEphysFile(LFPfilePath+'/'+LFPfileName,swp)        
+
+print('===================================')
+print('======Exporting LFPs complete =====')
+print('===================================')
+#%% Step 2: Now carry out step 2 by running EphysSparseTWFbatchAnalysis.m in Matlab
+    # once that's done:
+    
+    
+    
+    
+#%% Step 3: Read in cleaned data files and compute TWF with regression:
+#ff=0
+##%%  
+#    
+#filePath=filepaths[ff]   
+#ff+=1
+#%
+    
+from pathlib import PurePath
+plotLFPs=False
+# create column names and row names 
+
+
+sc=cond[-3:]
+rowName = [sc+' 0',sc+' 1',sc+' 2',sc+' 3']
+
+
+# processing 
+
+     
+for filePath in filepaths:  
+#%
+    #            
+    filename = os.path.splitext(filePath)[0]
+    fname = ntpath.basename(filename)
+    LFPfilePath, LFPfileName = os.path.split(filePath)
+    LFPfileName='CLEAN_'+LFPfileName[:-4]
+    
+    
+    #% step 3.1, read in cleaned data
+    print("===   ===   ===")        
+    print('Processing data from ' + LFPfileName)
+
+    
+    addPath = PurePath(LFPfilePath)       
+    readInPath = PurePath(PATH).joinpath(addPath)
+    
+    swpPath = PurePath(readInPath).joinpath(LFPfileName)
+    swpFilename = swpPath.with_suffix('.ephys').as_posix()
+    
+    stimPath = PurePath(readInPath).joinpath(fname)
+    stimFilename = stimPath.with_suffix('.csv').as_posix()
+
+
+## read in .ephys file
+#    swp, ignore = ep.readEphysFile(PATH+'/'+LFPfilePath+ '/'+LFPfileName)    
+#    stim = ep.readEphysStimTable((PATH+'/'+filePath)[:-4])   
+
+
+#%
+#    swp, ignore = ep.readEphysFile(swpFilename)  ## this causes error, no stim file found. The CLEAN dataset does not have a corresponding .csv file.
+    
+    swp =np.array([])      
+    with open(swpFilename,'rb') as fileObject:
+        swpN=ep.MCsweep(fileObject)
+        while not swpN.signal is None:
+            swp=np.append(swp,swpN)
+            swpN=ep.MCsweep(fileObject)
+    fileObject.close()
+            
+    # read in corresponding stim param table  
+    stim = ep.readEphysStimTable(stimFilename)   
+ 
+    
+        
+    if stim.shape[0] > len(swp):
+    
+        # last sweep may have been dropped. If so, drop last stim too
+        stim.drop(stim.index[-1])
+    print('Reading in data done!')
+   
+#%
+    meanLFP=copy.copy(swp[0].signal)
+    nswp=1
+    for ss in range(1,len(swp)):
+        meanLFP+=swp[ss].signal
+        nswp+=1
+    meanLFP=meanLFP/nswp
+    
+    if plotLFPs:
+     ## plotand save  meanLFP    
+        plt.figure(1)
+        plt.clf()
+        dchan=1
+        taxis=np.array(range(meanLFP.shape[0]))/swp[0].sampleRate
+        for chan in range(swp[0].signal.shape[1]): # loop through channels
+            plt.subplot(8,8,dchan)
+            if dchan==4:
+                plt.title(filename)
+            dchan+=1
+            if dchan in blankChans:
+                dchan += 1
+            plt.plot(taxis,meanLFP[:,chan])
+            plt.yticks([])
+    
+        plt.pause(0.05)
+        
+        if not os.path.exists(analysisPath+'/meanLFPfigures/'):
+            os.makedirs(analysisPath+'/meanLFPfigures/')
+        plt.savefig(analysisPath+'/meanLFPfigures/'+fname+'.svg',format='svg')
+    
+        #% step 3.2, compute responses 
+        print('Generate base line and response channel data.')
+        if 'preSilence (s)' in stim.columns : 
+            stimStart=stim['preSilence (s)'][0]
+        else:
+            stimStart=0
+
+#%
+    respWin=stimStart+np.array([0.001,0.1]) # resp window 
+#    baselineWin=stimStart+np.array([0.65,-1])  
+   
+    
+    # create np array header
+    resp=[]
+#    baseline=[]
+    for cc in range(swp[0].signal.shape[1]):
+        resp.append('RespChan{}'.format(cc+1))
+#        baseline.append('BaselineChan{}'.format(cc+1))
+    respM=np.array(resp)
+#    baselineM=np.array(baseline)
+    
+    for ss in range(swp.shape[0]):
+#        resp= swp[ss].response(respWin)
+#        resp= swp[ss].sigSnip(respWin).std(axis=0)
+        resp=np.sqrt(np.mean(swp[ss].sigSnip(respWin)**2,axis=0))
+        respM=np.vstack([respM,resp])
+#        baseline= swp[ss].response(baselineWin)
+#        baselineM=np.vstack([baselineM,baseline])
+    respD=pd.DataFrame(respM[1:,:],columns=respM[0,:])
+#    baselineD=pd.DataFrame(baselineM[1:,:],columns=baselineM[0,:])
+    #%
+    allD=pd.merge(stim,respD,left_index=True,right_index=True)
+    for cc in range(swp[0].signal.shape[1]):
+        allD['RespChan{}'.format(cc+1)]=allD['RespChan{}'.format(cc+1)].astype(float)
+        
+    #% save all LPF response data to csv
+    
+    saveFile = fname + '_LFPresponses.csv'
+    allD.to_csv(os.path.join(savePath, saveFile))  
+    print('All dataset saved!')
+    
+    
+    
+
+print('====================================')
+print('====== Analysis batch complete =====')
+print('====================================')        
+
+#%% step 3a - Let's look through the data to find an example with a big difference in LFP amplitude for all pos and all neg ITDs.
+# for illustrative purposes
+from scipy.stats import ranksums
+sc=cond[-3:]
+rowName = [sc+' 0',sc+' 1',sc+' 2',sc+' 3']
+filePath=filepaths[11]
+# filePath=filepaths[3]
+filename = os.path.splitext(filePath)[0]
+fname = ntpath.basename(filename)    
+if fname.find("900")<0:
+    clickRate=300
+else:
+    clickRate=900
+loadFile = fname + '_LFPresponses.csv'
+print('loading LFP responses from ',loadFile)
+allD=pd.read_csv(os.path.join(savePath, loadFile))
+sc=cond[-3:]
+
+for chan in [3, 6, 9, 12 ,15, 18, 21, 24, 27, 30]:
+    posITDs=allD[(allD["ITD 0"]>0) & (allD["ITD 1"]>0) & (allD["ITD 2"]>0) & (allD["ITD 3"]>0)]
+    negITDs=allD[(allD["ITD 0"]<0) & (allD["ITD 1"]<0) & (allD["ITD 2"]<0) & (allD["ITD 3"]<0)]
+    mu1=posITDs["RespChan{}".format(chan)].mean()
+    mu2=negITDs["RespChan{}".format(chan)].mean()
+    s1=posITDs["RespChan{}".format(chan)].std()
+    s2=negITDs["RespChan{}".format(chan)].std()
+    print((mu1-mu2)/np.sqrt(0.5*(s1**2+s2**2)))
+
+chan=29
+# chan=3
+plt.clf()
+plt.boxplot([posITDs["RespChan{}".format(chan)], negITDs["RespChan{}".format(chan)]])
+ranksums(posITDs["RespChan{}".format(chan)], negITDs["RespChan{}".format(chan)])
+#%% so with trial and error we found that filePath 11 channel 29 gives reasonably distinct LFP amplitudes for the all pos vs all neg ITDs on average.
+#  Let's check whether these would make for a good figure to illustrate the logic of the regression analysis
+chan=29
+filePath=filepaths[11]
+LFPfilePath, LFPfileName = os.path.split(filePath)
+LFPfileName='CLEAN_'+LFPfileName[:-4]
+filename = os.path.splitext(filePath)[0]
+fname = ntpath.basename(filename)    
+if fname.find("900")<0:
+    clickRate=300
+else:
+    clickRate=900
+loadFile = fname + '_LFPresponses.csv'
+print('loading LFP responses from ',loadFile)
+allD=pd.read_csv(os.path.join(savePath, loadFile))
+swpFilename = PATH+'/'+LFPfilePath+'/'+LFPfileName+'.ephys'   
+stimPath = PATH+'/'+LFPfilePath+'/'+fname+'.csv'
+
+
+## read in .ephys file
+#    swp, ignore = ep.readEphysFile(PATH+'/'+LFPfilePath+ '/'+LFPfileName)    
+#    stim = ep.readEphysStimTable((PATH+'/'+filePath)[:-4])   
+
+
+#%
+#    swp, ignore = ep.readEphysFile(swpFilename)  ## this causes error, no stim file found. The CLEAN dataset does not have a corresponding .csv file.
+posIdx=np.where((allD["ITD 0"]>0) & (allD["ITD 1"]>0) & (allD["ITD 2"]>0) & (allD["ITD 3"]>0))[0]
+negIdx=np.where((allD["ITD 0"]<0) & (allD["ITD 1"]<0) & (allD["ITD 2"]<0) & (allD["ITD 3"]<0))[0]
+pSwp = []
+nSwp = []
+swpIdx=0     
+with open(swpFilename,'rb') as fileObject:
+    swpN=ep.MCsweep(fileObject)
+    while not swpN.signal is None:
+        if swpIdx in posIdx:
+            pSwp.append(swpN.signal[:,chan])
+        if swpIdx in negIdx:
+            nSwp.append(swpN.signal[:,chan])
+        swpIdx+=1
+        swpN=ep.MCsweep(fileObject)
+fileObject.close()
+pSwp=np.array(pSwp)
+nSwp=np.array(nSwp)
+#%% now plot
+sampleRate= 6103.515625
+posMean=pSwp.mean(axis=0)*1000
+negMean=nSwp.mean(axis=0)*1000
+posSEM=pSwp.std(axis=0)/np.sqrt(pSwp.shape[0])*1000
+negSEM=nSwp.std(axis=0)/np.sqrt(nSwp.shape[0])*1000
+taxis=np.arange(pSwp.shape[1])/sampleRate*1000
+fig=plt.figure(4, figsize=(4,4))
+plt.clf()
+plt.fill_between(taxis, posMean - posSEM, posMean + posSEM,  color='red', alpha=0.5)
+plt.fill_between(taxis, negMean - negSEM, negMean + negSEM,  color='blue', alpha=0.5)
+plt.xlim([0,30])
+plt.xlabel('Time (ms)')
+plt.ylabel('mV')
+plt.legend(['all ITDs negative','all ITDs positive'])
+plt.savefig('Figure4.tif',dpi=300)
+plt.savefig('Figure4.svg')
+#%%   step 4 fit regression model  and make figures
+        
+#    allD=pd.merge(allD,baselineD,left_index=True,right_index=True)
+#    for cc in range(swp[0].signal.shape[1]):
+#        allD['BaselineChan{}'.format(cc+1)]=allD['BaselineChan{}'.format(cc+1)].astype(float)   
+    
+sc=cond[-3:]
+rowName = [sc+' 0',sc+' 1',sc+' 2',sc+' 3']
+
+anlResults = pd.DataFrame([ ])
+chAll = list(range(1,65))
+chBlank = {8,57, 64} 
+chList = [ele for ele in chAll if ele not in chBlank] 
+chList.reverse() # reversing the list order means that later chan 1 will be plotted last, which is good for sticking labels on panels.
+colNameP = []
+colNameB = []
+colNameS = []
+for chN in range(1,65): # naming loop through channels                   
+    chP = 'P_chan' + str(chN)  
+    colNameP.append(chP)
+    
+    chB = 'Beta_chan' + str(chN) 
+    colNameB.append(chB)
+    
+    chS = 'Sig_chan' + str(chN)
+    colNameS.append(chS)
+# we also create a CSV file that keeps all the betas and pvals for summary statistics
+betaFile=open(analysisPath +'/betas_and_pvals.csv','w')
+betaFile.write("filename,clickrate,channel,beta1,beta2,beta3,beta4,pval1,pval2,pval3,pval4\n")
+for filePath in filepaths:  
+    filename = os.path.splitext(filePath)[0]
+    fname = ntpath.basename(filename)    
+    if fname.find("900")<0:
+        clickRate=300
+    else:
+        clickRate=900
+    loadFile = fname + '_LFPresponses.csv'
+    print('loading LFP responses from ',loadFile)
+    allD=pd.read_csv(os.path.join(savePath, loadFile))
+    sc=cond[-3:]
+    chan=1
+    
+    #% Ordinary least squares regression
+    # build a chan x ITD array of pvals, a chan x ITD array of regression coeffs,  and a chan x ITD array of significance outcomes
+    pvals=np.zeros((4,64))
+    betas=np.zeros((4,64))
+    sigs = np.zeros((4,64))
+    for chan in chList: # loop through channels
+        responses=np.array(allD['RespChan'+str(chan+1)])
+        responses=np.log(responses)
+        responses=zscore(responses)
+        # find indeces of non-outliers (keep)
+        criterion=np.median(responses)+np.std(responses)*3
+        keep=np.where(responses<criterion)
+        responses=responses[keep]
+        model = sm.OLS(responses, sm.add_constant(allD[[sc+' 0',sc+' 1',sc+' 2',sc+' 3']].loc[keep])).fit()
+        # model = sm.OLS(responses, sm.add_constant(allD[[sc+' 0',sc+' 1',sc+' 2',sc+' 3']])).fit()
+        #model = sm.OLS(allD['RespChan7'], sm.add_constant(allD[['ILD 0','ILD 1','ILD 2','ILD 3']])).fit()
+        pvals[:,chan]=model.pvalues[1:]
+        #betas[:,chan]= -model.params[[sc+' 0',sc+' 1',sc+' 2',sc+' 3']]
+        betas[:,chan]= -model.params[1:]
+        betaFile.write("{},{},{}".format(fname,clickRate,chan))
+        for clck in range(4):
+            betaFile.write(",{}".format(betas[clck,chan]))
+        for clck in range(4):
+            betaFile.write(",{}".format(pvals[clck,chan]))
+        betaFile.write("\n")
+
+
+    ##  save regression results
+    # turn significant information into dataframe, 0 = no significance, 1 = significance level P < 0.05
+        
+    for par in range(pvals.shape[0]):
+        if pvals[par,chan] < 0.05:
+           sigs[par,chan] = 1
+           
+    sdf = pd.DataFrame(sigs, columns = colNameS)
+    
+    # turn betas and pvalues with column names into dataframe        
+    bvdf = pd.DataFrame(betas, columns = colNameB)
+    pvdf = pd.DataFrame(pvals, columns = colNameP)
+    
+    # concat betas, pvalues and significan indicates into a big dataframe
+    
+    results = pd.concat([bvdf,pvdf,sdf], axis=1) 
+    
+    # add rowname and the processed data filename as column          
+    results.insert(0,'Condition', fname)
+    results.insert(1, 'Click_N', rowName)
+    
+    # concat results from different datafiles  
+    anlResults = pd.concat([anlResults, results], ignore_index=True)
+      
+       
+  
+    
+  ## plot and save TWF weight figures  
+    plt.figure(2, figsize=(6,6))
+    plt.clf();  
+    plt.rc('font',size=9)
+    ylim=v=np.max(np.abs(betas))*1.05
+    for dchan in chList: # loop through channels         
+        plt.subplot(8,8,dchan)
+        plt.plot([0,3],[0,0],'k:')
+        if dchan==4:
+            plt.title(fname)
+        #%
+        plt.plot(betas[:,dchan])
+        plt.ylim([-ylim,ylim]) 
+        #plt.xlabel('Click number')
+        #plt.ylabel('Beta value')
+        plt.xticks(np.arange(4), ('1', '2', '3', '4'))
+        
+        # unshow the x and y tick labels in channels other than chan1
+        if not (dchan in [1,9,17,25,33,41,49]):
+           plt.yticks([])
+        if dchan < 57:
+            if not(dchan in [49,56]):
+                plt.xticks([]) 
+            
+        # highlight significant values in red
+        for par in range(pvals.shape[0]):
+            if pvals[par,dchan] < 0.05:
+                plt.plot(par,betas[par,dchan],'r*')
+    #%plt.title('Channel '+str(chan))
+        if dchan==25:
+            plt.ylabel('Weight β (1/ms)')
+        if dchan==60:
+            plt.xlabel('Click #')
+            # if dchan==1:
+            #     plt.text(-4,2,letter[index],fontsize=16)    
+
+    plt.pause(0.05)
+    if not os.path.exists(analysisPath+'/TWFweightFigures/'):
+        os.makedirs(analysisPath+'/TWFweightFigures/')
+    plt.savefig(analysisPath+'/TWFweightFigures/'+fname+'.svg',format='svg')
+    # select figures to save as panels for Fig 4
+    if fname=='trained_twf_1801_ITD_900Hz_leftCortex':
+        plt.text(-4,ylim*2,'A',fontsize=16)
+        plt.savefig('/home/jan/jan/document/science/cityu/TWF_ITD_ILD_NormalHearingRats/figures/Fig5A.svg')
+    if fname=='trained_twf_1803_ITD_900Hz':
+        plt.text(-4,ylim*2,'B',fontsize=16)
+        plt.savefig('/home/jan/jan/document/science/cityu/TWF_ITD_ILD_NormalHearingRats/figures/Fig5B.svg')
+    if fname=='trained_twf_1804_ITD_900Hz_leftCortex':
+        plt.text(-4,ylim*2,'C',fontsize=16)
+        plt.savefig('/home/jan/jan/document/science/cityu/TWF_ITD_ILD_NormalHearingRats/figures/Fig5C.svg')
+    if fname=='trained_twf_1803_ITD_300Hz_leftCortex':
+        plt.text(-4,ylim*2,'D',fontsize=16)
+        plt.savefig('/home/jan/jan/document/science/cityu/TWF_ITD_ILD_NormalHearingRats/figures/Fig5D.svg')
+                    
+## save resuls dataframe as csv file
+anlResults.to_csv(analysisPath + "/" + cond + "_cleanedLFP_TWFregressionResults.csv", index = False)
+betaFile.close()
+print('=======================================')
+print('====== Creating OLS TWFs complete =====')
+print('=======================================')    
+
+#%% Let's read in the regression results and make boxplots for Figure 6
+
+
+
+alpha=0.05
+bNp=np.genfromtxt('/twinkle/ephys/2/Cecilia_TWF/Analysis/Analysis_AC_ECoG_TWF_ITD_SparseData/betas_and_pvals.csv', delimiter=',')
+bNp2=np.genfromtxt('/twinkle/ephys/2/Cecilia_TWF/Analysis/Analysis_AC_ECoG_TWF_trained_ITD_SparseData/betas_and_pvals.csv', delimiter=',')
+allBeta=np.vstack((bNp[1:,3:7],bNp2[1:,3:7]))
+allP=np.vstack((bNp[1:,7:12],bNp2[1:,7:12]))
+cRate=np.hstack((bNp[1:,1],bNp2[1:,1]))
+del bNp, bNp2
+#%%
+def myBoxPlot(X):
+    bp=plt.boxplot(X)
+    plt.ylim([-0.01,1.2])
+    for ii in range(len(bp['medians'])):
+        bp['medians'][ii].set_color([0.4,0.4,0.4])
+        bp['medians'][ii].set_linewidth(2)
+        bp['fliers'][ii].set_markersize(2)
+        bp['fliers'][ii].set_markeredgecolor([0.6,0.6,0.6])
+    return bp
+
+idx300=np.where(cRate==300)
+beta300=allBeta[idx300,:].squeeze()
+pv300=allP[idx300,:].squeeze()
+idx900=np.where(cRate==900)
+beta900=allBeta[idx900,:].squeeze()
+pv900=allP[idx900,:].squeeze()
+plt.figure(5,figsize=(6,6))
+plt.clf()
+ax=plt.subplot(2,2,1)
+myBoxPlot(np.abs(beta300))
+plt.ylabel('Weight β (1/ms))')
+plt.xticks([])
+plt.title('All β, 300 Hz')
+ax=plt.subplot(2,2,2)
+myBoxPlot(np.abs(beta900))
+plt.xticks([])
+plt.yticks([])
+plt.title('All β, 900 Hz')
+sig300=(pv300<0.05)
+sig900=(pv900<0.05)
+sigWeights300=[]
+sigWeights900=[]
+for clkIdx in range(4):
+    sigWeights300.append(beta300[np.where(sig300[:,clkIdx]),clkIdx].squeeze())
+    sigWeights900.append(beta900[np.where(sig900[:,clkIdx]),clkIdx].squeeze())
+ax=plt.subplot(2,2,3)
+myBoxPlot(np.abs(sigWeights300))
+plt.xlabel('Click #')
+plt.ylabel('Weight β (1/ms))')
+plt.title('Significant β, 300 Hz')
+
+ax=plt.subplot(2,2,4)
+bp=myBoxPlot(np.abs(sigWeights900))
+plt.yticks([])
+plt.xlabel('Click #')
+plt.title('Significant β, 900 Hz')
+
+#%%
+plt.savefig('/home/jan/jan/document/science/cityu/TWF_ITD_ILD_NormalHearingRats/figures/BoxplotsFig6.svg')
+plt.savefig('/home/jan/jan/document/science/cityu/TWF_ITD_ILD_NormalHearingRats/figures/BoxplotsFig6.png')

ECoG/code/ecog_decoding_figure_allpanels.m

@@ -0,0 +1,454 @@
+% root='/ephys/2/Cecilia_TWF/Analysis'
+% if ~exist(root,'dir')
+%     root=['/home/colliculus/' root]
+% end
+% cd([root '/Analysis_AC_ECoG_TWF_ITD_SparseData/decoding'])
+% addpath(root)
+% fsz = 14; % font size %10
+%%
+% working on Cecilia's desktop
+cd Z:/ephys/2/Cecilia_TWF/Analysis/Analysis_AC_ECoG_TWF_ITD_SparseData/decoding
+addpath('Z:/ephys/2/Cecilia_TWF/Analysis')
+
+%% plot time series of decoding (based on transients / a sliding time window)
+
+temp=dir('*decoding_30ms.mat');  
+
+no_iter=1000;
+
+all_decoding=[];
+all_decoding_perm=[];
+for i=1:length(temp)
+    load(temp(i).name);
+    try
+        all_decoding(i,:,:)=mean(dist_all,3);
+        for n=1:no_iter
+            all_decoding_perm(i,n,:,:)=mean(dist_all.*reshape(randsample([-1 1],size(dist_all,1)*size(dist_all,2)*size(dist_all,3),'true'),size(dist_all)),3);
+        end
+    catch
+        all_decoding(i,:,1:43)=mean(dist_all,3);
+        for n=1:no_iter
+            all_decoding_perm(i,n,:,1:43)=mean(dist_all.*reshape(randsample([-1 1],size(dist_all,1)*size(dist_all,2)*size(dist_all,3),'true'),size(dist_all)),3);
+        end
+    end
+end
+
+
+confint_99=squeeze(mean(prctile(all_decoding_perm,[.5 99.5],2),1));
+confint_95=squeeze(mean(prctile(all_decoding_perm,[2.5 97.5],2),1));
+
+all_decoding(:,:,[1:3 end-3:end])=NaN;
+confint_99(:,:,[1:3 end-3:end])=NaN;
+confint_95(:,:,[1:3 end-3:end])=NaN;
+
+clear ps_fluc
+for i=1:size(all_decoding,2)
+    for j=4:size(all_decoding,3)-4
+%         [h,ps_fluc(i,j),ci,stats]=ttest(all_decoding(:,i,j));
+        ps_fluc(i,j)=signrank(all_decoding(:,i,j));
+    end
+end
+ps_fluc(:,1:3)=NaN;
+
+% for i=1:12
+%     for j=1:4
+%         for k=1:44
+%             all_decoding_perm_sorted(i,:,j,k)=sort(all_decoding_perm(i,:,j,k));
+%         end
+%     end
+% end
+% 
+% clear ps
+% for i=1:size(all_decoding,2)
+%     for j=4:size(all_decoding,3)-4
+%         ps(i,j)=length(find(mean(all_decoding_perm_sorted(:,:,i,j),1)>mean(all_decoding(:,i,j))))/1000;
+%     end
+% end
+% ps(:,1:3)=NaN;
+
+clear pfdr
+for i=1:4
+    psorted=sort(ps_fluc(i,:)); 
+    psorted=psorted(~isnan(psorted));
+    for j=1:length(psorted)
+        if psorted(j)<j*.005/length(psorted)
+            pfdr(j)=1;
+        else
+            pfdr(j)=0;
+        end
+    end
+    try
+        pfdrthresh(i)=psorted(max(find(pfdr==1)));  
+    catch
+        pfdrthresh(i)=0;
+    end
+end
+
+%% compute average 1-30ms trained
+
+temp=dir('trained*twin_perm.mat');
+% temp=dir('trained*900Hz*decoding_twin.mat');
+all_decoding_avg=[];
+all_decoding_perm=[];
+for i=1:length(temp)
+    load(temp(i).name);
+    
+    all_decoding_avg(i,:)=mean(dist_all,2);
+
+    for n=1:no_iter
+        all_decoding_perm(i,n,:)=mean(dist_all.*reshape(randsample([-1 1],size(dist_all,1)*size(dist_all,2),'true'),size(dist_all)),2);
+    end
+end
+
+confint_99_avg=squeeze(mean(prctile(all_decoding_perm,[.5 99.5],2),1));
+
+clear ps
+for i=1:4
+    [ps(i) h stats]=signrank(all_decoding_avg(:,i),0,'method','approximate');
+    zs(i) = stats.zval;
+end
+
+for i=1:4
+    for j=1:4
+        ps2(i,j)=signrank(all_decoding_avg(:,i),all_decoding_avg(:,j));
+    end
+end
+
+all_decoding_avg_trained = all_decoding_avg;
+ps_trained = ps;
+
+%% compute average 1-30ms naive
+
+temp=dir('TWF*twin_perm.mat');
+% temp=dir('trained*900Hz*decoding_twin.mat');
+all_decoding_avg=[];
+all_decoding_perm=[];
+for i=1:length(temp)
+    load(temp(i).name);
+    
+    all_decoding_avg(i,:)=mean(dist_all,2);
+
+    for n=1:no_iter
+        all_decoding_perm(i,n,:)=mean(dist_all.*reshape(randsample([-1 1],size(dist_all,1)*size(dist_all,2),'true'),size(dist_all)),2);
+    end
+end
+
+confint_99_avg=squeeze(mean(prctile(all_decoding_perm,[.5 99.5],2),1));
+
+clear ps
+for i=1:4
+    [ps(i) h stats]=signrank(all_decoding_avg(:,i),0,'method','approximate');
+    zs(i) = stats.zval;
+end
+
+for i=1:4
+    for j=1:4
+        ps2(i,j)=signrank(all_decoding_avg(:,i),all_decoding_avg(:,j));
+    end
+end
+
+all_decoding_avg_naive = all_decoding_avg;
+ps_naive = ps;
+
+%% Compute 300 hz 
+
+temp=dir('*300Hz*twin_perm.mat');
+all_decoding_300=[];
+all_decoding_perm=[];
+for i=1:length(temp)
+    load(temp(i).name);
+    
+    all_decoding_300(i,:)=mean(dist_all,2);
+
+    for n=1:no_iter
+        all_decoding_perm(i,n,:)=mean(dist_all.*reshape(randsample([-1 1],size(dist_all,1)*size(dist_all,2),'true'),size(dist_all)),2);
+    end
+end
+
+confint_99_300=squeeze(mean(prctile(all_decoding_perm,[.5 99.5],2),1)); 
+
+clear ps300
+for i=1:4
+    [ps300(i) h stats]=signrank(all_decoding_300(:,i),0,'method','approximate');
+    zs(i) = stats.zval;
+end
+
+for i=1:4
+    for j=1:4
+        ps300_2(i,j)=signrank(all_decoding_300(:,i),all_decoding_300(:,j));
+    end
+end
+
+%% Compute 900 Hz 
+
+temp=dir('*900Hz*twin_perm.mat');
+all_decoding_900=[];
+all_decoding_perm=[];
+for i=1:length(temp)
+    load(temp(i).name);
+    
+    all_decoding_900(i,:)=mean(dist_all,2);
+
+    for n=1:no_iter
+        all_decoding_perm(i,n,:)=mean(dist_all.*reshape(randsample([-1 1],size(dist_all,1)*size(dist_all,2),'true'),size(dist_all)),2);
+    end
+end
+
+confint_99_900=squeeze(mean(prctile(all_decoding_perm,[.5 99.5],2),1));
+
+clear ps900
+for i=1:4
+    [ps900(i) h stats]=signrank(all_decoding_900(:,i),0,'method','approximate');
+    zs(i) = stats.zval;
+end
+
+for i=1:4
+    for j=1:4
+        ps900_2(i,j)=signrank(all_decoding_900(:,i),all_decoding_900(:,j));
+    end
+end
+
+%% compute "good" example 
+
+temp=dir('trained_twf_1801_ITD_900Hz_leftCortex_decoding_snr_twin_perm.mat');
+all_decoding_good=[];
+all_decoding_perm=[];
+for i=1:length(temp)
+    load(temp(i).name);
+    
+    all_decoding_good(i,:)=mean(dist_all,2);
+
+    for n=1:no_iter
+        all_decoding_perm(i,n,:)=mean(dist_all.*reshape(randsample([-1 1],size(dist_all,1)*size(dist_all,2),'true'),size(dist_all)),2);
+    end
+end
+
+confint_99_good=squeeze(mean(prctile(all_decoding_perm,[.5 99.5],2),1));
+
+clear ps_good
+for i=1:4
+    [ps_good(i) h stats]=signrank(all_decoding_good(:,i),0,'method','approximate');
+    zs(i) = stats.zval;
+end
+
+for i=1:4
+    for j=1:4
+        ps_good_2(i,j)=signrank(all_decoding_good(:,i),all_decoding_good(:,j));
+    end
+end
+% compute "bad" example 
+
+temp=dir('trained_twf_1803_ITD_900Hz_decoding_snr_twin_perm.mat');
+all_decoding_bad=[];
+all_decoding_perm=[];
+for i=1:length(temp)
+    load(temp(i).name);
+    
+    all_decoding_bad(i,:)=mean(dist_all,2);
+
+    for n=1:no_iter
+        all_decoding_perm(i,n,:)=mean(dist_all.*reshape(randsample([-1 1],size(dist_all,1)*size(dist_all,2),'true'),size(dist_all)),2);
+    end
+end
+
+confint_99_bad=squeeze(mean(prctile(all_decoding_perm,[.5 99.5],2),1));
+
+clear ps_bad
+for i=1:4
+    [ps_bad(i) h stats]=signrank(all_decoding_bad(:,i),0,'method','approximate');
+    zs(i) = stats.zval;
+end
+
+for i=1:4
+    for j=1:4
+        ps_bad_2(i,j)=signrank(all_decoding_bad(:,i),all_decoding_bad(:,j));
+    end
+end
+%% Now plot
+fsz = 14;
+figure('units','centimeters','position',[1,1,15,13]); clf;
+ax2=subplot(3,4,[3 4 7 8])
+set(ax2,'fontsize',fsz)
+fsz=12;%8
+cols={'b' 'r' 'm' 'k'};
+hold on
+for i=1:4
+    shadedErrorBar(taxis_resampled,squeeze(mean(all_decoding(:,i,:),1)),squeeze(std(all_decoding(:,i,:),[],1))/sqrt(14),'lineprops',cols{i});  %all_decoding to plot all click rates
+%     shadedErrorBar(taxis_resampled,squeeze(mean(all_decoding(:,i,:),1)),squeeze(confint_95(2,i,:)),'lineprops',cols{i});  %all_decoding to plot all click rates
+end
+xlim([0 .2])
+% ylabel('decoding (a.u.)','fontsize',fsz)
+xlabel('time (s)','fontsize',fsz)
+t2=title('Decoding activity fluctuations','fontsize',fsz)
+p=get(t2,'pos'); p(2)=p(2)*1.07; set(t2,'pos',p)
+line([0 .2],[0 0],'linestyle','--','color',[.5 .5 .5])
+try
+    scatter(taxis_resampled(find(ps_fluc(1,:)<=pfdrthresh(1))),ones(1,length(find(ps_fluc(1,:)<=pfdrthresh(1))))*-2*10e-5,15,cols{1},'filled')
+catch
+end
+try
+    scatter(taxis_resampled(find(ps_fluc(2,:)<=pfdrthresh(2))),ones(1,length(find(ps_fluc(2,:)<=pfdrthresh(2))))*-4*10e-5,15,cols{2},'filled')
+catch
+end
+try
+    scatter(taxis_resampled(find(ps_fluc(3,:)<=pfdrthresh(3))),ones(1,length(find(ps_fluc(3,:)<=pfdrthresh(3))))*-6*10e-5,15,cols{3},'filled')
+catch
+end
+try
+    scatter(taxis_resampled(find(ps_fluc(4,:)<=pfdrthresh(4))),ones(1,length(find(ps_fluc(4,:)<=pfdrthresh(4))))*-8*10e-5,15,cols{4},'filled')
+catch
+end
+L=legend({'click 1' 'click 2' 'click 3' 'click 4'},'fontsize',fsz-1)
+
+ylim([-8/7 8]*10e-4)
+p=get(ax2,'pos'); p(2)=p(2)*1.1; p(4)=p(4)*.9; set(ax2,'pos',p)
+xl=xlim();xw=xl(2)-xl(1); lx=xl(1)-xw*0.10;
+yl=ylim();yw=yl(2)-yl(1); ly=yl(2)+yw*0.1;
+lab2=text(lx, ly, 'C', 'fontsize',fsz+3);
+%
+ax1=subplot(3,4,[1 5])
+set(ax1,'fontsize',fsz)
+
+cols={'b' 'r' 'm' [.5 .5 .5]};
+hold on
+for i=1:4
+    bar(i,mean(all_decoding_avg_trained(:,i),1),'facecolor',cols{i},'linestyle','none')
+    line([i i],[mean(all_decoding_avg_trained(:,i),1)-std(all_decoding_avg_trained(:,i),1)/sqrt(14) mean(all_decoding_avg_trained(:,i),1)+std(all_decoding_avg_trained(:,i),1)/sqrt(14)],'color','k')
+end
+xlim([0 5])
+ylabel('decoding (a.u.)','fontsize',fsz)
+xlabel('click','fontsize',fsz)
+% xticklabels({' ','1','2','3','4',' '}) 
+xlab = {' ','1','2','3','4',' '};
+set(gca, 'XTick', linspace(0,5,length(xlab)), 'XTickLabels', xlab)
+
+ylim([-2 14]*10e-4)
+scatter(find(ps_trained<.0125),mean(all_decoding_avg_trained(:,find(ps_trained<.0125)),1)*1.36,14,'b*')
+t1=title({'Decoding avg. activity (0-30 ms)' 'trained rats'},'fontsize',fsz)
+p=get(t1,'pos'); p(2)=p(2)*1.07; set(t1,'pos',p)
+p=get(ax1,'pos'); p(2)=p(2)*1.1; p(4)=p(4)*.9; set(ax1,'pos',p)
+xl=xlim();xw=xl(2)-xl(1); lx=xl(1)-xw*0.10;
+yl=ylim();yw=yl(2)-yl(1); ly=yl(2)+yw*0.1;
+lab1=text(lx, ly, 'A', 'fontsize',fsz+3);
+
+%
+ax1=subplot(3,4,[2 6])
+set(ax1,'fontsize',fsz)
+
+cols={'b' 'r' 'm' [.5 .5 .5]};
+hold on
+for i=1:4
+    bar(i,mean(all_decoding_avg_naive(:,i),1),'facecolor',cols{i},'linestyle','none')
+    line([i i],[mean(all_decoding_avg_naive(:,i),1)-std(all_decoding_avg_naive(:,i),1)/sqrt(14) mean(all_decoding_avg_naive(:,i),1)+std(all_decoding_avg_naive(:,i),1)/sqrt(14)],'color','k')
+end
+xlim([0 5])
+ylabel('decoding (a.u.)','fontsize',fsz)
+xlabel('click','fontsize',fsz)
+% xticklabels({' ','1','2','3','4',' '})
+set(gca, 'XTick', linspace(0,5,length(xlab)), 'XTickLabels', xlab)
+ylim([-2 14]*10e-4)
+scatter(find(ps_naive<.0125),mean(all_decoding_avg_naive(:,find(ps_naive<.0125)),1)*1.5,14,'b*')
+t1=title({'Decoding avg. activity (0-30 ms)' 'naive rats'},'fontsize',fsz)
+p=get(t1,'pos'); p(2)=p(2)*1.07; set(t1,'pos',p)
+p=get(ax1,'pos'); p(2)=p(2)*1.1; p(4)=p(4)*.9; set(ax1,'pos',p)
+xl=xlim();xw=xl(2)-xl(1); lx=xl(1)-xw*0.10;
+yl=ylim();yw=yl(2)-yl(1); ly=yl(2)+yw*0.1;
+lab1=text(lx, ly, 'B', 'fontsize',fsz+3);
+
+% plot 300 Hz
+ax3=subplot(3,4,9)
+set(ax3,'fontsize',fsz)
+cols={'b' 'r' 'm' [.5 .5 .5]};
+hold on
+for i=1:4
+    bar(i,mean(all_decoding_300(:,i),1),'facecolor',cols{i},'linestyle','none')
+    line([i i],[mean(all_decoding_300(:,i),1)-std(all_decoding_300(:,i),1)/sqrt(14) mean(all_decoding_300(:,i),1)+std(all_decoding_300(:,i),1)/sqrt(14)],'color','k')
+end
+xlim([0 5])
+ylabel('decoding (a.u.)','fontsize',fsz, 'fontweight','normal')
+xlabel('click','fontsize',fsz, 'fontweight','normal')
+xticks(1:4)
+xticklabels({'1','2','3','4'})
+ylim([-2 14]*10e-4)
+scatter(find(ps300<.05),mean(all_decoding_300(:,find(ps300<.05)),1)*1.5,14,'b*')
+t3=title('300 Hz','fontsize',fsz, 'fontweight','bold');
+p=get(t3,'pos'); p(2)=p(2)*1.15; set(t3,'pos',p)
+xl=xlim();xw=xl(2)-xl(1); lx=xl(1)-xw*0.15;
+yl=ylim();yw=yl(2)-yl(1); ly=yl(2)+yw*0.15;
+lab3=text(lx, ly, 'D', 'fontsize',fsz+3);
+% PLOT 900 HZ
+
+ax4=subplot(3,4,10)
+set(ax4,'fontsize',fsz)
+
+cols={'b' 'r' 'm' [.5 .5 .5]};
+hold on
+for i=1:4
+    bar(i,mean(all_decoding_900(:,i),1),'facecolor',cols{i},'linestyle','none')
+    line([i i],[mean(all_decoding_900(:,i),1)-std(all_decoding_900(:,i),1)/sqrt(14) mean(all_decoding_900(:,i),1)+std(all_decoding_900(:,i),1)/sqrt(14)],'color','k')
+end
+xlim([0 5])
+%ylabel('decoding')
+xlabel('click','fontsize',fsz, 'fontweight','normal')
+xticks(1:4)
+xticklabels({'1','2','3','4'})
+ylim([-2 14]*10e-4)
+scatter(find(ps900<.05),mean(all_decoding_900(:,find(ps900<.05)),1)*1.45,14,'b*')
+t4=title('900 Hz','fontsize',fsz, 'fontweight','bold');
+p=get(t4,'pos'); p(2)=p(2)*1.15; set(t4,'pos',p)
+xl=xlim();xw=xl(2)-xl(1); lx=xl(1)-xw*0.15;
+yl=ylim();yw=yl(2)-yl(1); ly=yl(2)+yw*0.15;
+lab3=text(lx, ly, 'E', 'fontsize',fsz+3);
+
+% plot "good" example
+ax5=subplot(3,4,11)
+set(ax5,'fontsize',fsz)
+
+cols={'b' 'r' 'm' [.5 .5 .5]};
+hold on
+for i=1:4
+    bar(i,all_decoding_good(i),'facecolor',cols{i},'linestyle','none')
+%     line([i i],[mean(all_decoding_good(:,i),1)-std(all_decoding_good(:,i),1)/sqrt(14) mean(all_decoding_good(:,i),1)+std(all_decoding_good(:,i),1)/sqrt(14)],'color','k')
+    line([i i],[all_decoding_good(i)+confint_99_good(1,i) all_decoding_good(i)+confint_99_good(2,i)],'color','k')
+end
+xlim([0 5])
+%ylabel('decoding')
+xlabel('click','fontsize',fsz, 'fontweight','normal')
+xticks(1:4)
+xticklabels({'1','2','3','4'})
+ylim([-7 49]*10e-4)
+% scatter(find(ps_good<.05),mean(all_decoding_good(:,find(ps<.05)),1)*1.4,10,'k*')
+t5=title('Example 1','fontsize',fsz, 'fontweight','bold');
+p=get(t5,'pos'); p(2)=p(2)*1.15; set(t5,'pos',p)
+xl=xlim();xw=xl(2)-xl(1); lx=xl(1)-xw*0.15;
+yl=ylim();yw=yl(2)-yl(1); ly=yl(2)+yw*0.15;
+lab5=text(lx, ly, 'F', 'fontsize',fsz+3);
+% plot "bad" example 
+
+
+ax6=subplot(3,4,12)
+set(ax6,'fontsize',fsz)
+cols={'b' 'r' 'm' [.5 .5 .5]};
+hold on
+for i=1:4
+    bar(i,all_decoding_bad(i),'facecolor',cols{i},'linestyle','none')
+%     line([i i],[mean(all_decoding_bad(:,i),1)-std(all_decoding_bad(:,i),1)/sqrt(14) mean(all_decoding_bad(:,i),1)+std(all_decoding_bad(:,i),1)/sqrt(14)],'color','k')
+    line([i i],[all_decoding_bad(i)+confint_99_bad(1,i) all_decoding_bad(i)+confint_99_bad(2,i)],'color','k')
+end
+xlim([0 5])
+%ylabel('decoding')
+xlabel('click','fontsize',fsz, 'fontweight','normal')
+xticks(1:4)
+xticklabels({'1','2','3','4'})
+ylim([-7 49]*10e-4)
+t6=title('Example 2','fontsize',fsz, 'fontweight','bold')
+% scatter(find(ps_bad<.05),mean(all_decoding_bad(:,find(ps_bad<.05)),1)*1.4,10,'k*')
+p=get(t6,'pos'); p(2)=p(2)*1.15; set(t6,'pos',p)
+xl=xlim();xw=xl(2)-xl(1); lx=xl(1)-xw*0.15;
+yl=ylim();yw=yl(2)-yl(1); ly=yl(2)+yw*0.15;
+lab6=text(lx, ly, 'G', 'fontsize',fsz+3);
+%%
+set(gcf,'paperpositionmode','auto')
+% print -dsvg '/home/jan/jan/document/science/cityu/TWF_ITD_ILD_NormalHearingRats/figures/ecog_decoding_figure_allpanels_trained_naive.svg'
+saveas(gcf, 'C:/Users/kongyanli2/OneDrive - City University of Hong Kong/BMS/TWF/BMC_reviewed/ecog_decoding_figure_allpanels_trained_naive','fig')
+print -dsvg 'C:/Users/kongyanli2/OneDrive - City University of Hong Kong/BMS/TWF/BMC_reviewed/ecog_decoding_figure_allpanels_trained_naive.svg'

ECoG/code/ecog_decoding_plot_left_right.m

@@ -0,0 +1,124 @@
+root='/ephys/2/Cecilia_TWF/Analysis'
+if ~exist(root,'dir')
+    root=['/home/colliculus/' root]
+end
+cd([root '/Analysis_AC_ECoG_TWF_ITD_SparseData/decoding'])
+addpath(root)
+fsz = 10; % font size
+
+% working on Cecilia's desktop
+% cd Z:/ephys/2/Cecilia_TWF/Analysis/Analysis_AC_ECoG_TWF_ITD_SparseData/decoding
+% addpath('Z:/ephys/2/Cecilia_TWF/Analysis')
+
+%% compute average 1-30ms  for left cortex 
+
+temp=dir('*leftCortex_decoding_snr_twin_perm.mat');
+all_decoding_avg=[];
+all_decoding_perm=[];
+for i=1:length(temp)
+    load(temp(i).name);
+    
+    all_decoding_avg(i,:)=mean(dist_all,2);
+
+    for n=1:no_iter
+        all_decoding_perm(i,n,:)=mean(dist_all.*reshape(randsample([-1 1],size(dist_all,1)*size(dist_all,2),'true'),size(dist_all)),2);
+    end
+end
+
+confint_99_avg=squeeze(mean(prctile(all_decoding_perm,[.5 99.5],2),1));
+
+clear ps
+for i=1:4
+    [ps(i) h stats]=signrank(all_decoding_avg(:,i),0,'method','approximate');
+    zs(i) = stats.zval;
+end
+
+for i=1:4
+    for j=1:4
+        ps2(i,j)=signrank(all_decoding_avg(:,i),all_decoding_avg(:,j));
+    end
+end
+
+all_decoding_avg_left = all_decoding_avg;
+ps_left = ps;
+
+
+%% compute average 1-30ms  for right cortex
+
+temp=dir('*Hz_decoding_snr_twin_perm.mat');
+all_decoding_avg=[];
+all_decoding_perm=[];
+for i=1:length(temp)
+    load(temp(i).name);
+    
+    all_decoding_avg(i,:)=mean(dist_all,2);
+
+    for n=1:no_iter
+        all_decoding_perm(i,n,:)=mean(dist_all.*reshape(randsample([-1 1],size(dist_all,1)*size(dist_all,2),'true'),size(dist_all)),2);
+    end
+end
+
+confint_99_avg=squeeze(mean(prctile(all_decoding_perm,[.5 99.5],2),1));
+
+clear ps
+for i=1:4
+    [ps(i) h stats]=signrank(all_decoding_avg(:,i),0,'method','approximate');
+    zs(i) = stats.zval;
+end
+
+for i=1:4
+    for j=1:4
+        ps2(i,j)=signrank(all_decoding_avg(:,i),all_decoding_avg(:,j));
+    end
+end
+
+all_decoding_avg_right = all_decoding_avg;
+ps_right = ps;
+
+%% Now plot
+figure('units','centimeters','position',[1,1,15,13]); clf;
+
+%
+ax1=subplot(1,2,1)
+set(ax1,'fontsize',fsz)
+
+cols={'b' 'r' 'm' [.5 .5 .5]};
+hold on
+for i=1:4
+    bar(i,mean(all_decoding_avg_left(:,i),1),'facecolor',cols{i},'linestyle','none')
+    line([i i],[mean(all_decoding_avg_left(:,i),1)-std(all_decoding_avg_left(:,i),1)/sqrt(14) mean(all_decoding_avg_left(:,i),1)+std(all_decoding_avg_left(:,i),1)/sqrt(14)],'color','k')
+end
+xlim([0 5])
+ylabel('decoding (a.u.)','fontsize',fsz)
+xlabel('click','fontsize',fsz)
+xticklabels({' ','1','2','3','4',' '})
+ylim([-2 14]*10e-4)
+scatter(find(ps_left<.05),mean(all_decoding_avg_left(:,find(ps_left<.05)),1)*1.4,10,'k*')
+t1=title('left cortex','fontsize',fsz)
+p=get(t1,'pos'); p(2)=p(2)*1.07; set(t1,'pos',p)
+p=get(ax1,'pos'); p(2)=p(2)*1.1; p(4)=p(4)*.9; set(ax1,'pos',p)
+xl=xlim();xw=xl(2)-xl(1); lx=xl(1)-xw*0.10;
+yl=ylim();yw=yl(2)-yl(1); ly=yl(2)+yw*0.1;
+lab1=text(lx, ly, 'A', 'fontsize',fsz+3);
+
+ax1=subplot(1,2,2)
+set(ax1,'fontsize',fsz)
+
+cols={'b' 'r' 'm' [.5 .5 .5]};
+hold on
+for i=1:4
+    bar(i,mean(all_decoding_avg_right(:,i),1),'facecolor',cols{i},'linestyle','none')
+    line([i i],[mean(all_decoding_avg_right(:,i),1)-std(all_decoding_avg_right(:,i),1)/sqrt(14) mean(all_decoding_avg_right(:,i),1)+std(all_decoding_avg_right(:,i),1)/sqrt(14)],'color','k')
+end
+xlim([0 5])
+ylabel('decoding (a.u.)','fontsize',fsz)
+xlabel('click','fontsize',fsz)
+xticklabels({' ','1','2','3','4',' '})
+ylim([-2 14]*10e-4)
+scatter(find(ps_right<.05),mean(all_decoding_avg_right(:,find(ps_right<.05)),1)*1.4,10,'k*')
+t1=title('right cortex','fontsize',fsz)
+p=get(t1,'pos'); p(2)=p(2)*1.07; set(t1,'pos',p)
+p=get(ax1,'pos'); p(2)=p(2)*1.1; p(4)=p(4)*.9; set(ax1,'pos',p)
+xl=xlim();xw=xl(2)-xl(1); lx=xl(1)-xw*0.10;
+yl=ylim();yw=yl(2)-yl(1); ly=yl(2)+yw*0.1;
+lab1=text(lx, ly, 'A', 'fontsize',fsz+3);

ECoG/code/ecog_decoding_plot_trained_naive.m

@@ -0,0 +1,124 @@
+root='/ephys/2/Cecilia_TWF/Analysis'
+if ~exist(root,'dir')
+    root=['/home/colliculus/' root]
+end
+cd([root '/Analysis_AC_ECoG_TWF_ITD_SparseData/decoding'])
+addpath(root)
+fsz = 10; % font size
+
+% working on Cecilia's desktop
+% cd Z:/ephys/2/Cecilia_TWF/Analysis/Analysis_AC_ECoG_TWF_ITD_SparseData/decoding
+% addpath('Z:/ephys/2/Cecilia_TWF/Analysis')
+
+%% compute average 1-30ms  for naive rats
+
+temp=dir('TWF*twin_perm.mat');
+all_decoding_avg=[];
+all_decoding_perm=[];
+for i=1:length(temp)
+    load(temp(i).name);
+    
+    all_decoding_avg(i,:)=mean(dist_all,2);
+
+    for n=1:no_iter
+        all_decoding_perm(i,n,:)=mean(dist_all.*reshape(randsample([-1 1],size(dist_all,1)*size(dist_all,2),'true'),size(dist_all)),2);
+    end
+end
+
+confint_99_avg=squeeze(mean(prctile(all_decoding_perm,[.5 99.5],2),1));
+
+clear ps
+for i=1:4
+    [ps(i) h stats]=signrank(all_decoding_avg(:,i),0,'method','approximate');
+    zs(i) = stats.zval;
+end
+
+for i=1:4
+    for j=1:4
+        ps2(i,j)=signrank(all_decoding_avg(:,i),all_decoding_avg(:,j));
+    end
+end
+
+all_decoding_avg_naive = all_decoding_avg;
+ps_naive = ps;
+
+
+%% compute average 1-30ms  for naive rats
+
+temp=dir('trained*twin_perm.mat');
+all_decoding_avg=[];
+all_decoding_perm=[];
+for i=1:length(temp)
+    load(temp(i).name);
+    
+    all_decoding_avg(i,:)=mean(dist_all,2);
+
+    for n=1:no_iter
+        all_decoding_perm(i,n,:)=mean(dist_all.*reshape(randsample([-1 1],size(dist_all,1)*size(dist_all,2),'true'),size(dist_all)),2);
+    end
+end
+
+confint_99_avg=squeeze(mean(prctile(all_decoding_perm,[.5 99.5],2),1));
+
+clear ps
+for i=1:4
+    [ps(i) h stats]=signrank(all_decoding_avg(:,i),0,'method','approximate');
+    zs(i) = stats.zval;
+end
+
+for i=1:4
+    for j=1:4
+        ps2(i,j)=signrank(all_decoding_avg(:,i),all_decoding_avg(:,j));
+    end
+end
+
+all_decoding_avg_trained = all_decoding_avg;
+ps_trained = ps;
+
+%% Now plot
+figure('units','centimeters','position',[1,1,15,13]); clf;
+
+%
+ax1=subplot(1,2,1)
+set(ax1,'fontsize',fsz)
+
+cols={'b' 'r' 'm' [.5 .5 .5]};
+hold on
+for i=1:4
+    bar(i,mean(all_decoding_avg_naive(:,i),1),'facecolor',cols{i},'linestyle','none')
+    line([i i],[mean(all_decoding_avg_naive(:,i),1)-std(all_decoding_avg_naive(:,i),1)/sqrt(14) mean(all_decoding_avg_naive(:,i),1)+std(all_decoding_avg_naive(:,i),1)/sqrt(14)],'color','k')
+end
+xlim([0 5])
+ylabel('decoding (a.u.)','fontsize',fsz)
+xlabel('click','fontsize',fsz)
+xticklabels({' ','1','2','3','4',' '})
+ylim([-2 14]*10e-4)
+scatter(find(ps_naive<.05),mean(all_decoding_avg_naive(:,find(ps_naive<.05)),1)*1.4,10,'k*')
+t1=title('naive rats','fontsize',fsz)
+p=get(t1,'pos'); p(2)=p(2)*1.07; set(t1,'pos',p)
+p=get(ax1,'pos'); p(2)=p(2)*1.1; p(4)=p(4)*.9; set(ax1,'pos',p)
+xl=xlim();xw=xl(2)-xl(1); lx=xl(1)-xw*0.10;
+yl=ylim();yw=yl(2)-yl(1); ly=yl(2)+yw*0.1;
+lab1=text(lx, ly, 'A', 'fontsize',fsz+3);
+
+ax1=subplot(1,2,2)
+set(ax1,'fontsize',fsz)
+
+cols={'b' 'r' 'm' [.5 .5 .5]};
+hold on
+for i=1:4
+    bar(i,mean(all_decoding_avg_trained(:,i),1),'facecolor',cols{i},'linestyle','none')
+    line([i i],[mean(all_decoding_avg_trained(:,i),1)-std(all_decoding_avg_trained(:,i),1)/sqrt(14) mean(all_decoding_avg_trained(:,i),1)+std(all_decoding_avg_trained(:,i),1)/sqrt(14)],'color','k')
+end
+xlim([0 5])
+ylabel('decoding (a.u.)','fontsize',fsz)
+xlabel('click','fontsize',fsz)
+xticklabels({' ','1','2','3','4',' '})
+ylim([-2 14]*10e-4)
+scatter(find(ps_trained<.05),mean(all_decoding_avg_trained(:,find(ps_trained<.05)),1)*1.4,10,'k*')
+t1=title('trained rats','fontsize',fsz)
+p=get(t1,'pos'); p(2)=p(2)*1.07; set(t1,'pos',p)
+p=get(ax1,'pos'); p(2)=p(2)*1.1; p(4)=p(4)*.9; set(ax1,'pos',p)
+xl=xlim();xw=xl(2)-xl(1); lx=xl(1)-xw*0.10;
+yl=ylim();yw=yl(2)-yl(1); ly=yl(2)+yw*0.1;
+lab1=text(lx, ly, 'A', 'fontsize',fsz+3);