# Characterization of Nb films for superconducting qubits using phase boundary
measurements

Kevin M. Ryan  Northwestern University, Department of Physics and Astronomy,
Evanston, Illinois 60208, USA  Carlos G. Torres-Castanedo  Northwestern
University, Department of Materials Science and Engineering, Evanston,
Illinois 60208, USA  Dominic P. Goronzy  Northwestern University, Department
of Materials Science and Engineering, Evanston, Illinois 60208, USA
Northwestern University, International Institute for Nanotechnology, Evanston,
Illinois 60208, USA  David A. Garcia Wetter  Northwestern University,
Department of Materials Science and Engineering, Evanston, Illinois 60208, USA
Matthew J Reagor  Rigetti Computing, Berkeley, California 94710, USA  Mark
Field  Rigetti Computing, Berkeley, California 94710, USA  Cameron J Kopas
Rigetti Computing, Berkeley, California 94710, USA  Jayss Marshall  Rigetti
Computing, Berkeley, California 94710, USA  Michael J. Bedzyk  Northwestern
University, Department of Materials Science and Engineering, Evanston,
Illinois 60208, USA  Northwestern University, Department of Physics and
Astronomy, Evanston, Illinois 60208, USA  Mark C. Hersam  Northwestern
University, Department of Materials Science and Engineering, Evanston,
Illinois 60208, USA  Northwestern University, Department of Chemistry,
Evanston, Illinois 60208, USA  Northwestern University, Department of
Electrical and Computer Engineering, Evanston, Illinois 60208, USA  Venkat
Chandrasekhar  Northwestern University, Department of Physics and Astronomy,
Evanston, Illinois 60208, USA<EMAIL_ADDRESS>
###### Abstract

Continued advances in superconducting qubit performance require more detailed
understandings of the many sources of decoherence. Within these devices, two-
level systems arise due to defects, interfaces, and grain boundaries, and are
thought to be a major source of qubit decoherence at millikelvin temperatures.
In addition to Al, Nb is a commonly used metalization layer for
superconducting qubits. Consequently, a significant effort is required to
develop and qualify processes that mitigate defects in Nb films. As the
fabrication of complete superconducting qubits and their characterization at
millikelvin temperatures is a time and resource intensive process, it is
desirable to have measurement tools that can rapidly characterize the
properties of films and evaluate different treatments. Here we show that
measurements of the variation of the superconducting critical temperature
$T_{c}$ with an applied external magnetic field $H$ (of the phase boundary
$T_{c}-H$) performed with very high resolution show features that are directly
correlated with the structure of the Nb films. In combination with x-ray
diffraction measurements, we show that one can even distinguish variations
quality and crystal orientation of the grains in a Nb film by small but
reproducible changes in the measured superconducting phase boundary.

††preprint: AIP/123-QED

Two-level systems (TLSs) in superconducting devices form significant sources
of qubit decoherence at millikelvin temperatures,Premkumar  _et al._ (2021)
and are an important factor limiting qubit performance. TLSs may exist at the
interfaces between Nb films and the substrate,Murthy  _et al._ (2022a) in the
oxide layers at the surface of Nb films,Romanenko  _et al._ (2020); Verjauw
_et al._ (2021) at grain boundaries in Nb films,Premkumar  _et al._ (2021) and
at defects in the Nb film itself.Zarea, Ueki, and Sauls (2022) As common
fabrication processes preclude the use of epitaxial single crystal Nb, a
number of techniques are being explored to mitigate the effects of interfaces,
grain boundaries, and native surface oxides as sources of TLSs using wafer
scale processing.

The ultimate test of these mitigation strategies would be to fabricate
superconducting qubits with processed Nb films and measure their coherence
times at millikelvin temperatures. However, this is a time consuming and
expensive process which severely limits the throughput of materials
development. Alternatively, materials science characterization tools such as
scanning transmission electron microscopy (STEM),Altoé  _et al._ (2020);
Murthy  _et al._ (2022b) secondary ion mass spectroscopy (SIMS)Murthy  _et
al._ (2022a) or atom probe tomography (APT) can also be used to quantify the
structural quality of the Nb films. While such probes give detailed
information about the disorder of films, oxides/interfaces, and concentration
of impurities, the information is usually obtained for microscopic cross-
sections of the device and the relationships between these observations and
the superconducting properties important for qubit operation are not
immediately evident.

Conversely, macroscopic electrical transport measurements are well suited for
evaluation of process changes to superconducting films. For example, simple
temperature dependent resistivity measurements can determine $T_{c}$, and the
low temperature sheet resistance and residual resistivity ratio (RRR) provide
information about the amount of disorder in the film. More sophisticated
measurements can provide important microscopic parameters of the
superconductivity itself.

Here, we focus on one such technique which measures the superconducting phase
boundary $T_{c}-H$ rapidly and with high precision. From Ginzburg-Landau (GL)
theory, $T_{c}$ of a superconducting film with an external field $H$ applied
perpendicular to the plane of the film varies with $H$ as

$T_{c}=T_{c0}\left(1-\frac{2\pi\xi_{0}^{2}}{\Phi_{0}}H\right)$ (1)

where $\xi_{0}$ is the zero-temperature GL coherence length, $T_{c0}$ is the
zero-field transition temperature, and $\Phi_{0}=h/2e$ is the superconducting
flux quantum.Tinkham (2004) Thus, $T_{c}$ is a linear function of $H$, and its
slope combined with the measured value of $T_{c0}$ directly gives a measure of
$\xi_{0}$.

When applied to the study of arbitrary superconducting films, this analysis is
most frequently performed by tracing the resistive transition of the film at
various fixed fields, defining $T_{c}$ as $T$ at a fixed fraction of the
normal state resistance, and plotting this $T_{c}$ as a function of $H$. This
approach suffers from several drawbacks. First, in cryogenic probes used for
such measurements, there is always a thermal lag between the temperature of
the thermometer and that of the sample. Thus it is important to sweep slowly
both up and down in temperature through the resistive transition so that the
hysteresis in the traces is minimized. This is time consuming and incurs
additional cryogen and operational costs. Second, the discrete nature of the
magnetic fields at which $T_{c}$ is determined may result in features of the
phase boundary being missed. To avoid these issues, we have used a feedback
technique that allows for the measurement of $T_{c}$ as a continuous function
of $H$. Our results show non-linear deviations from GL theory in the phase
boundary that, via comparison with XRD measurements, appear correlated with
the crystal grain structure and disorder in the films.

##  I Sample preparation and measurement

The Nb films examined in this work were deposited via High Power Impulse
Magnetron Sputtering (HiPIMS) at a nominal base pressure $\leq 10^{-9}$ Torr.
To study the role of the film/substrate interface 40 nm or 155 nm films were
deposited on three different substrates: 1) intrinsically doped (001) oriented
Si; 2) a-plane Al2O3 with (110) termination; and 3) c-plane Al2O3 with (006)
termination. These substrates were chosen to qualify and improve upon current
fabrication techniques, as (001) oriented Si is the substrate used for
commercial production of many superconducting qubits,Oliver and Welander
(2013) while Al2O3 can be used to control the epitaxial growth and orientation
of Nb filmsWildes, Mayer, and Theis-Bröhl (2001) X-ray studies show that Nb
films grow as (110) oriented grains on all substrates used in this work, and
show partial epitaxial growth on Al2O3 substrates (see Supp. S1).

Substrate | d [nm] |  $T_{c0}$ [K] |  $\Delta T_{c}$ [mK] | RRR
---|---|---|---|---
Si (001) | 155 | 9.054 | 5 | 4.91
Si (001) | 40 | 8.720 | 8 | 3.96
a-Al2O3 | 40 | 8.789 | 11 | 5.23
c-Al2O3 | 40 | 8.672 | 12 | 3.29
UHV a-Al2O3 | 40 | 7.010† | 39 | 2.59
Table 1: Zero field transition properties of HiPIMS deposited films presented
in this work. Values are reported under feedback. RRR is derived from
resistance at 300K and 10K. †Values reported along the $[1\overline{1}0]$
direction, see Fig. 6.

In addition, separate 40 nm thick films deposited on a-Al2O3 substrates were
vacuum annealed at 1000∘ C for 30 minutes at 10-10 Torr. This treatment
resulted in enhanced epitaxy of the Nb, with the $[1\bar{1}0]$ Nb plane
oriented along the$[\overline{1}00]$ direction on the $(110)$ Al2O3 surface.
X-ray diffraction (XRD) rocking curve data show that this treatment results in
the average grain size of the Nb increasing from 7 nm to 124 nm, with columnar
growth extending from the substrates to the surface oxide (see Supp. S2).
Figure 1 shows Atomic Force Microscopy (AFM) maps of Nb films on a-Al2O3
before and after this annealing process. After deposition the films were
patterned using standard photolithography techniques into
$1.2\text{\,}\mathrm{mm}$ by $30\text{\,}\mathrm{\SIUnitSymbolMicro m}$ Hall
bars, then defined using an Ar ion mill with the photoresist serving as an
etch mask.

Figure 1: AFM micrographs of pristine (left) and UHV annealed (right) 40 nm Nb
films on a-Al2O3 substrates. Both scale bars are 250 nm.

The devices were measured in a custom probe inserted into a He dewar with a 3
T maximum field perpendicular to the plane of the films. In the probe, the
devices were mounted in vacuum and thermally anchored to a variable
temperature stage that is weakly coupled to the 4K flange of the probe. By
controlling the power through a heater mounted on this stage, the temperature
of the stage could be varied with a short response time. The accuracy of the
thermometer was verified by measuring the transition of a 99.99% pure Nb wire,
which showed a $T_{c}$ of 9.22 K. The resistance of the devices were measured
by custom modified Adler-Jackson ac resistance bridges,Adler and Jackson
(1966) with ac excitations of $\sim 100$ nA at frequencies
$\leq$200\text{\,}\mathrm{Hz}$$. Table 1 summarizes details of the transitions
of the films included in this study. The transition temperatures of the
majority of the films we measured were in the range of 8.5-9 K with residual
resistance ratios (RRR, the ratio of the room temperature resistance to
$R_{n}$, the resistance just above the superconducting transition) of $\sim$5.

Figure 2: $R(T)$ for a 40 nm Nb film on a a-Al2O3 substrate, normalized to the
normal state resistance $R_{n}$ measured at various externally applied
perpendicular fields obtained at a temperature sweep rate of $<5$K/hr,
demonstrating $\sim$20 mK of thermal hysteresis. Red/blue curves indicate data
taken during warming/cooling.

To initially qualify the superconducting transition of the films, the power to
the stage heater was ramped slowly under computer control. Fig. 2 shows
examples of a superconducting transitions measured on a 40 nm thick film on an
a-Al2O3 substrate ramping both up and down in temperature at various magnetic
fields. For most films studied, the transition widths were $\lessapprox$ 20
mK. Figure 2 also illustrates the problem of using the resistive transition to
accurately measure the phase boundary: even with a slow ramp rate of 5 K/hr,
the traces corresponding to sweeps up and down in temperature show hysteresis
comparable to the transition width.

To accurately measure $T_{c}-H$, we use a feedback technique to maintain the
resistance of the sample at a specified value along the superconducting
transition by adjusting the set-point resistance of the bridge to the desired
value as an applied magnetic field is swept. The error signal from the bridge
is then used as the input to a homemade analog proportional-integral-
differential (PID) controller which controls the power to the stage heater to
maintain the sample resistance. By fixing the bridge resistance at $R_{n}/2$,
the temperature as read from the thermometer directly gives $T_{c}$ when the
error signal is minimized by the PID under feedback. To define the transition
width, one can additionally measure $T-H$ at other values of resistance, such
as $5\%$ $R_{n}$ and $95\%$ $R_{n}$ to measure the transition width (see Fig.
5). This technique has been used to measure coherence lengths and anomalous
Little-Parks oscillations in mesoscopic Al structures;Vloeberghs  _et al._
(1992) superconductivity in star shaped mesoscopic devices;Dikin  _et al._
(2003) and hysteretic superconducting phase boundaries arising from underlying
magnetism in two-dimensional LAO/STO heterostructures;Mehta (2015) but has not
previously been applied to the assessment of superconducting materials for
quantum applications.

