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



    


# 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/008C0665790F0763/ComboPanNeuronalGCaMP6/862_863/180315/863ss2onc250_862ss2onc250zx2cregcregc862_863dFF20spsfkf221Smith0_4_60TS.mat

















In [3]:



    


filename = "/media/test11/ComboPanNeuronalGCaMP6/862_863/180315/863ss2onc250_862ss2onc250zx2cregcregc862_863dFF20spsfkf221Smith0_4_60TS.mat"



    











In [4]:



    


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



    


















    
Out[4]:








(5957, 221)

















In [6]:



    


# 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/sophie/008C0665790F0763/ComboPanNeuronalGCaMP6/862_863/180315/863ss2onc250_862ss2onc250zx2cregcregc862_863dFF20spsfkf221Smith0_4_60IC.nii

















In [7]:



    


filename2 ="/media/test11/ComboPanNeuronalGCaMP6/862_863/180315/863ss2onc250_862ss2onc250zx2cregcregc862_863dFF20spsfkf221Smith0_4_60IC.nii"



    











In [8]:



    


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



    


















    
Out[8]:








(174, 93, 12, 221)

Z-score





In [9]:



    


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



    


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



    











In [11]:



    


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



    











In [ ]:



    


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



    


filename ="/media/test11/ComboPanNeuronalGCaMP6/862_863/XkStim.mat"



    











In [13]:



    


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



    











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



    


















    






/media/test11/ComboPanNeuronalGCaMP6/862_863/862_863registration/Template4alignement2anatomypsf.nii

















In [ ]:



    






    











In [15]:



    


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



    


















    
Out[15]:








(174, 93, 12)

Fit odor and light onsets





In [16]:



    


%%javascript
IPython.OutputArea.auto_scroll_threshold =4000;



    


















    



























In [17]:



    


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

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

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



    


















    






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










    
Out[17]:








<matplotlib.image.AxesImage at 0x7f03d8d09290>










    
























In [18]:



    


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)

algorithm = linear_model.LinearRegression()

Sxk=Xk.shape

Sxk

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

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

plt.plot(X[:,0])
plt.plot(X[:,1])



    


















    
Out[18]:








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

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





In [19]:



    


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

X.shape

DT.shape

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

plt.plot(DT[3850:4000,6])
plt.plot(X[3850:3900,1])



    


















    
Out[19]:








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










    
























In [20]:



    


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

Sx=X.shape

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

plt.plot(Betas[0,:])



    


















    
Out[20]:








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










    
























In [21]:



    


import random



    











In [22]:



    


del Final_map
del Fmaps



    


















    






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

NameError: name 'Final_map' is not defined
















In [23]:



    


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



    


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



    


















    






16
[ 0.03566413]
18
[ 0.03321012]
24
[ 0.1061012]
29
[ 0.02085529]
30
[ 0.01454197]
35
[ 0.01383488]
42
[ 0.02876138]
45
[ 0.0262224]
46
[ 0.09780365]
48
[ 0.01576255]
49
[ 0.00318861]
53
[ 0.04224128]
55
[ 0.04522537]
58
[ 0.0160867]
60
[ 0.08700904]
61
[ 0.02788882]
64
[ 0.00624137]
67
[ 0.04327817]
68
[ 0.03128859]
74
[ 0.0804405]
79
[ 0.0151209]
89
[ 0.0429691]
90
[ 0.03182817]
91
[ 0.01017075]
96
[ 0.00685366]
102
[ 0.04388732]
106
[ 0.00734878]
108
[ 0.04550682]
113
[ 0.01847885]
116
[ 0.0526088]
119
[ 0.04219017]
122
[ 0.02006807]
123
[ 0.04631941]
125
[ 0.03147026]
129
[ 0.05887153]
131
[ 0.04070635]
133
[ 0.03544915]
136
[ 0.01499375]
138
[ 0.01563543]
139
[ 0.02207778]
146
[ 0.02114218]
148
[ 0.02149185]
153
[ 0.05033412]
154
[ 0.04130495]
155
[ 0.01910912]
159
[ 0.04466967]
160
[ 0.04637926]
162
[ 0.05023123]
167
[ 0.04389098]
169
[ 0.04670044]
170
[ 0.07089345]
174
[ 0.03719359]
190
[ 0.02899371]
213
[ 0.02709707]

















In [25]:



    


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/200
    #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 [26]:



    


C=np.zeros((S[3],3))
i=0
l=0
Betas2=Betas
for j in range(S[3]):  
    if Betas2[1,j]>0.1*np.max(Betas2[1,:]) and abs(Betas2[0,j])<0.1*np.max(Betas2[0,:]):
    #if 1>0.1:
        #C[j,:]=C1[i%6][:]
        C[j,0]=1
        C[j,1]=Betas2[1,j]/np.max(Betas2[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.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)
        #if l==2:
         #   break
NumOdor=i
print('Number of odor components')
print(i)



    


















    






1
2
3
4
6
7
8
9
10
11
15
19
20
21
22
23
26
40
57
76
78
87
130
Number of odor components
23

















In [27]:



    


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/300
    #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 [28]:



    


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=1
l=0
Betas2=Betas
for j in range(2,S[3]):  
    if Betas2[0,j]>0.8*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)
        #if l==2:
            #break
    

C=np.zeros((S[3],3))
i=1
l=0
Betas2=Betas
for j in range(2,S[3]):  
    if Betas2[1,j]>0.8*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)
        #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)



    


















    






24










    
























In [29]:



    


Xthresh=abs(np.diff(X))-6
Xthresh[Xthresh<0]=0
Xthresh[Xthresh>0]=1
plt.plot(Xthresh)



    


















    
Out[29]:








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










    
























In [31]:



    


Xthresh.shape

AvLightIC=np.zeros(DT.shape[1])        
for j in range(S[3]):
    k=0
    for i in range(3000):
        if Xthresh[i,0]==1:
            AvLightIC[j]=AvLightIC[j]+(np.mean(DT[i:i+200,j])-np.mean(DT[i-200:i,j]))/(2*np.std(DT[i-200:i,j]))
            k=k+1
    AvLightIC[j]=AvLightIC[j]/k

AvOdorIC=np.zeros(DT.shape[1])        
for j in range(S[3]):
    k=0
    for i in range(3000,5500):
        if Xthresh[i,0]==1:
            AvOdorIC[j]=AvOdorIC[j]+(np.mean(DT[i:i+200,j])-np.mean(DT[i-200:i,j]))/(2*np.std(DT[i-200:i,j]))
            k=k+1
    AvOdorIC[j]=AvOdorIC[j]/k



    











In [33]:



    


filenameM='/home/sophie/LargeRegionList'
with open(filenameM) as f:
    content = f.readlines()
Names=[Line.replace('\n','').split(' ') for Line in content]
RegionName=[Names[i][1] for i in range(12)]
Num=[int(Names[i][0]) for i in range(12)]



    











In [34]:



    


C=np.zeros((S[3],3))
i=0
l=0

LightNuminRegionAv=np.zeros(12)
for j in range(S[3]):  
    #if AvLightIC[j]>1 and GoodICAnat[j]:
    if AvLightIC[j]>1 and AvOdorIC[j]>1:
        #C[j,:]=C1[i%6][:]
        C[j,2]=1
        C[j,1]=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.5*D2[:,:,:,j]*C[j,k]/M
        Final_map=Final_map+Fmaps
        #Betas[0,j]=0
        #print(Indexo[j])
        print(j+1)
        #print(CompMainName[j]) 
        print(RegionName[int(CompMainName[j])-1])
        
        i=i+1
        l=l+1
        LightNuminRegionAv[int(CompMainName[j])-1]=LightNuminRegionAv[int(CompMainName[j])-1]+1
        #if l==2:
            #break



    


















    






8










    






---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-34-c6e3f9b9a1f8> in <module>()
     19         print(j+1)
     20         #print(CompMainName[j])
---> 21         print(RegionName[int(CompMainName[j])-1])
     22 
     23         i=i+1

NameError: name 'CompMainName' is not defined
















In [ ]:

Content source: sophie63/FlyLFM

Similar notebooks:

notebook.community | gallery | about