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 [117]:
# 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)
In [118]:
Ua=sio.loadmat(filename)
Ua
Out[118]:
In [119]:
DT=Ua['TS']
In [120]:
DT.shape
Out[120]:
In [121]:
S1=DT.shape
In [122]:
DTmean=np.zeros(S1)
DTvar=np.zeros(S1)
Var=np.zeros(S1[1])
In [123]:
for i in range(S1[1]):
DTmean[:,i]=DT[:,i]-np.mean(DT[:,i],0)
In [124]:
for i in range(S1[1]):
Var[i]=np.sqrt(np.var(DTmean[:,i]))
DTvar[:,i]=DTmean[:,i]/Var[i]
In [125]:
DTvar.shape
Out[125]:
open maps
In [126]:
import nibabel as nb
In [13]:
# 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)
In [14]:
img1 = nb.load(filename2)
In [15]:
data = img1.get_data()
In [16]:
S=data.shape
In [17]:
S
Out[17]:
In [127]:
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 [128]:
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 [129]:
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 [130]:
datao=data
Dmapso=Dmaps
In [131]:
plt.plot(Var)
Out[131]:
In [132]:
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 [133]:
Dtemp=data[:,:,:,0]
In [134]:
%%javascript
IPython.OutputArea.auto_scroll_threshold =4000;
In [135]:
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 [136]:
DTvar.shape
Out[136]:
In [137]:
S
Out[137]:
In [29]:
# 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)
In [30]:
Ua=sio.loadmat(filename)
Xk=Ua['Xk']
#Xk[1,:]=Ua['Walk']
In [54]:
Xk.shape
Out[54]:
In [138]:
# 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
In [139]:
plt.imshow(Vmean,cmap=plt.cm.gray)
Out[139]:
In [33]:
Xk=Xk.T
In [34]:
Xk.shape
Out[34]:
In [144]:
Label_ICs=[]
In [104]:
Xksmoothed=np.zeros(Xk.shape)
for i in range(Xk.shape[1]):
Xksmoothed[:,i]=np.mean(Xk[:,max(0,i-999):min(Xk.shape[1],i+1000)],1)
Xkdff=Xk[range(3,6),:]-Xksmoothed[range(3,6),:]
In [141]:
plt.plot(Xksmoothed.T)
plt.show()
In [107]:
plt.plot(Xkdff.T)
Out[107]:
In [142]:
TimeCCMax=np.zeros(S[3])
In [149]:
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])
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)
S1=DT.shape
C=np.zeros((S[3],3))
i=1
l=0
for j in range(S[3]):
if Label_ICs[j]=='a' and TimeCCMax[j]>0 and TimeCCMax[j]<0.3:
#if 1>0.1:
#C[j,:]=C1[i%6][:]
C[j,0]=10*TimeCCMax[j]
C[j,1]=1-10*TimeCCMax[j]
#C[j,2]=1
for k in range(3):
M=np.max(np.squeeze(np.reshape(D2[:,:,:,j],S[0]*S[1]*5)))
Fmaps[:,:,:,k]=0.2*D2[:,:,:,j]*C[j,k]/M
Final_map=Final_map+Fmaps
#Betas2[1,j]=0
#print(Indexo[j])
i=i+1
l=l+1
#print(j+1)
#print(C[j,:])
#if l==2:
# break
for j in range(S[3]):
if Label_ICs[j]=='a' and TimeCCMax[j]<0 and TimeCCMax[j]>-0.3:
#if 1>0.1:
#C[j,:]=C1[i%6][:]
#print(j+1)
C[j,2]=-6*TimeCCMax[j]
C[j,1]=1+6*TimeCCMax[j]
#C[j,2]=1
for k in range(3):
M=np.max(np.squeeze(np.reshape(D2[:,:,:,j],S[0]*S[1]*5)))
Fmaps[:,:,:,k]=0.5*D2[:,:,:,j]*C[j,k]/M
Final_map=Final_map+Fmaps
#Betas[0,j]=0
#print(Indexo[j])
i=i+1
l=l+1
#print(j+1)
#print(C[j,:])
#if l==2:
#break
C=np.zeros((S[3],3))
i=1
l=0
pylab.rcParams['figure.figsize'] = (15, 6)
C2=np.zeros(3)
Df=np.zeros([S[0],S[1],5,3])
for i in range(3):
Df[:,:,:,i]=Final_map[:,:,:,i]+Dmean/90
#Df=Df/(np.max(np.max(np.max(Df),3)))
Df[Df>1]=1
Df[Df<0]=0
if S[2]>5:
N=Nstack
else:
N=S[2]
for i in range(N):
#if Good_ICs[j]:
plt.subplot(1,N,i+1)
plt.imshow(Df[:,:,i],cmap=plt.cm.gray)
plt.imshow(Df[:,:,i,:],cmap=my_cmap,interpolation='none')
frame1 = plt.gca()
frame1.axes.get_xaxis().set_visible(False)
frame1.axes.get_yaxis().set_visible(False)
plt.tight_layout(pad=0,w_pad=0,h_pad=0)
In [59]:
#zip(range(S[3]),Label_ICs)
In [ ]:
CCry=np.correlate(DT[:,5],DT[:,11],'full')
In [60]:
set(Label_ICs)
Out[60]:
In [ ]:
Label_ICs[94]='M'
In [ ]:
Label_ICs[104]='alpha'
Label_ICs[100]='alpha'
Label_ICs[138]='betap'
Label_ICs[94]='betap'
Label_ICs[90]='betap'
Label_ICs[67]='gamma'
Label_ICs[60]='gamma'
Label_ICs[24]='gamma'
Label_ICs[54]='gamma'
Label_ICs[33]='gamma'
In [85]:
Xk.shape
Out[85]:
In [66]:
Rsq=np.zeros((1,S[3]))
Betas=np.zeros((3,S[3]))
In [67]:
from sklearn import linear_model
In [68]:
algorithm = linear_model.LinearRegression()
In [70]:
for j in range(S[3]):
model = algorithm.fit(Xkdff.T, DT[:,j])
Betas[:,j] = model.coef_
Rsq[:,j] = model.score(Xkdff.T,DT[:,j])
In [71]:
Xkdff.shape
Out[71]:
In [ ]:
Label_ICsb=Label_ICs[i for i in range(4,256)]
In [ ]:
len(Label_ICsb)
In [72]:
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([16,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)
C1[13]=(0.2,0.5,0.5)
C1[14]=(0.5,0.2,0.5)
C1[15]=(0.5,0.5,0.2)
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[:]]
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%12+1][l]/M
Datasort[:,:,:,j,l]=Datasort[:,:,:,j,l]+Dmaps[:,:,:,Neworder[i]]*C1[i%15+1][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%12+1][:])
print(Neworder[i])
print(Rsq[:,Neworder[i]])
print(Betas[:,Neworder[i]])
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)
print(Neworder)
In [ ]:
plt.plot(Xkdff.T)
In [148]:
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])
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)
S1=DT.shape
C=np.zeros((S[3],3))
i=1
l=0
for j in range(S[3]):
if Label_ICs[j]=='a' and TimeCCMax[j]>0 and TimeCCMax[j]<0.3:
#if 1>0.1:
#C[j,:]=C1[i%6][:]
C[j,0]=1
C[j,1]=1-6*TimeCCMax[j]
#C[j,2]=1
for k in range(3):
M=np.max(np.squeeze(np.reshape(D2[:,:,:,j],S[0]*S[1]*5)))
Fmaps[:,:,:,k]=0.4*D2[:,:,:,j]*C[j,k]/M
Final_map=Final_map+Fmaps
#Betas2[1,j]=0
#print(Indexo[j])
i=i+1
l=l+1
print(j+1)
print(C[j,:])
#if l==2:
# break
for j in range(S[3]):
if Label_ICs[j]=='a' and TimeCCMax[j]<0 and TimeCCMax[j]>-0.3:
#if 1>0.1:
#C[j,:]=C1[i%6][:]
#print(j+1)
C[j,2]=1
C[j,1]=1+3*TimeCCMax[j]
#C[j,2]=1
for k in range(3):
M=np.max(np.squeeze(np.reshape(D2[:,:,:,j],S[0]*S[1]*5)))
Fmaps[:,:,:,k]=0.4*D2[:,:,:,j]*C[j,k]/M
Final_map=Final_map+Fmaps
#Betas[0,j]=0
#print(Indexo[j])
i=i+1
l=l+1
print(j+1)
print(C[j,:])
#if l==2:
#break
C=np.zeros((S[3],3))
i=1
l=0
pylab.rcParams['figure.figsize'] = (15, 6)
C2=np.zeros(3)
Df=np.zeros([S[0],S[1],5,3])
for i in range(3):
Df[:,:,:,i]=Final_map[:,:,:,i]+Dmean/70
#Df=Df/(np.max(np.max(np.max(Df),3)))
if S[2]>5:
N=Nstack
else:
N=S[2]
for i in range(N):
#if Good_ICs[j]:
plt.subplot(1,N,i+1)
plt.imshow(Df[:,:,i],cmap=plt.cm.gray)
plt.imshow(Df[:,:,i,:],cmap=my_cmap,interpolation='none')
frame1 = plt.gca()
frame1.axes.get_xaxis().set_visible(False)
frame1.axes.get_yaxis().set_visible(False)
plt.tight_layout(pad=0,w_pad=0,h_pad=0)
In [ ]: