<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Cell Types Database
This notebook demonstrates most of the features of the AllenSDK that help manipulate data in the Cell Types Database.  The main entry point will be through the CellTypesCache class.
CellTypesCache is responsible for downloading Cell Types Database data to a standard directory structure on your hard drive.  If you use this class, you will not have to keep track of where your data lives, other than a root directory.
Step1:  The data_set variable is an NwbDataSet instance, which has some methods we can use to access the injected current stimulus waveform and the voltage response waveform for all experimental sweeps.  Let's pull one sweep out and plot it.
Step2:  Cell Morphology Reconstructions
The Cell Types Database also contains 3D reconstructions of neuronal morphologies. The data are presented in the SWC format. We'll download a particular cell's reconstrution here.
Step3:  The AllenSDK contains a module that makes it easier to work with the SWC files. We'll see how the data is contained in the file by looking at the first node.
Step4:  Note that the type field refers to the type of neuronal compartment. The values can be 1 for the soma, 2 for the axon, 3 for dendrites, and 4 for apical dendrites (if present).
We can use this data to draw lines between each node and all its children to get a drawing of the cell. We'll do it looking at it from the front and from the side.
Step5:  Electrophysiology Features
The Cell Types Database contains a set of features that have already been computed, which could serve as good starting points for analysis. We can query the database using the SDK to get these features.
Step6:  That's how to get all the ephys features for a given specimen - what if we want a particular feature for all cells?
Step7:  Let's use numpy to fit a regression line to these data and plot it.
Step8:  It looks like there may be roughly two clusters in the data above. Maybe they relate to whether the cells are presumably excitatory (spiny) cells or inhibitory (aspiny) cells. Let's query the API and split up the two sets to see.
Step9:  Morphology Features
The Cell Types Database contains a set of precomputed morphological features for cells that have reconstructions.  You can access morphology features by themselves, or combined with the electrophysiology features.
Step10:  Computing Electrophysiology Features
The AllenSDK contains the code used to compute the electrophysiology features you accessed above.  You can run it yourself like this.
Step11:  A list comprehension is an easy way to pull out the spike times.
<ASSISTANT_TASK:>
Python Code:
from allensdk.core.cell_types_cache import CellTypesCache
# Instantiate the CellTypesCache instance.  The manifest_file argument
# tells it where to store the manifest, which is a JSON file that tracks
# file paths.  If you supply a relative path (like this), it will go
# into your current working directory
ctc = CellTypesCache(manifest_file='cell_types/manifest.json')
# this saves the NWB file to 'cell_types/specimen_464212183/ephys.nwb'
data_set = ctc.get_ephys_data(464212183)
Explanation: Cell Types Database
This notebook demonstrates most of the features of the AllenSDK that help manipulate data in the Cell Types Database.  The main entry point will be through the CellTypesCache class.
CellTypesCache is responsible for downloading Cell Types Database data to a standard directory structure on your hard drive.  If you use this class, you will not have to keep track of where your data lives, other than a root directory.
End of explanation
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
sweep_number = 30
sweep_data = data_set.get_sweep(sweep_number)
index_range = sweep_data["index_range"]
i = sweep_data["stimulus"][0:index_range[1]+1] # in A
v = sweep_data["response"][0:index_range[1]+1] # in V
i *= 1e12 # to pA
v *= 1e3 # to mV
sampling_rate = sweep_data["sampling_rate"] # in Hz
t = np.arange(0, len(v)) * (1.0 / sampling_rate)
plt.style.use('ggplot')
fig, axes = plt.subplots(2, 1, sharex=True)
axes[0].plot(t, v, color='black')
axes[1].plot(t, i, color='gray')
axes[0].set_ylabel("mV")
axes[1].set_ylabel("pA")
axes[1].set_xlabel("seconds")
Explanation: The data_set variable is an NwbDataSet instance, which has some methods we can use to access the injected current stimulus waveform and the voltage response waveform for all experimental sweeps.  Let's pull one sweep out and plot it.
End of explanation
from allensdk.core.cell_types_cache import CellTypesCache
ctc = CellTypesCache()
# this downloads metadata for all cells with reconstructions
cells = ctc.get_cells(require_reconstruction=True)
print "Cells with reconstructions: ", len(cells)
# download and open an SWC file
cell_id = 464212183
morphology = ctc.get_reconstruction(cell_id) 
print morphology
Explanation: Cell Morphology Reconstructions
The Cell Types Database also contains 3D reconstructions of neuronal morphologies. The data are presented in the SWC format. We'll download a particular cell's reconstrution here.
End of explanation
compartment_list = morphology.compartment_list
print compartment_list[0]
Explanation: The AllenSDK contains a module that makes it easier to work with the SWC files. We'll see how the data is contained in the file by looking at the first node.
End of explanation
fig, axes = plt.subplots(1, 2, sharey=True, sharex=True)
axes[0].set_aspect('equal', 'box-forced')
axes[1].set_aspect('equal', 'box-forced')
# Make a line drawing of x-y and y-z views
for n in morphology.compartment_list:
    child_nodes = [c for c in morphology.compartment_list if c['id'] in n['children']] 
    for c in child_nodes:
        axes[0].plot([n['x'], c['x']], [n['y'], c['y']], color='black')
        axes[1].plot([n['z'], c['z']], [n['y'], c['y']], color='black')
