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/test15/ArclightCombo/100632series/100633_100634/100633cregcdFF20spsfkf-50toendconcat202Smith0_4_60TS.mat

















In [3]:



    


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



    


















    
Out[3]:








(12079, 202)

















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/test15/ArclightCombo/100632series/100633_100634/100633cregcdFF20spsfkf-50toendconcat202Smith0_4_60IC.nii

















In [5]:



    


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



    


















    
Out[5]:








(78, 55, 10, 202)

















In [6]:



    


S=data.shape
S



    


















    
Out[6]:








(78, 55, 10, 202)

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



    


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



    


filenameM ="/media/test15/ArclightCombo/100632series/100632seriesRegistration/JFRCTransformedLargefullpsftrimmed.nii"



    











In [13]:



    


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



    











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








(78, 55, 10, 202)

















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








[0]

















In [22]:



    


# 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/test15/ArclightCombo/100632series/100633_100634Behavior.mat

















In [23]:



    


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



    











In [27]:



    


Xk.shape



    


















    
Out[27]:








(3, 12079)

















In [31]:



    


Xk=Xk.T



    











In [32]:



    


plt.plot(Xk[:,0]/np.std(Xk[:,0])+2,color=(1,0,0))   
plt.plot(Xk[:,1]/np.std(Xk[:,1])+3,color=(0.5,0.5,0))
plt.plot(Xk[:,2]/np.std(Xk[:,2])+6,color=(1,0,1))



    


















    
Out[32]:








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










    
























In [39]:



    


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

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<50:
            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)
            
    plt.plot(Xk[:,0]/np.std(Xk[:,0])+2,color=(1,0,0))   
    plt.plot(Xk[:,1]/np.std(Xk[:,1])+3,color=(0.5,0.5,0))
    plt.plot(Xk[:,2]/np.std(Xk[:,2])+6,color=(1,0,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)



    


















    






25
46
70
84
96
100
117
191
196
197
198
201
OL










    

















    

















    






28
174
VLNP










    

















    

















    






20
40
45
122
VMNP










    

















    

















    






4
9
26
85
89
98
173
AL










    

















    

















    






14
15
18
31
62
139
MB










    

















    

















    






21
38
53
54
55
77
79
124
LH










    

















    

















    






16
27
66
142
159
165
192
SNP










    

















    

















    






88
93
144
148
150
169
179
CX










    

















    

















    






2
6
32
143
INP










    

















    

















    






7
68
PENP










    

















    

















    






10
11
23
24
29
33
42
47
49
50
56
57
59
61
64
65
69
74
75
80
86
87
90
91
92
94
97
99
102
103
106
107
108
110
111
114
115
118
120
121
125
127
130
133
134
135
137
138
141
146
149
154
158
160
162
176
182
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 [24]:



    


BadICs=[]



    











In [25]:



    


for idx in BadICs:
    GoodICAnat[idx] = 0.0



    











In [28]:



    


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

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



    


















    






6
8
23
25
33
40
54
82
114
132
133
153
194
221
OL










    

















    

















    






20
46
92
96
200
VLNP










    

















    

















    






41
57
58
66
70
76
79
91
99
100
146
151
203
VMNP










    

















    

















    






12
14
44
53
81
98
188
AL










    

















    

















    






1
5
7
9
31
42
MB










    

















    

















    






13
21
22
24
30
45
56
59
63
85
104
111
LH










    

















    

















    






19
26
28
61
72
87
97
112
127
137
138
148
168
173
179
212
SNP










    

















    

















    






15
109
120
126
163
164
169
210
CX










    

















    

















    






2
3
18
36
51
139
165
INP










    

















    

















    






4
PENP










    

















    

















    
























In [27]:



    


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



    











In [ ]:



    






    











In [ ]:

Content source: sophie63/FlyLFM

Similar notebooks:

notebook.community | gallery | about