Unbiased Data Analysis for the Parameterization of Fast Translocation Events through Nanopores.

Single-molecule nanopore electrophysiology is an emerging technique for the detection of analytes in aqueous solutions with high sensitivity. These detectors have proven applicable for the enzyme-assisted sequencing of oligonucleotides. There has recently been an increased interest in the use of nanopores for the fingerprinting of peptides and proteins, referred to as single-molecule nanopore spectrometry. However, the analysis of the resulting electrophysiology traces remains complicated due to the fast unassisted translocation of such analytes, usually in the order of micro- to milliseconds, and the small ion current signal produced (in the picoampere range). Here, we present the application of a generalized normal distribution function (gNDF) for the characterization of short-lived ion current signals (blockades). We show that the gNDF can be used to determine if the observed blockades have adequate time to reach their maximum current plateau while also providing a description of each blockade based on the open pore current (I O), the difference caused by the pore blockade (ΔI B), the position in time (μ), the standard deviation (σ), and a shape parameter (β), leaving only the noise component. In addition, this method allows the estimation of an ideal range of low-pass filter frequencies that contains maximum information with minimal noise. In summary, we show a parameter-free and generalized method for the analysis of short-lived ion current blockades, which facilitates single-molecule nanopore spectrometry with minimal user bias.

Introduction

Single-molecule nanopore electrophysiology is an emerging field with the most notable application in label- and amplification-free nucleic acid sequencing, enabling the identification of long stretches of nucleic acids as they translocate through a nanopore.[1−4] The working principle is similar to nanometer-sized Coulter counters in allowing the differentiation between nucleotides (and other molecules) due to differences in ion displacement (Figure A).[5] The current generation of nanopores is now capable of detecting individual amino acids,[6] opening a new frontier for the sequence identification of proteins. However, whereas nucleic acid sequencing typically needs to differentiate between four bases, proteins can consist of 20 or more amino acids with widely diverging sizes and charges. A helicase attached to the top of the nanopore can be used to aid the translocation of uniformly charged nucleic acids, whereas due to the variability of the protein samples, alternative strategies are required. Even though proof-of-concept studies have shown that motor proteins can be used to control translocation of peptides across nanopores,[7−10] they are currently not able to resolve functional sequence information. Alternatively, proteins can be identified “bottom-up”, based on their fingerprint following a trypsin digest, as commonly done in mass spectrometry experiments.[11−13]

Figure 1

Nanopore experiments, dwell time skewness, and simulated pulse-dilated ideal pulses. (A) Schematic representation of nanopore electrophysiology. A measurement potential is applied across a nonconductive membrane with a nanometer-sized aperture inserted. A stable current can be measured across the open pore, and, upon entry of an analyte, the current is reduced relative to the ion exclusion caused by the analyte. (B) Raw current trace of 10 μM (added to cis) of a pentapeptide (YAGFL) measured in 1 M KCl, pH 4.2 buffered using 15 mM citric acid, and bis-tris-propane under an applied potential of negative 100 mV (cis–trans) at a sampling frequency of 500 kHz with a recording bandwidth of 100 kHz and filtered between 25 Hz and 5 kHz, each represented with a separate line using a digital Gaussian filter using fragaceatoxin C modified with a tryptophan at position 13 (G13W-FraC). (C) Excluded current set against the dwell time for a measurement performed as in part A. For (A,B), experiments were performed using an Axon 200B amplified (molecular devices) coupled to a Digidata 1550B (Molecular Devices). Event characteristics were determined using the single-channel search method as implemented in ClampFit 10. The events were exported to an Excel sheet and plotted using matplotlib (Python 3.9). (D) A 10 μs pulse is shown together with the resulting pulses after 50 and 10 kHz 4-pole Bessel filtering. (E) A 20 μs pulse is shown together with the resulting pulses after 50 and 10 kHz 4-pole Bessel filtering. Simulation was performed using a continuous-time linear system (lsim2) as implemented in SciPy on ideal pulses sampled at 100 MHz.

When a nanopore has a voltage applied across it, an ionic current can be measured, which directly relates to the physicochemical properties of that specific pore type, that is, its size and the charge distribution on its surface. A reduction in this measured ionic current, attributed to the entry of a molecule into the nanopore, is termed an “event” (Figure B). Fast translocation events have been historically dealt with by modifying nanopores chemically,[14] and to build fast and accurate current amplifiers,[15−18] events resulting from a tryptic digest of a protein should correspond to the entry and translocation of individual peptides through the pore. In order to identify events generated by the entry of molecules into a nanopore, the following analysis steps are generally required: (1) preprocessing, involving, for example, smoothing of the data, and removal of instrument noise, (2) event localization, that is, “when does an event occur”, which is commonly performed using a threshold search algorithm, and (3) event characterization, involving the extraction of event characteristics based on the fluctuations observed during an event. It is hypothesized that additional information from short-lived events contained within individual noise components, such as frequency fingerprinting or a Fourier transform, could allow individual characterization of molecules.[19−22]

The potential for obtaining additional information about the sample from short events occurring near the noise floor of the raw data creates a need for critical preprocessing and characterization. Optimal preprocessing is an intricate problem as the residence time of peptides inside the nanopore is generally in the order of micro- to milliseconds, and the change in signal induced by the current blockade is in the picoampere range. The root-mean-square noise of commonly used amplifiers is close to the same order of magnitude, creating a relatively low signal-to-noise ratio.[23,24] Thus, a preprocessing filter is nearly always a necessity. However, if the effective filter frequency of the event reaches the same order of magnitude as the residence time of an event, we observe that the blocked current appears reduced, while the residence time appears elongated. While the determination of the event localization can also be considered a preprocessing step in many applications, it does not transform the data and therefore requires less accuracy.

Here, we present a universally applicable method for the characterization of short-lived events, where the required sampling frequency to correctly capture the amplitude of the current blockade for each event can be calculated and verified. Herein, events are characterized based on a generalized normal distribution function (gNDF), which allows the description of events based on the open pore current (IO), difference caused by the blocked pore (ΔIB), position in time (μ), standard deviation (σ), a shape parameter (β), and the noise component. Additionally, this method allows the estimation of an ideal range of low-pass filter frequencies that contain maximum information with minimal noise. Taken together, this may provide a direct correlation between different data sets, independent of the measurement setup used.

Results and Discussion

Assumptions for Nanopore Event Detection

Generally, events are described as the relative blocked current [i.e., the residual current (Ires) or its complement to one: the excluded current (Iex)] and residence time (or dwell time). The description of events based on the residual current, as a direct result of analytes translocating pores, is subject to several assumptions:

Assumption 1

The diameters of the pore and analyte, conductivity of the buffer solution, and length of the pore remain constant (Supporting Information 1).[25]

Assumption 2

Because the nanopore acts as a resistor–capacitor circuit (RC circuit), we assume that the residence time of peptides inside the nanopore must be larger than the required time to charge the RC circuit (expressed as the RC time constant, or τsystem). We typically observe a time constant in the microsecond range. This can be corrected under assumption 1 and expressed as eq where Ires (t) is the residual current at time point t and Ires is the residual current, equal to the fraction of the difference in current caused by the block ΔIb and the open pore current Io. It is important to note that one can only measure the residual current at a given time point t. If this assumption is not satisfied, an increase in the residual current should be observed. However, the dwell time of an event should remain (roughly) constant.

Assumption 3

The current blockade does not have any bandwidth limitation. Current blockades are, however, always filtered as we observe a physical translocation through the nanopore, which is not instant. Thus, we describe short events as ideal pulses, which, in the process of the measurement, are dilated by filtering effects. We utilize the term filtering effects here to describe all phenomena causing divergence from an ideal pulse. Additionally, it is important to note that the effective sampling frequency is not always equal to the sampling frequency set during acquisition as events may be dilated by physical events, such as aliasing, analog-to-digital converter limitations, and the capacitance induced by the system. When the event duration reaches or is briefer than the effective sampling period, the event is broadened, and the dwell time, and therefore the excluded current, becomes probabilistic. Effectively, this means that the observed current is lower than the actual current block caused by the analyte, while the observed dwell time is increased.

We compare three representative algorithms to understand several artifacts that arise from these assumptions. The first algorithm is the Fetchan algorithm as implemented in pCLAMP (Molecular Devices), which considers assumptions 1, 2, and 3 to hold. For long-lived events, we notice that all assumptions hold. However, when the residence time of peptides is short, assumption 3 is not fulfilled. This is observed as an excluded current “tail”, where the excluded current seems to decrease while the dwell time is increased, which is biased toward the minimal observable residence time and therefore seems centered (Figure C). We visualize this effect using simulated data by application of a commonly used 4-pole (digital) Bessel filter.[26] When a 10 μs pulse (100 kHz) is filtered to 50 kHz, we notice that the event area remains equal, while the dwell time increases and the excluded current decreases (Figure D). When a 20 μs pulse (50 kHz) with half the excluded current undergoes the same 50 kHz filter, the excluded current is not decreased (Figure E).

The adaptive time-series analysis (ADEPT) by Balijepalli et al. considers assumption 1 and corrects for the capacitance in the system.[25] In ADEPT, the excluded current is represented as a Heaviside step function multiplied by the rise and fall time based on the RC time constant and considers that multiple pulses may be present in each blockade. Notably, the Heaviside step implicates that assumption 3 holds; however, the authors consider that the Heaviside step function could be replaced in future.

Lastly, we investigated the second-order-differential-based calibration (DBC) method,[27] as proposed by the Long group, which is able to accommodate any shape of an event[27] and, by extension, seems independent of assumptions 1, 2, and 3. The method attempts to correct the dwell time of a blockade based on the area under the event, which is (for the most part) not dilated. The DBC method assumes that the event has reached the full height. However, it is difficult to validate this assumption, as we exemplified for 10 μs and 20 μs pulses, which appear equal when they are filtered to 10 kHz (Figure D,E); therefore, the accurate determination of the dwell time and excluded current becomes virtually impossible.

Event Characterization

We hypothesized that a method combining ADEPT with DBC would allow event characterization relying merely on assumption 1. We sought to implement a well-defined function that assumes the shape of any event. To allow the full description of events with five parameters leaving only the noise component, we propose the use of a gNDF, f(x; IO,ΔIB,σ,β,μ), for the optimization of events (eq )with x, the dependent variable (e.g., the time), IO, the open pore current, ΔIB, the difference in current of the event relative to IO, β, the shape parameter, σ, the standard deviation, and μ, the location. In essence, the probability function is multiplied by the blocked pore current (ΔIB) to signify the observed event. The probability density function (PDF) of the gNDF is well characterized and shown in eq

The CDF is given by eq :with Γ, the gamma function, and γ, the unnormalized incomplete lower gamma function. The quantile function (inverse of the CDF) is given by eq with p, the probability. Finally, the full width at half-maximum (fwhm) is given by eq

The resulting function assumes the shape of an ideal pulse, equal to the ideal pulse of a Heaviside step function when the exponent β reaches infinity while displaying a Gaussian profile when this parameter equals 2. In previous contributions, we utilized a similar function for the detection and exclusion of events based on their shape;[12,13] however, it is important to note that the gNDF in this contribution follows a slightly modified version.

The length of events resulting from protein or peptide translocation is one of the most important features that can be extracted from a signal. Of major concern, for short-lived events, is the skewness in the estimated dwell time due to filtering effects. An interesting approach taken by Long and co-workers utilizes the area under the signal as this is almost fully independent of the effective sampling period.[27] However, the total area encompassed by an ideal event is assumed to be equal to the event length multiplied by the event height. As the dwell time approaches the effective sampling period, pulse dilation is expected such that the observed event height is reduced and the observed event length is increased.[26,28] In order to derive the characteristics of an event, we first consider that it reaches the true amplitude of the pulse. In this scenario, we can determine the dwell time of an event by reconstructing the ideal pulse. We know that the area under the event from the ideal pulse should be equal to the area under the fitted event, as expressed in eq and visualized in Figure A.where Δt is the total event time, h is the event height, and f(t) is the probability density of the gNDF at any given time t.

Figure 2

Conceptual representation of events described in this article. Each figure shows a gNDF as a continuous solid black line with β = 3, σ = 1, and μ = 0. The black dashed lines represent the reconstructed ideal pulse event. The gray dashed line represents the baseline at y = 0. (A) Back-mapped area where the gray patches (2) represent the same area. The arrow indicates that the area after the event has the equal area as the missing area of the event’s ideal pulse. (1) indicates the area under the observed event. (B) shows that the dwell time of an event must encompass the localization of the fitted function. The solid vertical line indicates the event localization (μ). (C) shows that the back-mapped curve must overlap with the forward event. The red dotted line represents the inverse of the gNDF shifted by the dwell time, and the solid red part of this line represents the overlap. (D) shows the manual event frequency, where the top red lines represent the probability thresholds (p = 0.001) at the top and bottom of the event. The indicator line with the two solid vertical lines shows the range of the curve where the difference in the cumulative distribution function (CDF) with respect to the localization is equal to the probability threshold.

By utilizing the maximum ordinate of the gNDF, we can derive the total event time (Supporting Information 2), resulting in eq where Δt is the total event time, β is the shape parameter, σ is the standard deviation, and Γ denotes the incomplete gamma distribution.

Importantly, we postulate that eq is only quantitative if the observed event encompassed the localization (μ, Figure B), which can be easily validated (Supporting Information 3). Moreover, we postulate that eq is also only valid if the event area can be explained completely by the back-mapped area. We can validate this by comparing the probability at half-width at half-maximum with the complement to one equidistant probability after t1 (Supporting Information 4, Figure C). Interestingly, this provides a singular, nonanalytically solvable, probability threshold for the detection of t0. However, the true event height must be validated differentlywhere Δt is the total event time, μ is the localization, and Q(p) denotes the quantile function at probability p (eq ).

If there exists a value for p that satisfies eq , we can determine the quantitative limit of detection. As a result, we can determine the time an event is significantly occupying its current plateau (maximum ordinate). Thus, we calculate the effective frequency as the minimal observed change using eq (Supporting Information 5, Figure D).where Fs,event is the effective frequency as the minimal observed change. It is important to note that the Fs,event is not necessarily related to the dwell time. Rather, it is a function of the shape of the observed event and is not defined when β reaches infinity. The gain (or loss) of the event amplitude can be calculated from the required sampling frequency of events (i.e., Fs,event);[29] however, as a rule of thumb, it is easiest to state that the sampling frequency must be equal or larger than the required sampling frequency of the event. Importantly, the second-order DBC method,[27] as proposed by the Long group, is able to assume any shape of an event.[26] However, it is difficult to validate the assumption that the maximum current blockade of an event is reached. Fortunately, Fs,event can be used to validate the DBC method and a possible conjunctive pipeline could be constructed based on a combination of the two to correct for (slight) variations in the shape of dilated events.

Additionally, usage of the fwhm (eq ) is also justified for shorter events as the maximum difference between the integrated back-mapped dwell time is less than 1% when β reaches 2. These deviations are typically within the error margin of the fit to the data and therefore become virtually indistinguishable. We show the difference between the fwhm and back-mapped dwell time in Figure S1. Ultimately, this agrees with previous contributions, where the event dwell time is defined as starting on the vertical rising edge and ending on the descending edge of the signal pulse. The method presented in this contribution states that this assumption is correct if the shape parameter is large. It also explains the tail-like widening dwell time distribution seen for short pulses analyzed using “standard” single-channel search as the uncertainty in the estimated dwell time is increased. This effect is especially notable when used in, for example, pCLAMP’s Fetchan (Figure B) or the “rapid event detection” (RED) method, which we introduced in an earlier contribution.[30]

Functional Implementation

We implemented the gNDF for fitting using a custom Python library (see code availability) suitable for the analysis of Axon Binary Files; however, other data loaders can be easily added. We implemented RED as described in the previous contribution[12,13,30] and utilized it for event localization. Subsequently, we fit the gNDF around each observed event, after which the event features are determined as described in this contribution. All events are stored in a SQL file containing the start and end times of each event, the extracted features, and utilized fitting functions. The optimized results are stored in a separate table from RED, allowing validation of the results.

Effect of Over- and Under-Filtering

The over-filtering of events causes the observed excluded current to be lower than the expected value, while the observed dwell time is increased (Figures and 3A). To minimize this artifact, we can localize events using a low frequency cutoff and characterize events with minimal filters. It may seem as if we can directly utilize unfiltered data; however, under-filtering may cause the observed excluded current to also be lower than the expected value as the baseline becomes indistinguishable from an event due to noise.[31] We observe this effect when comparing the residual current histograms at different filter frequencies (Figure ). Therefore, we postulate that there is a low-pass filter that can gain maximum information while minimizing noise. We pose that eq holds under the assumptions that the origin of events can be described by immediate (ideal) pulses, all filter effects cumulatively result in a Gaussian-like shape, eq is satisfied, and the event sampling frequency (eq ) is equal or larger than the effective sampling frequency. We observe that Fs,event as calculated by eq appears larger for signals that have been filtered at or under Fs,event and remains stable with higher frequency filters. If the filter frequency is lower than Fs,event, we observe that the estimated dwell time is elongated (Figure A). When the effective filter frequency is (nearly) equal to or larger than Fs,event, we observe the true dwell time (Figure B–E). This explains why a 5 kHz filter on the presented data results in a correct excluded current spectrum and dwell time, while the Fs,event appears larger. It is worth noting that for under-filtered events, while the excluded current changes due to baseline noise, the dwell time stays as expected.

Figure 3

Optimized event detection using different cutoff filters. All panels represent the same data where a 10 μM (added to cis) of penta-peptide (YAGFL) was measured in 1 M KCl, pH 4.2 buffered using 15 mM citric acid, and bis-tris-propane under an applied potential of negative 100 mV (cis–trans) at a sampling frequency of 500 kHz with a recording bandwidth of 100 kHz using fragaceatoxin C modified with a tryptophan at position 13 (G13W-FraC). Experiments were performed using an Axon Axopatch 200B amplifier (Molecular Devices) coupled to an Axon Digidata 1550B (Molecular Devices). All events were localized using a 5 kHz Gaussian filter and subsequently characterized using a Gaussian filter at the filter frequency as described above each graph. (A–E) Top graph: the excluded current (Iex) of each observed event (Fs,event ≤ 500 kHz) binned in 200 residual current bins evenly distributed between 0 and 1 Iex. Middle graph: minimum required sampling frequency of events (Fs,event) set against the Iex for different Gaussian filter frequencies (see the title in each panel). The dotted line represents the filter frequency. Bottom graph: Iex set against the integrated dwell time as estimated by eq of this article for different Gaussian filter frequencies (see the title in each panel).

Conclusions

In this work, we presented a generally applicable function that can accurately describe data acquired by nanopore sensors with only five parameters with noise. We also show a novel method for the determination of the minimal required sampling frequency, which allows a reproducible data analysis pipeline, with minimal user bias.

For this, we use a gNDF. While the gNDF is a continuous function, we demonstrated that the events are finite. We showed that the gNDF can be used to determine the minimum required sampling frequency for an event to reach its maximum current plateau. Therefore, we have shown that the dwell time of events can be accurately determined using this method under the assumptions that the diameter of the pore and analyte, conductivity of the buffer solution, and length of the pore remain constant. While other approaches may be beneficial depending on experimental conditions, we presented a unified approach for the determination of nanopore event characteristics.

Materials and Methods

Data Recording

Recordings of ionic currents were obtained using an Axopatch 200B (Axon Instruments) combined with a Digidata 1550B A/D converter (Axon instruments), similar to the preceding work.[12,13] The sampling frequency was set at 500 kHz for analyte recordings, the analog Bessel filter was set at 100 kHz, respectively. Data was recorded using Clampex 10 (Molecular Devices).

Data Analysis

Data was analyzed using Jupyter Notebook (version 5.5.0) running with Python 3.6.5 (64-bit), both within the Anaconda (version 5.2.0) environment. Events were localized using a RED algorithm.[12,13,30]

Data Availability

The authors declare that the data and code supporting the findings of this study are available within the article and its Supporting Information or from the corresponding authors upon reasonable request.

Code Availability

The code is available under DOI: https://doi.org/10.5281/zenodo.6546320.