\documentclass{bioinfo}
\copyrightyear{2021}
\pubyear{2021}

\usepackage{graphicx}
\usepackage{hyperref}
\usepackage{url}
\usepackage{amsmath}
\usepackage[ruled,vlined]{algorithm2e}
\newcommand\mycommfont[1]{\footnotesize\rmfamily{\it #1}}
\SetCommentSty{mycommfont}
\SetKwComment{Comment}{$\triangleright$\ }{}

\usepackage{natbib}
\bibliographystyle{apalike}

\DeclareMathOperator*{\argmax}{argmax}

\begin{document}
\firstpage{1}

\title[Improvements to minimap2]{New strategies to improve minimap2 alignment accuracy}
\author[Li]{Heng Li$^{1,2}$}
\address{$^1$Dana-Farber Cancer Institute, 450 Brookline Ave, Boston, MA 02215, USA,
$^2$Harvard Medical School, 10 Shattuck St, Boston, MA 02215, USA}

\maketitle

\begin{abstract}

\section{Summary:} We present several recent improvements to minimap2, a
versatile pairwise aligner for nucleotide sequences. Now minimap2 v2.22 can
more accurately map long reads to highly repetitive regions and align through
insertions or deletions up to 100kb by default, addressing major weakness in
minimap2 v2.18 or earlier.

\section{Availability and implementation:}
\href{https://github.com/lh3/minimap2}{https://github.com/lh3/minimap2}

\section{Contact:} hli@ds.dfci.harvard.edu
\end{abstract}

\section{Introduction}
Minimap2~\citep{Li:2018ab} is a widely used aligner for maping long sequence
reads and assembly contigs. \citet{Jain:2020aa} found minimap2 occasionally
misaligned reads from highly repetitive regions as minimap2 ignored seeds of
high occurrence. They also noticed minimap2 may misplace reads with structural
variations (SVs) in such regions~\citep{Jain2020.11.01.363887}. These
misalignments have become a pressing issue in the advent of
temolere-to-telomore human assembly~\citep{Miga:2020aa}. Meanwhile, minimap2
was unable to efficiently align long insertions/deletions (INDELs) and often
breaks an alignment around variable-number tandem repeats (VNTRs). This has
inspired new chaining algorithms~\citep{Li:2020aa,Ren:2021aa} which are not
integrated into minimap2. Here we will describe recent improvements implemented
in v2.19 through v2.22.

\begin{methods}
\section{Methods}

\subsection{Rescuing high-occurrence $k$-mers}
Minimap2 keeps all $k$-mer minimizers during indexing. Its original
implementation only selected low-occurrence minimizers during mapping. The
cutoff is a few hundred for mapping long reads against a human genome. If a
read habors only a few or even no low-occurrence minimizers, it will fail
chaining due to insufficient anchors.

To resolve this issue, we implemented a new heuristic to add additional
minimizers. Suppose we are looking at two adjacent low-occurence $k$-mers
located at position $x_1$ and $x_2$, respectively. If $|x_1-x_2|\ge500$, the
new minimap2 adds $\lfloor|x_1-x_2|/500\rfloor$ minimizers among
high-occurrence minimizers between $x_1$ and $x_2$. We use a binary heap data
structure to select minimizers of the lowest occurrence in this interval.
This strategy adds necessary anchors at the cost of increasing total alignment
time by a few percent on real data.

\subsection{Aligning through longer INDELs}
The original minimap2 may fail to align long INDELs due to its chaining
heuristics. Briefly, minimap2 applies dynamic programming (DP) to chain
minimizer anchors. It is a quadratic algorithm, which is slow for chaining
contigs. For acceptable performance, the original minimap2 uses a 500bp band by
default. If there is an INDEL longer than 500bp and the two chains around it
have no overlaps on either the query or the reference sequence, minimap2 may
join the two short chains later as a postprocessing step. We call it the
long-join heuristic. This heuristic may fail around VNTRs because short chains
often have overlaps in VNTRs. More subtly, minimap2 may escape the inner DP
loop early, again for performance, if the chaining result is not improved for
50 iterations. When there is a copy number change in a long segmental
duplication, the early escape may break around the event even if users
specify a large band.

In minigraph~\citep{Li:2020aa}, we developed a new chaining algorithm that
finds short INDELs with DP-based chaining and goes through long INDELs with a
subquadratic algorithm~\citep{DBLP:conf/wabi/AbouelhodaO03}. We ported the same
algorithm to minimap2 for contig mapping. For read mapping, the minigraph
algorithm is slower. The updated minimap2 still uses the DP-based algorithm to
find short chains and then uses the minigraph algorithm to rechain anchors in
these short chains. The rechaining steps achieves the same goal as long-join
but is more reliable as it can resolve overlaps between short chains. The old
long-join heuristic has since been removed.

\subsection{Properly mapping long reads with SVs}
The original minimap2 ranks an alignment by its Smith-Waterman score and
outputs the best scoring alignment. However, when there are SVs on the read,
the best scoring alignment is sometimes not the correct alignment.
\citet{Jain2020.11.01.363887} resolved this dilemma by altering the mapping
algorithm.

In our view, this problem is rooted in impropriate scoring: affine-gap penalty
over-penalizes a long INDEL that was often evolutionarily created in one event.
We should not penalize a SV linearly in its length. The new minimap2 rescores
an alignment with the following scoring function. Suppose an alignment consists
of $M$ matching bases, $N$ substitutions and $G$ gap opens, we empirically
score the alignment with
$$
M-\frac{N+G}{2d}-\sum_{i=1}^G\log_2(1+g_i)
$$
where $g_i\ge1$ is the length of the $i$-th gap and
$$
d=\max\left\{\frac{N+G}{M+N+G},0.02\right\}
$$
Here $d$ approximates per-base sequence divergence with the smallest value set
to 2\%. As an analogy to affine-gap scoring, the matching score in our scheme
is 1, the mismatch and gap open penalties are both $1/2d$ and the gap extension
penalty is a logarithm function of the gap length. Our scoring gives a long SV
a much milder penalty. In terms of time complexity, scoring an alignment is
linear in the length of the alignment. Time spent on rescoring is negligible in
practice.

%If we assume sequences evolve under a duplication-mutation model, we may have a
%better way to choose the best alignment. If a long read can be mapped to $n$
%loci, we can take the read as the template and build a
%pseudo-multi-sequence-alignment (pMSA) of $n+1$ sequences. In this pMSA, we say
%a site on the read is informative if the $n$ reference subsequences differ at
%the position.

\end{methods}

\section{Results}

\begin{table}
\processtable{Evaluation of minimap2 v2.22}
{\footnotesize\begin{tabular}{p{4.2cm}rrrr}
\toprule
$[$Benchmark$]$ Metric & v2.22 & v2.18 & Winno & lra \\
\midrule
$[$sim-map$]$ \% mapped reads at Q10      & 97.9 & 97.6 & {\bf 99.0} & 97.3 \\
$[$sim-map$]$ err. rate at Q10 (phredQ)   & {\bf 52}   & {\bf 52}   & 38   & 24 \\
$[$winno-cmp$]$  rate of diff. (phredQ)   & {\bf 41}   & 37   & N/A  & 18 \\
$[$sim-sv$]$  \% false negative rate      & {\bf 0.5}  & 2.0  & {\bf 0.5}  & 1.4  \\
$[$sim-sv$]$  \% false discovery rate     & {\bf 0.0}  & 0.1  & {\bf 0.0}  & 0.1  \\
$[$real-sv-1k$]$ \% false negative rate   & {\bf 7.3}  & 20.0 & 13.0 & N/A \\
$[$real-sv-1k$]$ \% false discovery rate  & 2.7  & {\bf 2.4}  & 2.7  & N/A \\
\botrule
\end{tabular}}
{In $[$sim-map$]$, 152,713 reads were simulated from the CHM13 telomere-to-telomere assembly v1.1
(AC: GCA\_009914755.3) with pbsim2~\citep{Ono:2021aa}: ``pbsim2 -{}-hmm\_model R94.model -{}-length-min
5000 -{}-length-mean 20000 -{}-accuracy-mean 0.95''. Alignments of mapping quality
10 or higher were evaluated by ``paftools.js mapeval''. The mapping error rate
is measured in the phred scale: if the error rate is $e$, $-10\log_{10}e$ is
reported in the table. In $[$winno-cmp$]$, 1.39 million CHM13 HiFi reads from
SRR11292121 were mapped against CHM13. 99.3\% of them were mapped by Winnowmap2
at mapping quality 10 or higher and were taken as ground truth to evaluate
minimap2 and lra with ``paftools.js pafcmp''. $[$sim-sv$]$ simulated 1,000
50bp to 1000bp INDELs from chr8 in CHM13 using SURVIVOR~\citep{Jeffares:2017aa} and simulated Nanopore
reads at 30 folds with the same pbsim2 command line. SVs were called with
``sniffles -q 10''~\citep{Sedlazeck:2018ab} and compared to the simulated truth with ``SURVIVOR eval
call.vcf truth.bed 50''. In $[$real-sv-1k$]$, small and long variants were
called by dipcall-0.3 for HG002 assemblies (AC: GCA\_018852605.1 and
GCA\_018852615.1) and compared to the GIAB truth~\citep{Zook:2020aa} using ``truvari -r 2000 -s
1000 -S 400 -{}-multimatch -{}-passonly'' which sets the minimum INDEL size to 1kb in evaluation. }
\end{table}

\section*{Acknowledgements}

\paragraph{Funding\textcolon} NHGRI R01HG010040

\bibliography{minimap2}

\end{document}