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/960/960ss1_500regcdFF20spsfkf192Smith0_4_60TS.mat

















In [3]:



    


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



    


















    
Out[3]:








(5310, 192)

















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/960/960ss1_500regcdFF20spsfkf192Smith0_4_60IC.nii

















In [5]:



    


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



    


















    
Out[5]:








(88, 77, 10, 192)

















In [6]:



    


S=data.shape
S



    


















    
Out[6]:








(88, 77, 10, 192)

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



    


# 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/960/960registration/JFRCTransformedLargefullpsftrimmed.nii

















In [11]:



    


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



    


RegionName



    


















    
Out[12]:








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

Average in masks to sort components by brain region





In [13]:



    


Dmaps.shape



    


















    
Out[13]:








(88, 77, 10, 192)

















In [14]:



    


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



    











In [15]:



    


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



    


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



    


J



    


















    
Out[17]:








[8]

















In [27]:



    


pylab.rcParams['figure.figsize'] = (13, 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
            
                    
    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)



    


















    






OL










    

















    

















    






VLNP










    

















    

















    






VMNP










    

















    

















    






AL










    

















    

















    






MB










    

















    

















    






SNP










    

















    

















    






CX










    

















    

















    






INP










    

















    

















    






PENP










    

















    

















    






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



    


BadICs=[178,115,87,138,5,18,58,125,154,22,21,10,26,134,136,23,50,7,108]



    











In [23]:



    


for idx in BadICs:
    GoodICAnat[idx] = 0.0



    











In [24]:



    


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

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



    


















    






66
VLNP










    

















    

















    






49
61
70
79
88
121
141
147
VMNP










    

















    

















    






2
4
13
15
16
17
20
24
25
27
30
33
37
43
46
51
52
54
55
60
62
63
65
68
75
77
78
94
98
102
107
AL










    

















    

















    






0
3
11
12
14
31
41
42
48
57
59
69
76
151
166
MB










    

















    

















    






8
9
19
29
32
34
40
47
56
64
71
80
82
83
86
93
99
100
101
103
105
110
112
113
114
116
122
128
129
135
142
153
159
180
181
SNP










    

















    

















    






73
85
89
90
91
92
95
104
111
117
119
120
126
131
132
137
140
144
149
157
161
163
165
172
176
177
189
CX










    

















    

















    






1
6
28
36
45
67
81
156
179
INP










    

















    

















    






44
109
127
152
PENP










    

















    

















    
























In [26]:



    


# Output number of component per region
np.savetxt('/'.join(filename.split('/')[:-1])+'/NumberInLargeRegionsV2.txt',NumberInLargeRegion)



    











In [ ]:

Content source: sophie63/FlyLFM

Similar notebooks:

notebook.community | gallery | about