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

    
# 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/sophie/db554c18-e3eb-41e2-afad-7de1c92bf4a5/THDDCGCaMP62/100597series/1000601ss1oncregdFF2000pointspsfkfconcat186Smith0_4_60TS.mat



In [5]:

    
Ua=sio.loadmat(filename)



In [6]:

    
DT=Ua['TSo']



In [7]:

    
DT.shape









    Out[7]:





(10444, 186)



In [8]:

    
S1=DT.shape



In [9]:

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



In [10]:

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



In [11]:

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



In [12]:

    
DTvar.shape









    Out[12]:





(10444, 186)

open maps



In [13]:

    
import nibabel as nb



In [14]:

    
# 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/sophie/db554c18-e3eb-41e2-afad-7de1c92bf4a5/THDDCGCaMP62/100597series/1000601ss1oncregdFF2000pointspsfkfconcat186Smith0_4_60IC.nii



In [15]:

    
img1 = nb.load(filename2)



In [16]:

    
data = img1.get_data()



In [17]:

    
S=data.shape



In [18]:

    
S









    Out[18]:





(85, 73, 11, 186)

Zscore maps



In [19]:

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

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

    
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



In [77]:

    
# 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/sophie/db554c18-e3eb-41e2-afad-7de1c92bf4a5/THDDCGCaMP62/100597series/100601/MAX_1000601ss1oncregdFF2000pointspsfkfconcat186Smith0_4_60ICtemp.nii






    Out[77]:





<matplotlib.image.AxesImage at 0x7f6e9aa31890>

Order ICs by variance



In [78]:

    
datao=data
Dmapso=Dmaps



In [79]:

    
plt.plot(Var)









    Out[79]:





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

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



In [80]:

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

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



In [82]:

    
%%javascript
IPython.OutputArea.auto_scroll_threshold =4000;



In [83]:

    
DTvar.shape









    Out[83]:





(10444, 186)



In [84]:

    
S









    Out[84]:





(85, 73, 11, 186)



In [64]:

    
# 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/sophie/db554c18-e3eb-41e2-afad-7de1c92bf4a5/THDDCGCaMP62/100597series/100601/100601X.mat



In [65]:

    
Ua=sio.loadmat(filename)
Xk=Ua['X']



In [66]:

    
Xk.shape









    Out[66]:





(4, 5206)



In [70]:

    
DTvarnew=DTvar[range(1,5206),:]



In [71]:

    
DTvarnew.shape









    Out[71]:





(5205, 186)



In [85]:

    
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(DTvarnew[:,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 [86]:

    
for j in range(S[3]):

    if Good_ICs[j]==1:       
        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(DTvarnew[:,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()



In [147]:

    
List1









    Out[147]:





[('', 0),
 ('', 1),
 ('', 2),
 ('', 3),
 ('', 4),
 ('', 5),
 ('', 6),
 ('', 7),
 ('', 8),
 ('', 9),
 ('', 10),
 ('', 11),
 ('', 12),
 ('', 13),
 ('DCB', 14),
 ('SOG', 15),
 ('', 16),
 ('', 17),
 ('DCB', 18),
 ('', 19),
 ('', 20),
 ('DCB', 21),
 ('', 22),
 ('DCB', 23),
 ('', 24),
 ('', 25),
 ('DCB', 26),
 ('', 27),
 ('', 28),
 ('', 29),
 ('', 30),
 ('DCB', 31),
 ('MB', 32),
 ('', 33),
 ('OL', 34),
 ('Cup', 35),
 ('', 36),
 ('SOG', 37),
 ('', 38),
 ('', 39),
 ('', 40),
 ('', 41),
 ('', 42),
 ('', 43),
 ('', 44),
 ('AL', 45),
 ('', 46),
 ('', 47),
 ('', 48),
 ('', 49),
 ('DCB', 50),
 ('', 51),
 ('SOG', 52),
 ('', 53),
 ('', 54),
 ('DCB', 55),
 ('MB', 56),
 ('', 57),
 ('', 58),
 ('', 59),
 ('', 60),
 ('', 61),
 ('', 62),
 ('', 63),
 ('', 64),
 ('', 65),
 ('', 66),
 ('', 67),
 ('MB', 68),
 ('', 69),
 ('', 70),
 ('', 71),
 ('SOG', 72),
 ('', 73),
 ('', 74),
 ('', 75),
 ('AL', 76),
 ('MB', 77),
 ('', 78),
 ('', 79),
 ('', 80),
 ('DCB', 81),
 ('', 82),
 ('', 83),
 ('', 84),
 ('SOG', 85),
 ('', 86),
 ('', 87),
 ('', 88),
 ('', 89),
 ('', 90),
 ('', 91),
 ('', 92),
 ('DCB', 93),
 ('Lat', 94),
 ('', 95),
 ('', 96),
 ('MB', 97),
 ('MB', 98),
 ('', 99),
 ('', 100),
 ('MB', 101),
 ('', 102),
 ('DCB', 103),
 ('MB', 104),
 ('', 105),
 ('', 106),
 ('Lat', 107),
 ('', 108),
 ('MB', 109),
 ('', 110),
 ('', 111),
 ('', 112),
 ('DCB', 113),
 ('', 114),
 ('', 115),
 ('', 116),
 ('', 117),
 ('', 118),
 ('', 119),
 ('', 120),
 ('', 121),
 ('MB', 122),
 ('', 123),
 ('', 124),
 ('', 125),
 ('', 126),
 ('SOG', 127),
 ('', 128),
 ('', 129),
 ('Dorsal', 130),
 ('', 131),
 ('MB', 132),
 ('', 133),
 ('', 134),
 ('', 135),
 ('', 136),
 ('', 137),
 ('', 138),
 ('', 139),
 ('', 140),
 ('MB', 141),
 ('DCB', 142),
 ('AL', 143),
 ('', 144),
 ('FB', 145),
 ('', 146),
 ('', 147),
 ('LAT', 148),
 ('', 149),
 ('', 150),
 ('', 151),
 ('', 152),
 ('', 153),
 ('MB', 154),
 ('', 155),
 ('', 156),
 ('', 157),
 ('DCB', 158),
 ('MB', 159),
 ('', 160),
 ('', 161),
 ('', 162),
 ('', 163),
 ('', 164),
 ('', 165),
 ('', 166),
 ('', 167),
 ('', 168),
 ('', 169),
 ('', 170),
 ('', 171),
 ('', 172),
 ('', 173),
 ('', 174),
 ('', 175),
 ('FSB', 176),
 ('', 177),
 ('', 178),
 ('', 179),
 ('', 180),
 ('', 181),
 ('', 182),
 ('', 183),
 ('', 184),
 ('', 185)]



In [153]:

    
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([13,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)
C1[6][:]=(1,0.5,0)
C1[7][:]=(0,1,0.5)
C1[8][:]=(0.5,0,1)
C1[9][:]=(0.8,0.8,0.5)
C1[10][:]=(0.5,1,1)
C1[11][:]=(1,0.5,1)
C1[12]=(0.5,0.5,0.5)
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.6*np.max(DMapsordered[:,:,:,i],2)*C1[i%13][l]/M
        Datasort[:,:,:,j,l]=Datasort[:,:,:,j,l]+Dmaps[:,:,:,Neworder[i]]*C1[i%13][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,0:5000]+h*m,color=C1[i%13][:])
    
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 [ ]:

    
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)

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[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=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,0:5000]+h*m,color=C1[i%6][:])



In [129]:

    
#Mapsordered=data[:,:,:,Good_ICs==1]
mapfilename=filename2+'GoodIC.nii'
img = nb.Nifti1Image(DMapsordered, np.eye(4))
#Mapsordered2=np.transpose(Mapsordered,(3,2,1,0))
nb.save(img,mapfilename)



In [130]:

    
DT.shape









    Out[130]:





(10444, 186)



In [131]:

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



In [ ]: