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
#from nifti import NiftiImage
import nibabel as nb
from scipy.interpolate import interp1d
from scipy import ndimage

Open data





In [2]:



    


# 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/40x/100660series/ForOdors/100659ss1regdFF20sconcatkf219Smith0_4_60TS.mat

















In [3]:



    


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



    


















    
Out[3]:








(7302, 219)

















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



    


















    






/media/sophie/008C0665790F0763/ComboPanNeuronalGCaMP6/40x/100660series/ForOdors/100659ss1regdFF20sconcatkf219Smith0_4_60IC.nii

















In [5]:



    


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



    


















    
Out[5]:








(104, 103, 36, 219)

















In [6]:



    


S=data.shape
S



    


















    
Out[6]:








(104, 103, 36, 219)

Z-score





In [7]:



    


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



    


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



    











In [9]:



    


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

Open Masks





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
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/sophie/008C0665790F0763/ComboPanNeuronalGCaMP6/40x/100660series/100660Registration/JFRC100660TransformedseparateLargetrimmed.nii

















In [14]:



    


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



    


RegionName



    


















    
Out[15]:








['OL',
 'VLNP',
 'VMNP',
 'AL',
 'MB',
 'LH',
 'SNP',
 'CX',
 'LX',
 'INP',
 'PENP',
 'GNG']

Average in masks to sort components by brain region





In [16]:



    


Dmaps.shape



    


















    
Out[16]:








(104, 103, 36, 219)

















In [17]:



    


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



    











In [18]:



    


for i in range(S[3]):
    Mapmean[i]=np.mean(np.mean(np.mean(Dmaps[:,:,:,i])))
    for j in range(12):
        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])



    











In [19]:



    


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



    


















    






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

















In [20]:



    


J



    


















    
Out[20]:








[11]

















In [21]:



    


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

h=5
tot=0
GoodICAnat=np.zeros(S[3])

for l in range(12):
    Final_maps=np.zeros((S[0],S[1],3))
    Fmap=np.zeros((S[0],S[1],3))
    C=np.zeros(3)

    n=0
    for i in range(len(CompMainName)):                    
        Dmmv=np.mean(data[:,:,:,i],2) 
        Dmmv[Dmmv<0.2*np.max(np.max(np.max(Dmmv)))]=0
        C=np.squeeze(np.random.rand(3,1))
        labeled, nrobject=ndimage.label(Dmmv>0)
        
        if CompMainName[i]==Names[l][0] and (sum(CompNameAdd[i,:])<5) and nrobject<200:
            n=n+1            
            
            for k in range(3):
                Fmap[:,:,k]=0.7*Dmmv*C[k]/np.max(C)
            Final_maps=Final_maps+Fmap
            #plt.plot(Time_fluoICA.T,(DT[:,i]/np.sqrt(np.var(DT[:,i]))-h*n+2),color=C/2)
            plt.plot((DT[:,i]/np.sqrt(np.var(DT[:,i]))-h*n+2),color=C/2)
            tot=tot+1
            GoodICAnat[i]=1
            print(i+1)
            
                    
    if n!=0:
        print(RegionName[l])
        plt.show()
        FM=Final_maps/np.max(np.max(Final_maps))
        FM[FM<0.1]=0
        plt.imshow(FM,interpolation='none')
        plt.show()
        frame1 = plt.gca()
        frame1.axes.get_xaxis().set_visible(False)
        frame1.axes.get_yaxis().set_visible(False)



    


















    






2
3
7
11
13
14
19
21
22
29
32
38
42
53
55
56
60
63
65
68
82
89
90
92
123
125
135
141
166
175
OL










    

















    

















    






17
34
118
140
162
194
VLNP










    

















    

















    






52
85
87
93
102
114
147
158
VMNP










    

















    

















    






27
41
48
57
59
66
69
77
84
99
100
104
127
130
131
145
160
167
169
179
187
192
195
214
AL










    

















    

















    






26
35
62
74
76
98
113
142
165
MB










    

















    

















    






16
40
44
72
101
106
170
184
LH










    

















    

















    






1
80
83
95
109
111
117
129
134
143
150
153
156
168
174
203
SNP










    

















    

















    






9
20
28
45
54
71
73
79
103
186
CX










    

















    

















    






50
159
190
INP










    

















    

















    






24
25
30
37
64
67
88
96
97
112
128
154
163
191
PENP










    

















    

















    






4
5
6
8
10
12
15
18
23
31
33
36
39
43
46
47
49
51
58
61
70
75
78
81
91
94
105
107
108
110
115
119
124
126
137
138
139
151
178
GNG
Looked at the components maps and time series and remove all the components which are localized on the edge of the brain and with activity unlike GCaMP6 transients.




In [19]:



    


BadICs=[178,115,87,138,5,18,58,125,154,22,21,10,26,134]



    











In [20]:



    


for idx in BadICs:
    GoodICAnat[idx] = 0.0



    











In [21]:



    


pylab.rcParams['figure.figsize'] = (13, 3)

h=5
tot=0
NumberInLargeRegion=np.zeros(13)

for l in range(12):
    Final_maps=np.zeros((S[0],S[1],3))
    Fmap=np.zeros((S[0],S[1],3))
    C=np.zeros(3)

    n=0
    for i in range(len(CompMainName)):                    
        Dmmv=np.mean(data[:,:,:,i],2) 
        Dmmv[Dmmv<0.2*np.max(np.max(np.max(Dmmv)))]=0
        C=np.squeeze(np.random.rand(3,1))
        labeled, nrobject=ndimage.label(Dmmv>0)
        
        if CompMainName[i]==Names[l][0] and (sum(CompNameAdd[i,:])<5) and nrobject<200 and GoodICAnat[i]==1:
            n=n+1            
            
            for k in range(3):
                Fmap[:,:,k]=0.7*Dmmv*(C[k]+0.2)/np.max(C+0.2)
            Final_maps=Final_maps+Fmap
            #plt.plot(Time_fluoICA.T,(DT[:,i]/np.sqrt(np.var(DT[:,i]))-h*n+2),color=C/2)
            plt.plot((DT[:,i]/np.sqrt(np.var(DT[:,i]))-h*n+2),color=C/2)
            tot=tot+1
            GoodICAnat[i]=1
            print(i)
                    
    if n!=0:
        print(RegionName[l])
        plt.show()
        FM=Final_maps/np.max(np.max(Final_maps))
        FM[FM<0.1]=0
        plt.imshow(FM,interpolation='none')
        plt.show()
        frame1 = plt.gca()
        frame1.axes.get_xaxis().set_visible(False)
        frame1.axes.get_yaxis().set_visible(False)
                
    NumberInLargeRegion[l]=n



    


















    






29
38
54
59
61
69
81
86
88
89
91
93
95
97
98
101
104
106
108
114
116
122
123
126
130
147
148
150
151
156
158
169
171
172
173
177
179
182
190
196
206
212
OL










    

















    

















    






9
14
33
37
44
67
110
128
192
205
VLNP










    

















    

















    






120
129
145
175
176
210
216
VMNP










    

















    

















    






0
1
2
4
6
11
12
17
32
35
39
41
53
AL










    

















    

















    






8
15
16
19
25
27
28
48
57
60
68
84
131
133
159
168
181
186
MB










    

















    

















    






3
13
23
51
56
71
73
90
92
94
121
135
195
227
LH










    

















    

















    






64
100
112
124
139
144
146
166
191
193
194
198
200
213
219
221
229
SNP










    

















    

















    






46
47
62
65
75
76
78
79
109
111
118
127
142
143
165
180
207
220
CX










    

















    

















    






107
167
211
LX










    

















    

















    






36
49
136
201
225
INP










    

















    

















    






7
20
63
70
99
140
149
184
PENP










    

















    

















    






24
31
34
45
77
83
105
117
119
141
152
160
161
162
164
170
174
188
197
199
208
223
226
GNG










    

















    

















    
























In [22]:



    


InterestingICs=[44,72,106,41,48,59,69,77,145]



    











In [23]:



    


pylab.rcParams['figure.figsize'] = (13, 3)

h=5
tot=0
NumberInLargeRegion=np.zeros(13)

for l in range(12):
    Final_maps=np.zeros((S[0],S[1],3))
    Fmap=np.zeros((S[0],S[1],3))
    C=np.zeros(3)

    n=0
    for i in range(len(CompMainName)):                    
        Dmmv=np.mean(data[:,:,:,i],2) 
        Dmmv[Dmmv<0.2*np.max(np.max(np.max(Dmmv)))]=0
        C=np.squeeze(np.random.rand(3,1))
        labeled, nrobject=ndimage.label(Dmmv>0)
        
        if i in InterestingICs and CompMainName[i]==Names[l][0]:
            n=n+1            
            
            for k in range(3):
                Fmap[:,:,k]=0.7*Dmmv*(C[k]+0.2)/np.max(C+0.2)
            Final_maps=Final_maps+Fmap
            #plt.plot(Time_fluoICA.T,(DT[:,i]/np.sqrt(np.var(DT[:,i]))-h*n+2),color=C/2)
            plt.plot((DT[:,i]/np.sqrt(np.var(DT[:,i]))-h*n+2),color=C/2)
            tot=tot+1
            GoodICAnat[i]=1
            print(i)
                    
    if n!=0:
        print(RegionName[l])
        plt.show()
        FM=Final_maps/np.max(np.max(Final_maps))
        FM[FM<0.1]=0
        plt.imshow(FM,interpolation='none')
        plt.show()
        frame1 = plt.gca()
        frame1.axes.get_xaxis().set_visible(False)
        frame1.axes.get_yaxis().set_visible(False)
                
    NumberInLargeRegion[l]=n



    


















    






