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]:

    
clear all



In [2]:

    
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 [3]:

    
import scipy.io as sio



In [4]:

    
# 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)









    



/home/sophie/Downloads/792ss1onc500regcU7sMpsfkf100Smith2_3TS.mat



In [5]:

    
Ua=sio.loadmat(filename)



In [6]:

    
Ua









    Out[6]:





{'TSo': array([[ -6.60401090e-02,   1.32515761e-01,  -9.14888277e-02, ...,
           1.40896898e-01,   1.22901536e-03,   1.66143481e-02],
        [ -1.67052622e-01,  -5.24866220e-02,  -6.81334714e-02, ...,
           4.02301106e-02,  -1.09442161e-01,   3.38636518e-02],
        [  7.99838190e-02,  -8.23296734e-03,   4.50064856e-02, ...,
          -5.36471952e-02,  -8.78289940e-03,   1.13714355e-02],
        ..., 
        [  3.36594174e-02,   7.28882965e-03,  -1.58677024e-02, ...,
           1.24172442e-02,  -2.52255829e-02,  -5.59877411e-03],
        [ -2.43206711e-02,   1.81880143e-02,  -3.10555841e-02, ...,
           2.04137304e-02,  -2.33141023e-02,   6.43387618e-05],
        [ -1.18028934e-02,   8.54281643e-03,  -1.42022739e-02, ...,
           1.03086964e-02,  -1.17417295e-02,  -1.44086452e-03]]),
 '__globals__': [],
 '__header__': 'MATLAB 5.0 MAT-file, Platform: GLNXA64, Created on: Fri Jul 15 17:03:57 2016',
 '__version__': '1.0'}



In [7]:

    
DT=Ua['u']









    



---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
<ipython-input-7-954f5e7f0b22> in <module>()
----> 1 DT=Ua['u']

KeyError: 'u'



In [8]:

    
DT=Ua['TSo']



In [9]:

    
DT.shape









    Out[9]:





(5616, 100)



In [10]:

    
S1=DT.shape



In [11]:

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



In [12]:

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



In [13]:

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



In [14]:

    
DTvar.shape









    Out[14]:





(5616, 100)

open maps



In [15]:

    
from nifti import NiftiImage



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)









    



/home/sophie/Downloads/792ss1onc500regcU7sMpsfkf100Smith2_3IC.nii



In [18]:

    
img1 = nb.load(filename2)



In [19]:

    
data = img1.get_data()



In [20]:

    
S=data.shape



In [21]:

    
S









    Out[21]:





(99, 47, 11, 100)

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 0x48d1650>]

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]:

    
Tk().withdraw() 
filenamet = askopenfilename() 
print(filenamet)
nimt=NiftiImage(filenamet)
Dtemp=np.squeeze(nimt.data.T)
Dtemp.shape









    



()






    



---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-28-d6de7db42465> in <module>()
      2 filenamet = askopenfilename()
      3 print(filenamet)
----> 4 nimt=NiftiImage(filenamet)
      5 Dtemp=np.squeeze(nimt.data.T)
      6 Dtemp.shape