Figure 3: $T_{c}-H$ phase boundary measured for a 40 nm Nb film on Si (001),
taken over repeated field sweeps. The dashed lines indicate linear trends
below 1 kOe (blue) and above 3 kOe (red).

##  II Experimental results

Figure 3 shows an example $T_{c}-H$ measurement for a 40 nm thick Nb film on
Si(001), ramping in both directions in field. As can be seen, the curves
corresponding to the different sweep directions cannot be distinguished, with
the reproducibility being better than 1 mK except at the field extrema where
the sweep direction changes (see Supp. S3). As noted, GL theory predicts
$T_{c}$ to be a linear function of $H$, with a slope determined by $\xi_{0}$.
The traces in Fig. 3 show a clear nonlinear structure, with the derivative
$|dT_{c}/dH|$ largest near zero field and decreasing asymptotically to a
constant value at larger fields. This positive curvature of the
superconducting phase boundary is seen to an extent in all as-deposited films
measured. Before we explore the origin of this curvature, we first discuss the
overall differences between Nb films grown on different substrates which can
be seen by this method.

Figure 4: $T_{c}-H$ phase boundaries for Nb films of 40 nm and 155 nm
thickness on Si (001). Both measurements demonstrate varying degrees of
positive curvature.

Fig. 4 shows $T_{c}-H$ for two Nb films grown on Si (001) substrates. The
$T_{c0}$ of the 40 nm film is lower than that of the corresponding 155 nm
film, in line with existing results that show a sharp drop in $T_{c0}$’s and
RRRs of Nb films below a thickness of 200 nm.Pinto  _et al._ (2018) The
corresponding slopes of the $T_{c}-H$ curves for these films give equivalent
coherence lengths $\xi_{0}\leq$13 nm (See Supp. S4), far less than the
$\geq$30 nm seen in bulk NbWilliamson (1970) but comparable to reports for
thin deposited films.Pinto  _et al._ (2018); Bose  _et al._ (2006) Notably,
the non-linearity is more evident in the thinner 40 nm film, appearing as a
decrease in slope between 2-3 kOe. This transition implies substantial
inhomogeneity throughout the film, possibly due to greater influence of the
Nb/substrate interface.

Figure 5: Temperature measured at fixed resistance values of $5\%$, $50\%$,
and $95\%$ of the normal state resistance for 40 nm Nb on a-Al2O3 (top) and
c-Al2O3 (bottom) substrates.

Fig. 5 shows equivalent data for 40 nm thick Nb films grown on a- and c-plane
Al2O3 substrates, where we now additionally plot the $T-H$ curves taken under
feedback at 5% and 95% of the normal state resistance, to indicate the full
width of the transition. These $T_{c0}$ values are comparable to 40 nm Nb
films grown on Si(001) (Fig. 4), but a clear $T_{c}$ and $\xi_{0}$ enhancement
can be seen in the a-Al2O3 film. We attribute this to improved epitaxy between
the Nb(110) termination and the termination a-Al2O3 of the substrate
(discussed in Supp. S1). The shape of the $T_{c}-H$ curves for the two
substrate orientations is characteristically different, and the curves for the
5% $R_{n}$ and 95% $R_{n}$ diverge more at larger fields for the a-Al2O3 film
than for the c-Al2O3 film. Thus there appears to be a correspondence between
greater transition width broadening, positive curvature of the phase boundary,
higher $T_{c0}$, and crystalline epitaxy.

Figure 6: Raw (top) and normalized (bottom) $T_{c}-H$ phase boundaries for UHV
annealed 40 nm thick Nb films on a-Al2O3 substrates.

These correlated features allow us to consider origins for the unexpected
nonlinearities in $T_{c}-H$. As the effect appears most strongly in thinner
films, and varies in magnitude for films of similar thickness on different
substrates, it is likely that film granularity and inhomogeneity at the Nb-
substrate interface contribute substantially. Evidence for this can be seen in
our measurements of the more epitaxial UHV annealed Nb films, when studied
along different crystal orientations with respect to the electronic current
direction.

As noted before, Nb films grow with a (110) termination on a-Al2O3, with the
$[1\bar{1}0]$ direction of the Nb films oriented along the $[100]$ surface
direction on the a-Al2O3. XRD measurements (Supp. S1) show that UHV annealing
yields columnar grains reaching from substrate to surface with lateral sizes
as large at $\sim$120 nm, much larger the $\xi_{0}$. To explore the dependence
of the superconducting properties as a function of the angle between the in-
plane crystal direction and the current direction, Hall bars oriented with
their lengths at specific angles to the Nb $[1\bar{1}0]$ orientation were
fabricated from the film. $T_{c}-H$ curves for these samples are shown in Fig.
6, and are significantly different from the samples discussed earlier. First,
$T_{c0}$ for these films is considerably reduced in comparison to unannealed
Nb films on all substrates. However this is consistent with the reduced RRR,
indicative of greater disorder in spite of the improved crystallinity,
possibly due to voids generated in the film after relaxation of the crystal or
diffusion of surface oxides into the bulk. Second, $T_{c0}$ appears to vary
systematically with the direction of the current with respect to the
crystalline axes, being a maximum when the current is aligned along the Nb
$[1\bar{1}0]$ plane, and minimum close to the Nb $[001]$ direction (Fig. 6).
Third, and most striking in view of the data on Nb films on other substrates,
$T_{c}$ is almost a perfectly linear function of $H$ except for a small field
regime near $H=0$ where one now sees a _negative_ curvature. Normalizing these
phase boundaries by $T_{c0}$, such that their slope is only defined by
$\xi^{2}$, collapses them to a single curve (Fig. 6(b)), implying that in
spite of the variation of $T_{c0}$ with angle, $\xi_{0}$ is the same for all
directions. Given the lack of correlation between $T_{c0}$ and $\xi$ after
annealing, we believe that this can occur when the grain size is much larger
than $\xi_{0}$.

These results enable us to propose a simple model to explain the
nonlinearities in $T_{c}-H$ observed in polycrystalline Nb films with small,
randomly oriented grains, wherein the Nb grains have a distribution of
$T_{c0}$ and $\xi_{0}$ with higher $T_{c0}$ corresponding to longer $\xi_{0}$.
This is possible when the GL parameter $\kappa$, the ratio of the London
penetration depth $\lambda$ and $\xi_{0}$ increases as a function of impurity
concentration, as is widely reported in Nb.Narlikar and Dew-Hughes (1966);
Kozhevnikov  _et al._ (2018); Finnemore, Stromberg, and Swenson (1966) Near
$H=0$, the grains with the largest $T_{c0}$ contribute to $T_{c}-H$. As $H$
increases, the $T_{c}$ of these grains is more rapidly reduced due to their
longer $\xi_{0}$, at which point grains with smaller $T_{c0}$ and shorter
$\xi_{0}$ dominate the transport, leading to a reduction in the slope of the
$T_{c}-H$ curves. For uniform epitaxial films with a single predominant
crystal axis, one observes a linear behavior in $T_{c}-H$ as in Fig. 6. For
films with grains with two distinct crystalline axes with respect the current,
one would observe two distinct slopes, as seen in Fig. 5. For films with more
complicated crystal grain structures, one might observe multiple slopes or
even a continual change in slope at low fields.

In summary, we observe unexpected nonlinearities in the superconducting phase
boundaries of thin polycrystalline Nb films fabricated on different
substrates. These nonlinearities are due the distributions of critical
temperatures and coherence lengths of crystal grains in the films. Careful
measurements of the superconducting phase boundary as demonstrated here, thus
enable a rapid means of determining the quality of Nb films used in the
fabrication of superconducting qubits.

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, National Quantum Information Science Research Centers,
Superconducting Quantum Materials and Systems Center (SQMS) under Contract No.
DEAC02-07CH11359. This work made use of the NUFAB facility of Northwestern
University’s NUANCE Center, which has received support from the SHyNE Resource
(NSF ECCS-2025633), the IIN, and Northwestern’s MRSEC program (NSF
DMR-1720139). This work made use of the Jerome B. Cohen X-Ray Diffraction
Facility supported by the MRSEC program of the National Science Foundation
(DMR-1720139) at the Materials Research Center of Northwestern University and
the Soft and Hybrid Nanotechnology Experimental (SHyNE) Resource (NSF
ECCS-2025633). Additional support and equipment was provided by DURIP grant
W911NF-20-1-0066.

## References

  * Premkumar  _et al._ (2021) A. Premkumar, C. Weiland, S. Hwang, B. Jäck, A. P. M. Place, I. Waluyo, A. Hunt, V. Bisogni, J. Pelliciari, A. Barbour, M. S. Miller, P. Russo, F. Camino, K. Kisslinger, X. Tong, M. S. Hybertsen, A. A. Houck, and I. Jarrige, “Microscopic relaxation channels in materials for superconducting qubits,” Communications Materials 2, 72 (2021).
  * Murthy  _et al._ (2022a) A. A. Murthy, J. Lee, C. Kopas, M. J. Reagor, A. P. McFadden, D. P. Pappas, M. Checchin, A. Grassellino, and A. Romanenko, “TOF-SIMS analysis of decoherence sources in superconducting qubits,” Applied Physics Letters 120, 044002 (2022a).
  * Romanenko  _et al._ (2020) A. Romanenko, R. Pilipenko, S. Zorzetti, D. Frolov, M. Awida, S. Belomestnykh, S. Posen, and A. Grassellino, “Three-Dimensional Superconducting Resonators at T < 20 mK with Photon Lifetimes up to $\tau$ = 2 s,” Physical Review Applied 13, 034032 (2020).
  * Verjauw  _et al._ (2021) J. Verjauw, A. Potočnik, M. Mongillo, R. Acharya, F. Mohiyaddin, G. Simion, A. Pacco, T. Ivanov, D. Wan, A. Vanleenhove, L. Souriau, J. Jussot, A. Thiam, J. Swerts, X. Piao, S. Couet, M. Heyns, B. Govoreanu, and I. Radu, “Investigation of Microwave Loss Induced by Oxide Regrowth in High- Q Niobium Resonators,” Physical Review Applied 16, 014018 (2021).
  * Zarea, Ueki, and Sauls (2022) M. Zarea, H. Ueki, and J. A. Sauls, “Effects of anisotropy and disorder on the superconducting properties of Niobium,” arXiv:2201.07403 [cond-mat] (2022), arXiv: 2201.07403.
  * Altoé  _et al._ (2020) M. V. P. Altoé, A. Banerjee, C. Berk, A. Hajr, A. Schwartzberg, C. Song, M. A. Ghadeer, S. Aloni, M. J. Elowson, J. M. Kreikebaum, E. K. Wong, S. Griffin, S. Rao, A. Weber-Bargioni, A. M. Minor, D. I. Santiago, S. Cabrini, I. Siddiqi, and D. F. Ogletree, “Localization and reduction of superconducting quantum coherent circuit losses,” arXiv:2012.07604 [cond-mat, physics:quant-ph] (2020), arXiv: 2012.07604.
  * Murthy  _et al._ (2022b) A. A. Murthy, P. M. Das, S. M. Ribet, C. Kopas, J. Lee, M. J. Reagor, L. Zhou, M. J. Kramer, M. C. Hersam, M. Checchin, A. Grassellino, R. d. Reis, V. P. Dravid, and A. Romanenko, “Insights into Improving Performance of Niobium Superconducting Transmon Qubit Devices through Advanced Microscopy,” Tech. Rep. arXiv:2203.08710 (arXiv, 2022) arXiv:2203.08710 [cond-mat, physics:quant-ph] type: article.
  * Tinkham (2004) M. Tinkham,  _Introduction to Superconductivity_, 2nd ed. (Dover Publications, 2004) published: Paperback.
  * Oliver and Welander (2013) W. D. Oliver and P. B. Welander, “Materials in superconducting quantum bits,” MRS Bulletin 38, 816–825 (2013).
  * Wildes, Mayer, and Theis-Bröhl (2001) A. Wildes, J. Mayer, and K. Theis-Bröhl, “The growth and structure of epitaxial niobium on sapphire,” Thin Solid Films 401, 7–34 (2001).
  * Adler and Jackson (1966) J. G. Adler and J. E. Jackson, “System for Observing Small Nonlinearities in Tunnel Junctions,” Review of Scientific Instruments 37, 1049–1054 (1966).
  * Vloeberghs  _et al._ (1992) H. Vloeberghs, V. V. Moshchalkov, C. Van Haesendonck, R. Jonckheere, and Y. Bruynseraede, “Anomalous Little-Parks oscillations in mesoscopic loops,” Phys. Rev. Lett. 69, 1268–1271 (1992), publisher: American Physical Society.
  * Dikin  _et al._ (2003) D. A. Dikin, V. Chandrasekhar, V. R. Misko, V. M. Fomin, and J. T. Devreese, “Nucleation of superconductivity in mesoscopic star-shaped superconductors,” The European Physical Journal B - Condensed Matter 34, 231–235 (2003).
  * Mehta (2015) M. Mehta,  _Interplay between superconductivity and ferromagnetism at the LaAlO3/SrTiO3 interface_ , PhD Thesis (2015).
  * Pinto  _et al._ (2018) N. Pinto, S. J. Rezvani, A. Perali, L. Flammia, M. V. Milošević, M. Fretto, C. Cassiago, and N. De Leo, “Dimensional crossover and incipient quantum size effects in superconducting niobium nanofilms,” Scientific Reports 8, 4710 (2018).
  * Williamson (1970) S. J. Williamson, “Bulk Upper Critical Field of Clean Type-II Superconductors: V and Nb,” Physical Review B 2, 3545–3556 (1970).
  * Bose  _et al._ (2006) S. Bose, P. Raychaudhuri, R. Banerjee, and P. Ayyub, “Upper critical field in nanostructured Nb: Competing effects of the reduction in density of states and the mean free path,” Physical Review B 74, 224502 (2006).
  * Narlikar and Dew-Hughes (1966) A. V. Narlikar and D. Dew-Hughes, “Superconductivity in deformed niobium alloys,” Journal of Materials Science 1, 317–335 (1966).
  * Kozhevnikov  _et al._ (2018) V. Kozhevnikov, A.-M. Valente-Feliciano, P. J. Curran, G. Richter, A. Volodin, A. Suter, S. J. Bending, and C. Van Haesendonck, “Equilibrium Properties of the Mixed State in Superconducting Niobium in a Transverse Magnetic Field: Experiment and Theoretical Model,” Journal of Superconductivity and Novel Magnetism 31, 3433–3444 (2018).
  * Finnemore, Stromberg, and Swenson (1966) D. K. Finnemore, T. F. Stromberg, and C. A. Swenson, “Superconducting Properties of High-Purity Niobium,” Physical Review 149, 231–243 (1966).

## Supplementary material for: Characterization of Nb Films for
Superconducting Qubits using Phase Boundary Measurements

##  S1 Epitaxy of Al2O3 deposited Nb films

Figure S1: In-plane XRD of the Nb($110$)/Al2O3($110$) film before (left) and
after the UHV annealing (right) at $1000^{\circ}$C for 30 minutes.

Epitaxial growth of Nb on Al2O3 was confirmed via X-ray characterization using
a Rigaku Smartlab 9kW diffractometer. Specular X-ray diffraction and X-ray
reflectivity (XRR) measurements were performed at
$\lambda=$1.5406\text{\,}\mathrm{\SIUnitSymbolAngstrom}$$ (Cu K$\alpha_{1}$)
with a Ge (220)x2 monochromator. In the case of in-plane X-ray diffraction
measurements, an incident and receiving $0.5^{\circ}$ Soller slits were
utilized with no monochromator (Cu K$\alpha$).

Epitaxy of the Al2O3(110) grown Nb films was studied by X-ray diffraction
before and after the UHV annealing (Fig. S1). In-plane XRD measurements show
2-fold symmetry in the Nb [$1\overline{1}0$] in-plane Bragg peak (Fig. S1) for
both the as-deposited and UHV-annealed samples. However, an additional 2-fold
symmetry domain is observed in the as-deposited film suggesting an impartial
epitaxy for the Nb[$1\overline{1}0$] || Al2O3[$1\overline{1}04]$ orientation.
The UHV annealed sample shows an overall improvement in the epitaxial
orientation of the films, with two sharp Nb [$1\overline{1}0$] peaks, 180°
apart, as expected for a uniform epitaxially aligned film.

##  S2 X-ray properties of UHV annealed Nb on Al2O3 (110)

Figure S2: Specular XRD of the Nb(110)/Al2O3(110) film before and after the
UHV annealing at 1000°C for 30 min.

The ultra-high vacuum (UHV) annealing was performed in a custom-built chamber
with a base pressure of
$\approx$2\text{\times}{10}^{-10}\text{\,}\mathrm{T}\mathrm{o}\mathrm{r}\mathrm{r}$$.
The annealing was monitored with a pyrometer. The final temperature of
$1000^{\circ}\,$C was held for 30 minutes at a base pressure of
$\approx$2\text{\times}{10}^{-9}\text{\,}\mathrm{T}\mathrm{o}\mathrm{r}\mathrm{r}$$.
UHV conditions were necessary to increase the crystallinity of the films
without oxidizing the Nb film. The morphology of the film was analyzed by
atomic force microscopy in an Asylum Cypher (Oxford Instruments) using tapping
mode, and is shown in Fig. 1 of the main text.

Nb/a-Al2O3 (110) | Pristine | UHV Annealed
---|---|---
Nb (110) | $2.649\text{\,}\mathrm{\SIUnitSymbolAngstrom}^{-1}$ | $2.687\text{\,}\mathrm{\SIUnitSymbolAngstrom}^{-1}$
$a$ | $3.354\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ | $3.307\text{\,}\mathrm{\SIUnitSymbolAngstrom}$
Lattice Expansion | 1.45% | 0.02%
$L_{z}$ | - | $327\text{\,}\mathrm{\SIUnitSymbolAngstrom}$
$L_{x}$ | - | $867\text{\,}\mathrm{\SIUnitSymbolAngstrom}$
$d$ | $320\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ | $333\text{\,}\mathrm{\SIUnitSymbolAngstrom}$
Nb Roughness | $7\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ | $5\text{\,}\mathrm{\SIUnitSymbolAngstrom}$
Table S1: XRD/XRR comparison of pristine and UHV annealed Nb films on a-Al2O3,
where $a$ is the BCC lattice parameter, $L_{z}$ is the measured vertical grain
size, $L_{x}$ is the measured horizontal grain size, and $d$ is the measured
thickness of the Nb layer.

The as-deposited film shows a wide Bragg Nb ($110$) peak suggesting low
crystallinity and a large offset with the Nb(110) bulk value as seen in Fig.
S2. After annealing, the peak shifts to the expected bulk Nb position. The
sharpening of the Nb (110) Bragg peak and the appearance of thickness fringes
is evidence of improved crystallinity of the film after annealing. Table S1
summarizes these findings. The vertical grain size of the annealed film
extracted from XRD rocking curve measurements reaches a value close to the
predicted thickness by XRR. This suggests that the grains form a columnar
structure up to the surface of the film ( 35 nm), which is ideal for
epitaxially grown Nb. While the as-received sample has a contraction in the
lattice (-1.28%), expected from the expansion in the specular directions,
while the annealed sample shows a slight expansion (0.17%), which may have
created voids across the film after relaxation.

##  S3 Phase boundary reproducibility and errors

Figure S3: Left: Measured $R$ and $T_{c}$ for a 40 nm Nb film on (001) Si over
two full traces over the field range presented in Fig. 4. Right: Raw $T-H$
data for the same 40 nm Nb film 4

With sufficient tuning of the PID controller, typically the resistance of the
sample can be fixed with mK precision at fixed field. As the midpoint of the
transition changes under applied fields, small resistance offsets on the order
of 1% of the normal state resistance can occur as the sample is heated or
cooled. This “thermal lag”, for samples with sharp transitions, yields
negligible $T_{c}$ hysteresis on the scale of the change in $T_{c}$ in field.
This can be seen by examining raw time series data presented in Fig. S3

The more prominent regions where error occours are at the field extrema and at
$H=0$, where the time rate of change $dT_{c}/dt$ changes sign, causing a rapid
change in the error signal provided to the PID controller, which requires a
finite integration time to return to the midpoint of the transition. By
monitoring the sample resistance/error signal during a phase boundary
measurement, one obtains a direct detection for periods where the sample
temperature is not at $T_{c}$ at these extrema. For the results presented in
the main text, we have suppressed these errors in two ways. First, by sweeping
past the field regime of interest (typically by 500 Oe or more), the field
reversal error can be safely ignored as the PID controller brings the sample
back to $T_{c}$ during this overshoot. Second, for errors near $H=0$, data-
points which have large error signals can be excluded by a cutoff criteria.
Again, for the results presented in the main text, data points where the
resistance/error signal exceeds $3\sigma$ have been filtered out before
interpolation/plotting. By definition, these points cannot be measurements of
$T_{c}$ under the $50\%$ $R_{n}$ criteria. Raw data, without filtering or
supression of the field range for a 40 nm film is also shown in Fig. S3.

##  S4 Effective Coherence Length

Figure S4: Effective coherence length for Nb films on (001) Si substrates.

For the examination of real phase boundary data, one may chose to redefine
$\xi_{0}$ in terms of the derivative of $T_{c0}(H)$ from Eq. 1 as

$\xi_{0}\equiv\sqrt{\left|\frac{dT_{c}}{dH}\frac{\Phi_{0}}{2\pi
T_{c0}}\right|}$ (2)

In the case of small non-linearity, this treatment can reveal subtle curvature
of the dataset, which can then be treated in terms of variation of the
effective GL coherence length. Considering the data presented in Fig. 4, when
plotted in this form, a sharp inflection in the curvature of the phase
boundary can be seen to be associated with a $\approx 1$ nm decrease the in
effective coherence length for the 40 nm Nb film on (001) Si.

Numerical differentiation was performed with a second order Savitzky-Golay
filter, with a window size of 155 Oe, which causes the dip at zero field where
${dT_{c}}/{dH}$ reverses sign, and at the extrema of the data set where the
range is restricted.


Department of Computer Science, University of Pisa,
<EMAIL_ADDRESS>of Tennessee Health Science Center,
<EMAIL_ADDRESS>Department of Computer Science, University of Pisa,
Italynadia.pisanti@unipi.ithttps://orcid.org/0000-0003-3915-7665 Njagi Moses
Mwaniki and Erik Garrison and Nadia Pisanti [100]68W01, 68W32, 68Q25, 68Q17
This paper is part of a project that has received funding from the European
Union’s Horizon 2020 research and innovation programme under the Marie
Skłodowska-Curie grant agreement No 956229. Co-financed by the Connecting
Europe Facility of the European Union. John Q. Open and Joan R. Access 2 42nd
Conference on Very Important Topics (CVIT 2016) CVIT 2016 CVIT 2016 December
24–27, 2016 Little Whinging, United Kingdom 42 XX

# Fast Exact String to D-Texts Alignments

Njagi Moses Mwaniki  Erik Garrison  Nadia Pisanti

###### Abstract

In recent years, aligning a sequence to a pangenome has become a central
problem in genomics and pangenomics. A fast and accurate solution to this
problem can serve as a toolkit to many crucial tasks such as read-correction,
Multiple Sequences Alignment (MSA), genome assemblies, variant calling, just
to name a few. In this paper we propose a new, fast and exact method to align
a string to a D-string, the latter possibly representing an MSA, a pan-genome
or a partial assembly.
An implementation of our tool dsa is publicly available at
https://github.com/urbanslug/dsa.

###### keywords:

D-texts, Pangenomes, Alignments

###### category:

##  1 Introduction

In recent years, aligning a sequence to a pangenome has become a central
problem in genomics and pangenomics ([5]). This problem has often been
addressed with the purpose of aligning sequencing reads to a complex structure
such as a degenerate string or a more general genome graph structure (e.g.[9,
22, 2, 30, 29, 19]), and hence assuming the query string to be considerably
shorter than the pangenome (e.g. [9, 22, 3, 1, 2]). As sequencing data is
recently offering longer sequences with high accuracy, the attention has very
recently been moved towards the problem of aligning to a graph-like structure
a string which is approximately as long: a fast and accurate solution to this
problem would serve as a toolkit to many crucial tasks such as read-
correction, Multiple Sequences Alignment (MSA), genome assemblies, variant
calling, just to name a few. To this purpose, because the requirement is an
alignment that is at the same time exact and sensible to indels as well as
mismatches, then new computational challenges are faced: current exact methods
that perform exact alignments to long strings would have a quadratic
complexity with high memory requirements, and hence be prohibitively slow.
Consider the following MSA of three closely-related sequences (on the left)
and its compact representation as a D-string (on the right):
GCAATCGGGTATT GCAATCGGGAATT GCACGCTGGATTT
$\hat{T}=\texttt{GCA}[{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\texttt{AT}/\texttt{CG}}]\texttt{C}\left[{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\texttt{G}/\texttt{T}}\right]\texttt{GG}[{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\texttt{TA}/\texttt{AA}/\texttt{AT}}]\texttt{TT}$

The D-string (or D-text) $\hat{T}$ contains some deterministic (shown in
black) and some non-deterministic (shown in red) segments. Formally, $\hat{T}$
is a sequence of n sets of strings where the ith set contains strings of the
same length $\ell_{i}$ (possibly $=1$ in the deterministic case) but this
length can vary between different sets.
An indeterminate string over an alphabet $\Sigma$ is a sequence of subsets of
$\Sigma$, and it basically corresponds to the special case of D-string in
which $\ell_{i}=1$ for all $1\leq i\leq n$. A great deal of research has been
conducted in the bioinformatic literature on this type of uncertain sequence
(see [20, 31, 32, 24, 28, 6] and references therein), that are equivalent to a
sequence written in the IUPAC notation [21] to represent a position in a DNA
sequence that can have multiple possible alternatives. These are commonly used
to encode the consensus of a population of sequences [37, 5, 26, 27, 33] in a
multiple sequence alignment (MSA).
The more general notion of ED-strings (where within a degenerate position
variants can have different sizes), and over them the short read matching
problem elastic-degenerate string matching (EDSM) problem has attracted some
attention in the combinatorial pattern matching community. Since its
introduction in 2017 [7], a series of results have been published both for the
exact ([38, 14, 25, 15]) as well as for the approximate ([13, 12]) version of
the problem.
In this paper we propose a new, fast and exact method to align a string to a
D-string, the latter possibly representing an MSA, a pan-genome or a partial
assembly. Ours is a base-level heuristic-free alignment method that allows
affine gap penalty score functions and the use of arbitrary scores as well as
weights. D-strings are one out of many others possible pangenome
representation ([40]). Our algorithm exploits the assumption that the
similarity the alignment must detect is high enough to suitably adapt the
seminal idea of [10] to D-strings, combining at the same time _partial order
alignments_ ([4]) and the _wavefront_ paradigm ([39]) managing to work with
affine gap penalty function to score gaps, which is a desirable property when
the target of the alignment is to account for INDELs variant events.

This paper is organised as follows: Section 2 gives some preliminary
definitions on D-strings and notions of partial order and wavefront alignment,
Section 3 describes our algorithm, and Section 4 the experimental validation
of the resulting tool dsa we realized.
The implementation of DSA is publicly available at
https://github.com/urbanslug/dsa, and our D-strings random generator can be
found at https://github.com/urbanslug/simed/.

##  2 Preliminary Notions

###  2.1 D-strings

A string $X$ is a sequence of elements on an alphabet $\Sigma$, where the
alphabet $\Sigma$ is a non-empty finite set of letters of size $|\Sigma|$. The
set of all finite strings over an alphabet $\Sigma$, including the empty
string $\varepsilon$ of length $0$, is denoted by $\Sigma^{*}$, while
$\Sigma^{+}$ denotes the set $\Sigma^{*}\setminus\\{\epsilon\\}$. The set of
all strings of length $k\\!>\\!0$ over $\Sigma$ is denoted by $\Sigma^{k}$.
For any string $X$, we denote by $X[i,j]$ the substring of $X$ that starts at
position $i$ and ends at position $j$. In particular, $X[0,j]$ is the prefix
of $X$ that ends at position $j$, and $X[i,|X|]$ is the suffix of $X$ that
starts at position $i$, where $|X|$ denotes the length of $X$.

######  Definition 2.1.

A _degenerate string (D-string)_
$\hat{S}=\hat{S}\\{1\\}\hat{S}\\{2\\}\dots\hat{S}\\{n\\}$ of length $n$ over
an alphabet $\Sigma$ is a finite sequence of $n$ degenerate letters
$\hat{S}\\{i\\}$’s. Each degenerate letter $\hat{S}\\{i\\}$ has _width_
$\ell_{i}\\!>\\!0$ and is a finite non-empty set of $|\hat{S}\\{i\\}|$ strings
of the same length $\ell_{i}$ (i.e.
$\hat{S}\\{i\\}[1],\ldots,\hat{S}\\{i\\}[|\hat{S}\\{i\\}|]\in\Sigma^{\ell_{i}}$).

For any $1\leq i\leq n$, we remark that a degenerate letter $\hat{S}\\{i\\}$
with $|\hat{S}\\{i\\}|=\ell_{i}=1$ is just a simple letter of $\Sigma$.
Whenever this happens, we will say that $i$ is a _solid position_ of
$\hat{S}$. We now define some parameters that measure the degeneracy of a
D-string.

######  Definition 2.2.

The _total size_ $N$ and _total width_ $w(\hat{T})=W$ of a D-string $\hat{S}$
are respectively defined as $N=\sum_{i=1}^{n}|\hat{S}\\{i\\}|\cdot\ell_{i}$
and $W=\sum_{i=1}^{n}\ell_{i}$.

######  Example 2.3.

In the D-string $\hat{T}=\texttt{GCA[AT/CG]C[G/T]GG[TA/AA/AT]TT}$,
$\ell_{4}=\ell_{9}=2$ while all other $\ell_{i}$’s are $1$, and
$|\hat{S}\\{4\\}|=|\hat{S}\\{6\\}|=2$, $|\hat{S}\\{9\\}|=3$ while all the
other $|\hat{S}\\{i\\}|$’s are $1$; finally, the solid positions of $\hat{T}$
are $1,2,3,5,7,8,10,11$.

We remark that a D-string of width $W$ actually represents a set of linear
strings of length $W$, corresponding to all the strings that can be read by
making any choice at a degenerate positions. Formally, given a D-string
$\hat{T}$ of width $W$, and a string $T$ with $|T|=W$ being any string $T$ in
such set, we say that $T$ belongs to $\hat{T}$ ($T\in\hat{T}$).

######  Example 2.4.

The strings GCACGCTGGAATT and GCAATCTGGTATT are two of the twelve strings that
belong to the D-string $\hat{T}=\texttt{GCA[AT/CG]C[G/T]GG[TA/AA/AT]TT}$.

The $i^{th}$ position $\hat{S}\\{i\\}$ of a D-string $\hat{S}$ ($1\leq i\leq
n$) should not be confused with its $i^{th}$ width $\hat{S}[i]$ ($1\leq i\leq
W$):

######  Definition 2.5.

For any D-string $\hat{S}$, its $i^{th}$ width (for $1\leq i\leq W$)
$\hat{S}[i]$ is the letter, or the set of alternative letters, that can appear
in $T[i]$ for string $T\in\hat{T}$. For any D-string $\hat{S}$, we denote by
$\hat{S}[i_{1},i_{2}]$ the _D-substring_ of $\hat{S}$ that starts at width
$i_{1}$ and ends at width $i_{2}$ with $1\leq i_{1}\leq i_{2}\leq W$.

We will also use the notation $s_{i}$ for both $|\hat{S}\\{i\\}|$ and
$|\hat{S}[i]|$.

######  Example 2.6.

The D-string $\hat{T}=\texttt{GCA[AT/CG]C[G/T]GG[TA/AA/AT]TT}$ has length
$n\\!=\\!11$, size $N\\!=\\!20$, and width $W\\!=\\!13$. We have that at
width, say, $3$ there is a solid position $\texttt{A}=\hat{S}[3]$, while at
width $4$ there is $\texttt{[A/C]}=\hat{S}[4]$. Furthermore, the D-substring
$\hat{S}[3,7]$ is A[AT/CG]C[G/T]. Here below we show positions and widths for
the D-string CA[AT/CG]C[G/T]GG[TA/AA]T having length 9, width 11, and size 16.

D-string | $\tt{C}$ | $\tt{A}$ | [A | T | / | C | G] | C | [G | / | T] | G | G | [T | A | / | A | A] | T
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
width | 1 | 2 | 3 | 4 |  | 3 | 4 | 5 | 6 |  | 6 | 7 | 8 | 9 | 10 |  | 9 | 10 | 11
position | 1 | 2 |  |  | 3 |  |  | 4 |  | 5 |  | 6 | 7 |  |  | 8 |  |  | 9

###  2.2 Partial Order Alignments

The use of Partial Order (graphs) for sequence Alignments (POA) was introduced
in [4] with the purpose of improving the iterative step of progressive ([8])
multiple sequence alignments (MSA), that is the step where a new sequence is
added to an MSA (_starting MSA_). Traditionally, this was performed by first
reducing the starting MSA to a linear profile, and then aligning it to the new
sequence in order to obtain a new MSA (_resulting MSA_), with the downside
that the reduction of the starting MSA to a linear profile might carry and
propagate alignment errors. The idea behind POA is to perform a direct
pairwise dynamic programming alignment between the new sequence and a suitable
graph representation of the starting MSA (thus replacing the somewhat lossy
linear profile representation), with the outcome of a guarantee that, in the
resulting MSA, the optimal alignment of the new sequence with respect to each
one of the sequences contributing to the starting MSA is maintained. Such
graph representation, the _partial order graph_ , basically replaces the
linear MSA profile with a DAG whose edges represent a partial order between
the letters of the MSA to be used as an alternative to the traditional total
order of the rows of an alignment.

###  2.3 Wave Front Alignments with affine gap penalty

In [11], Gene Myers introduced the linear time and space dynamic programming
(DP) alignment for similar strings: when aligning two strings of size $N$
whose distance is known to be upper-bounded by $D$, rather than computing the
whole DP table of size $N^{2}$, only a stripe of $\cal{O}{D}$ diagonals are
computed, as their similarity guarantees their optimal alignment to lay
therein. The result is a $\cal{O}(DN)$ algorithm replacing the quadratic one.
The _Affine Gap Penalty_ function to score the cost of gaps in alignments
evaluates a series of $k$ consecutive gaps as $w(k)=o+k\cdot e$ (where $o$ is
the cost of opening the gap, and $e$ that of extending it), rather than the
sum of $k$ single gap’s costs. Using such gap penalty score is almost
mandatory in genomic sequences analysis as this forces the optimal alignment
to detect and highlight INDELs variations that typically involve several
consecutive nucleotides. In order to restore the optimal substructure of the
alignment problem when using affine gap penalty score funcion, the DP
algorithm requires three alignment tables instead of one: matrices $I,D,M$
that store the score of the best alignments ending - respectively - with an
Insertions, a Deletions, or a (Mis)Match.
In [39], the linear time and space optimization was extended to the exact
computation of an optimal pairwise alignment of strings using the affine gap
penalty score function. Roughly, for the alignment of two strings of size $N$
and $M$, the three matrices $I,D,M$ of size $N\times M$ are replaced by _wave
front_ records $I_{d,k},D_{d,k},M_{d,k}$ that store, per each score/distance
$d$ and diagonal $k$, the furthest-reaching offset in $k$ that scores $d$. The
dynamic programming recurrence is then operated on increasing values of $d$
and adjacent diagonals $k$, up to the final distance.

##  3 The DSA Algorithm

Let us first formally state the problem we solve:

String Alignment to D-strings (StoDS)
Input: A D-string $\hat{T}$ of width $W$ and size $N$, a pattern $P$ of length
$m$, and penalty scores $a,x,o,e$. Output: An optimal alignment between $P$
and $\hat{T}$ using scores: $a$ for match, $x$ for mismatch, $o$ for gap
opening, and $e$ for gap extension.

In the problem statement above, by optimal alignment we mean the alignment
between $P$ and a string $T\in\hat{T}$ (see Section 2.1) that minimizes the
distance computed by scoring $a$ for a match, $x$ for a mismatch, and gap
affine penalty function with $o$ for gap opening and $e$ for gap extension.

In this section we describe our algorithm DSA that optimally solves StoDS. We
start with Section 3.1 where we show how we adapted the _Partial Order
Alignment_ (POA) framework to work with alignment to D-strings. Then, in
Section 3.2 we describe our extension of _Wave Front Alignment_ (WFA) to
D-strings with and show how we merge the two techniques into our algorithm DSA
and analyse its complexity.

###  3.1 Partial Order Alignment with D-strings

In this section, we show how to use partial order alignment to be able to
perform base level alignments with D-strings. Without loss of generality, for
the sake of simplifying the notation, within this section we will not use the
gap affine penalty score (whose implementation with D-strings will be
described in the next section), but rather account a cost $g$ to any gap.
Given a D-string $\hat{T}$ of width $W$ and pattern $P$ of length $m$, we
build a $(W+1)\times(m+1)$ DP table $M$ having a row $i$ per each $i^{th}$
width of $\hat{T}$, and a column per each letter of $P$ (plus the usual first
row and first column for their empty prefixes), and that will store in
$M[i,j]$ the best score of aligning $P[1,j]$ to $\hat{T}[1,i]$ (an example is
shown in Figure 1). When a row falls within a non-solid position having $s$
strings of length $\ell$, then the row is not associated with a simple letter,
but rather with $s$ of them, and this will be the case for $\ell$ consecutive
rows. The entries of these rows will not contain a single score, but rather a
tuple of size $s$: if row $i$ has letters
$a_{i,1},\dots,a_{i,h},\dots,a_{i,s}$, then $M[i,j]$ is a tuple of $s$ values
$<M[i,j]^{1},\dots,M[i,j]^{h},\dots,M[i,j]^{s}>$ where, for each $1\leq h\leq
s$, a match is accounted if and only if $a_{i,h}=P[j]$.
On top of this tuple representation, a partial order is assumed for the $N$
letters of the D-string: within solid positions, the usual order of sequences
applies; in a degenerate position, instead, distinct letters in the same tuple
are not comparable, and nor is letter $a_{i,h}$ of tuple
$a_{i,1},\dots,a_{i,h},\dots,a_{i,s_{i}}$ at row $i$ comparable with the other
letters of all the $\ell$ adjacent rows composing the same degenerate
position, except for the $\ell$ other letters $a_{i^{\prime},h}$ that are also
the $h^{th}$ letter in their tuples. Within all comparable letters, the order
corresponds to that of the rows of $M$.

We use this partial order to drive the alignment along the D-string: wherever
the traditional dynamic programming alignment algorithm refers to the previous
row and/or previous column, here we refer only to entries that - instead - are
preceding according to the partial order defined above.

######  Example 3.1.

For the D-string AC[GC/AT]A the partial order is defined by the graph shown on
the letters at the rows of the DP table of Fogure 1.

We now formalize with dynamic programming recurrence relations the way we
apply the ideas sketched above. For any two letters $x,y\in\Sigma$, we define
a function $m(x,y)$ as equal to $a$ if $x=y$, and equal to $x$ otherwise.
Denoting with $s_{i}$ the size of the tuple at row $i$, we compute $M[i,j]$
for all $1\leq j\leq m$ distinguishing the following cases:

$s_{i-1}=s_{i}=1$.

In this case no degeneracy is encountered and the traditional dynamic
programming alignment framework applies.

$s_{i-1}=1$ and $s_{i}=s>1$.

In this case, row $i$ corresponds to the opening of a degenerate position of
$s_{i}$ strings of length $\ell$, this row (like the next $\ell-1$ rows)
represents $s_{i}$ alternative letters $a_{i,1},\dots,a_{i,s_{i}}$ and each
entry $M[i,j]$ will contain tuples of size $s_{i}$ computed as follows for
$1\leq h\leq s_{i}$:

$M[i,j]^{h}=min\begin{cases}M[i-1,j-1]+m(P[j],a_{i,h})\\\ M[i-1,j]+g\\\
M[i,j-1]^{h}+g\\\ \end{cases}$ (1)

$s_{i-1}=s_{i}=s>1$.

In this case, row $i$ is inside the same degenerate position of the previous
row111Should rows $i-1$ and $i$ fall into distinct degenerate positions that
have the same size only by chance, this case would be managed as in the last
case below., this row still represents $s$ alternative letters
$a_{i,1},\dots,a_{i,s_{i}}$ and each entry $M[i,j]$ will contain tuples of
size $s_{i}$ computed as follows for $1\leq h\leq s_{i}$:

$M[i,j]^{h}=min\begin{cases}M[i-1,j-1]^{h}+m(P[j],a_{i,h})\\\
M[i-1,j]^{h}+g\\\ M[i,j-1]^{h}+g\\\ \end{cases}$ (2)

$s_{i-1}=s>1$ and $s_{i}=1$.

In this case, row $i-1$ was the last letter of a degenerate position, and row
$i$ is a solid position representing a single letter $\hat{T}[i]$. The entry
$M[i,j]$ will contain a value computed as follows:

$M[i,j]=min\begin{cases}M[i-1,j-1]^{1}+m(P[j],\hat{T}[i])\\\ \dots\hskip
28.45274pt\dots\hskip 28.45274pt\dots\\\
M[i-1,j-1]^{s_{i-1}}+m(P[j],\hat{T}[i])\\\ M[i-1,j]^{1}+g\\\ \dots\hskip
28.45274pt\dots\\\ M[i-1,j]^{s_{i-1}}+g\\\ M[i,j-1]+g\\\ \end{cases}$ (3)

$s_{i-1}>1$ and $s_{i}=s>1$ but not necessarily $s_{i-1}=s_{i}$.

In this case, row $i-1$ was the last letter of a degenerate position and row
$i$ is the first of a new different degenerate position. This row represents
$s$ alternative letters $a_{i,1},\dots,a_{i,s_{i}}$ and each entry $M[i,j]$
will contain tuples of size $s_{i}$ computed as follows for $1\leq h\leq
s_{i}$:

$M[i,j]^{h}=min\begin{cases}M[i-1,j-1]^{1}+m(P[j],a_{i,h})\\\ \dots\hskip
28.45274pt\dots\hskip 28.45274pt\dots\\\
M[i-1,j-1]^{s_{i-1}}+m(P[j],a_{i,h})\\\ M[i-1,j]^{1}+g\\\ \dots\hskip
28.45274pt\dots\\\ M[i-1,j]^{s_{i-1}}+g\\\ M[i,j-1]^{h}+g\\\ \end{cases}$ (4)

For example, Figure 1 shows the DP table for the alignment of the D-string
AC[GC/AT]A against ACGTA.

Figure 1: Dynamic Programming table for the Alignment of the string ACGTA
against the D-string AC[GC/AT]A partially ordered as shown by oriented edges,
and using scores: 0 for match, 1 for mismatch, and 2 for gap.

###  3.2 Wave Front Alignment with D-strings

In this section we show how to perform Wave Front Alignment of Partially
Ordered D-Strings with affine gap penalty score in almost linear time using a
D-strings customization of Wave Front Alignment.

Table 1: Notation Symbol | Description
---|---
a | score of match
x | score of mismatch
o | score of gap opening
e | score of gap extending
d | current distance (it will start from $0$ and reach the solution of the StoDS problem)
P | query pattern: $P=P_{0}P_{1}...P_{m}$ ($m=\Theta(W)$)
$\hat{T}$ | D-string $\hat{T}=\hat{T}[0]\hat{T}[1]...\hat{T}[W]$
$WF_{d}$ | the wavefront $WF_{d}$ for $d$
$\widetilde{M}_{d},\widetilde{I}_{d},\widetilde{D}_{d}$ | components of $WF_{d}$ that will be defined for each diagonal $k$ in $WF_{d}$.
$\widetilde{M}_{d,k},\widetilde{I}_{d,k},\widetilde{D}_{d,k}$ | contain the values offsetd,k for each diagonal $k$ in $WF_{d}$.
$lo_{d},hi_{d}$ | define the width of $WF_{d}$
offsetd,k | defines for each diagonal $k$ in $WF_{d}$ the furthest reaching points with distance $d$.
abandonedd,k | boolean telling whether furthest point has been reached in diagonal $k$ for distance $d$.

Similarly to [39], we replace the three affine gap penalty dynamic programming
matrices $I,D,M$ of size $W\times m$ ($W$ being the width of $\hat{T}$ and $m$
the length of $P$) with _wave front_ records $\widetilde{I}_{d}$,
$\widetilde{D}_{d}$, $\widetilde{M}_{d}$ starting with an initial distance
$d=0$. The value of $d$ can thus only increase along the computation, and for
each partial score $d$, the records $\widetilde{I}_{d}$ (resp.
$\widetilde{D}_{d}$, $\widetilde{M}_{d}$) store the following values:
$lo_{d},hi_{d}$ and a tuple offsetd,k per each diagonal $k$ in the range
defined by $lo_{d}$ and $hi_{d}$. Indeed, the values $lo_{d},hi_{d}$ define
$WF_{d}$: the set of diagonals that allow to reach alignment score $d$. Per
each partial distance $d$ and diagonal $k$, offsetd,k stores the furthest-
reaching offset in diagonal $k$ that scores $d$. Initially, $lo_{d}=hi_{d}=0$
and the only single starting diagonal is the main diagonal $k=W-m$.

Table 1 summarizes the symbols and data (structure) we are using.

When using $WF_{d}$ with D-texts,
$\widetilde{I}_{d,k},\widetilde{D}_{d,k},\widetilde{M}_{d,k}$ are in general
tuples of size $s$ (the number of variants). Observe that ($d$, $k$, and)
offsetd,k, identify a specific diagonal and how far you can go down along it
keeping score $d$, and therefore $k$ and offsetd,k actually precisely define
an entry in what would have been the matrices $D,I$ and $M$. Let us name $u$
the row and $v$ the column of the entry defined by offsetd,k. Now, if
$\hat{T}[u]$ is a single/solid letter, then the offset of
$\widetilde{I}_{d},\widetilde{D}_{d},\widetilde{M}_{d}$ is a single value, and
else it is a tuple of size $s_{u}$: a wave front per each letter $\hat{T}[u]$.
Starting with initial distance $d=0$ and the sole central diagonal, we compute
the wavefront for growing values of $d$ up the final result. The value of $d$
is increased, and new offsets are computed, when the borders of the wave front
have been reached: no further match is possible, and either a mismatch score
$x$ must be accounted (and new offset computed along the same diagonal $k$ for
distance $d+x$), or a gap must be accounted and therefore new diagonals must
be explored (and new offset computed for distance that has increased according
to the gap penalty function).
Let us now formalize the Wave Front Alignment for D-strings with the
recurrence relations that show how to update, for each diagonal $k$, the
offsets information, We will use the generic notation offsetd,k as well as the
more specific notation offset${}^{D}_{d,k}$ (resp. offset${}^{I}_{d,k}$,
offset${}^{M}_{d,k}$) to mean the offsets in $\widetilde{D}_{d}$ (resp.
$\widetilde{I}_{d}$, $\widetilde{M}_{d}$).
The formulae below are basically those of [39] with the tuple information
added. When computing $\widetilde{I}_{d,k}$ we are assuming that the last
event has been an insertion, in which case the previous event could either
have been (i) a (mis)match, in which case a new gap is being opened, or (ii)
another insertion, in which case an existing gap is being extended. In case
(i) the new offset is that of $\widetilde{M}_{d-o-e}$ because the distance is
increased by $o+e$, while in case (ii) the new offset is that of
$\widetilde{I}_{d-e}$ because the distance is only increased by the gap
extension $e$. In both cases, the new offset has to be picked from the
preceding diagonal $k-1$. The following formula formalizes this, describing
how to compute the new offsets for $\widetilde{I}_{d,k}$, assuming the general
case in which the entry on which we do the recurrence contains a tuple of size
$s$:

$\begin{array}[]{ll}offset^{I}_{d,k}=1+max\begin{cases}\underset{\texttt{in
tuple of size s}}{\max}(offset^{M}_{d-o-e,k-1})&\text{(Open insertion)}\\\
\underset{\texttt{in tuple of size
s}}{\max}(offset^{I}_{d-e,k-1})&\text{(Extend
insertion)}\end{cases}\end{array}$ (5)

Dually, $\widetilde{D}_{d,k}$ is computed assuming that the last event has
been a deletion, in which case the previous event could either have been a
(mis)match, in which case a new gap is being opened, or another deletion, in
which case an existing gap is being extended. The new offset is then either in
$\widetilde{M}_{d-o-e}$ or in $\widetilde{D}_{d-e}$, and in both cases the new
offset has to be picked from the next diagonal $k-1$. The following formula
formalizes this, describing how to compute the new offsets for
$\widetilde{D}_{d,k}$, assuming the general case in which the entry on which
we do the recurrence contains a tuple of size $s$:

$\begin{array}[]{ll}offset^{D}_{d,k}=max\begin{cases}\underset{\texttt{in
tuple of size s}}{\max}(offset^{M}_{d-o-e,k+1})&\text{(Open deletion)}\\\
\underset{\texttt{in tuple of size
s}}{\max}(offset^{D}_{d-e,k+1})&\text{(Extend
deletion)}\end{cases}\end{array}$ (6)

Dually, $\widetilde{M}_{d,k}$ is computed assuming that the last event has
been a (mis)match, in which case we stay in the same diagonal $k$, and the
previous event could either have been another (mis)match, or an insertion or
deletion. The new offset is then either in $\widetilde{M}_{d-x}$ or in
$\widetilde{D}_{d}$ or $\widetilde{I}_{d}$. The following formula formalizes
this, describing how to compute the new offsets for $\widetilde{M}_{d,k}$,
assuming the general case in which the entry on which we do the recurrence
contains a tuple of size $s$:

$\begin{array}[]{ll}offset^{M}_{d,k}=max\begin{cases}offset^{I}_{d,k}&\text{(Insertion)}\\\
offset^{D}_{d,k}&\text{(Deletion)}\\\ \underset{\texttt{in tuple of size
s}}{\max}(offset^{M}_{d-x,k})+1&\text{(Mismatch)}\\\ \end{cases}\end{array}$
(7)

In the formula for $\widetilde{M}_{d,k}$, the case of a match is not taken
into account because offsets are recomputed only when a match is not
available.

Notice that, among possible values, the maximum is always sought here because
the distance is a fixed parameter $d$ in
$\widetilde{M}_{d,k},\widetilde{I}_{d,k},\widetilde{D}_{d,k}$, and with that
one wants to go as far as possible along $k$ maintaining that distance (that
can only increase otherwise), therefore eventually minimizing the final $d$
whose $WF_{d}$ will reach the bottom right final entry of the matrix.

In all formulae above, whenever incurring in a gap and increasing or
decreasing the diagonal number, then this has to be intended as updating
$lo_{d}$ or $hi_{d}$.

The pseudocode of DWF_ALIGN is shown in Algorithm 1 (structured as in [39] but
modified to be adapted to D-strings and their tuples based partial order
alignment realization). The function DWF_EXTEND (Algorithm 2) is called to
extend the wavefront extends $\widetilde{M}_{d,k}$ for all diagonals $k$
between $lo_{d}$ and $hi_{d}$, and it calls the function $\lambda$ (Algorithm
3), which is the core of DSA . This function is where the offsets are actually
increased. When this happens, then $\lambda$:
(i) Updates the width of the offset tuple with that of the current row to
account for the case in which we extend a wavefront from an index $i$ by 1
into an index $i+1$ such that $|\widehat{T}[i]|\neq|\widehat{T}[i+1]|$ (the
size of the tuple changes), that happens when the wavefront has been extended
into a degenerate letter or out of a degenerate letter. In this case, we need
to increase or reduce the length of offsets to correspond to the number of
variants at $\widehat{T}[i+1]$. We do this by selecting $max(\widehat{T}[i])$
and setting it as the offset for each variant at $\widehat{T}[i+1]$.
(ii) Picking the maximum from previous row does not incur into an artefact
match222With artefact match we mean to, say, align exactly pattern $P=ACGT$ to
$\widehat{T}=A\\{TC,GA\\}T$ matching the two letters $CG$ that belong to
distinct variants. thanks to the boolean vector $abandoned$ (if that part of a
variant cannot be extended any further, then we mark as abandoned the
extension of the wavefront for the given score) that will allow to know along
which variant we can extend the wavefront, and will also be used by the
traceback function to compute the correct alignment

Algorithm 1

function DWF_ALIGN($\widehat{T},P,\lambda$)

$W\leftarrow|w(\widehat{T})|$$\triangleright$ width of the D-Text

$A_{k}\leftarrow|\widehat{T}|-|P|$$\triangleright$ central diagonal

$end\leftarrow max\\{W,|P|\\}$ $\triangleright$ end of alignment index

$d\leftarrow 0$; $s_{0}\leftarrow|\widehat{T}[0]|$ $\triangleright$ initial
edit distance and width

offset${}^{M}_{0,A_{k}}\leftarrow<0,...,0>$ $\triangleright$ tuple of length
$s_{0}$ initialized to 0s

abandoned${}^{M}_{0,A_{k}}\leftarrow<false,...,false>$ $\triangleright$ tuple
of length $s_{0}$ initialized to false

while true do extend the M-wavefront

DWF_EXTEND($\widetilde{M},\lambda$, d)

$f\leftarrow\texttt{max in tuple }offsets^{M}_{d,A_{k}}$

if $f\geq end$ then break;

end if$\triangleright$ reached end of alignment

$d\leftarrow d+1$ $\triangleright$ increase current distance

DWF_NEXT($\widetilde{M},\widetilde{I},\widetilde{D}$, d) $\triangleright$
compute the bounds of the next wavefront

end while

$A\leftarrow\textsc{DWF\\_TRACEBACK}(\widetilde{M},\widetilde{I},\widetilde{D},d)$
return the edit distance and the alignment

return $(d,A)$

end function

Algorithm 2

function DWF_EXTEND($\widetilde{M},\lambda$, d)

for $k$ from $\widetilde{M}_{d}^{lo}$ to $\widetilde{M}_{d}^{hi}$ do

$ab\leftarrow abandoned^{M}_{d,k}$ $\triangleright$ tuple of abandoned
variants

$of\leftarrow offsets^{M}_{d,k}$ $\triangleright$ tuple of offsets

$f\leftarrow\texttt{max in tuple}of$ $\triangleright$ index in $\widehat{T}$

$u\leftarrow f$ $\triangleright$ index in $P$

$v\leftarrow f-k$ $\triangleright$ extend $WF_{d}$

while $\lambda$($v,u,of,ab$) do

end while

end for

end function

Algorithm 3

function $\lambda$(v, u, offsets, abandoned)

if $u<0$ or $v<0$ or $v\geq|P|$ or $u\geq W$ then return $false$

end if

$s_{prev}\leftarrow length(offsets)$ $\triangleright$ size of previous tuple

$s_{curr}\leftarrow|\widehat{T}[u]|$ $\triangleright$ size of current tuple we
are stepping into or out of a degenerate letter

if $s_{curr}\neq s_{prev}$ then furthest offset

$f\leftarrow\texttt{max in tuple }offsets$ create tuple of size $s_{curr}$
with furthest offset

$offsets\leftarrow<f,...,f>$ create tuple of $s_{curr}$ size with false

$abandoned\leftarrow<false,...,false>$

end if

$extend\leftarrow false$

$previously\\_extended\leftarrow false$ attempt to extend each variant by 1

for $s^{\prime}$ from $0$ to $s_{curr}$ do

if $\widehat{T}[h][s^{\prime}]=P[v]\textbf{and}abandoned[s^{\prime}]=false$
then

$offsets[s^{\prime}]\leftarrow offsets[s^{\prime}]+1$

$extend\leftarrow true$

if $previously\\_extended=false$ then

$v\leftarrow v+1$

$u\leftarrow h+1$

$previously\\_extended\leftarrow true$

end if

end if

if $\widehat{T}[u][s^{\prime}]\neq P[v]$ then

$abandoned[s^{\prime}]\leftarrow true$

end if

end for

return $extend$

end function

The final result we have obtained is an exact algorithm that computes the edit
distance (or any distance accounting for insertions, deletions, and
substitutions with any desired score, including affine gap penalty score),
between a string and a D-string in time and space proportional to $D\cdot N$,
where $N$ id the size of the D-string, and $D$ is the computed distance.

We remark that both the optimization problem we addressed and the dynamic
programming solution we suggested can naturally make use of weights to be
added to input data to be included in the distance function to be optimised.
These could be confidence value for bases as well as multiplicity values in
the D-string (accounting for variants that are more frequent than others).

We implemented our algorithm in a prototype tool dsa that computes the edit
distance and an optimal alignment (thus keeping all the information needed for
the tracing back the optimal path). In the next section we show some first
comparative experimental validation of the performances of DSA .

##  4 Experimental validation

###  4.1 State of the art

To the best of our knowledge, there is no published software tool implementing
an algorithm specifically designed for a global alignment of a linear string
and a D-text. Some algorithm are designed for _semi_ global alignment, that is
with the pattern being (substantially) shorter than the D-string (e.g. [9, 22,
3, 1, 2]); some of them and others are designed for more general purposes,
such as aligning a linear string to a graph structure that generalizes
D-strings (e.g.[9, 22, 2, 30, 29, 19]). Some of these and other tools are not
exact and hence do not guarantee to find the optimal alignment or solution
(e.g [29, 30, 19, 16, 17, 18]). Finally, some algorithms that basically solve
a very similar problem, actually do not compute any distance or alignment, but
rather output a consensus string that suitably merges the input strings (e.g.
[41, 34, 36, 35]).

The Variation Graph Toolkit vg suite [9, 22] contains a tool that can solve
StoDS aligning a string to a more general data structure than D-strings
(Variation Graphs indeed), but this is specifically designed to map reads that
are (much) shorter that the graph/text, rather than making a global alignment,
and therefore a comparison would not be fair. The same holds for sopang and
its extension sopang2 ([3, 1]): they are designed for _short_ reads mapping on
ED-strings, and moreover they only detect exact matches without allowing gaps
nor mismatches. Also the tool vargas ([2]) works with short reads.
We thus considered the tool abPOA ([41]), a C library tool to align a sequence
to a directed acyclic graph that also uses partial order alignment and
supports global alignment, but we could not run in on data of length
$100,000$b (accuracy comparison would not have been possible anyhow, as abPOA
does not output an alignment nor a distance, but rather a consensus string.
The same (could not run on length $100,000$b strings and output is a consensus
string) holds for the tool Astarix ([34, 36, 35]), which addresses the string
to graph alignment by solving a shortest path problem on a suitably designed
graph representation of input data.

We therefore compared the performances of DSA with those of GraphAligner ([29,
30]) and Minigraph ([19, 16, 17, 18]) for solving StoDS. Both of them are
designed to align strings on a more general graph structure that D-strings and
therefore the comparison we show below should be intended as a validation of
the performances of DSA in solving StoDS, and not as DSA being in general a
better tool than any of the other two.
GraphAligner is designed to align reads to genome graphs as a key to many
applications, including error correction, genome assembly, and genotyping of
variants in a pangenome graph. In aligning long reads to genome graphs,
GraphAligner was 13x faster and used 3x less memory when compared in [29] to
the state-of-the-art tools. As for accuracy, when employing GraphAligner for
error correction, in [29] it turned out more than twice as accurate and over
12x faster than extant tools. GraphAligner uses a seed and extend method and
the bitvector alignment extension algorithm of [30].
Minigraph is a sequence-to-graph mapper and graph constructor. It aligns a
query sequence against a sequence graph using a gap affine penalty score.
Minigraph is a well maintained and highly optimized software tool that uses
minimizers to find strong colinear chains as starting point to build the
alignment ([19, 16, 17, 18]).

###  4.2 Experiments setup

We randomly generated a D-string of width $W=100,000$b by first generating a
random string of length $W$ on alphabet $\\{A,C,G,T\\}$, and then
inserting333The insertion was done forcing the width to remain $W$ therein $g$
(input parameter to our random generator of D-strings given as a percentage of
$W$) degenerate non-solid positions as follows: for each such position $p$ we
pick at random a value for its size $s_{p}$ between $1$ and $S$ (another input
parameter), and a randomly chosen length $\ell_{p}$ between $1$ and $L$ (input
parameter again).
We generated D-strings using width $W=100,000$b in all tests, degeneracy
frequencies values $g=1\%,10\%,20\%$, maximum variance values $S=2,5$, and
maximum variant lengths $L=1,4$. As a consequence the width of tested
D-strings will always be $W=100,000$b, while its total size $N$ will be
greater than $W$ and depending from input parameters $g,S,L$.
To this purpose, we implemented the random D-strings generator simed, which is
publicly available at https://github.com/urbanslug/simed/.
From the obtained synthetic D-string $\hat{T}$, we extracted a ground truth
exact pattern $P_{0}$ of size $W$ (that is, a string $P\in\hat{T}$ that thus
matches $\hat{T}$ with distance $0$), and we (possibly) modified $P_{0}$ into
the actual input query $P$ with different possible divergences using real .vcf
files.
The divergences we tested were to insert no divergence at all, or $0,1\%$ or
$1\%$ SNPs, or $0,1\%$ INDELs (percentages are on $W$). Hence, the size
$m=|P|$ of the query string will be in $\Theta(W)$ and so will the distance
$d$ between $P$ and $\hat{T}$.

###  4.3 Preliminary Results

All tests were ran on a laptop (single threaded) Intel® Core™ i7-11800H × 16
with 16.0 GiB RAM. Space and time was measured using /usr/bin/time -f"%S\t%M"
to extract system time (seconds) and maximum resident set size (kbytes). In
all tests below we used alignment scores $a=0,x=1,o=2,e=1$.

In all experiments, for all tools time was below $0.1$ seconds and therefore
negligible and not reported case per case. We do report, instead the memory
peak (weakness of the current prototype of DSA : as we mention in further work
session, we plan to improve that), the number of mistakes in the alignment
found (only for the cases with no INDELs because there the alignment is
affected by differences in penalty scores), and finally the last column shows
the number of detected events on optimal alignment: number M of matches, X for
mismatches, I for insertions, and D for deletions. When the pattern divergence
is only SNPs, then for a correct alignment it must be $I=D=0$, $X$
approximately equal to the number of SNPs (except where by chance the
divergence will not change the DNA base), and $M=W-X$. For the experiments
involving INDELs divergence in the pattern, we also report the number $G$ of
gaps that are opened: with $W=100,000$b and $0.1\%INDELs$, in a correct
alignment it must be $G=100$.

Table 2 below shows results with a D-string of size $N=101,000$ generated with
lowly frequent low degeneracy (only $g=1\%$ of positions are degenerate, with
at most $S=2$ variants of length $L=1$), and with four different pattern
divergence: no divergence, $0.1\%$ of SNPs, $1\%$ of SNPs, and $0.1\%$ of
INDELs.

Table 2: Performances with D-string of size $N=101,000$ with degeneracy $g=1\%$, $S=2$, $L=1$. Tool | Pattern | Peak memory | Total Error | Alignment
---|---|---|---|---
| divergence | (kbytes) | (absolute) | events
DSA | none | 3162 | 0 | 100,000M 0X 0I 0D
GraphAligner | " | 6244 | 0 | 100,000M 0X 0I 0D
Minigraph | " | 5607 | 9 | 99,991M 9X 0I 0D
DSA | 0.1 % SNPs | 42209 | 0 | 99902M 98X 0I 0D
GraphAligner | " | 23972 | 1 | 99903M 97X 0I 0D
Minigraph | " | 8058 | 2 | 99891M 100X 0I 0D
DSA | 1 % SNPs | 474963 | 0 | 99003M 997X 0I 0D
GraphAligner | " | 23960 | 10 | 99001M 995X 4I 4D
Minigraph | " | 7928 | 3 | 98991M 1000X 0I 0D
DSA | 0.1 % INDELs | 99322 | 0 | 99865M 0X 106I 133D 100G
GraphAligner | " | 23973 | - | 99757M 111X 105I 132D 140G
Minigraph | " | 8162 | - | 99865M 0X 106I 133D 100G

Table 3 below shows results with a D-string of size $N=106,147$ generated with
lowly frequent high degeneracy (only $g=1\%$ of positions are degenerate, with
up to $S=5$ variants of length up to $L=4$), and with four different pattern
divergence: no divergence, $0.1\%$ of SNPs, $1\%$ of SNPs, and $0.1\%$ of
INDELs.

Table 3: Performances with D-string of size $N=106,147$ with degeneracy $g=1\%,S=5,L=4$. Tool | Pattern | Peak memory | Total Error | Alignment
---|---|---|---|---
| divergence | (kbytes) | (absolute) | events
DSA | none | 3188 | 0 | 100,000= 0X 0I 0D
GraphAligner | " | 7418 | 0 | 100,000= 0X 0I 0D
Minigraph | " | 5585 | 11 | 99989= 0X 0I 0D
DSA | 0.1 % SNPs | 44096 | 0 | 99900= 100X 0I 0D
GraphAligner | " | 26799 | 1 | 99901M 99X 0I 0D
Minigraph | " | 8038 | 11 | 99889M 100X 0I 0D
DSA | 1 % SNPs | 475668 | 0 | 99007= 993X 0I 0D
GraphAligner | " | 26682 | 6 | 99005M 993X 2I 2D
Minigraph | " | 7910 | 11 | 98989M 1000X 0I 0D
DSA | 0.1 % INDELs | 131311 | 0 | 99813M 0X 135I 174D 100G
GraphAligner | " | 26816 | - | 99717M 109X 135I 174D 145G
Minigraph | " | 8213 | - | 99813M 1X 136I 175D 100G

Table 4 below shows results with a D-string of size $N=110,000$ generated with
medium frequent low degeneracy ($g=10\%$ of positions are degenerate, with
only up to $S=2$ variants of length $L=1$), and with four different pattern
divergence: no divergence, $0.1\%$ of SNPs, $1\%$ of SNPs, and $0.1\%$ of
INDELs.

Table 4: Performances with D-string of size $N=110,000$ with degeneracy $g=10\%,S=2,L=1$. Tool | Pattern | Peak memory | Total Error | Alignment
---|---|---|---|---
| divergence | (kbytes) | (absolute) | events
DSA | none | 3174 | 0 | 100,000M 0X 0I 0D
GraphAligner | " | 16405 | 0 | 100,000M 0X 0I 0D
Minigraph | " | 5609 | 6 | 99994M 0X 0I 0D
DSA | 0.1 % SNPs | 46942 | 0 | 99903M 97X 0I 0D
GraphAligner | " | 50386 | 0 | 99903M 97X 0I 0D
Minigraph | " | 8078 | 12 | 99894M 100X 0I 0D
DSA | 1 % SNPs | 456316 | 0 | 99022M 978X 0I 0D
GraphAligner | " | 50250 | 12 | 99020M 976X 4I 4D
Minigraph | " | 7934 | 94 | 98947M 997X 0I 0D
DSA | 0.1 % (INDELS) | 105904 | 0 | 99806M 0X 171I 183D 100G
GraphAligner | " | 50223 | - | 99708M 109X 171I 183D 160G
Minigraph | " | 8145 | - | 99806M 0X 171I 183D 100G

Table 5 below shows results with a D-string of size $N=160,327$ generated with
medium frequent high degeneracy ($g=10\%$ of positions are degenerate, with up
to $S=5$ variants of length up to $L=4$), and with four different pattern
divergence: no divergence, $0.1\%$ of SNPs, $1\%$ of SNPs, and $0.1\%$ of
INDELs.

Table 5: Performances with D-string of size $N=160,327$ with degeneracy $g=10\%,S=5,L=4$. Tool | Pattern | Peak memory | Total Error | Alignment
---|---|---|---|---
| divergence | (kbytes) | (absolute) | events
DSA | none | 3134 | 0 | 100,000= 0X 0I 0D
GraphAligner | " | 27417 | 0 | 100,000= 0X 0I 0D
Minigraph | " | 5536 | 4 | 99996= 0X 0I 0D
DSA | 0.1 % SNPs | 66656 | 0 | 99905= 95X 0I 0D
GraphAligner | " | 76167 | 0 | 99905M 95X 0I 0D
Minigraph | " | 8046 | 14 | 99896M 100X 0I 0D
DSA | 1 % SNPs | 469600 | 0 | 99022= 978X 0I 0D
GraphAligner | " | 76129 | 13 | 99024M 973X 3I 3D
Minigraph | " | 7878 | 48 | 98804M 1000X 0I 0D
DSA | 0.1 % (INDELS) | 126051 | 0 | 99814M 0X 108I 186D 100G
GraphAligner | " | 76626 | - | 99665M 147X 110I 188D 170G
Minigraph | " | 8199 | - | 99804M 0X 108I 186D 100G

Table 6 below shows results with a D-string of size $N=120,000$ generated with
frequent but low degeneracy ($g=20\%$ of positions are degenerate, with only
up to $S=2$ variants of length $L=1$), and with four different pattern
divergence: no divergence, $0.1\%$ of SNPs, $1\%$ of SNPs, and $0.1\%$ of
INDELs.

Table 6: Performances with D-string of size $N=120,000$ with degeneracy $g=20\%,S=2$ and $L=1$. Tool | Pattern | Peak memory | Total Error | Alignment
---|---|---|---|---
| divergence | (kbytes) | (absolute) | events
DSA | none | 3166 | 0 | 100000M 0X 0I 0D
GraphAligner | " | 21224 | 0 | 100000M 0X 0I 0D
Minigraph | " | 5604 | 6 | 99994M 0X 0I 0D
DSA | 0.1 % SNPs | 53252 | 0 | 99901M 99X 0I 0D
GraphAligner | " | 76291 | 2 | 99902M 98X 0I 0D
Minigraph | " | 8078 | 8 | 99894M 100X 0I 0D
DSA | 1 % SNPs | 436404 | 0 | 99048M 952X 0I 0D
GraphAligner | " | 76002 | 6 | 99047M 951X 2I 2D
Minigraph | " | 7913 | 102 | 98994M 1000X 0I 0D
DSA | 0.1 % (INDELS) | 103666 | 0 | 99908M 0X 96I 88D 100G
GraphAligner | " | 76305 | - | 99813M 91X 150I 96D 141G
Minigraph | " | 8132 | - | 99908M 0X 96I 88D 100G

In all tables, the first experiment is for an exact match: no divergence is
introduced in the pattern, and an exact match is supposed to be found. As we
can see above, DSA always does find the exact solution with 100000 matches, 0
mismatches, and 0 gaps (like GraphAligner does) taking less memory than
Minigraph and much less than GraphAligner.
In the second and third experiments we introduced SNPs (in $0.1\%$ and $1\%$
of the positions, respectively): in all experiments, DSA is the only one which
is exact (except in two cases where for $0.1\%$ SNPs also GraphAligner is
correct but with higher memory usage) at a cost of the highest memory
consumption with $1\%$ SNPs, and only higher than Minigraph’s for $0.1\%$ SNPs
(but therein Minigraph is always the least accurate).
Finally, with INDELs, we can see that the memory usage of DSA must be
improved, but our prototype is always exact (and Minigraph almost is, as only
in one case it detects a spurious mismatch). GraphAligner is almost always
wrong which is expected given it is more general and does not aim to find an
exact alignment.

## Conclusions and further work

As mentioned above, at the moment our implementation of DSA is just a
promising prototype. We plan to improve its memory consumption using the ideas
of [23]. Also, we are working on improve its speed (as well as consequent
memory use) by removing hopeless diagonals suffixes that are currently being
kept in the range.
Also, we remark that with our dynamic programming method, it is very natural
to add weights to letters and to use a sum of weights modification of the
actual score as objective function, in order to account for possible useful
metadata such as confidence level in the query string or D-strings bases, or
abundance in the MSA represented by the D-string.

## References

  * [1] A.Cisłak and Szymon S.Grabowski.  SOPanG 2: online searching over a pan-genome without false positives.  arXiv:2004.03033 [cs], 2020.  URL: https://arxiv.org/abs/2004.03033.
  * [2] C.A.Darby, R.Gaddipati, M.C.Schatz, and B.Langmead.  Vargas: heuristic-free alignment for assessing linear and graph read aligners.  Bioinformatics, 36(12):3712–3718, 2020.
  * [3] Aleksander Cisłak, Szymon Grabowski, and Jan Holub.  SOPanG: online text searching over a pan-genome.  Bioinformatics, 34(24):4290–4292, 06 2018.  arXiv:https://academic.oup.com/bioinformatics/article-pdf/34/24/4290/27088655/bty506.pdf, doi:10.1093/bioinformatics/bty506.
  * [4] C.Lee, C.Grasso, and M.F.Sharlow.  Multiple Sequence Alignment Using Partial Order Graphs.  Bioinformatics, 18(3):452–464, 2002.
  * [5] The Computational Pan-Genomics Consortium.  Computational Pan-Genomics: Status, Promises and Challenges.  Briefings in Bioinformatics, 19(1):118–135, 2018.
  * [6] C.S.Iliopoulos and J.Radoszewski.  Truly Subquadratic-Time Extension Queries and Periodicity Detection in Strings with Uncertainties.  In 27th Annual Symposium on Combinatorial Pattern Matching (CPM), volume 54 of LIPIcs, pages 8:1–8:12, 2016.
  * [7] C.S.Iliopoulos, R.Kundu, and S.P.Pissis.  Efficient Pattern Matching in Elastic-Degenerate Texts.  In 11th International Conference on Language and Automata Theory and Applications (LATA), volume 10168 of Springer LNCS, pages 131–142, 2017.
  * [8] D.-F.Feng and R.F. Doolittle.  Progressive sequence alignment as a prerequisite to correct phylogenetic trees.  Journal of Molecular Evolution, 25(4):351–360, 1987.
  * [9] E.Garrison, J.Sirén, A.M.Novak, G.Hickey, J.M.Eizenga, E.T.Dawson, W.Jones, S.Garg, C.Markello, M.F.Lin, B.Paten, and R.Durbin.  Variation Graph Toolkit Improves Read Mapping by Representing Genetic Variation in the Reference.  Nature Biotechnology, 36(9):875–879, 2018.
  * [10] E.W.Myers.  AnO(ND) difference algorithm and its variations.  Algorithmica, 1(1):251–266, November 1986.  Number: 1.  doi:10.1007/BF01840446.
  * [11] E.W.Myers.  An O(ND) difference algorithm and its variations.  Algorithmica, 1(2):251–266, 1986.
  * [12] G.Bernardini, N.Pisanti, S.P.Pissis, and G.Rosone.  Pattern matching on elastic-degenerate text with errors.  In 24th International Symposium on String Processing and Information Retrieval (SPIRE), volume 10508 of Springer LNCS, pages 74–90, 2017.
  * [13] G.Bernardini, N.Pisanti, S.P.Pissis, and G.Rosone.  Approximate pattern matching on elastic-degenerate text.  Theoretical Computer Science, 812:109–122, April 2020.
  * [14] G.Bernardini, P.Gawrychowski, N.Pisanti, S.P.Pissis, and G.Rosone.  Even Faster Elastic-Degenerate String Matching via Fast Matrix Multiplication.  In 46th International Colloquium on Automata, Languages, and Programming (ICALP), volume 132 of LIPIcs, pages 21:1–21:15, 2019\.
  * [15] G.Bernardini, P.Gawrychowski, N.Pisanti, S.P.Pissis, and G.Rosone.  Elastic-Degenerate String Matching via Fast Matrix Multiplication.  SIAM Journal of Computing, 51(3):549–576, 2022.
  * [16] H.Li.  Minimap and Miniasm: Fast Mapping and de Novo Assembly for Noisy Long Sequences.  Bioinformatics, 32(14):2103–2110, 2016.
  * [17] H.Li.  Minimap2: Pairwise Alignment for Nucleotide Sequences.  Bioinformatics, 34(18):3094–3100, 2018.
  * [18] H.Li.  New strategies to improve minimap2 alignment accuracy.  Bioinformatics, 37(23):4572–4574, 2021.
  * [19] H.Li, X.Feng, and C.Chu.  The Design and Construction of Reference Pangenome Graphs with Minigraph.  Genome Biology, 21(265), 2020.
  * [20] H.Soldano, A.Viari, and M.Champesme.  Searching for flexible repeated patterns using a non-transitive similarity relation.  Pattern Recognition Letters, 16(3):233–246, 1995.
  * [21] IUPAC-IUB Commission on Biochemical Nomenclature.  Abbreviations and symbols for nucleic acids, polynucleotides, and their constituents.  Biochemistry, 9(20):4022–4027, 1970.
  * [22] J.M.Eizenga, A.M.Novak, E.Kobayashi, F.Villani, C.Cisar, S.Heumos, G.Hickey, V.Colonna, B.Paten, and E.Garrison.  Efficient dynamic variation graphs.  Bioinformatics, 36(21):5139–5144, 2021.
  * [23] J.M.Eizenga and B.Paten.  Improving the time and space complexity of the wfa algorithm and generalizing its scoring.  bioRxiv, 2022.  URL: https://www.biorxiv.org/content/early/2022/01/15/2022.01.12.476087.
  * [24] K.Abrahamson.  Generalized string matching.  SIAM Journal of Computing, 16(6):1039–1051, 1987.
  * [25] K.Aoyama, Y.Nakashima, T.I, S.Inenaga, H.Bannai, and M.Takeda.  Faster Online Elastic Degenerate String Matching.  In 29th Annual Symposium on Combinatorial Pattern Matching (CPM), volume 105 of LIPIcs, pages 9:1–9:10, 2018.
  * [26] M.Alzamel, L.A.K.Ayad, G.Bernardini, R.Grossi, C.S.Iliopoulos, N.Pisanti, S.P.Pissis, and G.Rosone.  Degenerate String Comparison and Applications.  In 18th International Workshop on Algorithms in Bioinformatics (WABI), volume 113 of LIPIcs, pages 21:1–21:14, 2018.
  * [27] M.Alzamel, L.A.K.Ayad, G.Bernardini, R.Grossi, C.S.Iliopoulos, N.Pisanti, S.P.Pissis, and G.Rosone.  Comparing Degenerate Strings.  Fundamenta Informaticae, 175(1-4):41–58, 2020.
  * [28] M.Crochemore, C.S.Iliopoulos, T.Kociumaka, J.Radoszewski, W.Rytter, and T.Walen.  Covering problems for partial words and for indeterminate strings.  Theoretical Computer Science, 698:25–39, 2017.
  * [29] M.Rautiainen and T.Marschall.  Graphaligner: rapid and versatile sequence-to-graph alignment.  Genome Biology, 21(253), 2020.
  * [30] M.Rautiainen, V.Mäkinen, and T.Marschall.  Bit-parallel sequence-to-graph alignment.  Bioinformatics, 35(19):3599–3607, 2019.
  * [31] N.Pisanti, H.Soldano, and M.Carpentier.  Incremental inference of relational motifs with a degenerate alphabet.  In 16th Annual Symposium on Combinatorial Pattern Matching (CPM), volume 3537 of Springer LNCS, pages 229–240, 2005.
  * [32] N.Pisanti, H.Soldano, M.Carpentier, and J.Pothier.  A relational extension of the notion of motifs: Application to the common 3d protein substructures searching problem.  Journal of Computational Biology, 16(12):1635–1660, 2009.
  * [33] P.Gawrychowski, S.Ghazawi, and G.M.Landau.  On indeterminate strings matching.  In 31st Annual Symposium on Combinatorial Pattern Matching (CPM), volume 161 of LIPIcs, pages 14:1–14:14, 2020.
  * [34] P.Ivanov, B.Bichsel, H.Mustafa, A.Kahles, G.Rätsch, and M.T.Vechev.  Astarix: Fast and optimal sequence-to-graph alignment.  In 24th Annual International Conference on Research in Computational Molecular Biology (RECOMB), volume 12074 of Springer LNCS, pages 104–119, 2020.
  * [35] P.Ivanov, B.Bichsel, and M.T.Vechev.  Fast and optimal sequence-to-graph alignment guided by seeds.  In 26th Annual International Conference on Research in Computational Molecular Biology (RECOMB), Springer LNCS In Press, 2022.
  * [36] P.Ivanov, B.Bichsel, and M.Vechev.  Fast and optimal sequence-to-graph alignment guided by seeds.  bioRxiv, 2022.  URL: https://www.biorxiv.org/content/10.1101/2021.11.05.467453v2.
  * [37] P.Peterlongo, N.Pisanti, F.Boyer, A.Pereira do Lago, and M.-F.Sagot.  Lossless filter for multiple repetitions with hamming distance.  Journal of Discrete Algorithms, 6(3):497–509, 2008.
  * [38] R.Grossi, C.S.Iliopoulos, C.Liu, N.Pisanti, S.P.Pissis, A.Retha, G.Rosone, F.Vayani, and L.Versari.  On-Line Pattern Matching on Similar Texts.  In 28th Annual Symposium on Combinatorial Pattern Matching (CPM), volume 78 of LIPIcs, pages 9:1–9:14, 2017.
  * [39] S.Marco-Sola, J.C.Moure, M.Moreto, and A.Espinosa.  Fast Gap-Affine Pairwise Alignment Using the Wavefront Algorithm.  Bioinformatics, 37(4):456–463, 2021.
  * [40] V.Carletti, P.Foggia, E.Garrison, L.Greco, P.Ritrovato, and M.Vento.  Graph-based representations for supporting genome data analysis and visualization: Opportunities and challenges.  In 12th International Workshop on Graph-Based Representations in Pattern Recognition (GbRPR), volume 11510 of Springer LNCS, pages 237–246, 2019.
  * [41] Y.Gao, Y.Liu, Y.Ma, B.Liu, Y.Wang, and Y.Xing.  abPOA: an SIMD-based C library for fast partial order alignment using adaptive band.  Bioinformatics, 37(15):2209–2211, 2021.