added assembly analysis + occupancy plots for the DFG proposal

analysis/assembly.py

@@ -0,0 +1,228 @@
+'''
+	Codes for PCA/ICA methods described in Detecting cell assemblies in large neuronal populations, Lopes-dos-Santos et al (2013).
+											https://doi.org/10.1016/j.jneumeth.2013.04.010
+	This implementation was written in Feb 2019.
+	Please e-mail me if you have comments, doubts, bug reports or criticism (Vítor, vtlsantos@gmail.com /  vitor.lopesdossantos@pharm.ox.ac.uk).
+'''
+
+
+from sklearn.decomposition import PCA
+from scipy import stats
+import numpy as np
+from numpy import matlib as mb
+np.matlib = mb
+
+__author__ = "Vítor Lopes dos Santos"
+__version__ = "2019.1"
+
+def toyExample(assemblies, nneurons = 10, nbins = 1000, rate = 1.):
+        
+        np.random.seed()
+    
+        actmat = np.random.poisson(rate,nneurons*nbins).reshape(nneurons,nbins)
+        assemblies.actbins = [None]*len(assemblies.membership)
+        for (ai,members) in enumerate(assemblies.membership):
+            
+                members = np.array(members)
+                nact = int(nbins*assemblies.actrate[ai])
+                actstrength_ = rate*assemblies.actstrength[ai]
+                
+                actbins = np.argsort(np.random.rand(nbins))[0:nact]
+                
+                actmat[members.reshape(-1,1),actbins] = np.ones((len(members),nact))+actstrength_
+                
+                assemblies.actbins[ai] = np.sort(actbins)
+
+        return actmat
+
+class toyassemblies:
+    
+        def __init__(self, membership, actrate, actstrength):
+
+                self.membership  = membership
+                self.actrate        = actrate
+                self.actstrength = actstrength
+
+def marcenkopastur(significance):
+    
+        nbins = significance.nbins
+        nneurons = significance.nneurons
+        tracywidom = significance.tracywidom
+    
+        # calculates statistical threshold from Marcenko-Pastur distribution
+        q = float(nbins)/float(nneurons) # note that silent neurons are counted too
+        lambdaMax = pow((1+np.sqrt(1/q)),2)
+        lambdaMax += tracywidom*pow(nneurons,-2./3) # Tracy-Widom correction 
+        
+        return lambdaMax
+    
+def getlambdacontrol(zactmat_):
+
+        significance_ = PCA()
+        significance_.fit(zactmat_.T)
+        lambdamax_ = np.max(significance_.explained_variance_)
+        
+        return lambdamax_
+    
+def binshuffling(zactmat,significance):
+    
+        np.random.seed()
+
+        lambdamax_ = np.zeros(significance.nshu)
+        for shui in range(significance.nshu):
+                zactmat_ = np.copy(zactmat)
+                for (neuroni,activity) in enumerate(zactmat_):
+                        randomorder = np.argsort(np.random.rand(significance.nbins))
+                        zactmat_[neuroni,:] = activity[randomorder]
+                lambdamax_[shui] = getlambdacontrol(zactmat_)
+
+        lambdaMax = np.percentile(lambdamax_,significance.percentile)
+        
+        return lambdaMax
+    
+def circshuffling(zactmat,significance):
+    
+        np.random.seed()
+
+        lambdamax_ = np.zeros(significance.nshu)
+        for shui in range(significance.nshu):
+                zactmat_ = np.copy(zactmat)
+                for (neuroni,activity) in enumerate(zactmat_):
+                        cut = int(np.random.randint(significance.nbins*2))
+                        zactmat_[neuroni,:] = np.roll(activity,cut)
+                lambdamax_[shui] = getlambdacontrol(zactmat_)
+
+        lambdaMax = np.percentile(lambdamax_,significance.percentile)
+        
+        return lambdaMax
+
+def runSignificance(zactmat,significance):
+    
+        if significance.nullhyp == 'mp':
+                lambdaMax = marcenkopastur(significance)
+        elif significance.nullhyp == 'bin':
+                lambdaMax = binshuffling(zactmat,significance)
+        elif significance.nullhyp == 'circ':
+                lambdaMax = circshuffling(zactmat,significance)
+        else: 
+                print('ERROR !')
+                print('    nyll hypothesis method '+str(nullhyp)+' not understood')
+                significance.nassemblies = np.nan
+                
+        nassemblies = np.sum(significance.explained_variance_>lambdaMax)
+        significance.nassemblies = nassemblies
+        significance.lambdaMax = lambdaMax
+        
+        return significance
+        
+def extractPatterns(actmat,significance,method):
+        nassemblies = significance.nassemblies
+    
+        if method == 'pca':
+                idxs = np.argsort(-significance.explained_variance_)[0:nassemblies]
+                patterns = significance.components_[idxs,:]
+        elif method == 'ica':
+                from sklearn.decomposition import FastICA
+                ica = FastICA(n_components=nassemblies)
+                ica.fit(actmat.T)
+                patterns = ica.components_
+        else:
+                print('ERROR !')
+                print('    assembly extraction method '+str(method)+' not understood')
+                patterns = np.nan
+                
+        
+                
+        if patterns is not np.nan:
+            
+                patterns = patterns.reshape(nassemblies,-1)
+                
+                # sets norm of assembly vectors to 1
+                norms = np.linalg.norm(patterns,axis=1)
+                patterns /= np.matlib.repmat(norms,np.size(patterns,1),1).T
+        
+        return patterns
+
+def runPatterns(actmat, method='ica', nullhyp = 'mp', nshu = 1000, percentile = 99, tracywidom = False):
+    
+        '''
+        INPUTS
+        
+            actmat:     activity matrix - numpy array (neurons, time bins) 
+            
+            nullhyp:    defines how to generate statistical threshold for assembly detection.
+                            'bin' - bin shuffling, will shuffle time bins of each neuron independently
+                            'circ' - circular shuffling, will shift time bins of each neuron independently
+                                                                obs: mantains (virtually) autocorrelations
+                            'mp' - Marcenko-Pastur distribution - analytical threshold
+                            
+            nshu:       defines how many shuffling controls will be done (n/a if nullhyp is 'mp')
+            
+            percentile: defines which percentile to be used use when shuffling methods are employed.
+                                                                        (n/a if nullhyp is 'mp')
+                                                                         
+            tracywidow: determines if Tracy-Widom is used. See Peyrache et al 2010.
+                                                    (n/a if nullhyp is NOT 'mp')
+                                                    
+        OUTPUTS
+            
+            patterns:     co-activation patterns (assemblies) - numpy array (assemblies, neurons)
+            significance: object containing general information about significance tests 
+            zactmat:      returns z-scored actmat
+        
+        '''
+    
+        nneurons = np.size(actmat,0)
+        nbins = np.size(actmat,1)
+        
+        silentneurons = np.var(actmat,axis=1)==0
+        actmat_ = actmat[~silentneurons,:]
+        
+        # z-scoring activity matrix
+        zactmat_ = stats.zscore(actmat_,axis=1)
+        
+        # running significance (estimating number of assemblies)
+        significance = PCA()
+        significance.fit(zactmat_.T)
+        significance.nneurons = nneurons
+        significance.nbins = nbins
+        significance.nshu = nshu
+        significance.percentile = percentile
+        significance.tracywidom = tracywidom
+        significance.nullhyp = nullhyp
+        significance = runSignificance(zactmat_,significance)
+        if np.isnan(significance.nassemblies):
+                return
+
+        if significance.nassemblies<1:
+                print('WARNING !')
+                print('    no assembly detecded!')
+                patterns = []
+        else:
+                # extracting co-activation patterns
+                patterns_ = extractPatterns(zactmat_,significance,method)
+                if patterns_ is np.nan:
+                	return
+		
+		# putting eventual silent neurons back (their assembly weights are defined as zero)
+                patterns = np.zeros((np.size(patterns_,0),nneurons))
+                patterns[:,~silentneurons] = patterns_
+        zactmat = np.copy(actmat)
+        zactmat[~silentneurons,:] = zactmat_
+		
+        return patterns,significance,zactmat
+    
+def computeAssemblyActivity(patterns,zactmat,zerodiag = True):
+
+        nassemblies = len(patterns)
+        nbins = np.size(zactmat,1)
+
+        assemblyAct = np.zeros((nassemblies,nbins))
+        for (assemblyi,pattern) in enumerate(patterns):
+                projMat = np.outer(pattern,pattern)
+                projMat -= zerodiag*np.diag(np.diag(projMat))
+                for bini in range(nbins):
+                        assemblyAct[assemblyi,bini] = \
+                                np.dot(np.dot(zactmat[:,bini],projMat),zactmat[:,bini])
+                        
+        return assemblyAct