/usr/lib/pymodules/python2.7/nifti/image.pyc in __init__(self, source, header, load, **kwargs)
     76         """
     77         # setup all nifti header related stuff
---> 78         NiftiFormat.__init__(self, source, header, **kwargs)
     79 
     80         # where the data will go to

/usr/lib/pymodules/python2.7/nifti/format.pyc in __init__(self, source, header, loadmeta)
     83             raise ValueError, \
     84                   "Unsupported source type. Only NumPy arrays and filename " \
---> 85                   + "string are supported."
     86 
     87 

ValueError: Unsupported source type. Only NumPy arrays and filename string are supported.



In [29]:

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



In [30]:

    
%%javascript
IPython.OutputArea.auto_scroll_threshold =4000;



In [31]:

    
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[:,:,:])



In [32]:

    
DTvar.shape









    Out[32]:





(5616, 100)



In [33]:

    
S









    Out[33]:





(99, 47, 11, 100)



In [34]:

    
Time=np.array(range(10468))*0.0025



In [35]:

    
Time.shape









    Out[35]:





(10468,)



In [36]:

    
10/0.005









    Out[36]:





2000.0



In [37]:

    
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.show()
    
    Label_ICs.append(raw_input())
    if Label_ICs[j]=='':
        Good_ICs[j]=0
    else:
        Good_ICs[j]=1









    



0






    












    












    



Ca
1






    












    












    



Ca
2






    












    












    



Ca
3






    












    












    



Ca
4






    












    












    



Ca
5






    












    












    



Ca
6






    












    












    



Ca
7






    












    












    



8






    












    












    



9






    












    












    



Ca
10






    












    












    



11






    












    












    



Ca
12






    












    












    



Ca
13






    



/usr/local/lib/python2.7/dist-packages/numpy/ma/core.py:4085: UserWarning: Warning: converting a masked element to nan.
  warnings.warn("Warning: converting a masked element to nan.")






    












    












    



Ca
14






    












    












    



Ca
15






    












    












    



Ca
16






    












    












    



Ca
17






    












    












    



LH
18






    












    












    



Ca
19






    












    












    



Ca
20






    












    












    



PN
21






    












    












    



22






    












    












    



KC
23






    












    












    



KC
24






    












    












    



25






    












    












    



26






    












    












    



27






    












    












    



28






    












    












    



29






    












    












    



LH
30






    












    












    



31






    












    












    



32






    












    












    



PN
33






    












    












    



34






    












    












    



35






    












    












    



36






    












    












    



PB
37






    












    












    



Ca
38






    












    












    



A
39






    












    












    



PN
40






    












    












    



41






    












    












    



KC
42






    












    












    



43






    












    












    



PN
44






    












    












    



PN
45






    












    












    



46






    












    












    



47






    












    












    



48






    












    












    



49






    












    












    



50






    












    












    



A
51






    












    












    



52






    












    












    



53






    












    












    



54






    












    












    



A
55






    












    












    



PB
56






    












    












    



57






    












    












    



58






    












    












    



A
59






    












    












    



60






    












    












    



61






    












    












    



62






    












    












    



A
63






    












    












    



No
64






    












    












    



65






    












    












    



No
66






    












    












    



67






    












    












    



68






    












    












    



69






    












    












    



70






    












    












    



71






    












    












    



72






    












    












    



73






    












    












    



74






    












    












    



75






    












    












    



76






    












    












    



77






    












    












    



78






    












    












    



79






    












    












    



80






    












    












    



81






    












    












    



82






    












    












    



83






    












    












    



84






    












    












    



85






    












    












    



86






    












    












    



87






    












    












    



88






    












    












    



PI
89






    












    












    



90






    












    












    



91






    












    












    



A
92






    












    












    



93






    












    












    



EB
94






    












    












    



95






    












    












    



96






    












    












    



97






    












    












    



98






    












    












    



99



In [39]:

    
D2.shape









    Out[39]:





(99, 47, 5, 100)



In [40]:

    
D2=Dmaps[:,:,:,:]



In [41]:

    
'a'









    Out[41]:





'a'



In [43]:

    
np.save('792',[Label_ICs, Good_ICs])



In [44]:

    
set(Label_ICs)









    Out[44]:





{'', 'A', 'Ca', 'EB', 'KC', 'LH', 'No', 'PB', 'PI', 'PN'}



In [45]:

    
if len(Label_ICs)<S[3]:
    for j in range(S[3]-len(Label_ICs)):
      Label_ICs.append('')



In [46]:

    
#Dict={'a':0,'AL':0,'Co':1,'co':1,'CO':1,'TM':1,'Tm':1,'T':1,'Mt':2,'ME':2,'O':2,'OL':2,'o':2,'OG':3,'A':5,'G':4,'PN':5,'PI':5,'CA':5,'Ca':5,'L':5,'LH':5,'lh':5,'l':5,'iPN':5,'KC':7,'kc':7,'M':7,'MB':8,'FB':9,'F':9,'EB':10,'eb':10,'CC':11,'No':9,'C':11,'c':11,'PB':11,'pb':11,'AMMC':13,'V':12,'S':13,'s':13,'PS':13,'SOG':13,'I':12,'i':12,'D':15,'CB':15,'DCB':15,'PCB':15,'':16,'mvt':16,'0':16}



In [47]:

    
Dict={'a':0,'AL':0,'A':0,'G':1,'PN':1,'pn':1,'CA':2,'Ca':2,'LH':3,'lh':3,'KC':4,'kc':4,'MBON':5,'FB':6,'EB':7,'No':7,'PB':8,'pb':8,'AMMC':9,'PS':10,'SOG':10,'PI':11,'OG':12,'L':12,'SLP':13,'SMP':13,'DCB':14,'CB':15,'PCB':15,'':16,'mvt':16,'0':16}



In [48]:

    
Translated=[Dict[X] for X in Label_ICs]



In [49]:

    
G=Good_ICs.tolist();



In [50]:

    
len(Good_ICs)









    Out[50]:





100



In [51]:

    
G.count(1)









    Out[51]:





39



In [52]:

    
where_are_NaNs = np.isnan(D2)

D2[where_are_NaNs] = 0



In [53]:

    
D3=np.zeros([S[0],S[1],5,S[3]])



In [54]:

    
if S[2]>5:
    for i in range(5):
        D3[:,:,i,:]=np.mean(D2[:,:,Indices[i],:],2)



In [96]:

    
if S[2]>5:
    Final_map=np.zeros([S[0],S[1],5,3])
    Fmaps=np.zeros([S[0],S[1],5,3])
else:
    Final_map=np.zeros([S[0],S[1],3])
    Fmaps=np.zeros([S[0],S[1],3])    
C=np.zeros([S[3],3])



In [97]:

    
for j in range(S[3]):    
    if Good_ICs[j]:
        C[j,:]=np.squeeze(np.random.rand(3,1))
        for k in range(3):
            if S[2]>5:
                M=np.max(np.squeeze(np.reshape(D3[:,:,:,j],S[0]*S[1]*5)))
            else:
                M=np.max(np.squeeze(np.reshape(Dmaps[:,:,:,j],S[0]*S[1]*S[2])))
            Fmaps[:,:,:,k]=0.8*D3[:,:,:,j]*C[j,k]/np.max(np.squeeze(np.reshape(D3[:,:,:,j]*np.max(C[j,:]),S[0]*S[1]*5)))
        Final_map=Final_map+Fmaps



In [98]:

    
pylab.rcParams['figure.figsize'] = (14, 5)

#Final_map[Final_map<0.1]=np.NaN

if S[2]>5:
    N=Nstack
else:
    N=S[2]
for i in range(N):
        plt.subplot(1,N,i+1)
        #plt.imshow(Dmean[:,:,i],cmap=plt.cm.gray)
        plt.imshow(Final_map[:,:,i],interpolation='none')
        #plt.imshow(Final_map[:,:,i]) 
        frame1 = plt.gca()
        frame1.axes.get_xaxis().set_visible(False)
        frame1.axes.get_yaxis().set_visible(False)



In [99]:

    
pylab.rcParams['figure.figsize'] = (8, 8)
h=0.4
i=0

for j in range(S[3]):
    if Good_ICs[j]:
        plt.plot((DT[:,j]+h*i),color=C[j,:])
        i=i+1
        
        #print jdel Final_map
axes=plt.gca()
#axes.set_ylim([-0.1,0.4])
#axes.set_xlim([2000,3000])
plt.show









    Out[99]:





<function matplotlib.pyplot.show>

'a'



In [ ]:



In [56]:

    
time=range(30)
tar=np.array(time)



In [455]:

    
tar.shape









    Out[455]:





(30,)



In [457]:

    
plt.plot(tar*0.005,(DT[2690:2720,11]),color=C[j,:])









    Out[457]:





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



In [390]:

    
List1=[(Translated[i],i) for i in range(S[3])]



In [391]:

    
Newlist=sorted(List1, key=lambda List1: List1[0])



In [392]:

    
Neworder=[Newlist[i][1] for i in range(S[3])  if Newlist[i][0] != 16]



In [393]:

    
Neworder









    Out[393]:





[112, 7, 9, 16, 25, 68, 11, 12, 13, 14, 28, 41, 40, 24, 29, 33, 46, 5, 21]



In [394]:

    
Order[Neworder]









    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-394-a8a72a14421a> in <module>()
----> 1 Order[Neworder]

NameError: name 'Order' is not defined



In [395]:

    
[Label_ICs[Neworder[i]] for i in range(len(Neworder))]









    Out[395]:





['a',
 'Ca',
 'Ca',
 'LH',
 'Ca',
 'LH',
 'KC',
 'KC',
 'MB',
 'MB',
 'MB',
 'MB',
 'EB',
 'PB',
 'PB',
 'PB',
 'PB',
 'PS',
 'D']



In [397]:

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



In [398]:

    
NewDT.shape









    Out[398]:





(19, 10468)



In [399]:

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



In [400]:

    
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)



In [401]:

    
S1=DT.shape
S1









    Out[401]:





(10468, 150)



In [402]:

    
h=8
pylab.rcParams['figure.figsize'] = (7, 14)
i=0
t=[j*0.005 for j in range(S1[0])]

for j in range(len(Neworder)):
        plt.plot(t,NewDT[j,:]+h*j,color=C1[i%6][:])
        i=i+1
plt.ylim([-3,h*j+5])

frame1=plt.gca()
frame1.axes.get_yaxis().set_ticks([])
matplotlib.rcParams.update({'font.size': 18})
plt.show()



In [412]:

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



In [413]:

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



In [414]:

    
[Translated[Neworder[i]] for i in range(len(Neworder))]









    Out[414]:





[0, 5, 5, 5, 5, 5, 7, 7, 8, 8, 8, 8, 10, 11, 11, 11, 11, 13, 15]



In [415]:

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



In [416]:

    
Regionname=[]



In [417]:

    
Newmaps.shape









    Out[417]:





(95, 48, 9, 19)



In [418]:

    
Nstack









    Out[418]:





5



In [419]:

    
Mapsordered=datao[:,:,:,Neworder[:]]
mapfilename=filename2+'mapsordered.nii'
Mapsordered2=np.transpose(Mapsordered,(3,2,1,0))
#nimap = NiftiImage(Mapsordered2)
#nimap.save(mapfilename)



In [420]:

    
DMapsordered=Dmapso[:,:,:,Neworder[:]]
mapfilename=filename2+'zscoredmapsordered.nii'
DMapsordered2=np.transpose(DMapsordered,(3,2,1,0))
nimap = NiftiImage(DMapsordered2)
nimap.save(mapfilename)









    



/usr/lib/pymodules/python2.7/nifti/image.py:231: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.
  return (not self._data == None)



In [422]:

    
j=0
i=0
k=Translated[Neworder[i]]
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=Translated[Neworder[i]]
        
    if k1 != k:
        j=j+1
        k=k1
        m=0
        Regionname.append(Label_ICs[Neworder[i]])



In [423]:

    
Regionname









    Out[423]:





['a', 'a', 'Ca', 'KC', 'MB', 'EB', 'PB', 'PS', 'D']



In [424]:

    
Datasort.shape









    Out[424]:





(95, 48, 9, 8, 3)



In [425]:

    
DMapscolor=np.sum(Datasort, axis=3)
DMapscolor.shape









    Out[425]:





(95, 48, 9, 3)



In [426]:

    
mapfilename=filename2+'zscoredmapscolor.nii'
DMapscolor2=np.transpose(DMapscolor,(3,2,1,0))



In [427]:

    
DMapscolor2.shape









    Out[427]:





(3, 9, 48, 95)



In [428]:

    
nimap = NiftiImage(DMapscolor2)
nimap.save(mapfilename)



In [429]:

    
pylab.rcParams['figure.figsize'] = (7, 14)
import scipy
from scipy import ndimage

for i in range(L):
        plt.subplot(L,1,L-i)
        Rotated_Plot = ndimage.rotate(Regionmaps[:,:,i], -90)
        IM=plt.imshow(Rotated_Plot) 
        frame1 = plt.gca()
        frame1.axes.get_xaxis().set_visible(False)
        frame1.axes.get_yaxis().set_visible(False)



In [97]:

    
Sr=Regionmaps.shape
RegionmapsT=np.zeros([Sr[1],Sr[0],Sr[2],Sr[3]])
for i in range(Sr[2]):
    for j in range(Sr[3]):
        RegionmapsT[:,:,i,j]=Regionmaps[:,:,i,j].T
RegionmapsT[RegionmapsT==0]=np.nan
Dmeanproj=np.max(Dmean,2)



In [ ]:

    
pylab.rcParams['figure.figsize'] = (7, 14)
import scipy
from scipy import ndimage
#from skimage import color

Dmeanprojrgb=np.dstack((Dmeanproj.T,Dmeanproj.T,Dmeanproj.T))

for i in range(L):
        plt.subplot(L,1,L-i)
        #plt.imshow(0.5*Dmeanproj.T,cmap=plt.cm.gray)
        plt.imshow(0.05*Dmeanprojrgb)
        plt.imshow(np.squeeze(RegionmapsT[:,:,i]),interpolation='none')
        #Rotated_Plot = ndimage.rotate(Regionmaps[:,:,i], -90)
        #IM=plt.imshow(Rotated_Plot)
        #Rotated_Plot2= plt.imshow(mp.max(Dmean,3),cmap=plt.cm.gray)
        frame1 = plt.gca()
        frame1.axes.get_xaxis().set_visible(False)
        frame1.axes.get_yaxis().set_visible(False)



In [232]:

    
import scipy.io as sio
sio.savemat(foldername+'NewDT.mat',{'NewDT':NewDT})



In [233]:

    
names=set([Label_ICs[Neworder[i]] for i in range(len(Neworder))])  
for i in range(L):
    regionfilename=foldername+Regionname[i]+'.nii'
    D3=np.transpose(Datasort[:,:,:,i,:],(2,3,1,0))
    nim = NiftiImage(D3)
    nim.save(regionfilename)



In [234]:

    
%notebook -e 862_150.ipynb



In [ ]:

    
%store