\chapter{Methodology and Research Methods}
\label{ch:research}

The aim of this chapter is to explore in detail a particular DNA nanopore signal processing pipeline developed at the Computational Biology Group at Comenius University Bratislava. The purpose of this pipeline is to produce an alignment of the raw signal data to a known reference genome. Such alignment can be then used for various other applications~\cite{boza_precise_2023}.

As seen in~\cref{fig:alignment-pipeline}, the processing is done in several steps, which will be discussed further in this chapter. First, raw signals from the nanopore sequencer have to be acquired and stored. The raw signals are then converted into a DNA sequence, which is then aligned with a reference genome. However, this obtains only an approximate alignment, which is then further improved by a custom alignment algorithm.

\begin{figure}
    \centering
    \input{figures/30-alignment-pipeline}
    \caption{A diagram of the nanopore signal alignment pipeline.}
    \label{fig:alignment-pipeline}
\end{figure}

The pipeline itself was developed for use with the previous generation of nanopore sequencers. As the new generation of nanopore sequencing technology became gradually available, many of the software tools developed for this generation became obsolete and were replaced by new ones.

Unfortunately, many bioinformatics pipelines (including this one) depended heavily on these tools. This makes it increasingly difficult (or impossible) to use them with new data obtained from next-generation sequencers.

\section{Principles of Nanopore Sequencing}
\label{sec:nanopore}
The first step of the pipeline is to obtain the data that will be processed. In this case, that means using a DNA sequencer---a device that can determine order of the nucleotide bases of a given DNA. A DNA strand consists of nucleotides, which are composed of one of four different bases: \emph{adenine}, \emph{cytosine}, \emph{guanine}, or \emph{thymine}\todo{citation needed}.

Currently, there are several next-generation sequencing methods in active use and development. One of the most notable methods is the nanopore sequencing developed by Oxford Nanopore Technologies, mostly thanks to its ability to produce very long reads of the strands, its speed, portability, and low cost\todo{citation needed}.

The DNA strand passes through a small barrel-shaped protein, called a nanopore, embedded in a membrane (\cref{fig:nanopore-sequencing}). During this process, a special motor protein is attached to the DNA, which regulates the speed of its translocation. By applying voltage across the membrane, it is possible to measure the ionic current passing through the nanopore, which changes as DNA passes through the nanopore~\cite{wang_nanopore_2021}.

\begin{figure}
    \centering
    \includegraphics[width=0.5\textwidth]{figures/30-nanopore.pdf}
    \caption[Basic principle of nanopore sequencing.]{Basic principle of nanopore sequencing used by Oxford Nanopore Technologies. A DNA strand is moving through the nanopore, which is used to measure electrical current passing through it.}
    \label{fig:nanopore-sequencing}
\end{figure}

A sequence of raw signal measurements (also known as \emph{squiggle}) can be then used to determine nucleotide bases through a process called \emph{basecalling} (also discussed in~\cref{sec:basecalling}). Another common practice is to analyze the raw signal data directly (or in combination with basecalling) in an attempt to utilize the additional information present in the signal, or to skip computationally intensive basecalling~\cite{boza_precise_2023}.

A big advantage of nanopore sequencing technology is that the same principles can be used for direct RNA sequencing, whereas other sequencing technologies have to rely on complementary DNA (cDNA) that is synthesized through reverse transcription of the RNA. Such approach has several disadvantages, for example inability to sequence non-coding RNA~\cite{wongsurawat_native_2022}.

\section{File Formats for Storage and Processing of the Nanopore Data}
\label{sec:formats}
During (not just) nanopore sequencing, a large amount of various data is acquired or produced. As such, there are many different data formats designed for different kinds of data, and any data processing pipeline will inherently need to utilize multiple such formats.

This section will describe some of the most used file formats that are relevant to our pipeline---FAST5 for storage of raw nanopore signals, FASTA for representing any nucleotide sequence, and finally, SAM for storing alignments of sequences.

\subsection{Raw Nanopore Signal}
\label{sub:formats-fast5}
The measurements produced by a nanopore sequencer are transmitted to a computer for further processing, either in real time or asynchronously. However, even if the processing is performed in real time, it is often useful to store the raw signals itself. This makes it possible to perform a different kind of analysis or to repeat the analysis later, with more advanced tools\todo{citation needed}. Furthermore, storing the raw data is necessary for development of new processing techniques.

The amount of data produced can be quite large, as even a single MinION flow cell contains 512 nanopore channels, each of which can simultaneously read a different strand at speeds up to 450 bases per second~\cite{ni_benchmarking_2023}. To store such amounts of data and to facilitate fast and efficient access, an appropriate data storage format has to be used.

Until recently, all Oxford Nanopore Technologies sequencers used the custom FAST5 format to store the raw signals of reads~\cite{gamaarachchi_fast_2022}. The FAST5 format is based on the Hierarchical Data Format (HDF5) designed for efficient storage of large datasets in a structured way.

The HDF5 file logical data model organizes the data using \emph{groups} and \emph{datasets}. A \emph{dataset} is a logical representation of the data itself, with attached metadata such as the data type, size and other properties. Each \emph{group} can store both datasets and other groups, effectively organizing the data into a tree-like structure~\cite{the_hdf_group_hdf5_2024}.

The FAST5 format is essentially just a definition of the structure of data contained in a HDF5 file. Originally, a single FAST5 file contained just a single read of the nanopore, but the specification was later changed to allow multiple reads in a single file. This is currently the only supported option, with the single-read files being officially deprecated and unsupported.

Similarly, the FAST5 format was originally designed to contain not only the raw data from nanopore sequencer, but also results of other analyses performed on this data, such as basecalling. The use of this option is now discouraged, and the ONT tools have removed the support to write the results to FAST5~\cite{oxford_nanopore_technologies_guppy_2018}.

In case of the pipeline described in this chapter, the FAST5 files produced by the sequencer are expected to contain only the raw nanopore signal data of multiple reads. These raw signals are stored in FAST5 as a series of \qty{16}{\bit} integers representing quantized values from the analog-to-digital converter of the sequencer. These values can be converted to actual current using the following relation~\cite{gamaarachchi_fast_2022}:
%
\begin{equation}
    I_j = (s_i + o) \frac{r}{d},
\end{equation}
%
where \(I_j\) is the electrical current of the \(j\)-th sample in \unit{\pico\ampere}, \(s_j\) is the quantized value of the \(j\)-th sample, \(o\) is the offset error of the analog-to-digital converter, \(r\) is the measurement range in \unit{\pico\ampere}, and \(d\) is the number of quantization levels used. The values of \(o\), \(r\), \(d\) are stored together with the raw signal data in the FAST5 file.

\subsection{Nucleotide Base Sequences}
\label{sub:formats-fasta}
FASTA is a simple plain text-based format for the storage of nucleotide sequences. It consists of two parts: a single line starting with the \texttt{>} character containing the description of the sequence and the sequence itself stored as a string of letters \texttt{A}, \texttt{C}, \texttt{G}, and \texttt{T} (or \texttt{U} in RNA) on multiple lines (typically aligned to 80 characters)~\cite{mount_bioinformatics_2001}.

\lstinputlisting[caption={An example of a sequence stored in the FASTA format.},label={lst:melanogaster-fasta}]{code/30-melanogaster.fasta}

Over the years, FASTA became the de facto standard for handling nucleotide sequences, whether as a storage format for reference genomes, input format for sequence alignment, or querying genome databases. Its popularity, simplicity, and lack of any formal standardization have naturally led to the creation of many modifications and variants~\cite{cock_sanger_2010}.

One such variant is the FASTQ format, which extends the basic FASTA format with a sequencing quality score of each base. This is achieved by including a \emph{quality line} for each line of the sequence, which contains the same number of characters. These characters represent different estimated probabilities of error, which serves as a quality score for each base~\cite{cock_sanger_2010}.

\subsection{Sequence Alignment}
\label{sub:formats-sam}
One of the most common tasks in bioinformatics is the alignment of sequences (discussed further in~\cref{sec:mapping}. The primary goal is to rearrange one or more sequences so that shared subsequences are aligned. One of the most used formats for storing results of alignment is the Sequence Alignment Map (SAM) format and its compressed binary equivalent BAM.

SAM is a plain text-based format that consists of header lines containing various (optional) metadata, and multiple alignment lines, each describing one linear alignment of a segment to reference. Each line consists of multiple tab-separated fields, the most important being the position of the leftmost mapping and a \emph{Compact Idiosyncratic Gapped Alignment Report} (CIGAR) string representing the alignment~\cite{the_sambam_format_specification_working_group_sequence_2024}.

\begin{table}[h]
    \centering
    \begin{tabular}{ccc}
        \toprule
        Op & Description \\
        \midrule
        M  & alignment match \\
        I  & insertion to the reference \\
        D  & deletion from the reference \\
        N  & skipped region from the reference \\
        \bottomrule
    \end{tabular}
    \caption{Overview of basic operations in CIGAR strings.}
    \label{tab:cigar-ops}
\end{table}

The CIGAR string is a sequence of operations represented by letters (see~\cref{tab:cigar-ops}), prefixed by a number indicating how many times the operation is repeated. For example, \texttt{6M2D3M} represents a segment with six bases exactly matched to the reference, two deletions from the reference, and another three exact matches~\cite{the_sambam_format_specification_working_group_sequence_2024}.

As mentioned earlier, the binary counterpart to SAM is the BAM format, which has the same internal structure as SAM, but with added gzip compression. To enable fast random access, a BAI file with an index can be created~\cite{the_sambam_format_specification_working_group_sequence_2024}.

\section{Basecalling of the Nanopore Signal Data}
\label{sec:basecalling}
As explained in~\cref{sec:nanopore}, the direct output of nanopore sequencing is a sequence of electrical currents measured across the nanopore, affected by the DNA or RNA strand passing through it. In order to obtain a sequence of nucleotide bases, the raw signal has to be processed first. This process is commonly referred to as \emph{basecalling}.

While sequenced bases are not the desired output of this pipeline, basecalling is used internally for alignment and mapping performed by the Oxford Nanopore Technologies tool \emph{Megalodon} as the next step of the pipeline, described in~\cref{sec:mapping}~\cite{oxford_nanopore_technologies_megalodon_2019}.

The nature of the raw signal data makes basecalling (or any other kind of processing) a rather complex task. In addition to the usual issues such as noise in the signal, there are numerous other ones caused by the sequencing process itself.

Notably, a single measurement does not correspond to just one nucleotide, but to a whole sequence of nucleotides (also called \emph{k-mer}) that are currently in the nanopore and affect the measured current~\cite{rang_squiggle_2018}. The size of the k-mer depends on the used nanopore chemistry and has slowly increased over the years, with 9-mer models currently being used for basecalling~\cite{oxford_nanopore_technologies_nanoporetechkmer_models_2024}.

Furthermore, the motor protein attached to the DNA cannot maintain a completely constant translocation speed, while the sampling frequency remains the same. As a result, signal segments of the same duration do not correspond to the same number of nucleotide bases~\cite{rang_squiggle_2018}.

Another source of complexity are chemical modifications of nucleotides, such as methylations, where a regular nucleotide is modified by attaching a methyl group. This does not change the resulting DNA sequence, but the different chemical properties of these nucleotides can affect the resulting signal~\cite{rang_squiggle_2018}. On the one hand, this makes basecalling much more difficult, but the additional information in the raw signal makes it possible to develop novel methods that can detect these modifications~\cite{boza_precise_2023, liu_dna_2021}.

Even though early basecallers used hidden Markov models~\cite{david_nanocall_2017}, the complexity of the task quickly led to the adoption of more powerful approaches such as deep convolutional and recurrent networks~\cite{boza_deepnano_2017,teng_chiron_2018}. Currently, the research focuses on transformer-based basecallers, which achieved even higher accuracy~\cite{oxford_nanopore_technologies_nanoporetechbonito_2024}.

\begin{figure}
    \centering
    \includegraphics[width=0.9\textwidth]{figures/30-flip-flop.pdf}
    \caption[The flip-flop basecaller architecture.]{The flip-flop neural network architecture used in the Guppy basecaller~\cite{wick_performance_2019}.}
    \label{fig:flip-flop}
\end{figure}

Although the architecture of models used in the Guppy basecaller was always based on convolution and recurrent neural networks, there have been multiple changes in the architecture over time. One of the last architectures used is the \emph{flip-flop} architecture shown in~\cref{fig:flip-flop}, which is an evolution of architectures used previously.

In the flip-flop architecture, the length of the input signal is first reduced by the convolutional layer. Then, a series of Gated Recurrent Unit (GRU) hidden layers is used to process the data. Finally, the probabilities of each base are obtained in the softmax layer, which are then used to create the final resulting sequence of nucleotides.

Despite the generally good accuracy of current basecallers, the performance is not always consistent across different types of experiments and genomes. The neural network models do not always generalize sufficiently and there are many cases when a model trained on a dataset from certain species performs significantly worse on a different species~\cite{wick_performance_2019}.

Another important aspect of the neural network basecallers is their speed and computational complexity. Because nanopore sequencing technology allows for real-time sequencing, a significant effort is put into making the basecaller capable of real-time processing~\cite{wick_performance_2019}. Due to the complexity of the models, this often requires high-performance hardware.

\section{Approximate Alignment of the Raw Nanopore Signal}
\label{sec:mapping}
In the next part of the pipeline, the Oxford Nanopore tool \emph{Megalodon} is used to obtain an approximate alignment of the raw nanopore signal to the bases of the reference genome. As the data obtained from nanopore sequencers can have a quite high error rate, the alignment requires the use of more sophisticated methods that take these errors into account.

The \emph{Megalodon} first performs a regular basecalling using the Oxford Nanopore basecaller \emph{Guppy}, which translates the raw signal data into a sequence of bases. This sequence is then aligned using the \emph{minimap2} aligner~\cite{oxford_nanopore_technologies_guppy_2018}. As the output of the \emph{Guppy} neural network is anchored to the input data, it is then possible to create a direct mapping from the raw signal to reference genome~\cite{oxford_nanopore_technologies_megalodon_2019}.

Most aligners, including \emph{minimap2}, are based on an approach called \emph{seed-and-extend}, first used in BLAST~\cite{altschul_basic_1990}. In this approach, two sequences of nucleotide bases are first divided into (usually overlapping) k-mers, which are then compared across the sequences, obtaining positions of k-mers that are present in both sequences, called \emph{seeds}. This comparison can be performed efficiently by saving all the k-mers and their positions from the first sequence into a hash map, which is then used to look up k-mers from the second sequence~\cite{li_minimap2_2018}.

\begin{figure}
    \centering
    \missingfigure{Seed-and-extend figure}
    \caption[Seed-and-extend alignment of sequences.]{Seed-and-extend alignment of sequences.}
    \label{fig:seed-and-extend}
\end{figure}

As can be seen in~\cref{fig:seed-and-extend}, the seeds (which represent subsequences where the sequences match exactly) are joined to other seeds by expanding the subsequences, thus allowing errors in the joined subsequences. Lastly, subsequences that cannot be expanded are joined by adding spaces to the sequences using dynamic programming~\cite{li_minimap2_2018}. 

The number of possible k-mers grows exponentially with \(k\), so the amount of memory required for the storage of the hash table can be very large for very long sequences and high values of \(k\). A technique called \emph{minimizers} is employed by \emph{minimap2} (and various other aligners) to reduce the number of stored k-mers~\cite{li_minimap2_2018}. This is done by choosing and saving a representative k-mer (usually lexicographically the smallest, or the smallest after applying a hashing function) out of a sliding window of \(w\) consecutive k-mers. The same approach is then used for processing of the second sequence~\cite{roberts_reducing_2004}.

The \emph{minimap2} algorithm also includes numerous improvements and optimizations of this basic process, such as heuristics for expansion and alignment~\cite{li_minimap2_2018}. However, these are out of the scope of this work.

Although this method of alignment is relatively straightforward and has mostly acceptable accuracy, the same kind of complexities that make basecalling a rather difficult task also complicate the alignment. For example, modified bases which often decrese the accuracy of basecalling are one such source of errors~\cite{boza_precise_2023}. Therefore, if a higher quality of alignment is required, it is necessary to use the additional information contained in the raw signal data to improve accuracy.

\section{Realigning the Signal Mappings}
\label{sec:improving-mapping}
The last part of the described signal processing pipeline attempts to further improve the signal mappings obtained in the previous section. The method used in this part of the pipeline is inspired by the approach used in the tool \emph{Nadavca} developed by the Computational Biology Group at Comenius University Bratislava~\cite{boza_precise_2023}.

The method itself is based on \emph{dynamic time warping} (DTW), a dynamic programming algorithm that attempts to align two time series in a way that minimizes the difference between their corresponding points by non-linearly shrinking or expanding (warping) the sequences, as can be seen in~\cref{fig:dtw-alignment}~\cite{muller_information_2007}.

\begin{figure}
    \centering
    \includegraphics[width=0.4\textwidth]{figures/30-dtw-alignment.pdf}
    \caption{Alignment of two sequences using dynamic time warping.}
    \label{fig:dtw-alignment}
\end{figure}

The alignment of the sequences \(s = (s_1, \ldots, s_n)\) and \(r = (r_1, \ldots, r_m)\) can be described by an alignment path \(P\):
%
\begin{equation}
    P = (p_1, \ldots, p_l),
\end{equation}
%
where each tuple \(p_i = (n_i, m_i) \in (1, n) \times (1, m)\) represents the alignment of element \(s_{n_i}\) to element \(r_{m_i}\). The alignment path \(P\) must also satisfy the following conditions\cite{muller_information_2007}:
%
\begin{itemize}
    \item boundary condition: \(p_1 = (1, 1)\) and \(p_l = (n, m)\),
    \item step-size condition: \(p_{i+1} - p_i \in \{(0, 1), (1, 0), (1, 1)\}\) for each successive \(p_i\) and \(p_{i+1}\).
\end{itemize}

The quality of the alignment path \(P\) is given as follows:
%
\begin{equation}
    c_P(s, r) = \sum_{i=1}^{|P|} c(s_{n_i}, r_{m_i}),
\end{equation}
%
where the value of the cost function \(c(s_{n_i}, r_{m_i})\) is small if the elements \(s_{n_i}\) and \(r_{m_i}\) are similar. Then, the optimal alignment path is a path \(P\), for which the value of \(c_P\) is the lowest possible. More formally, the cost of the optimal alignment path between sequences \(r\) and \(s\) is:
%
\begin{equation}
    \text{DTW}(r, s) = \min_{p \in A} \left(c_p\left(r, s\right)\right),
\end{equation}
%
where \(A\) is the set of all alignment paths between \(r\) and \(s\) that satisfy the conditions mentioned above~\cite{muller_information_2007}.

To avoid costly and inefficient computation of all possible alignment paths, a dynamic programming approach can be used. In essence, this solution works by constructing an accumulated cost matrix \(D\) with dimensions \((n, m)\). Each element \(D(i, j)\) has a value that corresponds to the cost of optimal alignment between subsequences \(s(1:i)\) and \(r(1:j)\)~\cite{muller_information_2007}:
%
\begin{equation}
    D(i, j) = \text{DTW}\left( s(1:i), r(1:j) \right).
\end{equation}

The matrix \(D\) can be computed iteratively using the following relation~\cite{muller_information_2007}:
%
\begin{equation}
    D(i, j) =
    \begin{cases}
        c(1, 1)    & i=1 \land j=1 \\
        \infty              & i=1 \lor j=1  \\
        \min(D(i-1, j), D(i, j-1), D(i-1, j-1)) + c(i, j).
    \end{cases} 
\end{equation}

In case of our pipeline, a modified dynamic time warping based algorithm is used to align the raw signal data to the expected signal constructed from the reference genome. First, the raw signal data are normalized with the normalization algorithm from the Oxford Nanopore Technologies research basecaller \emph{bonito}, which can be described by the following equations~\cite{oxford_nanopore_technologies_nanoporetechbonito_2024}:
%
\begin{align}
    \text{shift} &= \max(10, m_\text{shift} (Q_a(s) + Q_b(s)) \\
    \text{scale} &= \max(1, m_\text{scale} (Q_b(s) - Q_a(s))  \\
    s_i' &= \frac{s_i - \text{shift}}{\text{scale}},
\end{align}
%
where \(Q_a(s)\) and \(Q_b(s)\) denote the \(a\)-th and \(b\)-th quantiles of the raw data sequence \(s\). The default values of the parameters used by \emph{bonito} are \(a = 0.2\), \(b = 0.9\), \(m_\text{shift} = 0.51\), and \(m_\text{scale} = 0.53\)~\cite{oxford_nanopore_technologies_nanoporetechbonito_2024}.

The normalized raw data pre-aligned to the reference genome are then used to calculate medians for each 5-mer, which are used to construct expected signal of the reference genome (denoted \(r\)). The alignment itself is performed by a custom version of dynamic time warping, in which the most notable difference is the modified step-size condition: \(p_{i+1} - p_i \in \{(1, 0), (1, 1)\}\) for each successive \(p_i\) and \(p_{i+1}\). In other words, this means that the warping is performed only in one direction, by "shrinking" the raw signal \(s\). The accumulated cost matrix \(D'\) is then computed using this modified relation:
%
\begin{equation}
    D'(i, j) =
    \begin{cases}
        c(1, 1)    & i=1 \land j=1 \\
        D(i, j-1) + c(i, j) & i=1 \land j \neq 1 \\
        \infty              & i \neq 1 \land j=1  \\
        \min(D(i, j-1), D(i-1, j-1)) + c(i, j).
    \end{cases} 
\end{equation}

As was shown previously, the dynamic programming approach works by constructing an accumulated cost matrix with dimensions \((n, m)\). The length of the nanopore reads can be quite long\todo{citation needed}, which could cause the matrix to be too large to fit into memory. To avoid this, the signal is split into chunks of more manageable sizes, which are processed by DTW independently.

After the accumulated cost matrix is fully calculated, the final alignment is reconstructed iteratively by retracing the optimal choices from the lower right corner of the matrix, obtaining an alignment of the reference genome to the raw data. Finally, all outputs (squiggles, reference, and mapping) for a single read are saved for further use.