axes[0].set_ylabel('y')
axes[0].set_xlabel('x')
axes[1].set_xlabel('z')
Explanation: Note that the type field refers to the type of neuronal compartment. The values can be 1 for the soma, 2 for the axon, 3 for dendrites, and 4 for apical dendrites (if present).
We can use this data to draw lines between each node and all its children to get a drawing of the cell. We'll do it looking at it from the front and from the side.
End of explanation
# download all electrophysiology features for all cells
ephys_features = ctc.get_ephys_features()
# filter down to a specific cell
specimen_id = 464212183
cell_ephys_features = [f for f in ephys_features if f['specimen_id'] == specimen_id]
print cell_ephys_features
Explanation: Electrophysiology Features
The Cell Types Database contains a set of features that have already been computed, which could serve as good starting points for analysis. We can query the database using the SDK to get these features.
End of explanation
updown = np.array([f['upstroke_downstroke_ratio_short_square'] for f in ephys_features], dtype=float)
fasttrough = np.array([f['fast_trough_v_long_square'] for f in ephys_features], dtype=float)
plt.figure()
plt.scatter(fasttrough, updown, color='#2ca25f')
plt.ylabel("upstroke-downstroke ratio")
plt.xlabel("fast trough depth (mV)")
Explanation: That's how to get all the ephys features for a given specimen - what if we want a particular feature for all cells?
End of explanation
A = np.vstack([fasttrough, np.ones_like(updown)]).T
print "First 5 rows of A:"
print A[:5, :]
m, c = np.linalg.lstsq(A, updown)[0]
print "m", m, "c", c
plt.figure()
plt.scatter(fasttrough, updown, color='#2ca25f')
plt.plot(fasttrough, m * fasttrough + c, c='gray')
plt.ylabel("upstroke-downstroke ratio")
plt.xlabel("fast trough depth (mV)")
Explanation: Let's use numpy to fit a regression line to these data and plot it.
End of explanation
cells = ctc.get_cells()
cell_index = { c['id']: c for c in cells}
dendrite_types = ['spiny', 'aspiny']
data = {}
# group fast trough depth and upstroke downstroke ratio values by cell dendrite type
for dendrite_type in dendrite_types:
    type_features = [f for f in ephys_features if cell_index[f['specimen_id']]['dendrite_type'] == dendrite_type]
    data[dendrite_type] = {
        "fasttrough": [f['fast_trough_v_long_square'] for f in type_features],
        "updown": [f['upstroke_downstroke_ratio_short_square'] for f in type_features],
    }
    
plt.figure()
for a_type, color in zip(dendrite_types, ["#d95f02", "#7570b3"]):
    plt.scatter(data[a_type]['fasttrough'], data[a_type]['updown'], color=color, label=a_type)
plt.legend(loc='best')
plt.ylabel("upstroke-downstroke ratio")
plt.xlabel("fast trough depth (mV)")
Explanation: It looks like there may be roughly two clusters in the data above. Maybe they relate to whether the cells are presumably excitatory (spiny) cells or inhibitory (aspiny) cells. Let's query the API and split up the two sets to see.
End of explanation
# download all morphology features for cells with reconstructions
morphology_features = ctc.get_morphology_features()
# or download both morphology and ephys features
# this time we'll ask the cache to return a pandas dataframe
all_features = ctc.get_all_features(dataframe=True, require_reconstruction=True)
all_features
Explanation: Morphology Features
The Cell Types Database contains a set of precomputed morphological features for cells that have reconstructions.  You can access morphology features by themselves, or combined with the electrophysiology features.
End of explanation
from allensdk.ephys.feature_extractor import EphysFeatureExtractor
sweep_number = 35
sweep_data = data_set.get_sweep(sweep_number)
index_range = sweep_data["index_range"]
i = sweep_data["stimulus"][0:index_range[1]+1] # in A
v = sweep_data["response"][0:index_range[1]+1] # in V
i *= 1e12 # to pA
v *= 1e3 # to mV
sampling_rate = sweep_data["sampling_rate"] # in Hz
t = np.arange(0, len(v)) * (1.0 / sampling_rate)
fx = EphysFeatureExtractor()
stim_start = 1.0
stim_duration = 1.0
fx.process_instance("", v, i, t, stim_start, stim_duration, "")
feature_data = fx.feature_list[0].mean
print "Avg spike width: {:.2f} ms".format(feature_data['width'])
print "Avg spike threshold: {:.1f} mV".format(feature_data["threshold"])
import pprint
pp = pprint.PrettyPrinter(indent=2)
pp.pprint(feature_data["spikes"][0])
Explanation: Computing Electrophysiology Features
The AllenSDK contains the code used to compute the electrophysiology features you accessed above.  You can run it yourself like this.
End of explanation
spike_times = [s["t"] for s in feature_data["spikes"]]
print spike_times[:5] 
plt.figure()
plt.plot(t, v, color='black')
min_v = v.min()
min_v -= 5.0
plt.scatter(spike_times, np.ones(len(spike_times)) * min_v, c='r')
plt.xlim(0.9, 1.2)
Explanation: A list comprehension is an easy way to pull out the spike times.
End of explanation
<END_TASK>