41
59
OL










    

















    

















    






44
72
CX










    

















    

















    






48
69
77
106
145
GNG










    

















    

















    
























In [24]:



    


# 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)
Ua=sio.loadmat(filename)
Xk=Ua['Xk']



    


















    






()










    






---------------------------------------------------------------------------
IOError                                   Traceback (most recent call last)
<ipython-input-24-33e37e50fa51> in <module>()
      3 filename = askopenfilename() # show an "Open" dialog box and return the path to the selected file
      4 print(filename)
----> 5 Ua=sio.loadmat(filename)
      6 Xk=Ua['Xk']

/usr/local/lib/python2.7/dist-packages/scipy/io/matlab/mio.pyc in loadmat(file_name, mdict, appendmat, **kwargs)
    133     """
    134     variable_names = kwargs.pop('variable_names', None)
--> 135     MR = mat_reader_factory(file_name, appendmat, **kwargs)
    136     matfile_dict = MR.get_variables(variable_names)
    137     if mdict is not None:

/usr/local/lib/python2.7/dist-packages/scipy/io/matlab/mio.pyc in mat_reader_factory(file_name, appendmat, **kwargs)
     56        type detected in `filename`.
     57     """
---> 58     byte_stream = _open_file(file_name, appendmat)
     59     mjv, mnv = get_matfile_version(byte_stream)
     60     if mjv == 0:

/usr/local/lib/python2.7/dist-packages/scipy/io/matlab/mio.pyc in _open_file(file_like, appendmat)
     34         file_like.read(0)
     35     except AttributeError:
---> 36         raise IOError('Reader needs file name or open file-like object')
     37     return file_like
     38 

IOError: Reader needs file name or open file-like object
















In [83]:



    


# 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



    


















    






/media/sophie/008C0665790F0763/ComboPanNeuronalGCaMP6/40x/100670series/100670Registration/ICandAVG_100670regpsf.nii










    
Out[83]:








(104, 103, 10)

















In [84]:



    


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)



    


















    
Out[84]:








<matplotlib.image.AxesImage at 0x7f3739b1edd0>










    
























In [64]:



    


from sklearn import linear_model



    











In [36]:



    


my_cmap=plt.cm.jet
my_cmap.set_bad(alpha=0)



    











In [65]:



    


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








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










    
























In [75]:



    


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

X.shape

DT.shape
DT[DT<0]=0

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)
plt.plot(X)



    


















    
Out[75]:








[<matplotlib.lines.Line2D at 0x7fc227e27110>,
 <matplotlib.lines.Line2D at 0x7fc226fdcd90>]










    
























In [76]:



    


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








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










    
























In [31]:



    


import random



    











In [94]:



    


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



    


C=np.zeros((S[3],3))
i=0
l=0
#Betas2=Betas
LightNuminRegion=np.zeros(12)
for j in range(S[3]):  
    if j in InterestingICs:
    #if 1>0.1:
        C[j,:]=C1[i%6][:]
       
        
        #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.9*D2[:,:,:,j]*C[j,k]/M
        Final_map=Final_map+Fmaps
        #Betas[0,j]=0
        #print(Indexo[j])
        print(j+1)
        print(RegionName[int(CompMainName[j])-1])     
        LightNuminRegion[int(CompMainName[j])-1]=LightNuminRegion[int(CompMainName[j])-1]+1
        i=i+1
        l=l+1

        #if l==2:
            #break



    


















    






16
MB
22
LH
57
LH
95
LH

















In [96]:



    


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



    


C=np.zeros((S[3],3))
i=0
l=0
Betas2=Betas
OdorNuminRegion=np.zeros(12)

for j in range(S[3]):  
    if Betas2[1,j]>0.5*np.max(Betas2[1,:]) and Betas2[0,j]<0.5*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.25*D2[:,:,:,j]*C[j,k]/M
        Final_map=Final_map+Fmaps
        #Betas2[1,j]=0
        #print(Indexo[j])
        OdorNuminRegion[int(CompMainName[j])-1]=OdorNuminRegion[int(CompMainName[j])-1]+1
        print(RegionName[int(CompMainName[j])-1])
        i=i+1
        l=l+1
        print(j+1)
        #if l==2:
         #   break

NumOdor=i
print('Number of odor components')
print(i)



    


















    






AL
5
PENP
8
VLNP
15
MB
16
Number of odor components
4

















In [135]:



    


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



    


np.savetxt('/'.join(filename.split('/')[:-1])+'/OdorNumberInLargeRegions.txt',OdorNuminRegion)
np.savetxt('/'.join(filename.split('/')[:-1])+'/LightNumberInLargeRegions.txt',LightNuminRegion)



    











In [ ]:



    


plt.plot(OdorNuminRegion)
plt.plot(LightNuminRegion)



    











In [ ]:

Content source: sophie63/FlyLFM

Similar notebooks:

notebook.community | gallery | about