

In [3]:

    
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy import io
import scipy.io as sio
%matplotlib inline 
import pylab
import csv
from Tkinter import Tk
from tkFileDialog import askopenfilename
from tkFileDialog import askdirectory
import nibabel as nb
from scipy import io
import nibabel as nb
from scipy.interpolate import interp1d
from scipy import ndimage



In [4]:

    
from sklearn import linear_model

Open data



In [66]:

    
# 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/test5/FreeBehaviorPanNeuronalGCaMP6/100133/old/100133Final/100133ss2on250cregcdFF30sMpsfkfint599Smith0_4_60TS.mat



In [5]:

    
filename="/media/test5/FreeBehaviorPanNeuronalGCaMP6/100133/old/100133Final/100133ss2on250cregcdFF30sMpsfkfint599Smith0_4_60TS.mat"



In [6]:

    
Ua=sio.loadmat(filename)



In [7]:

    
Ua









    Out[7]:





{'TSo': array([[-0.02002288, -0.0005812 , -0.00710151, ..., -0.00947755,
          0.01904074,  0.01483885],
        [-0.01774609, -0.00212954, -0.00717034, ..., -0.0073894 ,
          0.0137778 ,  0.0079594 ],
        [-0.01546929, -0.00367788, -0.00723916, ..., -0.00530125,
          0.00851486,  0.00107995],
        ..., 
        [-0.00856585, -0.00344055, -0.00847336, ...,  0.00019336,
          0.00141308,  0.00258333],
        [-0.007443  , -0.00397578, -0.00766295, ...,  0.00028112,
          0.00121594,  0.00225143],
        [-0.00632016, -0.00451102, -0.00685254, ...,  0.00036887,
          0.0010188 ,  0.00191952]]),
 '__globals__': [],
 '__header__': 'MATLAB 5.0 MAT-file, Platform: GLNXA64, Created on: Sat Nov  5 17:21:57 2016',
 '__version__': '1.0'}



In [8]:

    
Ua=sio.loadmat(filename)
DT=Ua['TSo']
DT.shape









    Out[8]:





(20651, 599)



In [5]:

    
# 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
filename2 = askopenfilename() # show an "Open" dialog box and return the path to the selected file
print(filename2)









    



/media/test5/FreeBehaviorPanNeuronalGCaMP6/100133/old/100133Final/100133ss2on250cregcdFF30sMpsfkfint599Smith0_4_60IC.nii



In [9]:

    
filename2="/media/test5/FreeBehaviorPanNeuronalGCaMP6/100133/old/100133Final/100133ss2on250cregcdFF30sMpsfkfint599Smith0_4_60IC.nii"



In [10]:

    
img1 = nb.load(filename2)
data = img1.get_data()
S=data.shape
S









    Out[10]:





(180, 97, 10, 599)

Z-score



In [11]:

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



In [12]:

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



In [13]:

    
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-1
    Tvar[i]=np.var(DT[i,:])
Dmaps[Dmaps<0]=0



In [10]:

    
# 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/test5/FreeBehaviorPanNeuronalGCaMP6/100133/old/100133Final/100133Xk.mat



In [14]:

    
filename="/media/test5/FreeBehaviorPanNeuronalGCaMP6/100133/old/100133Final/100133Xk.mat"



In [15]:

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



In [ ]:

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



In [16]:

    
filenamet="/media/test5/FreeBehaviorPanNeuronalGCaMP6/100133/old/100133Final/AVG_100133ss2on250cregcpsf.nii"



In [17]:

    
nimt=nb.load(filenamet)
Dtemp=np.squeeze(nimt.get_data())
Dtemp.shape









    Out[17]:





(180, 97, 10)

Fit turns



In [18]:

    
%%javascript
IPython.OutputArea.auto_scroll_threshold =4000;



In [19]:

    
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=np.squeeze(data[:,:,range(Nstack),2])
    for i in range(Nstack):
        Vmean=np.mean(Dtemp[:,:,Indices[i]],2)
        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 [20]:

    
for j in range(S[3]):

    a=''
    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



In [21]:

    
plt.imshow(Dmean[:,:,1],cmap=plt.cm.gray)









    Out[21]:





<matplotlib.image.AxesImage at 0x7f03a067e490>



In [22]:

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



In [23]:

    
algorithm = linear_model.LinearRegression()



