This pipeline opens the result of ICAalamelodic.m, lets the user interactively label the components that look like neuronal activity (rather than movement artefacts or noise), sort them by label, plots a final summary for the chosen components, and save the reordered maps and time series.



In [1]:

    
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy import io
%matplotlib inline 
import pylab

Open time series



In [2]:

    
import scipy.io as sio



In [3]:

    
# from http://stackoverflow.com/questions/3579568/choosing-a-file-in-python-with-simple-dialog
from Tkinter import Tk
from tkFileDialog import askopenfilename

Tk().withdraw() # we don't want a full GUI, so keep the root window from appearing
filename = askopenfilename() # show an "Open" dialog box and return the path to the selected file
print(filename)









    



/media/test/100431/100431ss1regdFF20spsfkf90Smith0_4_60TS.mat



In [4]:

    
Ua=sio.loadmat(filename)



In [9]:

    
DT=Ua['TSmean']



In [10]:

    
DT.shape









    Out[10]:





(1095, 90)



In [11]:

    
S1=DT.shape



In [12]:

    
DTmean=np.zeros(S1)
DTvar=np.zeros(S1)
Var=np.zeros(S1[1])



In [13]:

    
for i in range(S1[1]):
    DTmean[:,i]=DT[:,i]-np.mean(DT[:,i],0)



In [14]:

    
for i in range(S1[1]):
    Var[i]=np.sqrt(np.var(DTmean[:,i]))
    DTvar[:,i]=DTmean[:,i]/Var[i]



In [15]:

    
DTvar.shape









    Out[15]:





(1095, 90)

open maps



In [16]:

    
import nibabel as nb



In [17]:

    
# from http://stackoverflow.com/questions/3579568/choosing-a-file-in-python-with-simple-dialog
from Tkinter import Tk
from tkFileDialog import askopenfilename

Tk().withdraw() # we don't want a full GUI, so keep the root window from appearing
filename2 = askopenfilename() # show an "Open" dialog box and return the path to the selected file
print(filename2)









    



/media/test/100431/100431ss1regdFF20spsfkf90Smith0_4_60IC.nii



In [18]:

    
img1 = nb.load(filename2)



In [19]:

    
data = img1.get_data()



In [20]:

    
S=data.shape



In [21]:

    
S









    Out[21]:





(106, 106, 8, 90)

Zscore maps



In [22]:

    
Demean=np.zeros(S)
Dmaps=np.zeros(S)
Dvar=np.zeros(S)
Var=np.zeros(S[3])
D2=np.zeros([S[0],S[1],5,S[3]])
Tvar=np.zeros(S[3])

Transform the maps to have zero mean



In [23]:

    
for i in range(S[3]):
    Demean[:,:,:,i]=data[:,:,:,i]-np.mean(np.mean(np.mean(data[:,:,:,i],0),0),0)

Transform the maps to have unit variance and zscore



In [24]:

    
for i in range(S[3]):
    Dsq=np.reshape(Demean[:,:,:,i],S[0]*S[1]*S[2])
    Var[i]=np.sqrt(np.var(Dsq))
    Dvar=Demean[:,:,:,i]/Var[i]
    Dmaps[:,:,:,i]=Dvar-2.5
Dmaps[Dmaps<0]=0

Order ICs by variance



In [25]:

    
datao=data
Dmapso=Dmaps



In [26]:

    
plt.plot(Var)









    Out[26]:





[<matplotlib.lines.Line2D at 0x7f8c4aa9d050>]

Separate maps in substacks, sort the independent components by brain regions



In [27]:

    
my_cmap=plt.cm.jet
my_cmap.set_bad(alpha=0)
Good_ICs=np.zeros(S[3])
Label_ICs=[]
pylab.rcParams['figure.figsize'] = (13, 2.5)



In [28]:

    
Dtemp=data[:,:,:,0]



In [29]:

    
%%javascript
IPython.OutputArea.auto_scroll_threshold =4000;



In [30]:

    
if S[2]>5:
    Nstack=5
    Int100=[(i+1)*100/Nstack for i in range(Nstack)]
    Percs=np.percentile(range(S[2]),Int100)
    Indices=np.split(range(S[2]),Percs)
    D1=np.zeros([S[0],S[1],Nstack])
    Dmean=Dtemp[:,:,range(Nstack)]
    for i in range(Nstack):
        Vmean=np.mean(Dtemp[:,:,Indices[i]],2)
        #Dmean[:,:,i]=np.max(Vmean,0)
        Dmean[:,:,i]=Vmean
else:
    Nstack=S[2]
    D1=np.zeros([S[0],S[1],S[2]])
    Dmean=data[:,:,range(S[2])]  
    Dmean=np.squeeze(Dtemp[:,:,:])









    



/usr/local/lib/python2.7/dist-packages/numpy/lib/shape_base.py:422: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
  sub_arys.append(_nx.swapaxes(sary[st:end], axis, 0))



In [31]:

    
DTvar.shape









    Out[31]:





(1095, 90)



In [32]:

    
S









    Out[32]:





(106, 106, 8, 90)



In [33]:

    
# from http://stackoverflow.com/questions/3579568/choosing-a-file-in-python-with-simple-dialog
Tk().withdraw() # we don't want a full GUI, so keep the root window from appearing
filename = askopenfilename() # show an "Open" dialog box and return the path to the selected file
print(filename)









    



/media/test/100431/MAX_100431ss1regdFF20spsfkf90Smith0_4_60IC.nii



In [34]:

    
Ua









    Out[34]:





{'TSmean': array([[  1.24743107e-01,   3.54995452e-02,  -2.54939253e-02, ...,
           7.73417170e-03,   1.99389138e-02,  -4.58816337e-03],
        [  1.74235784e-01,   1.54234967e-02,  -1.92229023e-02, ...,
          -2.51549515e-03,  -5.96453404e-04,   2.78127907e-03],
        [  2.04678666e-01,   1.09864731e-02,  -1.59863112e-02, ...,
          -4.96297832e-03,  -1.05652704e-02,   3.22859873e-03],
        ..., 
        [  2.20148348e-01,   2.20538451e-02,  -2.06020309e-02, ...,
          -8.43721813e-05,  -1.13540602e-03,   9.81601592e-05],
        [  2.63678384e-01,   2.69465535e-02,  -1.88091289e-02, ...,
          -3.75771262e-04,  -6.24219786e-04,   3.02683555e-04],
        [  2.41405248e-01,   2.50890927e-02,  -2.00133516e-02, ...,
          -2.64972104e-04,  -6.23886708e-04,   2.22896138e-04]]),
 '__globals__': [],
 '__header__': 'MATLAB 5.0 MAT-file, Platform: GLNXA64, Created on: Thu Apr  4 12:28:06 2019',
 '__version__': '1.0'}



In [35]:

    
Ua=sio.loadmat(filename)
Xk[0,:]=Ua['Groom']
#Xk[1,:]=Ua['Walk']









    



---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-35-a1a915206512> in <module>()
----> 1 Ua=sio.loadmat(filename)
      2 Xk[0,:]=Ua['Groom']
      3 #Xk[1,:]=Ua['Walk']

/usr/local/lib/python2.7/dist-packages/scipy/io/matlab/mio.pyc in loadmat(file_name, mdict, appendmat, **kwargs)
    134     variable_names = kwargs.pop('variable_names', None)
    135     MR = mat_reader_factory(file_name, appendmat, **kwargs)
--> 136     matfile_dict = MR.get_variables(variable_names)
    137     if mdict is not None:
    138         mdict.update(matfile_dict)

/usr/local/lib/python2.7/dist-packages/scipy/io/matlab/mio4.pyc in get_variables(self, variable_names)
    392         mdict = {}
    393         while not self.end_of_stream():
