Single-molecule protein sequencing through fingerprinting: computational assessment

We simulated 2000 read-outs, each for a different protein. The proteins are randomly picked from the database and thus contain random amino acids and fingerprint lengths. Next, to assess the robustness of the method against inaccuracies that are expected from actual experiments, errors are iteratively introduced for each read-out up to the error level we want to investigate.

We expect that actual data will be convoluted with poor dye-labeling, photoblinking and photobleaching of dyes, local structures of a substrate protein, non-uniform speed of substrate translocation, proximity between dyes etc. The poor labeling, photoblinking, and photobleaching of acceptor dyes will appear as deletion errors (figure 4). The non-uniform speed of translocation will introduce insertion and deletion errors to CK-dist fingerprinting. The proximity of acceptor dyes will bring deletion and transposition/substitution errors. If a donor dye is photobleached during a measurement, it will appear as a truncation error. We do not consider this error for fingerprinting analysis since donor photobleaching can be determined from SM time traces and thus can be easily excluded from further analysis. Other complications, such as aggregation of denatured proteins, may also be expected but are not considered in our analysis. See a pseudo-code for simulating these errors (supplementary information).

Zoom In Zoom Out Reset image size

In figure 3, we investigated one combination of errors (70% deletions, 20% insertions, 10% transpositions) for CK fingerprinting, in which we assigned the largest percentage to deletions since this error is the most likely to occur (poor labeling of acceptors due to incomplete denaturation of proteins, photoblinking/photobleaching of acceptors, and presence of consecutive identical acceptor fluorophores). For CK-dist fingerprinting analysis, we considered the same combination of errors as for CK fingerprinting but with errors at CK residues and errors of the distance between CK residues equally likely to occur. In supplementary figure 3, we expanded the error space that we explored and obtained trends nearly identical to that found in figure 3(a).

The 2000 simulated readouts are searched for in the database, and the numbers of true positives and the number of matches are recorded. To examine the performance variability of our algorithm in retrieving proteins using fingerprints, three independent repetitions are executed. In each repetition, detection precision (P) (figure 3) and detection recall (R) (supplementary figure 4) are calculated based on the outputs. P is defined as the number of true positives divided by the number of read-outs returned by the algorithm. R is the number of true positives divided by the number of conditional matches.

The inputs to our method are a reference database $R$ containing fingerprint representations of protein sequences, a query fingerprint $Q$ and an error level $\alpha .$ The alphabet is ${\rm{\Sigma }}=\{C,\;K\},$ since we only compare fingerprints of these two amino acids. Let ${L}_{Q}$ be the query length and ${R}_{x}\in R$ denote the $x$th reference sequence in the database $R$ with length ${L}_{x}^{R}.$ The distance $S({R}_{x},\;Q)$ between a reference fingerprint ${R}_{x}$ and a query $Q$ is the minimal number of steps required to transform $Q$ into ${R}_{x}$. Formally, given $Q,$ $R,$ and $\alpha ,$ the problem is to find all ${R}_{x}\in R$ for which $S({R}_{x},\;Q)$ is smaller than $k$ = $\alpha \;\times \;{L}_{Q}$.

Given the inputs, the algorithm takes two steps to retrieve matches: (1) a filtration strategy is applied to identify candidate sequences in $R;$ and (2) a verification method is employed to examine all candidates for possible matches.

Dynamic programming is computationally costly, prohibiting direct application on large databases in a high-throughput setting [13]. In order to reduce the running time without affecting sensitivity, we use filtration to remove those references that definitely cannot match the query fingerprint $Q$ with distance smaller than or equal to $k.$ Filtration exploits the fact that it is easier to tell a reference fingerprint that does not match a query fingerprint than to tell one that does match. Typically, it uses a simple and highly efficient filter criterion to analyze the reference sequences, leaving only a small number of ${R}_{x}'s$ for further (more expensive) analysis. We devised a new filtration method combining two existing algorithms, partial exact matching and $n-{\rm gram}$ counting.

In partial exact matching, the query fingerprint $Q$ is divided into $(k+1)$ pieces ${q}^{0},\;{q}^{1},\ldots ,\;{q}^{k},$ where $k$ equals $\alpha \times {L}_{Q}.$ For a match to be possible, there must be at least one piece that appears exactly in a reference sequence ${R}_{x}$ [15]. If this is not the case, ${R}_{x}$ is discarded.

A faster filtration method is $n-{\rm gram}$ counting, which compares the $n-{\rm grams}$ of two fingerprints. An $n-{\rm gram}$ [16] on the alphabet set ${\rm{\Sigma }}=\{C,\;K\}$ is any string in ${{\rm{\Sigma }}}^{n},$ where ${{\rm{\Sigma }}}^{n}$ is the set of all possible strings of length $n$ over ${\rm{\Sigma }}.$ For example, the $2-{\rm grams}$ for ${\rm{\Sigma }}=\{C,\;K\}$ are CC, CK, KC and KK. The $n-{\rm gram}$ distance is defined as the sum of the absolute differences between the numbers of occurrences of each $n-gram.$ If the $n-{\rm gram}$ distance exceeds $2nk,$ ${R}_{x}$ is discarded [16].

We combined the partial exact matching and $n-{\rm gram}$ counting approaches to decide whether there exists at least one piece in $Q$ that appears with a limited amount of errors as a piece of ${R}_{x}$ [17]. The distance function between two pieces of $Q$ and ${R}_{x},$ ${q}^{j}$ and ${r}_{x}^{j},$ based on their $n-{\rm grams}$ was defined as:

$$\begin{eqnarray*}{S}_{npm}\left({r}_{x}^{j},\;{q}^{j}\right) & = & \;\displaystyle \sum _{\nu \in {{\rm{\Sigma }}}^{n}}{\rm max}\left(G\left({q}^{j}\right)\left[\nu \right]\right.\\ & & \left.-G\left({r}_{x}^{j}\right)\left[\nu \right],\;0\right),\end{eqnarray*}$$

where $\left[\nu \right]$ is an $n-{\rm gram}$ and $G\left({q}^{j}\right)\left[\nu \right]$ $\text{and}\;G\left({r}_{x}^{j}\right)\left[\nu \right]$ denote the total number of times $\left[\nu \right]$ occurs in ${q}^{j}$ and ${r}_{x}^{j},$ respectively.

For each piece ${q}^{j}$ in the query, the corresponding piece ${r}_{x}^{j}$ contains the same letters in the reference sequence with an additional $k$ letters on both sides, as shown in figure 5. It is sufficient to compare the ${r}_{x}^{j}$ in the reference with the ${q}^{j}$ in the query to determine whether the piece ${q}^{j}$ appears in the reference ${R}_{x},$ since $k$ errors cannot alter more than $k$ positions. Since a query piece is searched in a limited range in the reference, it can discard more entries in the reference database than the partial exact matching method, in which the ${q}^{j}$ is compared with the entire reference sequence.

Zoom In Zoom Out Reset image size

The distance between a piece ${q}^{j}$ in query $Q$ and the corresponding piece ${r}_{x}^{j}$ in ${R}_{x}$ is computed to determine whether ${R}_{x}$ is a candidate match. For each ${q}^{j}$ and its corresponding ${r}_{x}^{j},$ we check whether any $n-{\rm gram}$ occurs more often in ${q}^{j}$ than in ${r}_{x}^{j}.$ If not, the ${S}_{npm}({r}_{x}^{j},\;{q}^{j})$ is zero, i.e. the n-grams in ${q}^{j}$ appear exactly in ${r}_{x}^{j}$ . Only if for at least one ${q}^{j},$ ${S}_{npm}({r}_{x}^{j},\;{q}^{j})\;$ is zero, ${R}_{x}$ is kept as a candidate.

The remaining candidate matches are examined by a global alignment dynamic programming approach considering a number of possible error types. In our analysis, four types of error may occur: deletion, insertion, mismatching an amino acid with another one (substitution), and swapping (transposition).