In [24]:

    
Sxk=Xk.shape



In [25]:

    
Sxk









    Out[25]:





(20651, 8)



In [26]:

    
X=np.zeros((Sxk[0],2))



In [27]:

    
X[:,0]=(Xk[:,0]-np.mean(Xk[:,0]))/np.std(Xk[:,0])
X[:,1]=(Xk[:,1]-np.mean(Xk[:,1]))/np.std(Xk[:,1])
#X[:,2]=(Xk[:,3]-np.mean(Xk[:,3]))/np.std(Xk[:,3])
#X[:,3]=(Xk[:,4]-np.mean(Xk[:,4]))/np.std(Xk[:,4])
#X[:,4]=(Xk[:,6]-np.mean(Xk[:,6]))/np.std(Xk[:,6])
#X[:,5]=(Xk[:,7]-np.mean(Xk[:,7]))/np.std(Xk[:,7])



In [28]:

    
Center=20651
CCcb=np.correlate(DT[:,13],DT[:,7],'full')
#CCry=np.correlate(DT[:,4],DT[:,11],'full')
print(np.argmax(CCcb)-Center)
#print(np.argmax(CCry)-Center)
plt.plot(CCcb[range(Center-200,Center+500)])
#plt.plot(CCry[range(Center-200,Center+500)])









    



15






    Out[28]:





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



In [29]:

    
DT.shape









    Out[29]:





(20651, 599)



In [30]:

    
plt.plot(np.array(range(20651))*0.03,X[:,0]/2)
plt.plot(np.array(range(20651))*0.03,X[:,1]/2)

plt.plot(np.array(range(20651))*0.03,DT[:,7]*2-2,'c')
plt.plot(np.array(range(20651))*0.03,DT[:,13]*2-3,'b')
plt.plot(np.array(range(20651))*0.03,DT[:,8]*2-4,'y')
#plt.plot(X[:,2])
#plt.plot(X[:,3])
#plt.plot(X[:,4])
#plt.plot(X[:,5])









    Out[30]:





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

Clearer with components than averaged brain regions.

Plot all components for turning left, right, walking, and grooming



In [31]:

    
Rsq=np.zeros((1,S[3]))
Betas=np.zeros((2,S[3]))



In [32]:

    
X.shape









    Out[32]:





(20651, 2)



In [33]:

    
DT.shape









    Out[33]:





(20651, 599)



In [34]:

    
for j in range(S[3]):
    model = algorithm.fit(X, DT[:,j])
    Betas[:,j] = model.coef_
    Rsq[:,j] = model.score(X,DT[:,j])



In [35]:

    
RsqUni=np.zeros((6,S[3]))
BetaUni=np.zeros((6,S[3]))



In [36]:

    
Sx=X.shape



In [37]:

    
for k in range(2):
    for j in range(S[3]):
        model = algorithm.fit(np.reshape(X[:,k],(Sx[0],1)), DT[:,j])
        BetaUni[k,j] = model.coef_
        RsqUni[k,j] = model.score(np.reshape(X[:,k],(Sx[0],1)),DT[:,j])



In [38]:

    
plt.plot(RsqUni[0,:])









    Out[38]:





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



In [39]:

    
import random



In [40]:

    
del Final_map
del Fmaps









    



-------------------------------------------------
NameError       Traceback (most recent call last)
<ipython-input-40-afdbaf60a348> in <module>()
----> 1 del Final_map
      2 del Fmaps

NameError: name 'Final_map' is not defined



In [44]:

    
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



In [45]:

    
C=np.zeros((S[3],3))
i=0
l=0
Betas2=Betas
for j in range(S[3]):  
    if Betas2[0,j]>0.87*np.max(Betas2[0,:]):
    #if 1>0.1:
        #C[j,:]=C1[i%6][:]
        C[j,2]=1
        C[j,1]=random.uniform(0,1)
        #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.8*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(Rsq[:,j])
        #if l==2:
            #break









    



5
[ 0.59363838]
8
[ 0.68349285]
14
[ 0.65424356]



In [46]:

    
C=np.zeros((S[3],3))
i=0
l=0
Betas2=Betas
for j in range(S[3]):  
    if Betas2[1,j]>0.78*np.max(Betas2[1,:]):
    #if 1>0.1:
        #C[j,:]=C1[i%6][:]
        C[j,0]=1.5
        C[j,1]=random.uniform(0,1)
        #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.8*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)
        #if l==2:
         #   break

5
9



In [47]:

    
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/30
    #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 [53]:

    
del Final_map
del Fmaps



In [55]:

    
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
Betas2=Betas
for j in range(6,S[3]):  
    if (Betas2[0,j]>0.8*np.max(Betas2[0,:])) and (j != 9):
    #if 1>0.1:
        #C[j,:]=C1[i%6][:]
        C[j,2]=1
        C[j,1]=i/2
        #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.75*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,:])
        print(RsqUni[:,j])
        #if l==2:
            #break
    

C=np.zeros((S[3],3))
i=1
l=0
Betas2=Betas
for j in range(6,S[3]):  
    if Betas2[1,j]>0.8*np.max(Betas2[1,:]) and j != 9:
    #if 1>0.1:
        #C[j,:]=C1[i%6][:]
        C[j,0]=1
        C[j,1]=1-i/2
        #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.75*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,:])
        print(RsqUni[:,j])
        #if l==2:
         #   break

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









    



8
[ 0.  0.  1.]
[ 0.61099972  0.05273261  0.          0.          0.          0.        ]
14
[ 0.  1.  1.]
[ 0.59279346  0.04363842  0.          0.          0.          0.        ]
9
[ 1.  1.  0.]
[ 0.21378623  0.3889573   0.          0.          0.          0.        ]



In [44]:

    
# 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
filenameM = askopenfilename() # show an "Open" dialog box and return the path to the selected file
print(filenameM)
img1 = nb.load(filenameM)
Masks = img1.get_data()
Sm=Masks.shape
Masks=np.array(Masks)









    



/media/test5/FreeBehaviorPanNeuronalGCaMP6/100133/100133Registration/JFRC100133Transformedfullpsftrimmed.nii



In [50]:

    
Masks.shape









    Out[50]:





(180, 97, 10, 87)



In [78]:

    
filenameM='/home/sophie/RegionList'
with open(filenameM) as f:
    content = f.readlines()
Names=[Line.split('\t') for Line in content]
RegionName=[Names[i][0] for i in range(75)]
Num=[int(Names[i][2]) for i in range(75)]

# Average in masks to sort components by brain region

M=np.zeros((S[3],86))
Mapmean=np.zeros(S[3])
MMasks=np.zeros(86)

for i in range(S[3]):
    Mapmean[i]=np.mean(np.mean(np.mean(Dmaps[:,:,:,i])))
    for j in range(86):
        MMasks[j]=np.mean(np.mean(np.mean(Masks[:,:,:,j])))
        if MMasks[j]:
            M[i,j]=np.mean(np.mean(np.mean(Masks[:,:,:,j]*Dmaps[:,:,:,i])))/(MMasks[j]*Mapmean[i])


NoList=[]
M[98,:]=0

CompMainName=S[3]*['']
CompNameAdd=np.zeros((S[3],86))
for i in range(S[3]):
    Max=np.max(M[i,:])
    I=np.argmax(M[i,:])+1
    if M[i,5]>0.5*Max:
        NoList.append(i)
    J=[l for l in range(74) if Num[l]==I]
    if J!=[]:
        CompMainName[i]=Names[np.array(J)][0]

NoList









    



/usr/local/lib/python2.7/dist-packages/ipykernel/__main__.py:34: VisibleDeprecationWarning: converting an array with ndim > 0 to an index will result in an error in the future






    Out[78]:





[3, 5, 14, 27, 50, 480, 529, 558]



In [88]:

    
del Final_map
del Fmaps

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=0
l=0
Betas2=Betas
for j in range(S[3]):  
    #if j in NoList:
    if j==572:
    #if 1>0.1:
        #C[j,:]=C1[i%6][:]
        C[j,2]=1
        C[j,1]=1-i/8
        #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(Rsq[:,j])
        print(Betas[:,j])
        print(RsqUni[:,j])
        print(BetaUni[:,j])
        #if l==2:

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









    



573
[ 0.03969343]
[-0.00040273  0.00717512]
[ 0.00044518  0.03956978  0.          0.          0.          0.        ]
[-0.0007632   0.00719535  0.          0.          0.          0.        ]



In [ ]: