abf_analysis

abf_analysis.py

"""
This file read a .abf file, find the peaks and troughs, and calculate delta values

Requirements:
pyabf
numpy
matplotlib
"""
import numpy as np
import sys
import pyabf
import matplotlib.pyplot as plt

def sample(y, num):
    N = len(y)
    assert N > 10*num
    step = N // num
    idx = np.arange(0, N-1, step=step)
    return idx

def smooth(y, ks):
    assert ks % 2 == 1
    N = len(y)
    new_y = np.zeros_like(y)
    for i in range(ks//2, N-ks//2-1):
        new_y[i] = np.mean(y[i-ks//2:i+ks//2+1])
    return new_y

def gradient(y, step):
    grads = np.zeros_like(y)
    N = len(y)
    for i in range(1, N-step):
        grads[i] = (y[i+step] - y[i-step]) / (step*2)
    return grads

if __name__ == '__main__':
    fname = sys.argv[1]
    abf = pyabf.ABF(fname)
    abf.setSweep(sweepNumber=0, channel=1)
    x = abf.sweepX
    y = abf.sweepY
    # print(f'original number: {len(y)}') 
    idx = sample(y, 1000)
    x, orig_y = x[idx], y[idx]
    y = smooth(orig_y, 11)
    # print(f'new number: {len(y)}')
    # plt.plot(y)

    grads = gradient(y, step=1)
    # plt.plot(grads)

    N = len(grads)
    gap = 5
    ends = []
    new_round = True
    for i in range(gap, N-gap, gap):
        if new_round and grads[i] < -0.5 and grads[i+gap] > 0:
            ends.append(i+gap)
            new_round = False
        if grads[i] > 1:
            new_round = True
    print(f'Found {len(ends)} troughs, x: {x[ends].astype(int)}')

    prev = 0
    starts = [[] for _ in ends]
    for i in range(len(ends)):
        for j in range(prev, ends[i], gap):
            if abs(grads[j]) < 1 and grads[j+gap] < -1:
                starts[i].append(j)
        # print(i, starts)
        prev = ends[i] + gap

    real_ends = []
    real_starts = []
    for i in range(len(ends)):
        if len(starts[i]) != 1:
            print(f'For {i}th trough, invalid {len(starts[i])} peaks, it will be removed.')
        else:
            real_ends.append(ends[i])
            real_starts.append(starts[i][0])

    ends = np.array(real_ends)
    starts = np.array(real_starts)

    for start, end in zip(starts, ends):
        high = np.max(orig_y[start:end])
        low = np.min(orig_y[start:end])
        print(high - low)

    plt.plot(x, y)
    plt.scatter(x[starts], y[starts], color='green')
    plt.scatter(x[ends], y[ends], color='red')
    plt.show()