SciPyProject.py

# SciPy project by Pavel Svehla

# this script processes csv files generated by a software called Neuronstdio, in which dendritic spines in confocal images were manually reconstructed. 
# the aim of the experiment is to see whether pathogenic antibodies isolated from patients affect neuronal morphology, namely post-synapses.
# mice were chronicly infused with these antibodies targeting NMDAR and cauing its downregulation. There are two groups: control antibody and NMDAR antibody.
# Data is organized in folder 'data' containing subfolders representing the experimental groups, each group contains subfolders representing a mouse.


import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import os
from os import listdir
from os.path import isfile, isdir, join
import scipy.stats as stats
from sklearn.decomposition import PCA
from sklearn import cluster

data = pd.DataFrame() # initialize a dataframe to append data to it

mypath = join(os.path.dirname(__file__), 'data') 

##### get the data
groups = [f for f in listdir(mypath)]
for group in groups:
    
    subfolders = [f for f in listdir(join(mypath, group))] # extraxt folder names
    for subfolder in subfolders: # extract subfolder names
        files = [f for f in listdir(join(mypath, group, subfolder))]
        for file in files: # extract files from subfolders
            data_tmp = pd.read_csv(join(mypath, group, subfolder, file))
            
            # add more columns to the ata
            ls = [1]
            ls[0] = group
            ls = ls*len(data_tmp)
            data_tmp['ExpGroup'] = ls
            
            ls = [1]
            ls[0] = subfolder
            ls = ls*len(data_tmp)
            data_tmp['MouseID'] = ls
            
            ls = [1]
            ls[0] = join(group, subfolder, file) # use this if filenames are not unique
            #ls[0] = file
            ls = ls*len(data_tmp)
            data_tmp['FileID'] = ls
            
            # append spine density to the data
            auto_ids = data_tmp.loc[:,'AUTO'] == 'yes' 
            SpineDensity = [1]
#            DendriticLength = np.array(data_tmp.loc[auto_ids,:]['SECTION-LENGTH'].unique()).sum()
            DendriticLength = data_tmp.loc[auto_ids,:]['SECTION-LENGTH'].unique().sum()
            nSpines = data_tmp.loc[:,:].nunique()['ID']
            SpineDensity[0] = nSpines / DendriticLength
            SpineDensity = SpineDensity * len(data_tmp)
            data_tmp['SpineDensity'] = SpineDensity
            
            # classify the spines
            NeckRatio = 1.1
            ThinRatio = 2.5
            HeadDiameterCrit = 0.35
            
            HeadDiameter = data_tmp.loc[:,'HEAD-DIAMETER']
            NeckDiameter = data_tmp.loc[:,'NECK-DIAMETER']
            
            MushroomCrit1 = HeadDiameter / NeckDiameter > NeckRatio
            MushroomCrit2 = HeadDiameter >= HeadDiameterCrit
            Mushroom = [ MushroomCrit1[i] and MushroomCrit2[i] for i in range(len(MushroomCrit1))] 
            
            ThinCrit1 = HeadDiameter / NeckDiameter > NeckRatio
            ThinCrit2 = HeadDiameter < HeadDiameterCrit
            Thin = [ ThinCrit1[i] and ThinCrit2[i] for i in range(len(ThinCrit1))]
            
            Stubby = np.isnan(NeckDiameter)

            SpineClass = np.array([np.nan] * len(ThinCrit1))
            SpineClass[Mushroom] = 3 
            SpineClass[Thin] = 2
            SpineClass[Stubby] = 1 
            
            data_tmp['TYPE'] = SpineClass
            
            data = data.append(data_tmp) # combine in one big dataframe


############## analysis
# print some output
print('total of %i spines were counted in %i images' %(len(data), len(data.groupby('FileID'))))

# spine density
auto_ids = data.loc[:,'AUTO'] == 'yes' # indexes to automatic spines, manual spines are problematic (incorrect data)
SpineDensity = [None]*len(groups)
for group_ind, group in enumerate(groups): 
    group_ids = data.loc[:,'ExpGroup'] == group
    # retrieve section lengths of individual images
    DendriticLength = np.array(data.loc[group_ids & auto_ids,:].groupby('FileID')['SECTION-LENGTH'].unique())
    for i in range(len(DendriticLength)):
        DendriticLength[i] = DendriticLength[i].sum() # sum each array withing the array of arrays
        
    # number of spines per image
    nSpines = np.array(data.loc[group_ids,:].groupby('FileID').nunique()['ID'])   
    # number of spines divided by dendritic length
    SpineDensity[group_ind] = nSpines / DendriticLength

# do a t test and print the result
t,p = stats.ttest_ind(SpineDensity[0], SpineDensity[1])
print('t-test p-value is %g' %(p))

############## plots
# spine morphology
data.loc[:,'NECK-DIAMETER'][np.isnan(data.loc[:,'NECK-DIAMETER']) == True] = 0 # get rid of nans

f, axs = plt.subplots(1,3)
nbins = 15
 
StuffToPlot = [data.loc[:,'HEAD-DIAMETER'][data.loc[:,'ExpGroup'] == groups[0]], data.loc[:,'HEAD-DIAMETER'][data.loc[:,'ExpGroup'] == groups[1] ]]
axs[0].hist(StuffToPlot, bins=nbins, density=True, cumulative=True)    
axs[0].set_title('Head diameter')
axs[0].set_xlabel('um')    
StuffToPlot = [data.loc[:,'NECK-DIAMETER'][data.loc[:,'ExpGroup'] == groups[0]], data.loc[:,'NECK-DIAMETER'][data.loc[:,'ExpGroup'] == groups[1] ]]
axs[1].hist(StuffToPlot, bins=nbins, density=True, cumulative=True)  
axs[1].set_title('Neck diameter')            
axs[1].set_xlabel('um')
StuffToPlot = [data.loc[:,'MAX-DTS'][data.loc[:,'ExpGroup'] == groups[0]], data.loc[:,'MAX-DTS'][data.loc[:,'ExpGroup'] == groups[1] ]]
axs[2].hist(StuffToPlot, label=groups, bins=nbins, density=True, cumulative=True)   
axs[2].set_title('Spine length')           
axs[2].set_xlabel('um')
f.legend()
plt.savefig('histograms.png')

# PCA and clustering analysis
f,ax = plt.subplots()
pca = PCA(n_components=2)
data_pca = pca.fit_transform(data.loc[:,('HEAD-DIAMETER','NECK-DIAMETER', 'MAX-DTS','SOMA-DISTANCE', 'XYPLANE-ANGLE')]) # subset meaningful data

colorset = np.array(['red', 'orange', 'green', 'blue', 'magenta', 'gray', 'brown', 'black']*10)

dbscan = cluster.DBSCAN(eps=2.5, min_samples=1000) 
dbscan.fit(data_pca) 
cids = dbscan.labels_ # get resulting cluster IDs from the kmeans object, one for each sample
colors = colorset[cids] # convert cluster IDs to colors

ax.scatter(data_pca[:, 0], data_pca[:, 1], color=colors) 
ax.set_title('PCA reduced data')
ax.set_xlabel('PC1')
ax.set_ylabel('PC2')
plt.savefig('PCA.png')


######### looking for outliers
idxs = data['SpineDensity'] < 1
data.loc[idxs,:].groupby('FileID').mean()

## mess goes here
autoIDs = data.loc[:, 'AUTO'] == 'yes' # use only auto spines for dendritic lengths 
DendriticLength = np.array(data.loc[autoIDs,:].groupby('FileID')['SECTION-LENGTH'].unique())

dataAuto = data.loc[autoIDs, :].groupby('FileID')