The dynamic programming algorithm is designed to provide the optimal gapped alignment between two sequences, i.e. an alignment with long regions of identical amino acid pairs and very few mismatches and gaps [14]. As the sequences become more dissimilar, more mismatched amino acid pairs and gaps should appear. To find the optimal alignment, a dynamic programming matrix $M$ first needs to be calculated. Each element ${M}_{i,j}$ represents the maximum score of aligning the substrings $Q[\mathrm{1..}.i]$ and ${R}_{x}[\mathrm{1..}.j].$ Let $c$ denote the scores of the four operations. The base cases, ${M}_{0,j}$ and ${M}_{i,0},$ are defined as (${c}_{{\rm{del}}}\times j)$ and (${c}_{{\rm{ins}}}\times i)$ for all $1\;\leqslant \;j\;\leqslant \;{L}_{x}^{R}$ (length of ${R}_{x})$ and $1\;\leqslant \;i\;\leqslant \;{L}_{Q}$ (length of $Q)$ respectively. Then, considering the four possibilities, $M$ is updated using the following recursive relation

$$\begin{eqnarray*}&&{M}_{i,j}=\;{\rm max}\left\{\begin{array}{l}\begin{array}{c}{M}_{i-1,j-1}+{c}_{{\rm{sub}}},\\ {M}_{i-1,j}+{c}_{{\rm{ins}}},\end{array}\\ \begin{array}{c}{M}_{i,j-1}+{c}_{{\rm{del}}},\\ {M}_{i-2,j-2}+{c}_{{\rm{swap}}}.\end{array}\end{array}\right.\end{eqnarray*}$$

The score for each operation is set based on the estimation of how likely each error is to occur in our measurements. Currently, deletions caused by low labeling efficiency are the dominating errors, followed by insertions, transpositions and substitutions (i.e. matching C to K or vice versa) (see error simulation). Hence we choose a relatively low penalty (negative) for deletions and higher penalties (negative) for transpositions and substitutions. For the matching positions, the score is positive (see supplementary table 2).

By memorizing the solutions to the subproblems for $1\;\leqslant \;j\;\leqslant \;{L}_{x}^{R}$ and $1\;\leqslant \;i\;\leqslant \;{L}_{Q}$ stored in the dynamic matrix, we can recursively compute the maximum score of aligning ${R}_{x}$ and $Q.$ Therefore we find the score of the optimal alignment of the two sequences starting from the maximum value in the last row or last column. We maintain a matrix of traceback pointers in the recursion, so that we remember which case was used to calculate every cell ${M}_{i,j},$ allowing to reconstruct the optimal alignment.

From this alignment the numbers of errors for different types as well as the total number of errors can be calculated. The distance between the query and the reference $S({R}_{x},Q)$ is defined as the total number of errors. If this distance is smaller than $k,$ the reference sequence ${R}_{x}$ is considered as a match. Otherwise, it is not a match of the query sequence within the error bound $k.$ A match is considered a true positive match when the match is the exact query protein. If a match has the same fingerprint but a different amino acid sequence, it is not considered to be a true positive match. In our analysis, this is determined by checking the protein accession codes.

Additional information, such as the distance between two read-outs, can be deduced from the measurements. This distance is the space between two labeled amino acids, which is the number of non-labeled amino acids in between, which show a different pattern in the measurement. For this to be estimated reliably, proteins will have to be sequenced at a relatively constant speed, an assumption which is not a priori valid. From the sequencing signals, we cannot easily determine the start or the end of proteins in the time trace if they do not correspond to a labeled amino acid. Thus, the starting and ending non-labeled amino acids are not included when we construct the fingerprint with distance information.

This distance information is added to the original CK fingerprints using an additional symbol (say, 'o'), occurring multiple times (representing the length of distance). Next, two distances between query and reference are calculated to examine whether a reference sequence is a match. One is the $S\left({R}_{x},Q\right)$ between fingerprints with distance information, the other the $S^{\prime} \left({R}_{x},Q\right)$ between CK fingerprints only. Reference sequence ${R}_{x}$ is considered a match if and only if $S^{\prime} \left({R}_{x},Q\right)$ is smaller than $k\text{'}\;=\;(\alpha \;\times \;{L}_{Q}^{\prime })$ and $\;S\left({R}_{x},Q\right)$ is smaller than $k\;=\;\alpha \;\times \;{L}_{Q},$ where ${L}_{Q}$ is the length of the query CK fingerprint, ${L}_{Q}^{\prime }$ is the length of the query fingerprint with distance information and $k^{\prime}$ and $k$ represent the numbers of errors allowed. Experimental error on the distance information is also taken into consideration.

Note: Supplementary information is available in the online version of the paper. An animated experimental scheme is available at https://youtube.com/watch?v=YpWCCWO5q10.