--> 394             hdr, next_position = self.read_var_header()
    395             name = asstr(hdr.name)
    396             if variable_names is not None and name not in variable_names:

/usr/local/lib/python2.7/dist-packages/scipy/io/matlab/mio4.pyc in read_var_header(self)
    348            position in stream of next variable
    349         '''
--> 350         hdr = self._matrix_reader.read_header()
    351         n = reduce(lambda x, y: x*y, hdr.dims, 1)  # fast product
    352         remaining_bytes = hdr.dtype.itemsize * n

/usr/local/lib/python2.7/dist-packages/scipy/io/matlab/mio4.pyc in read_header(self)
    119         O, rest = divmod(rest, 100)  # unused, should be 0
    120         if O != 0:
--> 121             raise ValueError('O in MOPT integer should be 0, wrong format?')
    122         P, rest = divmod(rest, 10)  # data type code e.g miDOUBLE (see above)
    123         T = rest  # matrix type code e.g. mxFULL_CLASS (see above)

ValueError: O in MOPT integer should be 0, wrong format?



In [37]:

    
# from http://stackoverflow.com/questions/3579568/choosing-a-file-in-python-with-simple-dialog
from Tkinter import Tk
from tkFileDialog import askopenfilename

Tk().withdraw() # we don't want a full GUI, so keep the root window from appearing
filenamet = askopenfilename() # show an "Open" dialog box and return the path to the selected file
print(filenamet)
nimt=nb.load(filenamet)
Dtemp=np.squeeze(nimt.get_data())
Dtemp.shape


if S[2]>5:
    Nstack=5
    Int100=[(i+1)*100/Nstack for i in range(Nstack)]
    Percs=np.percentile(range(S[2]),Int100)
    Indices=np.split(range(S[2]),Percs)
    Dmean=np.zeros([S[0],S[1],Nstack])
    #Dmean=np.squeeze(data[:,:,range(Nstack),2])
    for i in range(Nstack):
        Vmean=np.mean(Dtemp[:,:,Indices[i]],2)
        Dmean[:,:,i]=Vmean

plt.imshow(Vmean,cmap=plt.cm.gray)









    



/media/test/100431/MAX_100431ss1regdFF20spsfkf90Smith0_4_60IC.nii






    Out[37]:





<matplotlib.image.AxesImage at 0x7f8c4a90d590>



In [38]:

    
for j in range(S[3]):

    if S[2]>5:
        for i in range(Nstack):
            V=Dmaps[:,:,Indices[i],j]
            D1[:,:,i]=np.max(V,2)
        D2[:,:,:,j]=D1
        D1[D1==0]=np.nan
           
    else:
        for i in range(S[2]):
            V=Dmaps[:,:,i,j]
            D1[:,:,i]=V 
            

    print(j)
    for i in range(Nstack):
        plt.subplot(1,5,i+1)
        plt.imshow(Dmean[:,:,i],cmap=plt.cm.gray)
        plt.imshow(D1[:,:,i], cmap=my_cmap,interpolation='none')
        frame1 = plt.gca()
        frame1.axes.get_xaxis().set_visible(False)
        frame1.axes.get_yaxis().set_visible(False)
        
    plt.show()
    
   # plt.plot(TS_ROI[Order[j],:])
    plt.plot(DTvar[:,j])
    #plt.plot(Xk[0,:]/np.std(Xk[0,:])+0.5,color=(1,0,0))   
    #plt.plot(Xk[1,:]/np.std(Xk[1,:])+0.5,color=(0,1,0))
    #plt.plot(Xk[2,:]/np.std(Xk[1,:])+0.5,color=(0.5,0.5,0))    
    #plt.plot(Xk[3,:]/np.std(Xk[1,:])+0.5,color=(0,0.5,1))
    
    plt.show()
    
    Label_ICs.append(raw_input())
    if Label_ICs[j]=='':
        Good_ICs[j]=0
    else:
        Good_ICs[j]=1



In [39]:

    
List1=[(Label_ICs[i],i) for i in range(S[3])]
Newlist=sorted(List1, key=lambda List1: List1[0])

Neworder=[Newlist[i][1] for i in range(S[3])]

NewDT=DTvar[:,Neworder[:]].T

for j in range(len(Neworder)):
    A=NewDT[:,j]
    V=np.sqrt(np.var(A))
    NewDT[:,j]=A/V

C1=np.zeros([6,3])
C1[0][:]=(1,0,0)
C1[1][:]=(0,1,0)
C1[2][:]=(0,0,1)
C1[3][:]=(0.8,0.8,0)
C1[4][:]=(0,1,1)
C1[5][:]=(1,0,1)
h=3

Newmaps=Dmaps[:,:,:,Neworder[:]]

L=len(set([Label_ICs[Neworder[i]] for i in range(len(Neworder))]))

Regionmaps=np.zeros([S[0],S[1],L,3])
Datasort=np.zeros([S[0],S[1],S[2],L,3])

Regionname=[]

DMapsordered=Dmapso[:,:,:,Neworder[:]]

DMapsordered=Dmapso[:,:,:,Neworder[:]]

j=0
i=0
k=Label_ICs[Neworder[0]]
m=0
Regionname.append(Label_ICs[Neworder[i]])
for i in range(len(Neworder)):
    #C2=C1[i%6][:]
    for l in range(3):
        M=np.max(np.squeeze(np.reshape(Newmaps[:,:,:,i],S[0]*S[1]*S[2])))
        Regionmaps[:,:,j,l]=Regionmaps[:,:,j,l]+0.7*np.max(DMapsordered[:,:,:,i],2)*C1[i%6][l]/M
        Datasort[:,:,:,j,l]=Datasort[:,:,:,j,l]+Dmaps[:,:,:,Neworder[i]]*C1[i%6][l] 
    i=i+1
    m=m+1
    if i<len(Neworder):
        k1=Label_ICs[Neworder[i]]
        
    if k1 != k:
        j=j+1
        k=k1
        m=0
        Regionname.append(Label_ICs[Neworder[i]])

pylab.rcParams['figure.figsize'] = (14, 5)
import scipy
from scipy import ndimage
j=0
m=0
L=0
k=Label_ICs[Neworder[0]]
for i in range(len(Neworder)):
    m=m+1
    
    
    if i<len(Neworder):
        k1=Label_ICs[Neworder[i]]
        
    if k1 != k:
        
        k=k1
        m=0
        
        plt.show()
        plt.figure(2*j+1)
        Rotated_Plot = ndimage.rotate(Regionmaps[:,:,j], -90)
        IM=plt.imshow(Rotated_Plot) 
        frame1 = plt.gca()
        frame1.axes.get_xaxis().set_visible(False)
        frame1.axes.get_yaxis().set_visible(False)
        j=j+1
        plt.figure(2*j)
        #plt.plot(Xk[0,:]/np.std(Xk[0,:])+0.5,color=(1,0,0))   
        #plt.plot(Xk[1,:]/np.std(Xk[1,:])+0.5,color=(0,1,0))
        #plt.plot(Xk[2,:]/np.std(Xk[1,:])+0.5,color=(0.5,0.5,0))    
        #plt.plot(Xk[3,:]/np.std(Xk[1,:])+0.5,color=(0,0.5,1))
    plt.plot(NewDT[i,:]+h*m,color=C1[i%6][:])
    
plt.figure(2*j+1)
Rotated_Plot = ndimage.rotate(Regionmaps[:,:,j], -90)
IM=plt.imshow(Rotated_Plot) 
frame1 = plt.gca()
frame1.axes.get_xaxis().set_visible(False)
frame1.axes.get_yaxis().set_visible(False)



In [ ]: