notebook.community

Edit and run 





In [1]:



    


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



    


from sklearn import linear_model

Open data





In [8]:



    


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



    


filename="/media/test5/FreeBehaviorPanNeuronalGCaMP6/100106/100106Final/100106ss2on500cregcdFF20sMpsfkfint169Smith0_4_60TS.mat"



    



















In [4]:



    


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



    


















    
Out[4]:








(10608, 169)

























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/100106/100106Final/100106ss2on500cregcdFF20sMpsfkfint169Smith0_4_60IC.nii

























In [5]:



    


filename2="/media/test5/FreeBehaviorPanNeuronalGCaMP6/100106/100106Final/100106ss2on500cregcdFF20sMpsfkfint169Smith0_4_60IC.nii"



    



















In [6]:



    


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



    


















    
Out[6]:








(189, 125, 10, 169)

























Z-score
















In [70]:



    


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



    


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



    



















In [72]:



    


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/100106/100106Final/100106Xk.mat

























In [10]:



    


filename="/media/test5/FreeBehaviorPanNeuronalGCaMP6/100106/100106Final/100106Xk.mat"



    



















In [11]:



    


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



    


filenamet='/media/test5/FreeBehaviorPanNeuronalGCaMP6/100106/100106Final/AVG_100106ss2on500cregcpsf.nii'
nimt=nb.load(filenamet)
Dtemp=np.squeeze(nimt.get_data())
Dtemp.shape



    


















    
Out[12]:








(189, 125, 10)

























Fit turns
















In [47]:



    


%%javascript
IPython.OutputArea.auto_scroll_threshold =4000;



    


















    



































In [73]:



    


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



    



















In [74]:



    


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



    


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



    


















    
Out[75]:








<matplotlib.image.AxesImage at 0x7fb9ca9670d0>










    
































In [76]:



    


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



    


algorithm = linear_model.LinearRegression()



    



















In [53]:



    


Sxk=Xk.shape



    



















In [54]:



    


Sxk



    


















    
Out[54]:








(10608, 8)

























In [55]:



    


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



    



















In [56]:



    


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



    


Center=10608
CCcb=np.correlate(DT[:,4],DT[:,8],'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)])



    


















    






24










    
Out[57]:








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










    
































In [58]:



    


Time=np.array(range(6000))*0.01
plt.plot(Time,X[range(6000),0])
plt.plot(Time,X[range(6000),1])

plt.plot(Time,DT[range(6000),4]*2-2,'c')
plt.plot(Time,DT[range(6000),8]*2-3,'b')
plt.plot(Time,DT[range(6000),5]*2-4,'y')
plt.plot(Time,DT[range(6000),3]*2-5,'r')

#plt.plot(X[:,0])
#plt.plot(X[:,1])
#plt.plot(X[:,1]-X[:,0])
#plt.plot(X[:,2])
#plt.plot(X[:,3])
#plt.plot(X[:,4])
#plt.plot(X[:,5])
zero_crossings = np.where(np.diff(np.sign(X[:,1]-X[:,0])))[0]
print(zero_crossings.shape)



    


















    






(46,)










    
































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
















In [59]:



    


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



    



















In [60]:



    


X.shape



    


















    
Out[60]:








(10608, 2)

























In [61]:



    


DT.shape



    


















    
Out[61]:








(10608, 169)

























In [62]:



    


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



    



















In [63]:



    


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



    



















In [64]:



    


Sx=X.shape



    



















In [65]:



    


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



    


plt.plot(RsqUni[0,:])



    


















    
Out[66]:








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










    
































In [67]:



    


import random



    



















In [68]:



    


Betas[:,0]=0
Betas[:,1]=0
Betas[:,2]=0



    



















In [77]:



    


pylab.rcParams['figure.figsize'] = (11, 2.5)

#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 Betas2[0,j]>0.85*np.max(Betas2[0,:]):
    #if 1>0.1:
        #C[j,:]=C1[i%6][:]
        C[j,2]=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
        #Betas[0,j]=0
        #print(Indexo[j])
        i=i+1
        l=l+1
        print(j+1)
        print(RsqUni[:,j])
        print(C[j,:])
        #if l==2:
            #break
    

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



    


















    






4
[ 0.20785991  0.21197437  0.          0.          0.          0.        ]
[ 0.  1.  1.]
5
[ 0.31034004  0.0749746   0.          0.          0.          0.        ]
[ 0.  1.  1.]
9
[ 0.34584815  0.160656    0.          0.          0.          0.        ]
[ 0.  0.  1.]
4
[ 1.  0.  0.]
[ 0.20785991  0.21197437  0.          0.          0.          0.        ]
6
[ 1.  1.  0.]
[ 0.12823683  0.23013777  0.          0.          0.          0.        ]










    
































Open Masks
















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



    



















In [40]:



    


filenameM="/media/test5/FreeBehaviorPanNeuronalGCaMP6/100106/100106Registration/JFRC100106Transformedfullpsftrimmed.nii"



    



















In [41]:



    


img1 = nb.load(filenameM)
Masks = img1.get_data()
Sm=Masks.shape
Masks=np.array(Masks)



    



















In [42]:



    


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[63,:]=0
M[126,:]=0
M[166,:]=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]



    


















    






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

























In [43]:



    


NoList

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==20:
    #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)



    


















    






14
[ 0.31782934]
[ 0.15404119  0.12122238]
[ 0.18940895  0.11046115  0.          0.          0.          0.        ]
[ 0.14696898  0.11223549  0.          0.          0.          0.        ]
15
[ 0.17402635]
[ 0.11583206  0.08243694]
[ 0.11231383  0.05218713  0.          0.          0.          0.        ]
[ 0.11102262  0.07567921  0.          0.          0.          0.        ]
21
[ 0.19274535]
[ 0.1082719   0.09153523]
[ 0.10779157  0.07388482  0.          0.          0.          0.        ]
[ 0.10293166  0.08521856  0.          0.          0.          0.        ]
28
[ 0.25495697]
[ 0.09927398  0.08694054]
[ 0.13789818  0.10233036  0.          0.          0.          0.        ]
[ 0.09420181  0.08114882  0.          0.          0.          0.        ]
34
[ 0.43756993]
[ 0.12659812  0.09408965]
[ 0.27322138  0.14003589  0.          0.          0.          0.        ]
[ 0.12110886  0.08670382  0.          0.          0.          0.        ]
38
[ 0.47074982]
[ 0.10736708  0.08578229]
[ 0.27692872  0.16711748  0.          0.          0.          0.        ]
[ 0.10236247  0.07951841  0.          0.          0.          0.        ]
54
[ 0.37776397]
[ 0.06955123  0.05882518]
[ 0.21117816  0.14488991  0.          0.          0.          0.        ]
[ 0.06611932  0.0547675   0.          0.          0.          0.        ]
55
[ 0.37497031]
[ 0.07944375  0.03835702]
[ 0.30094307  0.05741358  0.          0.          0.          0.        ]
[ 0.07720597  0.03372221  0.          0.          0.          0.        ]
58
[ 0.31753909]
[ 0.06165133  0.04915502]
[ 0.18713383  0.11240195  0.          0.          0.          0.        ]
[ 0.05878359  0.04555823  0.          0.          0.          0.        ]
85
[ 0.14927033]
[-0.02491477 -0.03188899]
[ 0.05135288  0.08949908  0.          0.          0.          0.        ]
[-0.02305434 -0.03043544  0.          0.          0.          0.        ]
87
[ 0.34354948]
[ 0.03594086  0.04725595]
[ 0.11371854  0.21060438  0.          0.          0.          0.        ]
[ 0.03318391  0.04515914  0.          0.          0.          0.        ]
101
[ 0.05062846]
[ 0.01436464  0.01403884]
[ 0.02445233  0.0232233   0.          0.          0.          0.        ]
[ 0.0135456  0.0132008  0.         0.         0.         0.       ]
109
[ 0.19499037]
[ 0.02739957  0.0255771 ]
[ 0.0989122   0.08473251  0.          0.          0.          0.        ]
[ 0.02590738  0.02397859  0.          0.          0.          0.        ]
112
[ 0.01416953]
[-0.00914728 -0.00254148]
[ 0.01312616  0.00065343  0.          0.          0.          0.        ]
[-0.00899901 -0.00200782  0.          0.          0.          0.        ]
118
[ 0.00325414]
[ 0.00346081  0.0026867 ]
[ 0.00196155  0.00110939  0.          0.          0.          0.        ]
[ 0.00330407  0.00248479  0.          0.          0.          0.        ]
131
[ 0.44718314]
[-0.03772359 -0.02559449]
[ 0.29864877  0.12451202  0.          0.          0.          0.        ]
[-0.03623039 -0.02339367  0.          0.          0.          0.        ]
137
[ 0.08940748]
[-0.01596753 -0.00872887]
[ 0.06784757  0.01726241  0.          0.          0.          0.        ]
[-0.01545828 -0.00779731  0.          0.          0.          0.        ]
143
[ 0.10312955]
[ 0.01108781  0.01414498]
[ 0.03564474  0.06166346  0.          0.          0.          0.        ]
[ 0.01026258  0.01349811  0.          0.          0.          0.        ]
144
[ 0.14761188]
[ 0.01761711  0.01157534]
[ 0.10074315  0.03904803  0.          0.          0.          0.        ]
[ 0.0169418   0.01054754  0.          0.          0.          0.        ]
151
[ 0.00026218]
[ 0.00077649 -0.00019013]
[  2.47789364e-04   2.21421824e-05   0.00000000e+00   0.00000000e+00
   0.00000000e+00   0.00000000e+00]
[ 0.00078759 -0.00023543  0.          0.          0.          0.        ]
159
[ 0.11494946]
[ 0.01284015  0.00956057]
[ 0.0716714   0.03688733  0.          0.          0.          0.        ]
[ 0.01228238  0.00881147  0.          0.          0.          0.        ]
163
[ 0.02123251]
[ 0.00571172  0.00150443]
[ 0.01981912  0.00085953  0.          0.          0.          0.        ]
[ 0.00562395  0.0011712   0.          0.          0.          0.        ]
165
[ 0.02054149]
[-0.00296114 -0.00446851]
[ 0.00550885  0.01394021  0.          0.          0.          0.        ]
[-0.00270045 -0.00429576  0.          0.          0.          0.        ]










    
































In [ ]:



    






    
















    


    

        

            

            

                Content source: sophie63/FlyLFM
            

            
Similar notebooks:

            

            
                notebook.community |
                gallery |
                about