In [1]:
clear all
In [2]:
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
import scipy.ndimage
In [3]:
# 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
filename = askopenfilename() # show an "Open" dialog box and return the path to the selected file
print(filename)
In [4]:
Ua=sio.loadmat(filename)
DT=Ua['TSo']
DT.shape
Out[4]:
In [5]:
# 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
filename2 = askopenfilename() # show an "Open" dialog box and return the path to the selected file
print(filename2)
In [6]:
img1 = nb.load(filename2)
data = img1.get_data()
S=data.shape
S
Out[6]:
In [7]:
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 [8]:
#Old version from csv
#with open(filename) as inputfile:
# results = list(csv.reader(inputfile))
#Raw_timea=[float(results[i][2]) for i in (range(1,len(results)))]
#Raw_time=np.asarray(Raw_timea[0:len(results)])
#Raw_time.shape
#on=126
#off=6253
#Time_fluo=Raw_time[on:off]-Raw_time[on]
Ua=sio.loadmat(filename)
Time_fluo=Ua['TimeFluoOn']
Time_fluo.shape
Out[8]:
In [56]:
Time_fluoICA=Time_fluo[:,250:7213]
In [57]:
Time_fluoICA.shape
500
Out[57]:
In [58]:
#if stimulus
#Tstim=np.array(range(1,100*(int(Tmax)+1)))
#Tstim=Tstim/100
#Odor=np.zeros(len(Tstim*100))
#UV=np.zeros(len(Tstim*100))
#for i in range(0,7):
# UV[10*100+i*600:100*10+i*600+300]=1
#
#for i in range(0,5):
# Odor[(10+8*6+15)*100+i*600:(10+8*6+15)*100+i*600+200]=1
#plt.plot(Tstim,UV)
#plt.plot(Tstim,Odor)
In [80]:
# 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)
In [82]:
filenameM='/home/sophie/Desktop/Registration/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)]
In [552]:
M=np.zeros((S[3],86))
Mby2=np.zeros((S[3],86))
Mask of outside of the brain
In [553]:
MasksSum=np.sum(Masks[:,:,:,2:len(Masks)], axis=3)
In [554]:
MasksSum[MasksSum==0]=1000
MasksSum[MasksSum!=1000]=0
MasksSum[MasksSum==1000]=1
MaskOutside=MasksSum
In [555]:
MasksSumOL=np.zeros(MasksSum.shape)
for l in range(75):
if Names[l][0]=='AME_R' or Names[l][0]=='AME_L' or Names[l][0]=='ME_R' or Names[l][0]=='ME_L' or Names[l][0]=='LO_R' or Names[l][0]=='LO_L' or Names[l][0]=='LOP_R' or Names[l][0]=='LOP_L':
MasksSumOL=MasksSumOL+Masks[:,:,:,Num[l]]
In [391]:
MasksPBMean=np.mean(MasksSumOL,2)
plt.imshow(MasksPBMean,cmap=plt.cm.gray)
Out[391]:
In [556]:
Masks[:,:,:,0]=MaskOutside
In [557]:
del Masksby2
In [558]:
97/2
Out[558]:
In [559]:
Masksby2=np.zeros((np.floor(S[0]/2.0),np.ceil(S[1]/2.0),S[2],87))
Databy2=np.zeros((np.ceil(S[0]/2.0),np.ceil(S[1]/2.0),S[2],S[3]))
for i in range(S[2]):
for j in range(87):
Masksby2[:,:,i,j]=scipy.ndimage.zoom(Masks[:,:,i,j], 0.5, order=0)
for j in range(S[3]):
Databy2[:,:,i,j]=scipy.ndimage.zoom(data[:,:,i,j], 0.5, order=0)
In [560]:
for i in range(S[3]):
for j in range(1,86):
MMasks=np.mean(np.mean(np.mean(Masks[:,:,:,j])))
if MMasks:
M[i,j]=np.mean(np.mean(np.mean(Masks[:,:,:,j]*Dmaps[:,:,:,i])))/MMasks
In [561]:
plt.plot(M[20,:])
Out[561]:
In [562]:
np.std(M[350,:])
Out[562]:
In [563]:
CompName=S[3]*['']
for i in range(S[3]):
Max=np.max(Mby2[i,:])
I=np.argmax(Mby2[i,:])
if Max>1.2:
J=[l for l in range(75) if Num[l]==I]
CompName[i]=Names[np.array(J)][0]
In [564]:
h=1
tot=0
GoodICAnat=np.zeros(S[3])
for l in range(75):
Final_maps=np.zeros((S[0],S[1],3))
Fmap=np.zeros((S[0],S[1],3))
C=np.zeros(3)
print(Names[l])
n=0
for i in range(len(CompName)):
if CompName[i]==Names[l][0]:
n=n+1
Dm=np.mean(data[:,:,:,i],2)
C=np.squeeze(np.random.rand(3,1))
#C=(1,1,1)
for k in range(3):
#Maximum=np.max(np.squeeze(np.reshape(Dm[:,:],S[0]*S[1])))
#Dmm=Dm-m
V=np.var(np.reshape(Dm,S[0]*S[1]))
Dmmv=Dm/V
#Dmmv[Dmmv<2]=0
#Dmmv=Dmmv-2
#M[M==0]=1
# Fmaps[:,:,:,k]=0.7*Dmaps[:,:,:,j]*C[j,k]/np.max(np.squeeze(np.reshape(Dmaps[:,:,:,j]*np.max(C[j,:],S[0]*S[1]*5))
Fmap[:,:,k]=0.7*Dmmv*C[k]/np.max(C)
Final_maps=Final_maps+Fmap
plt.plot(Time_fluoICA.T,(DT[:,i]+h*n+2),color=C/2)
plt.plot(Tf,straightintbin*1.5)
plt.plot(Tf,leftintbin*1.5)
plt.plot(Tf,rightintbin*1.5)
plt.plot(Tf,FrontGroomint.T)
plt.plot(Tf,RearGroomint)
plt.plot(Tf,Walkint)
plt.plot(Tf,Weirdint)
tot=tot+1
GoodICAnat[i]=1
plt.show()
FM=Final_maps/np.max(np.max(Final_maps))
FM[FM<0.1]=0
plt.imshow(FM,interpolation='none')
plt.show()
In [488]:
tot
Out[488]:
In [204]:
np.max(np.max(Dm))
Out[204]:
In [191]:
C
Out[191]:
In [157]:
plt.imshow(np.mean(data[:,:,:,11],2),cmap=plt.cm.gray)
Out[157]:
In [114]:
np.array(Num).shape
Out[114]:
In [12]:
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)
B=sio.loadmat(filename)
In [13]:
B
Out[13]:
In [14]:
b=B['behavior']
In [15]:
FrontGroom=b[0,:]
RearGroom=b[1,:]
Rest=b[3,:]
Tbehavior=B['T']
Walk=b[2,:]
Weird=b[4,:]
In [16]:
Weird.shape
Out[16]:
In [17]:
#Threshold ball movements
In [19]:
left=B['Left'].T
right=B['Right'].T
straight=B['Straight'].T
In [20]:
plt.plot(left)
Out[20]:
In [21]:
max(left)
Out[21]:
In [22]:
#Fill gaps (1 if i 1 and i+4 1)
In [23]:
List=['left','right','straight','FrontGroom','RearGroom','Rest','Walk','Weird']
In [24]:
dataset=eval(List[0])
for i in range(len(List)):
del dataset
dataset=eval(List[i])
for j in range(len(dataset)-3):
if dataset[j]==1 and dataset[j+3]==1:
dataset[j+1]=1
dataset[j+2]=1
exec(List[i]+'=dataset')
In [25]:
Tb=np.squeeze(Tbehavior)
In [66]:
Tfi=Time_fluoICA.T
Tf=Tfi[range(len(Time_fluoICA.T)-5)]
f=interp1d(Tb,FrontGroom)
FrontGroomint=np.squeeze(f(Tf))
f=interp1d(Tb,np.squeeze(RearGroom))
RearGroomint=f(Tf)
f=interp1d(Tb,np.squeeze(Rest))
Restint=f(Tf)
f=interp1d(Tb,np.squeeze(Walk))
Walkint=f(Tf)
f=interp1d(Tb,np.squeeze(Weird))
Weirdint=f(Tf)
In [67]:
Tball=np.zeros(len(left))
for j in range(len(left)):
Tball[j]=(Tb[j]+Tb[j+1])/2
In [68]:
f=interp1d(Tball,left.T)
leftint=np.squeeze(f(Tf))
f=interp1d(Tball,np.squeeze(right))
rightint=f(Tf)
f=interp1d(Tball,np.squeeze(straight))
straightint=f(Tf)
In [69]:
leftintbin=leftint
lowvals=leftint<50
leftintbin[lowvals]=0
highvals=leftint>50
leftintbin[highvals]=1
rightintbin=rightint
lowvals=rightint<50
rightintbin[lowvals]=0
highvals=rightint>50
rightintbin[highvals]=1
straightintbin=straightint
lowvals=straightint<50
straightintbin[lowvals]=0
highvals=straightint>50
straightintbin[highvals]=1
In [70]:
#Align times (verify what frames have been removed)
In [71]:
#interpolate for ball movements
In [72]:
#check where the first image is compared to extrapolated full light
Tb=Tb
max(Tb)
Out[72]:
In [73]:
Tb.shape
Out[73]:
In [74]:
max(np.squeeze(Time_fluoICA))
Out[74]:
In [75]:
#falls between +/- 5ms. note that the offset for the video is about 6 ms later thanlast image
In [76]:
#Do regressions with delay
In [77]:
#Regression integrator straight
In [78]:
#scikit learn GUI
In [39]:
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 [40]:
for i in range(S[3]):
Demean[:,:,:,i]=data[:,:,:,i]-np.mean(np.mean(np.mean(data[:,:,:,i],0),0),0)
In [41]:
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
In [42]:
Order=np.argsort(Var)[::-1]
datao=data[:,:,:,Order[:]]
Dmapso=Dmaps[:,:,:,Order[:]]
Varor=Var[Order]
In [43]:
my_cmap=plt.cm.jet
my_cmap.set_bad(alpha=0)
Good_ICs=np.zeros(S[3])
Label_ICs=[]
pylab.rcParams['figure.figsize'] = (13, 2.5)
In [44]:
Dtemp=np.mean(data,3)
In [45]:
Dtemp.shape
Out[45]:
In [46]:
%%javascript
IPython.OutputArea.auto_scroll_threshold =4000;
In [47]:
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 [48]:
plt.imshow(Dmean[:,:,1],cmap=plt.cm.gray)
Out[48]:
In [49]:
Tb.shape
Out[49]:
In [50]:
S
Out[50]:
In [51]:
RearGroomint.shape
Out[51]:
In [64]:
straightintbin.shape
Out[64]:
In [65]:
Tf.shape
Out[65]:
In [504]:
for j in range(S[3]):
if S[2]>5:
for i in range(Nstack):
V=Dmaps[:,:,Indices[i],Order[j]]
D1[:,:,i]=np.max(V,2)
D2[:,:,:,j]=D1
D1[D1==0]=np.nan
else:
for i in range(S[2]):
V=Dmaps[:,:,i,Order[j]]
D1[:,:,i]=V
for j in range(S[3]):
print(j)
for i in range(Nstack):
plt.subplot(1,5,i+1)
plt.imshow(Dmean[:,:,i],cmap=plt.cm.gray)
plt.imshow(D1[:,:,i], cmap=my_cmap,interpolation='none')
frame1 = plt.gca()
frame1.axes.get_xaxis().set_visible(False)
frame1.axes.get_yaxis().set_visible(False)
plt.show()
plt.plot(Tf,straightintbin*1.5)
plt.plot(Tf,leftintbin*1.5)
plt.plot(Tf,rightintbin*1.5)
plt.plot(Tf,FrontGroomint.T)
plt.plot(Tf,RearGroomint)
plt.plot(Tf,Walkint)
plt.plot(Tf,Weirdint)
plt.plot(Time_fluoICA.T,DT[:,Order[j]]+2)
plt.show()
Label_ICs.append(raw_input())
if Label_ICs[j]=='':
Good_ICs[j]=0
else:
Good_ICs[j]=1
In [806]:
Tk().withdraw()
filenamet = askopenfilename()
print(filenamet)
nimt=NiftiImage(filenamet)
Dtemp=np.squeeze(nimt.data.T)
Dtemp.shape
In [298]:
Dmaps.shape
Out[298]:
In [299]:
set(Label_ICs)
Out[299]:
In [300]:
if len(Label_ICs)<S[3]:
for j in range(S[3]-len(Label_ICs)):
Label_ICs.append('')
In [301]:
Dict={'a':0,'AL':0,'Co':1,'CO':1,'O':1,'o':1,'OG':2,'AVLP':3,'A':0,'G':3,'PN':3,'CA':4,'L':4,'LH':4,'lh':4,'SMP':5,'SLP':5,'l':6,'KC':7,'kc':7,'M':7,'MB':7,'mb':7,'MBON':7,'FB':9,'EB':10,'eb':10,'C':10,'c':10,'PB':9,'pb':9,'AMMC':12,'ammc':12,'V':5,'S':13,'SOG':13,'PS':13,'s':13,'DCB':14,'I':12,'i':12,'D':5,'mvt':15,'':15,'0':15}
In [302]:
Translated=[Dict[X] for X in Label_ICs]
In [303]:
G=Good_ICs.tolist();
In [304]:
len(Good_ICs)
Out[304]:
In [305]:
G.count(1)
Out[305]:
In [505]:
where_are_NaNs = np.isnan(D2)
D2[where_are_NaNs] = 0
In [541]:
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=npLabel_ICs[Neworder[i]]
Fmaps=np.zeros([S[0],S[1],3])
C=np.zeros([S[3],3])
In [542]:
D2.shape
Out[542]:
Inverse of Order
In [565]:
[I,OrderInv]=np.sort([Order,range(len(Order))],axis=0)
In [566]:
for j in range(S[3]):
if GoodICAnat[j]:
#if Good_ICs[j]:
#if Translated[j]==12 or Translated[j]==13:
C[j,:]=np.squeeze(np.random.rand(3,1))
for k in range(3):
M=np.max(np.squeeze(np.reshape(D2[:,:,:,I[j]],S[0]*S[1]*5)))
#M[M==0]=1
# Fmaps[:,:,:,k]=0.7*Dmaps[:,:,:,j]*C[j,k]/np.max(np.squeeze(np.reshape(Dmaps[:,:,:,j]*np.max(C[j,:],S[0]*S[1]*5))
Fmaps[:,:,:,k]=0.7*D2[:,:,:,j]*C[j,k]/(M*np.max(C[j,:]))
Final_map=Final_map+Fmaps
In [567]:
pylab.rcParams['figure.figsize'] = (18, 5)
#Final_map[Final_map<0.1]=np.NaN
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(Dmean[:,:,i],cmap=plt.cm.gray)
plt.imshow(Final_map[:,:,i],interpolation='none')
#plt.imshow(Final_map[:,:,i])
frame1 = plt.gca()
frame1.axes.get_xaxis().set_visible(False)
frame1.axes.get_yaxis().set_visible(False)
In [568]:
pylab.rcParams['figure.figsize'] = (18, 10)
h=0.5
i=0
for j in range(S[3]):
if GoodICAnat[j]:
#if Good_ICs[j]:
#if Translated[j]==12 or Translated[j]==13:
#plt.plot(Time_fluoICA.T,DT[:,Order[j]]+2)
plt.plot(Time_fluoICA.T,(DT[:,Order[j]]+h*i+2),color=C[j,:])
plt.plot(Tf,straightintbin*1.5)
plt.plot(Tf,leftintbin*1.5)
plt.plot(Tf,rightintbin*1.5)
plt.plot(Tf,FrontGroomint.T)
plt.plot(Tf,RearGroomint)
i=i+1
#plt.xlim([10,17])
plt.show()
In [327]:
from sklearn import linear_model
import scipy
Extend linearregression class to get p values for the betas (from http://stackoverflow.com/questions/27928275/find-p-value-significance-in-scikit-learn-linearregression)
In [328]:
class LinearRegression(linear_model.LinearRegression):
def __init__(self,*args,**kwargs):
# *args is the list of arguments that might go into the LinearRegression object
# that we don't know about and don't want to have to deal with. Similarly, **kwargs
# is a dictionary of key words and values that might also need to go into the orginal
# LinearRegression object. We put *args and **kwargs so that we don't have to look
# these up and write them down explicitly here. Nice and easy.
if not "fit_intercept" in kwargs:
kwargs['fit_intercept'] = False
super(LinearRegression,self).__init__(*args,**kwargs)
# Adding in t-statistics for the coefficients.
def fit(self,x,y):
# This takes in numpy arrays (not matrices). Also assumes you are leaving out the column
# of constants.
# Not totally sure what 'super' does here and why you redefine self...
self = super(LinearRegression, self).fit(x,y)
n, k = x.shape
yHat = np.matrix(self.predict(x)).T
# Change X and Y into numpy matricies. x also has a column of ones added to it.
x = np.hstack((np.ones((n,1)),np.matrix(x)))
y = np.matrix(y).T
# Degrees of freedom.
df = float(n-k-1)
# Sample variance.
sse = np.sum(np.square(yHat - y),axis=0)
self.sampleVariance = sse/df
# Sample variance for x.
self.sampleVarianceX = x.T*x
# Covariance Matrix = [(s^2)(X'X)^-1]^0.5. (sqrtm = matrix square root. ugly)
self.covarianceMatrix = scipy.linalg.sqrtm(self.sampleVariance[0,0]*self.sampleVarianceX.I)
# Standard erros for the difference coefficients: the diagonal elements of the covariance matrix.
self.se = self.covarianceMatrix.diagonal()[1:]
# T statistic for each beta.
self.betasTStat = np.zeros(len(self.se))
for i in xrange(len(self.se)):
# self.betasTStat[i] = self.coef_[0,i]/self.se[i] sophie: I get an error here so change the line to the one below:
self.betasTStat[i] = self.coef_[i]/self.se[i]
# P-value for each beta. This is a two sided t-test, since the betas can be
# positive or negative.
self.betasPValue = 1 - scipy.stats.t.cdf(abs(self.betasTStat),df)
In [329]:
clf2=LinearRegression()
In [330]:
DT.shape
Out[330]:
In [331]:
X=np.squeeze(np.array([np.squeeze(straightint),np.squeeze(rightint),np.squeeze(leftint)]))
Nreg=X.shape[0]
Sig=np.zeros([S[3],Nreg])
beta=np.zeros([S[3],Nreg])
for j in range(S[3]):
if Good_ICs[j]:
Out=DT[range(len(Tf)),j]
clf2.fit(X.T,Out)
if clf2.betasPValue[2]<(0.005/(len(Good_ICs)*3))and clf2.coef_[2]>0:
Sig[j,2]=1
beta[j,2]=clf2.betasPValue[2]
if clf2.betasPValue[1]<(0.005/(len(Good_ICs)*3))and clf2.coef_[1]>0:
Sig[j,1]=1
beta[j,1]=clf2.betasPValue[1]
if clf2.betasPValue[0]<(0.005/(len(Good_ICs)*3))and clf2.coef_[0]>0:
Sig[j,0]=1
beta[j,0]=clf2.betasPValue[0]
In [332]:
Nreg
Out[332]:
In [333]:
plt.plot(Sig)
Out[333]:
In [334]:
del Final_map
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=npLabel_ICs[Neworder[i]]
Fmaps=np.zeros([S[0],S[1],3])
C=np.zeros([S[3],3])
In [335]:
for j in range(S[3]):
if ((Sig[j,2] or Sig[j,1] or Sig[j,0]) and Good_ICs[j]):
C[j,:]=Sig[j,:]
for k in range(3):
M=np.max(np.squeeze(np.reshape(D2[:,:,:,j],S[0]*S[1]*5)))
#M[M==0]=1
# Fmaps[:,:,:,k]=0.7*Dmaps[:,:,:,j]*C[j,k]/np.max(np.squeeze(np.reshape(Dmaps[:,:,:,j]*np.max(C[j,:],S[0]*S[1]*5))
Fmaps[:,:,:,k]=0.7*D2[:,:,:,j]*C[j,k]/M
Final_map=Final_map+Fmaps
In [336]:
pylab.rcParams['figure.figsize'] = (14, 5)
#Final_map[Final_map<0.1]=np.NaN
if S[2]>5:
N=Nstack
else:
N=S[2]
for i in range(N):
plt.subplot(1,N,i+1)
#plt.imshow(Dmean[:,:,i],cmap=plt.cm.gray)
plt.imshow(Final_map[:,:,i],interpolation='none')
#plt.imshow(Final_map[:,:,i])
frame1 = plt.gca()
frame1.axes.get_xaxis().set_visible(False)
frame1.axes.get_yaxis().set_visible(False)
In [337]:
pylab.rcParams['figure.figsize'] = (7, 10)
h=1
i=0
for j in range(S[3]):
if (Sig[j,2] or Sig[j,1] or Sig[j,0]) and Good_ICs[j]:
plt.plot(Tball,straightbin*1.5)
plt.plot(Tball,leftbin*1.5)
plt.plot(Tball,rightbin*1.5)
plt.plot(Tb,FrontGroom.T)
plt.plot(Tb,RearGroom.T)
#plt.plot(Time_fluoICA.T,DT[:,Order[j]]+2)
plt.plot(Time_fluoICA.T,(DT[:,Order[j]]+h*i+2),color=C[j,:]/2)
i=i+1
plt.show()
In [338]:
Groom=FrontGroomint+np.squeeze(RearGroomint)
Groom=np.ceil(Groom/2)
In [339]:
Walk=np.squeeze(straightint)+np.squeeze(rightint)
Walk=np.ceil(Walk/2)
Walk=Walk+leftint
Walk=np.ceil(Walk/2)
In [340]:
len(Good_ICs)
Out[340]:
In [341]:
X=np.squeeze(np.array([np.squeeze(Walk),np.squeeze(Groom),np.squeeze(Restint)]))
Nreg=X.shape[0]
Sig=np.zeros([S[3],Nreg])
for j in range(S[3]):
if Good_ICs[j]==1:
Out=DT[range(len(Tf)),j]
clf2.fit(X.T,Out)
if clf2.betasPValue[2]<(0.005/(len(Good_ICs)*3)) and clf2.coef_[2]>0:
Sig[j,2]=1
if clf2.betasPValue[1]<(0.005/(len(Good_ICs)*3))and clf2.coef_[1]>0:
Sig[j,1]=1
if clf2.betasPValue[0]<(0.005/(len(Good_ICs)*3))and clf2.coef_[0]>0:
Sig[j,0]=1
In [342]:
plt.plot(Sig)
Out[342]:
In [343]:
del Final_map
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=npLabel_ICs[Neworder[i]]
Fmaps=np.zeros([S[0],S[1],3])
C=np.zeros([S[3],3])
In [344]:
for j in range(S[3]):
if (Sig[j,2] or Sig[j,1] or Sig[j,0]) and Good_ICs[j]:
# C[j,:]=np.squeeze(np.random.rand(3,1))
C[j,:]=Sig[j,:]
C[j,2]=0
for k in range(3):
M=np.max(np.squeeze(np.reshape(D2[:,:,:,j],S[0]*S[1]*5)))
#M[M==0]=1
# Fmaps[:,:,:,k]=0.7*Dmaps[:,:,:,j]*C[j,k]/np.max(np.squeeze(np.reshape(Dmaps[:,:,:,j]*np.max(C[j,:],S[0]*S[1]*5))
Fmaps[:,:,:,k]=0.7*D2[:,:,:,j]*C[j,k]/M
Final_map=Final_map+Fmaps
In [345]:
pylab.rcParams['figure.figsize'] = (14, 5)
#Final_map[Final_map<0.1]=np.NaN
if S[2]>5:
N=Nstack
else:
N=S[2]
for i in range(N):
plt.subplot(1,N,i+1)
plt.imshow(Dmean[:,:,i],cmap=plt.cm.gray)
plt.imshow(Final_map[:,:,i],interpolation='none')
#plt.imshow(Final_map[:,:,i])
frame1 = plt.gca()
frame1.axes.get_xaxis().set_visible(False)
frame1.axes.get_yaxis().set_visible(False)
In [346]:
pylab.rcParams['figure.figsize'] = (7, 10)
h=1
i=0
for j in range(S[3]):
if (Sig[j,2] or Sig[j,1] or Sig[j,0]) and Good_ICs[j]:
plt.plot(Tball,straightbin*1.5)
plt.plot(Tball,leftbin*1.5)
plt.plot(Tball,rightbin*1.5)
plt.plot(Tb,FrontGroom.T)
plt.plot(Tb,RearGroom.T)
#plt.plot(Time_fluoICA.T,DT[:,Order[j]]+2)
plt.plot(Time_fluoICA.T,(DT[:,Order[j]]+h*i+2),color=C[j,:]/2)
i=i+1
plt.show()
In [347]:
plt.plot(range(600,800),Walk[range(600,800)])
Out[347]:
Detect onsets and offsets
In [352]:
OnSteps=[]
OffSteps=[]
PreStart=np.zeros(Walk.shape)
PreStop=np.zeros(Walk.shape)
for j in range(20,len(Walk)-20):
if (not(any(Walk[range(j-10,j)])) and all(Walk[range(j,j+10)])):
OnSteps.append(j)
PreStart[range(j-20,j)]=1
elif (all(Walk[range(j-10,j)]) and not(any(Walk[range(j,j+10)]))):
OffSteps.append(j)
PreStop[range(j-20,j)]=1
In [353]:
plt.plot(PreStart)
plt.plot(PreStop)
plt.plot(Walk)
plt.show()
In [354]:
X=np.squeeze(np.array([np.squeeze(PreStart),np.squeeze(PreStop)]))
Nreg=X.shape[0]
Sig=np.zeros([S[3],Nreg])
for j in range(S[3]):
if Good_ICs[j]==1:
Out=DT[range(len(Tf)),j]
clf2.fit(X.T,Out)
if clf2.betasPValue[1]<(0.0005/(len(Good_ICs)*2))and clf2.coef_[1]>0:
Sig[j,1]=1
if clf2.betasPValue[0]<(0.0005/(len(Good_ICs)*2))and clf2.coef_[0]>0:
Sig[j,0]=1
In [355]:
del Final_map
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=npLabel_ICs[Neworder[i]]
Fmaps=np.zeros([S[0],S[1],3])
C=np.zeros([S[3],3])
In [356]:
for j in range(S[3]):
if (Sig[j,1] or Sig[j,0]) and Good_ICs[j]:
# C[j,:]=np.squeeze(np.random.rand(3,1))
C[j,range(2)]=Sig[j,:]
for k in range(3):
M=np.max(np.squeeze(np.reshape(D2[:,:,:,j],S[0]*S[1]*5)))
#M[M==0]=1
# Fmaps[:,:,:,k]=0.7*Dmaps[:,:,:,j]*C[j,k]/np.max(np.squeeze(np.reshape(Dmaps[:,:,:,j]*np.max(C[j,:],S[0]*S[1]*5))
Fmaps[:,:,:,k]=0.7*D2[:,:,:,j]*C[j,k]/M
Final_map=Final_map+Fmaps
In [357]:
pylab.rcParams['figure.figsize'] = (14, 5)
#Final_map[Final_map<0.1]=np.NaN
if S[2]>5:
N=Nstack
else:
N=S[2]
for i in range(N):
plt.subplot(1,N,i+1)
plt.imshow(Dmean[:,:,i],cmap=plt.cm.gray)
plt.imshow(Final_map[:,:,i],interpolation='none')
#plt.imshow(Final_map[:,:,i])
frame1 = plt.gca()
frame1.axes.get_xaxis().set_visible(False)
frame1.axes.get_yaxis().set_visible(False)
In [380]:
# 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
filename2 = askopenfilename() # show an "Open" dialog box and return the path to the selected file
print(filename2)
In [389]:
img1 = nb.load(filename2)
InitData = np.array(img1.get_data())
S=InitData.shape
S
Out[389]:
In [390]:
InitDatap=np.zeros(S)
for i in range(S[0]):
for j in range(S[1]):
for k in range(S[2]):
InitDatap[i,j,k,:]=InitData[i,j,k,:]-np.min(InitData[i,j,k,:])
In [393]:
StraightAv=np.zeros([S[0],S[1],S[2]])
for i in range(S[3]-200):
StraightAv=StraightAv+InitDatap[:,:,:,i]*straightintbin[i]
LeftAv=np.zeros([S[0],S[1],S[2]])*len(OffSteps)
for i in range(S[3]-200):
LeftAv=LeftAv+InitDatap[:,:,:,i]*leftintbin[i]
rightAv=np.zeros([S[0],S[1],S[2]])
for i in range(S[3]-200):
rightAv=StraightAv+InitDatap[:,:,:,i]*rightintbin[i]
SLR=np.zeros([S[0],S[1],S[2],3])
SLR[:,:,:,0]=StraightAv/(S[3]-2)
SLR[:,:,:,1]=LeftAv/(S[3]-2)
SLR[:,:,:,2]=rightAv/(S[3]-2)
In [394]:
S2=StraightAv.shape
In [395]:
N=S2[2]/2
for i in range(N):
plt.subplot(1,N,i+1)
plt.imshow(1+SLR[:,:,i*2,:],interpolation='none')
#plt.imshow(Final_map[:,:,i])
frame1 = plt.gca()
frame1.axes.get_xaxis().set_visible(False)
frame1.axes.get_yaxis().set_visible(False)
In [397]:
WalkAv=np.mean(SLR,3)*3
GroomAv=np.zeros([S[0],S[1],S[2]])
for i in range(S[3]-200):
GroomAv=GroomAv+InitDatap[:,:,:,i]*(FrontGroomint[i]+RearGroomint[i])
WG=np.zeros([S[0],S[1],S[2],3])
WG[:,:,:,0]=WalkAv
WG[:,:,:,1]=GroomAv/(S[3]-2)
WG[:,:,:,2]=0
In [399]:
N=S2[2]/2
for i in range(N):
plt.subplot(1,N,i+1)
plt.imshow(WG[:,:,i*2,:],interpolation='none')
#plt.imshow(Final_map[:,:,i])
frame1 = plt.gca()
frame1.axes.get_xaxis().set_visible(False)
frame1.axes.get_yaxis().set_visible(False)
In [400]:
OnSteps
Out[400]:
In [401]:
OffSteps
Out[401]:
In [402]:
len(OnSteps)
Out[402]:
In [403]:
DataOnset=np.zeros([S[0],S[1],S[2],100,len(OnSteps)])
for i in range(len(OnSteps)):
j=OnSteps[i]
DataOnset[:,:,:,:,i]=InitData[:,:,:,range(j-50,j+50)]
In [404]:
DataOffset=np.zeros([S[0],S[1],S[2],100,len(OffSteps)])
for i in range(len(OffSteps)):
j=OffSteps[i]
DataOffset[:,:,:,:,i]=InitData[:,:,:,range(j-50,j+50)]
In [405]:
DataOnsetAv=np.mean(DataOnset,4)
DataOffsetAv=np.mean(DataOffset,4)
In [413]:
D3=np.transpose(DataOffsetAv[:,:,:,:],(3,2,1,0))
nim = NiftiImage(D3)
nim.save('/home/sophie/100106Off.nii')
In [414]:
D3=np.transpose(DataOnsetAv[:,:,:,:],(3,2,1,0))
nim = NiftiImage(D3)
nim.save('/home/sophie/100106On.nii')
In [415]:
S
Out[415]:
In [419]:
DataOnsetcat=np.zeros([S[0],S[1]*len(OnSteps),S[2],100])
for i in range(len(OnSteps)):
j=OnSteps[i]
DataOnsetcat[:,range(S[1]*i,S[1]*(i+1)),:,:]=InitData[:,:,:,range(j-50,j+50)]
DataOffsetcat=np.zeros([S[0],S[1]*len(OnSteps),S[2],100])
for i in range(len(OffSteps)):
j=OffSteps[i]
DataOffsetcat[:,range(S[1]*i,S[1]*(i+1)),:,:]=InitData[:,:,:,range(j-50,j+50)]
In [420]:
D3=np.transpose(DataOnsetcat[:,:,:,:],(3,2,1,0))
nim = NiftiImage(D3)
nim.save('/home/sophie/100106Oncat.nii')
D3=np.transpose(DataOffsetcat[:,:,:,:],(3,2,1,0))
nim = NiftiImage(D3)
nim.save('/home/sophie/100106Offcat.nii')
In [270]:
List1=[(Translated[i],i) for i in range(S[3])]
Newlist=sorted(List1, key=lambda List1: List1[0])
In [271]:
Neworder=[Newlist[i][1] for i in range(S[3]) if Newlist[i][0] != 15]
In [272]:
[Label_ICs[Neworder[i]] for i in range(len(Neworder))]
Out[272]:
In [273]:
NewDT=DT[:,Order[Neworder[:]]].T
for j in range(len(Neworder)):
A=NewDT[:,j]
V=np.sqrt(np.var(A))
NewDT[:,j]=A/V
In [274]:
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 [275]:
h=0.5
pylab.rcParams['figure.figsize'] = (7, 14)
i=0
t=[j*0.005 for j in range(S1[0])]
for j in range(len(Neworder)):
plt.plot(t,NewDT[j,:]+h*j,color=C1[i%6][:])
i=i+1
plt.ylim([-3,h*j+5])
frame1=plt.gca()
frame1.axes.get_yaxis().set_ticks([])
matplotlib.rcParams.update({'font.size': 18})
plt.show()
In [276]:
Newmaps=Dmaps[:,:,:,Order[Neworder[:]]]
In [277]:
L=len(set([Translated[Neworder[i]] for i in range(len(Neworder))]))
Regionmaps=np.zeros([S[0],S[1],L,3])
Datasort=np.zeros([S[0],S[1],S[2],L,3])
Mapsordered=datao[:,:,:,Neworder[:]]
mapfilename=foldername+'mapsordered.nii'
Mapsordered2=np.transpose(Mapsordered,(3,2,1,0))
#nimap = NiftiImage(Mapsordered2)
#nimap.save(mapfilename)
DMapsordered=Dmapso[:,:,:,Neworder[:]]
mapfilename=foldername+'zscoredmapsordered.nii'
DMapsordered2=np.transpose(DMapsordered,(3,2,1,0))
#nimap = NiftiImage(DMapsordered2)
#nimap.save(mapfilename)
In [278]:
j=0
i=0
k=Translated[Neworder[i]]
m=0
Regionname.append(Label_ICs[Neworder[i]])
for i in range(len(Neworder)):
#C2=C1[i%6][:]
for l in range(3):
M=np.max(np.squeeze(np.reshape(Newmaps[:,:,:,i],S[0]*S[1]*S[2])))
Regionmaps[:,:,j,l]=Regionmaps[:,:,j,l]+0.7*np.max(DMapsordered[:,:,:,i],2)*C1[i%6][l]/M
Datasort[:,:,:,j,l]=Datasort[:,:,:,j,l]+Dmaps[:,:,:,Order[Neworder[i]]]*C1[i%6][l]
i=i+1
m=m+1
if i<len(Neworder):
k1=Translated[Neworder[i]]
if k1 != k:
j=j+1
k=k1
m=0
Regionname.append(Label_ICs[Neworder[i]])
Datasort.shape
Out[278]:
In [279]:
Regionname
Out[279]:
In [281]:
DMapscolor=np.sum(Datasort, axis=3)
DMapscolor.shape
mapfilename=foldername+'zscoredmapscolor.nii'
DMapscolor2=np.transpose(DMapscolor,(3,2,1,0))
DMapscolor2.shape
nimap = NiftiImage(DMapscolor2)
nimap.save(mapfilename)
Dmeanproj=np.max(Dmean,2)
Regionmaps[Regionmaps==0]=np.nan
In [282]:
pylab.rcParams['figure.figsize'] = (7, 14)
import scipy
from scipy import ndimage
#from skimage import color
Dmeanprojrgb=np.dstack((Dmeanproj.T,Dmeanproj.T,Dmeanproj.T))
In [284]:
j=0
i=0
k=Translated[Neworder[i]]
m=0
Regionname.append(Label_ICs[Neworder[i]])
for i in range(len(Neworder)):
#C2=C1[i%6][:]
for l in range(3):
M=np.max(np.squeeze(np.reshape(Newmaps[:,:,:,i],S[0]*S[1]*S[2])))
Regionmaps[:,:,j,l]=Regionmaps[:,:,j,l]+0.7*np.max(DMapsordered[:,:,:,i],2)*C1[i%6][l]/M
Datasort[:,:,:,j,l]=Datasort[:,:,:,j,l]+Dmaps[:,:,:,Order[Neworder[i]]]*C1[i%6][l]
i=i+1
m=m+1
if i<len(Neworder):
k1=Translated[Neworder[i]]
if k1 != k:
j=j+1
k=k1
m=0
Regionname.append(Label_ICs[Neworder[i]])
In [286]:
Regionmaps=np.zeros([S[0],S[1],L,3])
Datasort=np.zeros([S[0],S[1],S[2],L,3])
In [287]:
j=0
i=0
k=Translated[Neworder[i]]
m=0
Regionname.append(Label_ICs[Neworder[i]])
for i in range(len(Neworder)):
#C2=C1[i%6][:]
for l in range(3):
M=np.max(np.squeeze(np.reshape(Newmaps[:,:,:,i],S[0]*S[1]*S[2])))
Regionmaps[:,:,j,l]=Regionmaps[:,:,j,l]+0.7*np.max(DMapsordered[:,:,:,i],2)*C1[i%6][l]/M
Datasort[:,:,:,j,l]=Datasort[:,:,:,j,l]+Dmaps[:,:,:,Order[Neworder[i]]]*C1[i%6][l]
i=i+1
m=m+1
if i<len(Neworder):
k1=Translated[Neworder[i]]
if k1 != k:
j=j+1
k=k1
m=0
Regionname.append(Label_ICs[Neworder[i]])
In [288]:
pylab.rcParams['figure.figsize'] = (7, 14)
import scipy
from scipy import ndimage
for i in range(L):
plt.subplot(L,1,L-i)
Rotated_Plot = ndimage.rotate(Regionmaps[:,:,i], -90)
IM=plt.imshow(Rotated_Plot)
frame1 = plt.gca()
frame1.axes.get_xaxis().set_visible(False)
frame1.axes.get_yaxis().set_visible(False)
In [156]:
import scipy.io as sio
sio.savemat(foldername+'NewDT.mat',{'NewDT':NewDT})
In [157]:
names=set([Label_ICs[Neworder[i]] for i in range(len(Neworder))])
for i in range(L):
regionfilename=foldername+Regionname[i]+'.nii'
D3=np.transpose(Datasort[:,:,:,i,:],(2,3,1,0))
nim = NiftiImage(D3)
nim.save(regionfilename)
In [ ]: