<jupyter_start><jupyter_text>pyABF-APisolation - Code Walkthrough pyABF-APisolation, is a simplistic module built to quickly extract action potentials and their raw features from abf files.   This module is built on top of swhardens excellent 'pyABF'. Features are defined by their values located in the Allen Institute's electrophysiology white papers The goal of this document is to show the exactly how the code goes about extracting the features.  Step 0: import the dependencies, import the ABF To begin with we need to import the dependencies. This code utilises numpy, and pandas for data storage and output. Additionally, we utilise the previously written abfderivative tools found here: https://github.com/smestern/abfderivative.<jupyter_code>import numpy as np
from numpy import genfromtxt
import matplotlib.pyplot as plt
from abfderivative import *
from nuactionpotential import *
import pyabf
from pyabf.tools import *
from pyabf import filter
import os
import pandas as pd
import statistics<jupyter_output>C:\Users\SMest\Anaconda3\lib\site-packages\pyabf\tools\__init__.py:8: UserWarning: All pyabf.tools modules are experimental (APIs may change)
  warnings.warn("All pyabf.tools modules are experimental (APIs may change)")<jupyter_text>Now we simply import our ABF file<jupyter_code>abf = pyabf.ABF('./processed/1971003-p.abf')<jupyter_output><empty_output><jupyter_text>After importing the ABF file, we can take a look to ensure its working.<jupyter_code>abf.setSweep(4) #Set to sweep number 4
plt.plot(abf.sweepX, abf.sweepY, 'b-') #plot our sweep<jupyter_output><empty_output><jupyter_text>Now, we can simply call apisolate. This is the function you will want to use most often. Ap isolate returns both a feature array and the raw traces for all the action potentials found in a file.<jupyter_code>aparray, featarray, abf = apisolate(abf, 0, 'tag', False, True, plot=1)
display(featarray.head())<jupyter_output>0
%5 threhold avg: 14373.779296876584
Ap count: 1
1
%5 threhold avg: 13122.558593751448
Ap count: 4
2
%5 threhold avg: 11298.285590278416
Ap count: 13
3
%5 threhold avg: 11102.294921872028
Ap count: 23
4
%5 threhold avg: 10453.444260817821
Ap count: 36
5
%5 threhold avg: 10149.002075196082
Ap count: 52
6
%5 threhold avg: 9719.04553865222
Ap count: 71
7
%5 threhold avg: 9431.384858631
Ap count: 92
8
%5 threhold avg: 10242.919921873214
Ap count: 117
9
%5 threhold avg: 8829.345703125266
Ap count: 142
feattag.csv saved<jupyter_text>As you can see, apisolate() returns a reasonably complete feature array. Now lets look at how it determines these values.  Step 1: Isolate the raw traces First, apisolate() calls appreprocess(). This function serves to isolate the 'raw' traces, as well as some feature data.  appreprocess begins by constructing arrays to fill. Then it calls thresholdavg() to find the 5% threshold avg.<jupyter_code>thresholdsl = (thresholdavg(abf, 4) * 0.05)
print(thresholdsl)<jupyter_output>10453.444260817821<jupyter_text>This function utilizes SWhardens pyabf.tools.ap.ap_points_currentSweep() function. The gist of this function is to look for points where the derivative of the sweep crosses 15 mV/mS.It somewhat acomplishes this by fitting a line and looking for crossings. Then centres on the dV/dT peak.<jupyter_code>slopex, slopey = derivative(abf,4)
index = pyabf.tools.ap.ap_points_currentSweep(abf)
plt.plot(abf.sweepX[0:700], slopey[0:700])
plt.plot(abf.sweepX[0:700], np.full(700, 20))
plt.plot(abf.sweepX[index], slopey[index], 'rx')
plt.xlim(0.020, 0.025)
plt.xlabel('time (S)')
plt.ylabel('dV/dT')<jupyter_output><empty_output><jupyter_text>thresholdavg() simply tabulates the mean of these returned indicies.  From here, appreprocess() takes each index and steps backwards until it reaches below 5% of the threshold. This is marked as our (single) action potential threshold, and the start of our ap.<jupyter_code># Now find the point where DVDT falls below the 5% threshold
apstrt = (int(index[0] - (abf.dataPointsPerMs * 2)))
thresholdslloc = index[0]
indexloc = np.nonzero(np.where(slopey[apstrt:thresholdslloc] < thresholdsl, 1, 0))[0]
if indexloc.size < 1:
    idx = apstrt
else:
    indexloc += apstrt
    idx = indexloc[-1]
apstrt = idx
plt.plot(abf.sweepX[0:700], slopey[0:700])
plt.plot(abf.sweepX[0:700], np.full(700, thresholdsl))
plt.axvspan(abf.sweepX[int(index[0] - (abf.dataPointsPerMs * 2))], abf.sweepX[thresholdslloc], facecolor='#2ca02c', alpha=0.25, label='search area')
plt.plot(abf.sweepX[apstrt], slopey[apstrt], 'gx', label='Intersection at threshold')
plt.xlim(0.020, 0.025)
plt.xlabel('time (S)')
plt.ylabel('dV/dT')
plt.legend()<jupyter_output><empty_output><jupyter_text>Next, the algorithm searches forwards for the absolute peak of the Action potential<jupyter_code>apstrt = (int(index[0] - (abf.dataPointsPerMs * 2)))
if apstrt < 0: 
    apstrt=0
apend = int(index[0]  + (abf.dataPointsPerMs * 3)) 
aploc = np.argmax(abf.sweepY[apstrt:apend]) + apstrt
plt.plot(abf.sweepX[0:700], abf.sweepY[0:700])
plt.axvspan(abf.sweepX[apstrt], abf.sweepX[apend], facecolor='#2ca02c', alpha=0.25, label='search area')
plt.plot(abf.sweepX[aploc], abf.sweepY[aploc], 'gx', label='Peak mV')
plt.xlim(0.020, 0.0275)
plt.xlabel('time (S)')
plt.ylabel('mV')
plt.legend()<jupyter_output><empty_output><jupyter_text>Now that we have an mostly established peak, we need to cap off the action potential. We cap at either 10ms post peak or the next ap, whatever is sooner.<jupyter_code>## Now we check to ensure the action potentials do not over lap
if((index[0] - index[1]) > (abf.dataPointsPerMs * 10)): ##if the next ap is over 10ms away then we simple cap off at 10ms
    apend = abs(int(aploc + abf.dataPointsPerMs * 10))
elif apend > index[1]:
    apend = index[1] #otherwise we cap the end at the next threshold
    aploc = np.argmax(abf.sweepY[apstrt:apend]) + apstrt #and re-find the peak
else:
        apend = abs(int(aploc + abf.dataPointsPerMs * 10)) #if this is the last ap in the sweep we cap at 10ms
k,  = abf.sweepY.shape
print(index)
if apend > k:
        apend = int(k) - 1
plt.plot(abf.sweepX[0:apend+200], abf.sweepY[0:apend+200])
plt.axvspan(abf.sweepX[apstrt], abf.sweepX[(abf.dataPointsPerMs * 10)+aploc], facecolor='#2ca02c', alpha=0.25, label='search area')
plt.plot(abf.sweepX[apend], abf.sweepY[apend], 'gx', label='Peak mV')
plt.xlim(0.020, 0.035)
plt.xlabel('time (S)')
plt.ylabel('mV')
plt.legend()<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup