# A Catalog of LAMOST Variable Sources Based on Time-domain Photometry of ZTF
111Released on ***, ***, ***

Tingting Xu Center For Astrophysics, Guangzhou University, Guangzhou,
Guangdong, China, 510006  Great Bay Center, National Astronomical Data Center,
Guangzhou, Guangdong, China, 510006  Astronomy Science and Technology Research
Laboratory of Department of Education of Guangdong Province, Guangzhou,
Guangdong, China, 510006  Key Laboratory for Astronomical Observation and
Technology of Guangzhou,Guangzhou, Guangdong, China, 510006  Chao Liu Key
Laboratory of Space Astronomy and Technology, National Astronomical
Observatories, CAS, Beijing, China,100101  University of Chinese Academy of
Sciences, Beijing, China,100049  Feng Wang Center For Astrophysics, Guangzhou
University, Guangzhou, Guangdong, China, 510006  Great Bay Center, National
Astronomical Data Center, Guangzhou, Guangdong, China, 510006  Astronomy
Science and Technology Research Laboratory of Department of Education of
Guangdong Province, Guangzhou, Guangdong, China, 510006  Key Laboratory for
Astronomical Observation and Technology of Guangzhou,Guangzhou, Guangdong,
China, 510006  University of Chinese Academy of Sciences, Beijing,
China,100049  Weirong Huang Center For Astrophysics, Guangzhou University,
Guangzhou, Guangdong, China, 510006  Great Bay Center, National Astronomical
Data Center, Guangzhou, Guangdong, China, 510006  Astronomy Science and
Technology Research Laboratory of Department of Education of Guangdong
Province, Guangzhou, Guangdong, China, 510006  Key Laboratory for Astronomical
Observation and Technology of Guangzhou,Guangzhou, Guangdong, China, 510006
Hui Deng Center For Astrophysics, Guangzhou University, Guangzhou, Guangdong,
China, 510006  Great Bay Center, National Astronomical Data Center, Guangzhou,
Guangdong, China, 510006  Astronomy Science and Technology Research Laboratory
of Department of Education of Guangdong Province, Guangzhou, Guangdong, China,
510006  Key Laboratory for Astronomical Observation and Technology of
Guangzhou,Guangzhou, Guangdong, China, 510006  Ying Mei Center For
Astrophysics, Guangzhou University, Guangzhou, Guangdong, China, 510006  Great
Bay Center, National Astronomical Data Center, Guangzhou, Guangdong, China,
510006  Astronomy Science and Technology Research Laboratory of Department of
Education of Guangdong Province, Guangzhou, Guangdong, China, 510006  Key
Laboratory for Astronomical Observation and Technology of Guangzhou,Guangzhou,
Guangdong, China, 510006  Zhong Cao  Center For Astrophysics, Guangzhou
University, Guangzhou, Guangdong, China, 510006  Great Bay Center, National
Astronomical Data Center, Guangzhou, Guangdong, China, 510006  Astronomy
Science and Technology Research Laboratory of Department of Education of
Guangdong Province, Guangzhou, Guangdong, China, 510006  Key Laboratory for
Astronomical Observation and Technology of Guangzhou,Guangzhou, Guangdong,
China, 510006

###### Abstract

The identification and analysis of different variable sources is a hot issue
in astrophysical research. The Large Sky Area Multi-Object Fiber Spectroscopic
Telescope (LAMOST) spectroscopic survey has accumulated massive spectral data
but contains no information about variable sources. Although a few related
studies present variable source catalogs for the LAMOST, the studies still
have a few deficiencies regarding the type and number of variable sources
identified. In this study, we presented a statistical modeling approach to
identify variable source candidates. We first crossed the Kepler, Sloan
Digital Sky Survey (SDSS), and Zwicky Transient Facility (ZTF) catalogs to
obtain light curves data of variable and non-variable sources. The data are
then modeled statistically using commonly used variability parameters,
respectively. And then, an optimal variable source identification model is
determined using the Receiver Operating Characteristic (ROC) curve and four
credible evaluation indices such as precision, accuracy, recall, and F1_score.
Based on this identification model, a catalog of LAMOST variable sources
(including 631,769 variable source candidates with a probability greater than
95$\%$ and so on) is obtained. To validate the correctness of the catalog, we
performed a two-by-two cross-comparison with the GAIA catalog and other
published variable source catalogs. We achieved the correct rate ranging from
50% to 100%. Among the 123,756 sources cross-matched, our variable source
catalog identifies 85,669 with a correct rate of 69%, which indicates that the
variable source catalog presented in this study is credible.

Variable sources; Light curves; Variability parameters; Statistical modeling;
Cross-identification

††journal: APJS

##  1 Introduction

The research of variable sources has been one of the frontier topics in
astronomy research and is also the core of many research topics in
astrophysics (Eyer & Mowlavi (2008)). For example, eruptive and episodic
systems can improve the understanding of accretion, mass loss, and stellar
birth(Crawford (1955)). Eclipsing systems constrain exoplanet demographics,
mass transfer, binary evolution, and the mass-radius-temperature relation of
stars. Pulsating sources are essential for probing stellar structure and
stellar evolution theory(Walkowicz et al. (2009)). In addition, some eclipsing
systems and many of the most common pulsating systems (e.g., RR Lyrae,
Cepheids, and Mira variables) are the fundamental means to measure precise
distances to clusters, to the local group of galaxies, and to relic streams of
disrupted satellites around the Milky Way(Richards et al. (2011),Kim et al.
(2003)).

Although variable sources have been studying for several hundred years, the
discovery and identification of many variable sources rely heavily on the
emergence of new detection tools. With the advent of the Charge-Coupled Device
(CCD), the Optical Gravitational Lensing Experiment(OGLE) detected more than
900,000 variables during its 20-year observation(Udalski et al. (1992),Udalski
et al. (2015)). The first all-sky variability survey was carried out by the
All-Sky Automated Survey (ASAS), which contained the eclipsing binaries and
periodic pulsating sources(Pojmanski et al. (2005)). Meanwhile, a series of
catalogs of variable sources have been published based on different types of
telescope(Drake et al. (2014), Drake et al. (2017), Cioni et al. (2011)).
Subsequently, the development of time-domain astronomy has further promoted
the identification of variable sources. A number of variable source catalogs
were published, including the Wide-field Infrared Survey Explorer(WISE)
catalog of periodic variable sources(Chen et al. (2018)), the variable
catalogs of Gaia DR2 (Clementini et al. (2019)) and DR3 (Brown et al. (2020)),
the Zwicky Transient Facility(ZTF) catalog of periodic variable sources(Chen
et al. (2020)), and the catalog of over ten million variable source candidates
based on ZTF DR1(Ofek et al. (2020)).

The Large Sky Area Multi-Object Fiber Spectroscopic Telescope(LAMOST) (Zhao et
al. (2012),Cui et al. (2012)) is a Schmidt telescope with an effective
aperture of 3.6 $-$ 4.9m with the field of view of about 5∘, which was
designed to collect 4000 spectra in a single exposure (spectral resolution R
$\sim$ 1800, limiting magnitude r $\sim$ 19 mag, wavelength coverage 3700 $-$
9000Å). LAMOST survey contains two main topics: the LAMOST ExtraGAlactic
Survey (LEGAS) and the LAMOST Experiment for Galactic Understanding and
Exploration (LEGUE) survey of the Milky Way stellar structure(Paunzen et al.
(2020)). At the end of September 2020, the LAMOST DR6 catalog was officially
released to the world. It includes the 4901 observation areas and 11.27
million spectra data. LAMOST provides valuable information on atmospheric
parameters, which are essential parameters for variable source studies and can
help make more accurate classifications and improve the classification
accuracy of variable sources in combination with light curves. However, there
are no variable sources information in the current LAMOST catalog, imposing
significant limitations on using LAMOST data for variable sources studies.

The issue of variable source identification for LAMOST catalogs has been
studied in several related studies, and the corresponding catalogs have been
published. Tian et al. (2020) presented a LAMOST radial velocity variable
sources catalog. They first selected the sources observed multiple times in
LAMOST and compared the observed radial velocity variation with the simulated
radial velocity variation to estimate the probability of the radial velocity
variable source. They finally obtained 80,702 radial velocity variable sources
with a probability greater than 60%, and 3138 sources have been classified by
crossing the data from other published catalogs. In addition, cross-matching
with other catalogs is also a very effective way to identify variable sources
for the LAMOST. Ofek et al. (2020) obtained 10 million variable source data by
standard deviation and period analysis methods, which contains several new
short-period (less than 90 minutes) candidates, about 60 new dwarf nova
candidates, two candidate eclipsing systems, and so on. They also provided
cross-matching information with catalogs such as GAIA, SDSS, and LAMOST.

These works have effectively advanced the identification of variable sources
for the LAMOST. However, these works still have some minor limitations. For
example, the identification method of Tian et al. (2020) is based on the
radial velocity. The identification requires a sufficient number of
observations of the variable source candidates. This results in a small number
of variables and a low identification rate for other variable source types.
The study of Ofek et al. (2020) is mainly based on the periodic analysis
method, which works well for periodic variable sources but somewhat limits the
discovery of other types of variable sources, especially non-periodic or long-
period variables. In addition, the variable source catalogs presented by Ofek
et al. (2020) and Chen et al. (2020) are limited in the number of sources that
can be crossed with the LAMOST catalog. Although the catalog of Ofek et al.
(2020) has 10 million variable sources, only 216,007 sources have been
observed by the LAMOST. Within 216,007 sources, 58% of the data with a signal-
to-noise ratio in g band(SNRg) less than 15, and 66% of the data with SNRg
less than 20. According to the data quality statement published on the LAMOST
official website, data with SNR less than 15 generally suffer from significant
uncertainties in either radial velocity or stellar parameter estimation.

In this study, we tried a statistical modeling approach for the identification
of the LAMOST variable sources. The specific implementation is described in
section 2 in detail. We further apply these models to the data sets of LAMOST
DR6 and ZTF DR2 and get the final catalog of LAMOST variable source candidates
in section 3. In section 4, the catalog of LAMOST variable source candidates
is evaluated by cross-identifying with the other variable sources catalogs
obtained from different survey projects. And then, we make a discussion for
models and results in section 5. Finally, we present the conclusion in section
6.

##  2 The Identification Approach of Variable Sources Based on Light Curves

A statistical modeling approach was used in the study. The main idea is to
first model the variability parameters by the light curves of variable and
non-variable sources. Then the optimal model for the variability parameters is
determined by model testing. Finally, the model is used to identify the
variable source candidates. A data-flow diagram introducing modeling
identification is shown in Fig. 1. The entire statistical modeling process is
divided into four main parts: searching trustful identification parameters,
generating a sample set of light curves, creating statistical modeling for
identification, and model testing.

Figure 1: The data-flow diagram of identifying variable sources for LAMOST
sources.

###  2.1 Variability Parameters

The variability parameters of light curves are the keys for statistical
modeling, ranging from basic statistical properties (e.g., mean, standard
deviation, and so on) to more complex time series characteristics (e.g., auto-
correlation function).

According to the characteristics of light curves, we select seven variability
parameters, i.e., Std, Iter-std, $C_{\nu}$, $\kappa$, $\gamma$, MAD, and Amp,
from previous studies (Nun et al. (2015),Nun et al. (2017)). These variability
parameters had been proven to work well for classifying variable sources
through a machine learning approach (Cabral et al. (2018), Coughlin et al.
(2021), van Roestel et al. (2021)). In addition, we introduce three
variability parameters (Q, Q1, and Q2). All parameters are described as
follows.

(1) Q value(Q):

$Q=\frac{\left|{m_{max}-m_{min}}\right|}{\sqrt{{\sigma_{max}^{2}+\sigma_{min}^{2}}}}$
(1)

where $m_{max}$ and $m_{min}$ are the maximum and minimum magnitude in the
light curves, respectively. $\sigma_{max}$ and $\sigma_{min}$ are their
magnitude measurement errors.

(2)Q1 value(Q1):

Q1 is the variant form of Q. After the maximum and minimum magnitude of the
light curves are removed, we re-calculated this parameter by the same
calculation method as Q.

(3)Q2 value(Q2):

Q2 is also the variant form of Q. After removing the maximum, minimum, sub-
maximum, and sub-minimum magnitudes of the light curves, re-calculated this
parameter by the same calculation method as Q.

(4) Standard deviation(Std):

$\sigma=\sqrt{{\frac{1}{N-1}}{\sum\limits_{i}(m_{i}-\overline{m})}}$ (2)

where $N$ is the detection times of light curves from the ZTF catalog, mi is
the magnitude of each observation in the light curves, and $\overline{m}$ is
the mean of magnitude.

(5) Iterative standard deviation(Iter-std):

After calculating the standard deviation of the light curve, the data other
than the median plus or minus twice the standard deviation are removed, and
the standard deviation is re-calculated. This process is repeated until the
resulting standard deviation converges to a stable value.

(6) Coefficient of variation(Cν):

$C_{\nu}=\frac{\sigma}{\overline{m}}$ (3)

The $C_{\nu}$ is a simple variability index and is defined as the ratio of the
standard deviation to the mean magnitude. If a light curve has substantial
variability, the $C_{\nu}$ of this light curve is generally significant.

(7)Small Kurtosis($\kappa$):

Small sample kurtosis of the magnitudes:

$\kappa=\frac{N(N+1)}{(N-1)(N-2)(N-3)}\sum\limits_{i=1}^{N}(\frac{m_{i}-\hat{m}}{\sigma})^{4}-\frac{3(N-1)^{2}}{(N-2)(N-3)}$
(4)

Where $\hat{m}$ is the median of magnitude. For a normal distribution, the
small Kurtosis should be zero.

(8)Skewness($\gamma$):

The skewness of a source is defined as follow:

$\gamma=\frac{N)}{(N-1)(N-2)}\sum\limits_{i=1}^{N}(\frac{m_{i}-\hat{m}}{\sigma})^{3}$
(5)

For a normal distribution, it should be equal to zero.

(9)Median absolute deviation(MAD):

The MAD is described as the median discrepancy of the data from the median
data:

$MAD=\hat{m}(|m_{i}-\hat{m}|)$ (6)

A normal distribution should have a value of about 0.675. The interquartile
ranges of a normal distribution can be used to illustrate this.

(10)Amplitude(Amp):

The amplitude is half of the difference between the median of the maximum and
minimum 5$\%$ magnitudes. The amplitude of a set of numbers from 0 to 1000
should be 475.5.

###  2.2 Data Set Preparation

As described in the previous section, the statistical modeling approach
seriously depends on data samples from sources. Constructing credible data set
from variable and non-variable sources is the key to modeling.

1\. Variable source data set

We build a variable source data set based on the variable sources explicitly
labeled in the Kepler-related catalogs. The Kepler Space Telescope makes long-
term continuous photometric observations in specific sky regions. A series of
variable source catalogs based on these photometric data are published
(Nielsen et al. (2013), Abdul-Masih et al. (2016), Reinhold et al.
(2013),Bowman et al. (2016)). The variable sources in the Kepler catalogs are
accurately labeled. We first selected 12,151 rotating variable sources
(Nielsen et al. (2013)) and 983 pulsating variable sources (Bowman et al.
(2016)). In addition, a catalog of 2926 eclipsing binaries was selected from
the Kepler website (http://keplerebs.villanova.edu/). A total of 3,752 common
sources were obtained by cross-matching these Kepler variable sources with the
ZTF catalog, and the light curves data for these sources were obtained from
the ZTF DR2.

2\. Non-variable source data set

To ensure the correctness of the non-variable sources, we finally select the
light curves of the standard stars as the non-variable source dataset.
Standard stars are used as benchmarks in astrophysics observations such as
photometry and spectral classification, and usually, the magnitude of light
curves is constant. Based on the multiple photometric observations of the
Sloan Digital Sky Survey(SDSS), Ivezić et al. (2007) provided a standard star
catalog that includes 1.01 million non-variable unresolved objects. All
sources are located near the equator. A total of 263,670 sources are obtained
by cross-matching the SDSS standard stars catalog with ZTF DR2 catalogs. We
also collect the light curves of these sources from ZTF DR2 to build the non-
variable data set.

And then, we randomly select SDSS standard stars equal to the number of Kepler
variable sources. The total number of the sample data set that we used in the
statistical modeling is 7,504. The distribution of detection times and mean
magnitude of light curves for these sources is shown in Fig. 2. It shows the
detection times of the standard stars are significantly less than that of the
variable sources, and most of the standard stars are faint.

Figure 2: The statistical histograms of detection times and mean magnitude vs.
variable sources/standard stars distributions.

To further determine the availability of the standard stars, we counted the
magnitude error distribution of the selected standard stars (see Fig. 3). From
the statistical histograms, the magnitude errors are within the 3 sigma
level(N_3) for basically all the standard stars, which indicates that the
light curves of the standard stars are essentially constant.

Figure 3: The statistical histograms of sigma level vs. variable
sources$/$standard stars distributions.

Based on the sample data set, we calculate the values of each of the ten
variability parameters. The results show significant differences between the
values of the ten parameters calculated by the variable sources and the non-
variable sources (see Fig. 4).

Figure 4: The statistical histograms of Q and Q2 vs. variable
sources$/$standard stars distributions.

Meanwhile, considering that the standard deviation may be affected by the
anomalous data in the light curves, we removed the anomalous values larger
than two times the standard deviation from the light curve data and then re-
calculated the standard deviation. After three iterations, the standard
deviation finally calculated was treated as Iter-std. The statistical
histograms for different values of Iter-std within three iterations for each
variable source are given in Figure 5. The distribution of the variable
sources gradually concentrates toward the direction where the value of Iter-
std becomes smaller with continuous iterations. Eventually, the position with
the highest number of variable sources is stabilized at about 0.014.

Figure 5: The distribution of standard deviation under three iterations for
Kepler variable sources.

###  2.3 Statistics Modeling Based on Variability Parameters

We take Q-value as an example and introduce the corresponding modeling
process. Total 80$\%$ of the data in the sample set is used to perform
statistical modeling, and the rest are used for model testing.

1\. Source labeling

The first step is labeling all sources based on whether the source is a
variable source or not. The label of the source $L_{i}$ is given by:

$L_{i}={\begin{cases}1&\text{ if }S_{i}\text{ is a variable source }\\\
0&\text{ if }S_{i}\text{ is a non-variable source }\end{cases},i\in[1,N]}$ (7)

where S is the sample set, $S_{i}$ is one source in the sample set. N is the
number of sources in the sample set.

2\. Variability parameter calculation

For each source $S_{i}$, the corresponding variability parameter $V_{i}$ is
given by:

$V={\left\\{V_{i}=\frac{\left|{m_{max}-m_{min}}\right|}{\sqrt{{\sigma_{max}^{2}+\sigma_{min}^{2}}}},i\in[1,N]\right\\}}$
(8)

where V is the data set of the variability parameters.

3\. Variability probability calculation

At first, $V_{i}$ is sorted in ascending order, and $L_{i}$ is expanded sort
based on $V_{i}$. Then the variability probabilities $P(V_{i})$ is calculated
through the statistics of $L_{i}$ within an interval. This interval is that
the 1$\%$ data of sample set adjacent to Vi. The $P(V_{i})$ is given by:

$P(V_{i})=\frac{\left|\left\\{L_{j}|\
L_{j}=1\right\\}\right|}{\left|\left\\{L_{j}|\
L_{j}=1\right\\}\right|+\left|\left\\{L_{j}|\
L_{j}=0\right\\}\right|},j\in[i-0.01*N,i+0.01*N],i\in[1,N]$ (9)

where, $P(V_{i})$ is the variability probabilities corresponding to $V_{i}$ .
And $\left|\left\\{L_{j}|\ L_{j}=1\right\\}\right|$ is the cardinality of the
set $\left\\{L_{j}|\ L_{j}=1\right\\}$ , $\left|\left\\{L_{j}|\
L_{j}=0\right\\}\right|$ is the cardinality of the set $\left\\{L_{j}|\
L_{j}=0\right\\}$, respectively.

After the above three steps on the sample data set, we get the probabilities
of the ten variability parameters at different values (see Fig. 6). These
diagrams reflect the relationship between the variability parameters and the
corresponding variability probabilities.

Figure 6: The variability probabilities at different variability parameters.

4\. Variability parameters evaluation

We take the variability parameter Q as an example to demonstrate the
evaluating procedures. The evaluation for other variability parameters is
consistent and will not be repeated here. Based on 20$\%$ of the data in the
sample set, we counted the number of true positive (TP), true negative (TN),
false positive (FP), and false-negative (FN) under 100 variability
probabilities ranging from 0 to 1. The accuracy, precision, recall rate, and
F1-score were also calculated (Fig. 7). Meanwhile, the ROC curve area of
variability parameter Q is 0.75.

Figure 7: The performance of variability parameter Q at different variability
probabilities.

The four credible evaluation indexes of the ten variability parameters models
under the different probabilities are presented in Table 1 and Table 2.
Moreover, the areas of ROC curves are shown in Fig 8. According to the
evaluation results, the Q2 parameters have better performance than the Q, Q1,
Iter-std, Std, $C_{\nu}$, $\kappa$, $\gamma$, MAD and Amp parameters.
Therefore, our final identification model for variable source candidates is
given by parameters Q2. Of course, the other nine variability parameters can
still be used as a reference.

Table 1: The evaluation indexes of ten variability parameters (probability $\geq$ 0.68$\%$ be recognized as variable sources)  [.5]0 Parameters | Accuracy | Recall | Precision | F1 score
---|---|---|---|---
Q | 0.69 | 0.46 | 0.84 | 0.55
Q1 | 0.77 | 0.62 | 0.89 | 0.69
Q2 | 0.83 | 0.72 | 0.92 | 0.77
Std | 0.77 | 0.71 | 0.80 | 0.74
Iter-std | 0.78 | 0.71 | 0.82 | 0.74
$C_{\nu}$ | 0.76 | 0.73 | 0.78 | 0.75
$\kappa$ | 0.54 | 0.10 | 0.80 | 0.17
$\gamma$ | 0.53 | 0.09 | 0.73 | 0.16
MAD | 0.79 | 0.73 | 0.83 | 0.76
Amp | 0.76 | 0.73 | 0.78 | 0.75
Table 2: The evaluation indexes of ten variability parameters(probability $\geq$ 0.95$\%$ be recognized as variable sources)  [.5]0 Parameters | Accuracy | Recall | Precision | F1 score
---|---|---|---|---
Q | 0.53 | 0.07 | 0.98 | 0.12
Q1 | 0.64 | 0.28 | 0.99 | 0.39
Q2 | 0.69 | 0.38 | 0.99 | 0.49
Std | 0.52 | 0.04 | 1.00 | 0.07
Iter-std | 0.52 | 0.04 | 0.94 | 0.07
$C_{\nu}$ | 0.50 | 0 | NaN | 0
$\kappa$ | 0.50 | 0 | NaN | 0
$\gamma$ | 0.50 | 0 | NaN | 0
MAD | 0.51 | 0.03 | 0.88 | 0.06
Amp | 0.54 | 0.07 | 1.00 | 0.12
Figure 8: The ROC curves of the ten variability parameters.

##  3 The Identification of LAMOST Sources

The models of variable sources identification are constructed based on the
variability parameters and statistical modeling. Therefore, in the follow-up
LAMOST sources identification, we calculated the variability parameters from
the light curves data. Moreover, we can obtain the final LAMOST variable
source candidates based on the identification model, respectively.

###  3.1 Data Preparation

We select the low-resolution data from LAMOST DR6 V1 with 9.91 million
spectra, consisting of stars, galaxies, quasars, and other unknown objects. To
ensure the quality of the spectral data, we limit the SNRg $>$ 20 by referring
to the related research of LAMOST data analysis (Liu et al. (2015),Wang et al.
(2019)). Finally, 4.68 million sources with high SNRg spectra and position
information are extracted from LAMOST data and saved into an IPAC table.

The light curves data of these 4.68 million sources are matched from ZTF DR2
to identify variable sources. The ZTF DR2 contains light curves data acquired
between 2018 March and 2019 June, covering a time span of around 470
days(Masci et al. (2019),Graham et al. (2019)). The ZTF DR2 includes more than
2 billion sources, about half of which have greater than 20 observations(Bellm
et al. (2019)). These light curves data can complement each other with the
spectral data of the LAMOST, and the limiting magnitude and observational
areas of the ZTF survey are close to those of the LAMOST survey. Referring to
the pre-processing approach of the (Ofek et al. (2020),Chen et al. (2020)), we
merely use light curves of g and r band from the ZTF DR2 data set and matching
with the IPAC table that has been generated previously. And the i band is not
considered at all because there is no target source in ZTF DR2 with an
detection times larger than 20.

At the same time, we attached two additional conditions to make the
identification process of variable sources easier and quicker. 1) low-quality
images and photometry data were excluded by adopting INFOBITS $<$ 33,554,432
and Catflags $\neq$ 32,768 (Masci et al. (2019)). 2) We mainly selected
sources in the ZTF DR2 data that were observed at least 50 times because the
period’s false-alarm probability (FAP) is relatively high based on less than
50 detections. After restricting the above two conditions and eliminating
duplicated sources, we finally obtained 2.66 million LAMOST data sources.

###  3.2 Sources Identification and Variable Sources Catalog Generation

We perform variable source identification on LAMOST sources. After calculating
the variability parameters obtained from light curves for each source, we
determine whether the source is variable in the different variability
parameters. Finally, we get the LAMOST variable source candidates based on the
model of variability parameter Q2. As a result, all the catalogs of LAMOST
variable source candidates in the different probabilities and conditions are
shown in Table 3.

Table 3: The catalogs of LAMOST variable source candidates($\alpha$ is the probability) Band Probability |  $\alpha$ $\geq$ 95$\%$ |  $\alpha$ $\geq$ 68$\%$
---|---|---
g band | 434,256 | 954,308
r band | 324,139 | 738,365
Union of g and r band | 631,769 | 1,280,096

We take the catalog that contains 631,769 LAMOST variable source candidates as
an example, and to list the detailed information are shown in Table 4. This
catalog contains the parameters provided by the LAMOST survey and the
variability parameters obtained from the ZTF light curves data and
corresponding probabilities.

Table 4: LAMOST variable source candidates catalog

R.A.(J2000) | Dec.(J2000) |  ZTF$\\_$oid(g/r) | … | Q2(g/r) | P(Q2) (g/r)
---|---|---|---|---|---
333.264311 | 1.835927 | 494101100004654/494201100007520 | … | 40.00/45.04 | 0.99/0.99
330.641667 | 1.239245 | 494102300003222/494202300005493 | … | 9.51/7.52 | 0.99/0.98
330.63137 | 1.835927 | 494102200004936/494202200009307 | … | 23.54/22.27 | 0.99/0.99
44.537916 | 0.952695 | —/453204400006321 | … | —/4.95 | —/0.98
10.08529 | 41.68036 | 695111200016244/695211200039474 | … | 9.50/12.27 | 0.99/1.00
10.214958 | 40.550835 | 695111300013377/695211300020792 | … | 6.38/5.69 | 1.00/0.95
81.4197442 | 29.6917296 | —/658201300022529 | … | —/4.82 | —/0.98
333.613593 | 29.758858 | 690101400006952/691204300010198 | … | 5.14/4.98 | 0.98/0.98
333.537014 | 30.355967 | 690101400008984/— | … | 5.64/— | 0.95/—
333.974481 | 29.998455 | 691104300005030/691204300016607 | … | 9.25/11.73 | 0.99/1.00
333.768941 | 29.761166 | 691104300007551/— | … | 4.71/— | 0.97/—
330.24329 | 30.66345 | 690102200008729/690202200013445 | … | 14.87/13.69 | 0.99/1.00
330.497222 | 31.192316 | 690102200001888/690202200002783 | … | 39.74/31.76 | 0.99/0.99
333.472737 | 32.11478 | 691208300004224/690105400002562 | … | 8.65/12.08 | 0.99/1.00
334.784969 | 31.591894 | —/691208400027388 | … | —/5.01 | —/0.98
330.698645 | 28.929285 | —/644215400001694 | … | —/12.45 | —/1.00
330.50609 | 29.093706 | 644115100020865/644215100014392 | … | 6.69/10.41 | 0.99/1.00
331.126246 | 28.929237 | 644114300001023/— | … | 5.85/— | 0.96/—
331.094036 | 29.315796 | 644114200006522/644214200010240 | … | 6.04/4.76 | 0.98/0.97
333.593922 | 32.639942 | 690105100006583/691208200010470 | … | 6.66/6.11 | 0.99/0.99
… | … | … | … | … | …
… | … | … | … | … | …
… | … | … | … | … | …

  * *

20 lines are shown here for guidance regarding its form and content.

##  4 Analysis and Evaluation for the Catalog

In this section, we cross-identify variable source catalogs previously
published to evaluate our catalog of 631,769 LAMOST variable source
candidates.

There are many catalogs of variable sources obtained from different survey
projects, such as Kepler(Kirk et al. (2016)), GAIA(Clementini et al. (2019),
Mowlavi et al. (2018), Roelens et al. (2018)), LAMOST(Tian et al. (2020)) and
ZTF(Chen et al. (2020),Ofek et al. (2020)). These catalogs contained variable
sources of different types, such as period variable sources, eclipsing
binaries, RV variable sources, and so on. Thus, we use cross-identification to
search for common sources between published variable source catalogs and
LAMOST variable sources identified in the study and verify their correctness.

###  4.1 Cross-match with GAIA Variable Sources

The Gaia delivers nearly simultaneous measurements in the three observational
domains on which most stellar astronomical studies are based: astrometry,
photometry, and spectroscopy(Brown et al. (2016)). The Gaia data releases
provide accurate astrometric measurements for an unprecedented number of
objects. In particular, trigonometric parallaxes carry invaluable information
since they are required to infer stellar luminosities, which form the basis of
understanding much of stellar astrophysics.

Based on these properties, the GAIA survey data can be complemented and
verified with other telescope projects. The catalogs of LAMOST variable source
candidates obtained in the study are also cross-match with the Gaia DR2
catalog in the Tool for OPerations on Catalogues And Tables(TOPCAT)
software(Taylor (2005)). The location distribution of LAMOST variable source
candidates is obtained through Colour-absolute Magnitude Diagram(CaMD) as well
(see Fig.9).

It can be seen from Fig.9 that the LAMOST variable source candidates in our
catalog are concentrated in the main sequence belt, followed by the giant star
branch. In contrast, the white dwarf star sequence has fewer sources, which is
consistent with the overall distribution of stars.

Figure 9: The CaMD of the LAMOST variable source candidates. The G_BP - G_RP
is color of these candidates in the G band of Gaia DR2. The absolute Gaia
magnitude in the G band for individual sources are estimated through by M_G =
G + 5 + 5*$\log_{10}(\varpi/1000)$ , with G is magnitude in the G band,
$\varpi$ is the parallax in milliarcseconds(Gaia Collaboration et al. (2018))

Gaia data include the Cepheid(Cep), RR Lyrae(RRL), long-period variable(LPV),
short-period variable(SPV) (Clementini et al. (2019), Mowlavi et al. (2018),
Roelens et al. (2018)). There are 2106 common sources between the catalog of
GAIA variable sources catalogs and LAMOST identification targets, and these
sources are recognized as variable sources in our catalogs. The detection rate
is 100$\%$ through our variability parameter models. Among them, these sources
is contain the 370 LPVs, 60 Ceps, 1667 RRLs and 8 SPVs.

###  4.2 Cross-match with Catalog of Chen et al. (2020)

ZTF DR2 represents a highly suitable database for the detection and
exploration of new variable source candidates. The catalog of 781,602 periodic
variables was published by Chen et al. (2020). Comparison with previously
published catalogs shows that 621,702 objects (79.5$\%$) are newly discovered
or newly classified, including $\sim$ 700 Cepheids, $\sim$ 5000 RR Lyrae
stars, $\sim$ 15,000 $\delta$ Scuti variables, $\sim$ 350,000 eclipsing
binaries, $\sim$ 100,000 long-period variables, and about 150,000 rotational
variables.

However, we only selected the sources crossed between ZTF DR2 and LAMOST DR6
V1 and imposed a strict limit on the number of observations. Therefore, the
data set used in our work is only a tiny fraction of the complete ZTF DR2.
Nonetheless, there are 17,711 common sources between the catalog of ZTF DR2
periodic variable sources(PVS) and LAMOST identification targets. Among them,
17,305 objects are recognized as PVS under a probability greater than 95$\%$.
Based on that, the detection rate of PVS is about 98$\%$ through our
variability parameter models. The distribution of the variable sources types
in our catalog is shown in Fig. 10 through the cross match Chen et al. (2020).
It is shown that our variability parameter model is valid for the most common
types of variable sources.

Figure 10: The types distribution of variable sources (The left is Chen et al.
(2020), and the right is Our catalog cross-match with Chen et al. (2020)).

###  4.3 Cross-match with Catalog of Ofek et al. (2020)

The catalog presented by Ofek et al. (2020) includes 10 million variable
sources, the largest published catalog of ZTF variable sources. However, as
mentioned in the previous section, only 216,007 sources can be cross-matched
with LAMOST. The distribution diagram of SNRg of these sources is shown in
Fig. 11. After removing data with SNR less than 20 and data points in the
light curves less than 50, there are 71,454 common sources. Our model
successfully identified 46,073 variable sources with a success rate of 64%. We
also conducted a cross-match of the three variable catalogs. There are 7534
common sources among Ofek et al. (2020), Chen et al. (2020) and our LAMOST
identification targets. 7383 sources are included in our catalog.

Figure 11: The SNRg distribution of LAMOST variable sources presented by Ofek
et al. (2020)

We further investigated the reasons that theses sources were not identified.
We manually analyzed the light curves of a random sample of these sources. We
note that the magnitude error of each photometric observation is given in the
ZTF catalog. The magnitude variations of these unidentified sources are
basically within the observed magnitude errors. Therefore, we realize that the
difference between our results and Ofek et al. (2020) is due to the definition
of variable sources. The variable sources identified in this study must ensure
two observations with magnitude variations exceeding a factor of 2 magnitude
error. While Ofek et al. (2020) selected variable candidates mainly give
specific thresholds under the robust std or classical periodogram.

To validate this analysis, we directly examined the variable source data given
by Ofek using our given variability parameter ($P(Q2)\geq 68\%$). We randomly
selected two-variable source sets from (Ofek et al., 2020). Table 2 shows that
our model successfully identified 91% objects, proving that the Q2 is robust
for variables identification.

Table 5: The performance of variability parameter Q2 for variable sources catalog of Ofek et al. (2020) Number of variable sources that randomly selected from Ofek et al. (2020) |  Number of variable sources identified by our model | Rate
---|---|---
10,029 | 9,087 | 0.91
474,186 | 431,428 | 0.91

###  4.4 Cross-match with Catalog of Kirk et al. (2016)

The Kepler Mission provided nearly continuous monitoring of 200,000 objects
with unprecedented photometric precision. The catalog of eclipsing binary
systems within the 105 deg∘ Kepler field of view was published by Kirk et al.
(2016). This catalog lists the KIC, ephemeris, morphology, principle
parameters, and so on. There are 859 common sources between the Kepler
Eclipsing Binaries(EBs) and the variable catalog proposed in the study. Among
them, 662 sources are identified as eclipsing binaries under a probability
greater than 95$\%$. The detection rate of EBs is about 77$\%$ through our
statistical model.

###  4.5 Cross-match with Catalog of Tian et al. (2020)

LAMOST spectroscopic survey has provided $\sim$ 4.7 million unique sources
were targeted and $\sim$ 1 million stars observed repeatedly. The
probabilities of stars being RV variables are estimated by comparing the
observed radial velocity variations (RVVs) with simulated ones. The catalog of
RVVs based on the LAMOST survey has been published by Tian et al. (2020). This
catalog collects 80702 variable sources, including 77 $\%$ binary systems and
7$\%$ pulsating stars as well as 16$\%$ pollution by single stars. There are
24,092 common sources between the LAMOST RVVs catalog and LAMOST
identification targets. Among them, 12,141 objects are recognized as RVVs
under a probability greater than 95$\%$. The total detection rate of RVVs is
about 50$\%$ through our variability parameter models.

At the same time, we also performed an additional cross-checked with other
catalogs. There are 12,141 common sources between our catalog with Tian et al.
(2020) with a rate of 1.9%, 1500 common variable sources between the Chen et
al. (2020) and Tian et al. (2020) with a rate of 0.19%, and 686 common
variable sources between Ofek et al. (2020) and Tian et al. (2020) with a rate
of 0.0064%.

###  4.6 Performance of the Catalog

A summary of the numbers of common sources between the published catalogs and
LAMOST identification targets is listed in Table 6. Note that some variable
sources are identified repeatedly in different published variable sources
catalogs. There are 123,756 common sources between our sources to be
identified and the published catalogs such as GAIA DR2 long-period variable
sources, LAMOST RV variable sources, KEPLER eclipsing binaries, and ZTF DR2
period variable sources. Total 85,669 sources are detected as variable sources
in our catalog under a probability greater than 95$\%$, and are classified as
different types. The detection rate of our catalog is 69$\%$ for the variable
sources published in the referred catalogs.

Table 6: The cross-identification between our catalog and published catalogs

Published Catalogs |  Common sources with LAMOST identification targets |  Identified by our model ($\alpha$ $\geq$ 0.95) |  Identified rate
---|---|---|---
Clementini et al. (2019),etc. | 2106 | 2105 | 1.00
Chen et al. (2020) | 17711 | 17305 | 0.98
Ofek et al. (2020) | 71454 | 46073 | 0.64
Chen et al. (2020) $\ast$ Ofek et al. (2020) | 7534 | 7383 | 0.98
Kirk et al. (2016) | 859 | 662 | 0.77
Tian et al. (2020) | 24092 | 12141 | 0.50
Total | 123756 | 85669 | 0.69

  * *

$\ast$ is indicate the cross-match between two variable source catalogs, and
$\alpha$ is the probability.

##  5 Discussion

###  5.1 Data Sets with Different Bias

The data sets containing the variable sources and standard stars are from
different catalogs in the statistical modeling, which may have some potential
biases due to the number of observations, types of source, distribution of
magnitude, and so on. However, some minor biases exist from the same data set
due to the influence of observation conditions, i.e., weather conditions and
poor seeing. Therefore, constructing a completely unbiased data set is very
difficult. The study obtained the corresponding light curves of these sources
collected from the ZTF DR2 by cross-matching other catalogs, but we only rely
on the ZTF light curve data.

In addition, our non-variable sources in the study were derived from the SDSS
catalog of standard stars. Some new research work suggests that some standard
stars from different survey projects may be less standard, which may affect
the final results of this study. For this reason, in the statistical modeling,
we only randomly selected standard stars with the same number of variable and
non-variable sources. In addition, we performed some experiments to replace
the samples in the original data set with randomly selected data (1/2, 1/3,
and 1/4) from the standard star catalog. We find a slight change in the
evaluation index, but the variability parameter Q2 still performs the best,
indicating that the results of this study are plausible.

###  5.2 Variable Source Identification Method

The identification method of variable sources is a very classical problem, and
there has been a lot of preliminary research work. These methods are generally
based on multiple photometric data, calculating the magnitude variance to
determine whether it is a variable source, or using period analysis to
determine whether it is a periodic variable source. The variable source
identification method presented in the study draws fully on these methods. In
the statistical modeling, we extract the corresponding targets in the ZTF by
using the variable sources confirmed in the Kepler catalog and statistically
modeling the light curves data of these targets. Compared with the classical
method, the calculation is fast and straightforward, well-suited for variable
source identification in large-scale catalogs. Meanwhile, the experimental
results show that the obtained accuracy is at least comparable to the
conventional method. In addition, the statistical modeling method in this
paper is a general identification method for different types of variable
sources, which facilitates the classification of variable sources.

Of course, there are some limitations to this statistical modeling approach.
For example, the number of observations of a source can directly affect the
correctness of the identification. At the same time, a large accidental error
during one photometry may lead to a final misjudgment. In addition, our method
does not sensitive for some special variable sources(outbursting stars,
cataclysmic variables and so on). For these types of variable sources,
specific methods may be needed for identification, and we will consider this
as a future work.

###  5.3 The Correlation Analysis of Variability Parameters

We performed the correlation analysis on the ten variability parameters using
for modeling, and the corresponding correlation coefficient matrix is shown in
(Fig.12). We can know that linear relationships between all features to some
extent, and most of them are positively correlated. Of course, these
variability parameters still have other functional relationships, such as
square relationships, logarithmic relationships, and so on.

Figure 12: The Correlation matrix ten variability parameters for variable
source candidates in two band of ZTF (Top: g band, Bottom: r band.)

###  5.4 The Limitations of Variability Parameter Models

In the process of variable source identification, ten variability parameters
(i.e., Q, Q1, Q2, Std, Iter-std, $C_{\nu}$, $\kappa$, $\gamma$, MAD and Amp)
are calculated based on the light curves data obtained from ZTF DR2. These
parameters are intrinsic statistical properties relating to the scale(Std,
Iter-std), morphology($\kappa$, $\gamma$, Amp), or other properties. These
parameters are highly explainable and robust against bias (Cabral et al.
(2018)). However, they may also have varying utility. For example, the lack of
some statistical attributes may affect the results of variable source
identification. For example, the time statistical properties(Period) can
reflect more information in the time series data. In future work, we must
consider introducing more parameters that contain more helpful information to
identify variable sources.

##  6 Conclusion

The LAMOST survey has accumulated massive amounts of spectral data, but the
research on variable sources is limited to a certain extent due to the lack of
photometric information. We combined the light curves data obtained from ZTF
time-domain survey to identify variable source candidates in this work. A
catalog of 631,769 LAMOST variable source candidates is constructed with a
probability greater than 95$\%$, and this process is based on the variability
parameters and statistical modeling for light curves data. This catalog is a
robust database of variable source candidates. Moreover, the methods in this
work are not only for the identification of periodic variables but also for
general variables, including transient sources. They will be beneficial in a
few future time-domain survey projects, such as WFST and MEPHESTO, in
construction in China.

Meanwhile, the cross-identification and classification are carried out by
matching with the published catalogs. Total 85,669 variable sources(the
probability greater than 95$\%$) in our catalog are identified and classified,
containing the cepheid, RR lyrae, Eclipsing binaries, long-period, short-
period, RV variable sources, and so on. Although recognized as variable
sources, most of our catalog objects are not classified based on data-match
with other catalogs.

In further work, we will use spectral data and photometric information from
LAMOST and other time-domain surveys to classify the catalog of variable
sources as a follow-up to this work. Meanwhile, We will use machine learning
methods for variable sources classification based on the calculated
variability parameter of light curves data. Besides, LAMOST has officially
entered the mid-resolution surveys and accumulating massive spectral data. At
the same time, ZTF DRs will cover both increased timespans and large numbers
of exposures. These aspects will help improve the determination of the long-
period variables and increase the completeness of the catalog. In addition to
the main types of variables discussed here, there are many cataclysmic
variables, low-amplitude pulsating stars, binaries with compact objects, and
stars with exoplanets. We will identify these variables to construct a large
sample for further study based on future ZTF DRs and LAMOST DRs.

We thank the anonymous referee for valuable and helpful comments and
suggestions. This work is supported by the National SKA Program of China No
2020SKA0110300, the Joint Research Fund in Astronomy (U1831204, U1931141)
under cooperative agreement between the National Natural Science Foundation of
China (NSFC) and the Chinese Academy of Sciences (CAS), National Science
Foundation for Young Scholars (11903009). Funds for International Cooperation
and Exchange of the National Natural Science Foundation of China
(11961141001). Fundamental and Application Research Project of
Guangzhou(202102020677). The Innovation Research for the Postgraduates of
Guangzhou University (2020GDJC-D20). This work is also supported by
Astronomical Big Data Joint Research Center, co-founded by National
Astronomical Observatories, Chinese Academy of Sciences and Alibaba Cloud.

## References

  * Abdul-Masih et al. (2016) Abdul-Masih, M., Prša, A., Conroy, K., et al. 2016, AJ, 151, 101, doi: 10.3847/0004-6256/151/4/101
  * Bellm et al. (2019) Bellm, E. C., Kulkarni, S. R., Barlow, T., et al. 2019, PASP, 131, 068003, doi: 10.1088/1538-3873/ab0c2a
  * Bowman et al. (2016) Bowman, D. M., Kurtz, D. W., Breger, M., Murphy, S. J., & Holdsworth, D. L. 2016, MNRAS, 460, 1970, doi: 10.1093/mnras/stw1153
  * Brown et al. (2020) Brown, A. G., Vallenari, A., Prusti, T., et al. 2020, arXiv preprint arXiv:2012.01533
  * Brown et al. (2016) —. 2016, Astronomy & Astrophysics, 595, A2
  * Cabral et al. (2018) Cabral, J. B., Sánchez, B., Ramos, F., et al. 2018, Astronomy and Computing, 25, 213, doi: 10.1016/j.ascom.2018.09.005
  * Chen et al. (2018) Chen, X., Wang, S., Deng, L., de Grijs, R., & Yang, M. 2018, ApJS, 237, 28, doi: 10.3847/1538-4365/aad32b
  * Chen et al. (2020) Chen, X., Wang, S., Deng, L., et al. 2020, ApJS, 249, 18, doi: 10.3847/1538-4365/ab9cae
  * Cioni et al. (2011) Cioni, M. R. L., Clementini, G., Girardi, L., et al. 2011, A&A, 527, A116, doi: 10.1051/0004-6361/201016137
  * Clementini et al. (2019) Clementini, G., Ripepi, V., Molinaro, R., et al. 2019, A&A, 622, A60, doi: 10.1051/0004-6361/201833374
  * Coughlin et al. (2021) Coughlin, M. W., Burdge, K., Duev, D. A., et al. 2021, MNRAS, 505, 2954, doi: 10.1093/mnras/stab1502
  * Crawford (1955) Crawford, J. A. 1955, ApJ, 121, 71, doi: 10.1086/145965
  * Cui et al. (2012) Cui, X.-Q., Zhao, Y.-H., Chu, Y.-Q., et al. 2012, Research in Astronomy and Astrophysics, 12, 1197, doi: 10.1088/1674-4527/12/9/003
  * Drake et al. (2014) Drake, A. J., Graham, M. J., Djorgovski, S. G., et al. 2014, ApJS, 213, 9, doi: 10.1088/0067-0049/213/1/9
  * Drake et al. (2017) Drake, A. J., Djorgovski, S. G., Catelan, M., et al. 2017, MNRAS, 469, 3688, doi: 10.1093/mnras/stx1085
  * Eyer & Mowlavi (2008) Eyer, L., & Mowlavi, N. 2008, in Journal of Physics Conference Series, Vol. 118, Journal of Physics Conference Series, 012010, doi: 10.1088/1742-6596/118/1/012010
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Babusiaux, C., van Leeuwen, F., et al. 2018, A&A, 616, A10, doi: 10.1051/0004-6361/201832843
  * Graham et al. (2019) Graham, M. J., Kulkarni, S. R., Bellm, E. C., et al. 2019, PASP, 131, 078001, doi: 10.1088/1538-3873/ab006c
  * Ivezić et al. (2007) Ivezić, Ž., Smith, J. A., Miknaitis, G., et al. 2007, AJ, 134, 973, doi: 10.1086/519976
  * Kim et al. (2003) Kim, S. L., Lee, J. W., Kwon, S. G., et al. 2003, A&A, 405, 231, doi: 10.1051/0004-6361:20030630
  * Kirk et al. (2016) Kirk, B., Conroy, K., Prša, A., et al. 2016, AJ, 151, 68, doi: 10.3847/0004-6256/151/3/68
  * Liu et al. (2015) Liu, C., Cui, W.-Y., Zhang, B., et al. 2015, Research in Astronomy and Astrophysics, 15, 1137, doi: 10.1088/1674-4527/15/8/004
  * Masci et al. (2019) Masci, F. J., Laher, R. R., Rusholme, B., et al. 2019, PASP, 131, 018003, doi: 10.1088/1538-3873/aae8ac
  * Mowlavi et al. (2018) Mowlavi, N., Lecoeur-Taïbi, I., Lebzelter, T., et al. 2018, A&A, 618, A58, doi: 10.1051/0004-6361/201833366
  * Nielsen et al. (2013) Nielsen, M. B., Gizon, L., Schunker, H., & Karoff, C. 2013, A&A, 557, L10, doi: 10.1051/0004-6361/201321912
  * Nun et al. (2015) Nun, I., Protopapas, P., Sim, B., et al. 2015, arXiv e-prints, arXiv:1506.00010.  https://arxiv.org/abs/1506.00010
  * Nun et al. (2017) —. 2017, FATS: Feature Analysis for Time Series.  http://ascl.net/1711.017
  * Ofek et al. (2020) Ofek, E. O., Soumagnac, M., Nir, G., et al. 2020, MNRAS, 499, 5782, doi: 10.1093/mnras/staa2814
  * Paunzen et al. (2020) Paunzen, E., Huemmerich, S., & Bernhard, K. 2020, arXiv e-prints, arXiv:2010.03280.  https://arxiv.org/abs/2010.03280
  * Pojmanski et al. (2005) Pojmanski, G., Pilecki, B., & Szczygiel, D. 2005, arXiv preprint astro-ph/0508017
  * Reinhold et al. (2013) Reinhold, T., Reiners, A., & Basri, G. 2013, A&A, 560, A4, doi: 10.1051/0004-6361/201321970
  * Richards et al. (2011) Richards, J. W., Starr, D. L., Butler, N. R., et al. 2011, ApJ, 733, 10, doi: 10.1088/0004-637X/733/1/10
  * Roelens et al. (2018) Roelens, M., Eyer, L., Mowlavi, N., et al. 2018, A&A, 620, A197, doi: 10.1051/0004-6361/201833357
  * Taylor (2005) Taylor, M. B. 2005, in Astronomical Society of the Pacific Conference Series, Vol. 347, Astronomical Data Analysis Software and Systems XIV, ed. P. Shopbell, M. Britton, & R. Ebert, 29
  * Tian et al. (2020) Tian, Z., Liu, X., Yuan, H., et al. 2020, ApJS, 249, 22, doi: 10.3847/1538-4365/ab9904
  * Udalski et al. (1992) Udalski, A., Szymanski, M., Kaluzny, J., Kubiak, M., & Mateo, M. 1992, Acta Astronomica, 42, 253
  * Udalski et al. (2015) Udalski, A., Szymański, M. K., & Szymański, G. 2015, Acta Astron., 65, 1\.  https://arxiv.org/abs/1504.05966
  * van Roestel et al. (2021) van Roestel, J., Duev, D. A., Mahabal, A. A., et al. 2021, AJ, 161, 267, doi: 10.3847/1538-3881/abe853
  * Walkowicz et al. (2009) Walkowicz, L. M., Becker, A. C., Anderson, S. F., et al. 2009, in astro2010: The Astronomy and Astrophysics Decadal Survey, Vol. 2010, 307.  https://arxiv.org/abs/0902.3981
  * Wang et al. (2019) Wang, R., Luo, A.-l., Zhang, S., et al. 2019, PASP, 131, 024505, doi: 10.1088/1538-3873/aaf25f
  * Zhao et al. (2012) Zhao, G., Zhao, Y.-H., Chu, Y.-Q., Jing, Y.-P., & Deng, L.-C. 2012, Research in Astronomy and Astrophysics, 12, 723, doi: 10.1088/1674-4527/12/7/002


# Scaling a Variant Calling Genomics Pipeline with FaaS

Aitor Arjona  Universitat Rovira i VirgiliTarragonaSpain<EMAIL_ADDRESS>,  Arnau Gabriel-Atienza  Universitat Rovira i VirgiliTarragonaSpain
<EMAIL_ADDRESS>,  Sara Lanuza-Orna  Universitat Rovira i
VirgiliTarragonaSpain<EMAIL_ADDRESS>,  Xavier Roca-Canals  Universitat
Rovira i VirgiliTarragonaSpain<EMAIL_ADDRESS>,  Ayman Bourramouss
Universitat Rovira i VirgiliTarragonaSpain<EMAIL_ADDRESS>,  Tyler
K. Chafin  Biomathematics and Statistics ScotlandEdinburghUnited Kingdom
<EMAIL_ADDRESS>,  Lucio Marcello  Biomathematics and Statistics
ScotlandEdinburghUnited Kingdom<EMAIL_ADDRESS>,  Paolo Ribeca
Biomathematics and Statistics ScotlandEdinburghUnited Kingdom
<EMAIL_ADDRESS>and  Pedro García-López  Universitat Rovira i
VirgiliTarragonaSpain<EMAIL_ADDRESS>
(2023)

###### Abstract.

With the escalating complexity and volume of genomic data, the capacity of
biology institutions’ HPC faces limitations. While the Cloud presents a viable
solution for short-term elasticity, its intricacies pose challenges for
bioinformatics users. Alternatively, serverless computing allows for workload
scalability with minimal developer burden. However, porting a scientific
application to serverless is not a straightforward process. In this article,
we present a Variant Calling genomics pipeline migrated from single-node HPC
to a serverless architecture. We describe the inherent challenges of this
approach and the engineering efforts required to achieve scalability. We
contribute by open-sourcing the pipeline for future systems research and as a
scalable user-friendly tool for the bioinformatics community.

serverless, FaaS, workflow, genomics

††copyright: acmcopyright††journalyear: 2023††conference: 9th International
Workshop on Serverless Computing; December 11–15, 2023; Bologna,
Italy††booktitle: 9th International Workshop on Serverless Computing (WoSC
’23), December 11–15, 2023, Bologna, Italy††price: 15.00††doi:
10.1145/3631295.3631403††isbn: 979-8-4007-0455-0/23/12††ccs: Computer systems
organization Cloud computing††ccs: Applied computing Genomics

##  1\. Introduction

Genome sequence analysis is a compute- and data-intensive task. Current
sequencing techniques generate an overwhelming volume of data, with
repositories accumulating several petabytes each year. The European Nucleotide
Archive alone surpassed 40 petabytes of data in 2021 (Burgin et al., 2023),
and this trend is expected to continue in the future (Archive, 2022).
Consequently, institutions struggle to meet the ever-growing demands for
genomics workloads. Although High Performance Computing (HPC) installations
are commonplace in biology departments, they may fall short in providing
sufficient capacity to handle peak demands, reanalyze extensive public
datasets generated by large consortia, or minimize the run time needed for
specific analyses.

Faced with this situation, many institutions often turn to public Clouds in
seek of compute and storage resource elasticity. The Cloud allows to meet
exceptional capacity needs on demand, while avoiding the burden of acquiring
and maintaining physical infrastructure on a permanent basis. Although there
are genomics-centered distributed computing frameworks available (e.g., Next-
flow (Di Tommaso et al., 2017) or GATK for Apache Spark), they present
numerous challenges for bio-informatics users. Evaluating, configuring,
deploying and scaling the necessary resources and services for each
application becomes a daunting task due to the multitude of heterogeneous
Cloud providers and the complexities of their services.

Figure 1. Big picture of the serverless Variant Calling workflow.

Function-as-a-Service (FaaS) offers a compelling alternative for highly
parallel data-intensive applications (Jonas et al., 2017; Pu et al., 2019;
Shankar et al., 2018), with minimal configuration burden with the ability to
instantly scale up and down to zero, making it an ideal choice for less
experienced Cloud users to deploy workloads. Recent advancements in FaaS
services, such as Amazon Lambda’s increased resource
allocation111https://aws.amazon.com/about-aws/whats-new/2020/12/aws-lambda-
supports-10gb-memory-6-vcpu-cores-lambda-functions/, open up new possibilities
for tackling demanding workloads with serverless. However, adopting the
serverless paradigm, despite its cost-effectiveness and scalability benefits,
introduces its own set of challenges (Hellerstein et al., 2018). While some
efforts have been made to ease this transition (Spillner et al., 2018), it
remains a complex task, often necessitating substantial application re-
architecting.

Here we present a case study involving the porting of a genomics Variant
Calling workflow from HPC to a serverless architecture. Our focus lies on
demonstrating the key challenges that arise from this approach, particularly
those concerning concurrency, input data partitioning, and stateful data
dependencies. The objective of our work is to harness the massive and instant
parallelism of FaaS (Barcelona-Pons and García-López, 2021) to achieve high
scalability in a genomic variant calling pipeline, all while ensuring a
streamlined experience with minimal burden for the end user. Despite the large
data movements involved, by decomposing the workflow into fine-grained tasks
and leveraging parallelism in multi-CPU functions, we can ultimately achieve
high performance and scalability in a cost-effective way.

Our primary contribution lies in the insights gained from the engineering
efforts undertaken during the migration. In summary, serverless enables to
massively scale scientific workloads, but adapting them from HPC requires
thorough adjustment and complex workflow restructuring. We’ve open-sourced
this work, on GitHub222https://github.com/CLOUDLAB-URV/serverless-genomics-
variant-calling, both for the evaluation of serverless research systems and as
an open-source tool for scalable bioinformatics variant calling in the Cloud.

##  2\. Serverless Variant Calling pipeline

Variant calling is a key analysis step in genomics. It typically involves
“alignment” (i.e., a string similarity search) of sequencing reads stored as
FASTQ files to a reference genome stored as a FASTA file. Alignment mismatches
are filtered and form the basis for establishing (“calling”) any mutations
(“variants”) in the samples analyzed. In general, alignment is a stateful
computation requiring the storage into memory of large data structures
(“indices”) allowing for fast similarity searches in the reference genome,
which in this scenario might be as large as the human genome, or $\approx 3$
billion DNA characters, resulting in indices as large as 5-50 GB depending on
the alignment algorithm chosen.

In order to scale, we must distribute the workload, but requiring each
serverless function to store a full index for the reference genome is
unrealistic due to the memory constraints. However, one can split the FASTA
reference genome into smaller chunks and align reads to all of them, at the
price of making the computation less efficient and having to introduce a
reduce step to collect and score all alignments for each read later on. As
each sequencing read is aligned to the genome independently, one can also
split the original FASTQ file into smaller chunks in order to trivially
increase parallelism. The resulting architecture is shown in Figure 1.

The original HPC implementation involved bash scripts and CLI tools that
utilized the local file system and pipe commands to pass data between
processes. To transition this pipeline into a serverless architecture, we
employed Python code to wrap the genomic CLI tools and introduce new logic for
data partitioning and load distribution across fine-grained serverless tasks.
The implementation required 5588 LOC on top of the genomic tools of combined
Python, AWK and Bash scripts. This serverless architecture follows PyWren’s
(Jonas et al., 2017) serverless map-reduce model, featuring multiple parallel
stages with object storage for input and intermediate data passing. We
leveraged the Lithops framework (Sampe et al., 2021) for deploying functions
on AWS Lambda and used AWS S3 for storage. More in detail, there are three
main pipeline phases:

Pre-processing phase: This phase prepares the data so that it can be consumed
in a parallel and distributed manner. First, indices are generated to allow
random access read to FASTA files (more details in Section 3). Secondly, an
indexed version of the genome partition is generated by the GEM3 Indexer
(Marco-Sola et al., 2012) to speed up the following alignment of the
sequencing reads.

Map phase: This phase consists of aligning the sampled read to the reference
genome. Having partitioned both, this results in a cartesian product between
the FASTA and FASTQ partition sets. All FASTA× FASTQ chunk combinations are
aligned using GEM3 mapper (Marco-Sola et al., 2012). The indexed FASTA chunk
is used as the reference for alignment of a given FASTQ chunk, generating an
alignment file in .map format. Because any given alignment arises from a
partial scan of the genome (a given FASTA chunk), we must add an extra step
(index correction) to discard suboptimal alignments arising from the mapping
of any given read to multiple genome chunks. Alignments are converted to the
mpileup format, which collects information from every read to provide a
snapshot of the cumulative alignment results for each chromosome position,
providing the number of reads supporting each base call detected across all
overlapping alignments.

Reduce phase: This phase involves a shuffle of all alignment results for
variant calling. To distribute the load, the map output is distributed over
multiple parallel reduce functions. Data processing involves merging mpileup
data across sets of functions that share the same reference FASTA chunk, i.e.
those files that share the same chromosome positions. The merged mpileups are
used as input for the SiNPle variant caller (Ferretti et al., 2019). The
partial output from each reduce function is then concatenated to produce the
final output file. This can be used for any subsequent nucleotide mutation
analysis workflow.

Porting this pipeline to a serverless architecture for scalability poses
several technical challenges, including the need for effective and efficient
FASTA and FASTQ file partitioning with appropriate partition sizes, dealing
with blocking code due to data dependencies, and managing stateful data
shuffling. In the following sections, we describe how we addressed these
challenges.

##  3\. Unstructured data partitioning

Effective data partitioning is crucial to fully exploit parallelism on
serverless workloads. When decomposing a workload into fine-grained serverless
tasks, we must provide the partitioned data with fine granularity to avoid
bottlenecks in input data consumption. However, unstructured data formats
present challenges when accessed in parallel. As it contains sequences of
variable length, accessing an arbitrary portion of a FASTA file presents the
challenge of establishing the correct sequence identifiers and relevant
offsets.

Splitting a large FASTA file into many smaller files is not ideal as it
implies reading, partitioning, and writing back to storage the entire file,
which can become costly and complex for large files. In addition, a fixed
partition size is not ideal, as it will depend on the requirements of each
workload. As an alternative, we propose FASTA indexing to enable data
partitioning directly from Cloud Object Storage. Our approach involves
performing pre-processing to generate indexes for each FASTA file required.
This pre-processing is read-only, meaning that data is only read from storage,
and only the generated index is written to storage, leaving the original file
intact. These indexes are reusable, meaning that they can be used for
different partition sizes for different workflows. The generated indexes only
consume a minimal portion of storage, significantly reducing the cost and
complexity of data partitioning.

We implement an inverse index for FASTA files that allows to lookup contextual
genome information of any arbitrary partition, regardless of its size or
offset. The index is based on Samtools faidx
333https://www.htslib.org/doc/samtools-faidx.html. The read-only indexing
process leverages parallel serverless functions to pre-process large FASTA
files which are kept as-is in object storage. When partitioning a FASTA file,
we leverage byte-range GET requests and look up the inverse index to annotate
each partition with the required genome metadata. Thanks to this approach, we
avoid downloading the entire sequence and storing duplicated partitioned data
on object storage.

##  4\. Data dependencies and synchronous code

In HPC, a fixed number of parallel processes is commonly maintained throughout
the workflow (e.g. rank size in MPI). This, however, clashes with many FaaS
services where only concurrency is guaranteed instead of parallelism,
potentially causing related functions to not execute simultaneously. This
limitation presents challenges when porting stateful, synchronous blocking HPC
code to a serverless architecture, mainly due to task synchronization and
blocking function calls (e.g. all gather in MPI). Blocking serverless
functions can be risky due to concurrency limits: the function responsible for
releasing resources may not be allocated or may fail, potentially resulting in
a deadlock. Additionally, FaaS systems lack direct function communication,
necessitating a rendezvous service for synchronization, like a Redis server on
a Virtual Machine, defeating the whole purpose of a serverless architecture.

In our variant calling workflow, the Alignment tasks involved blocking calls
to wait for the Index Correction tasks to be run in background. An initial
implementation with blocking functions, utilizing Redis for rendezvous and
temporary data exchange proved impractical due to scalability limitations.
Primarily, large data inputs led to an excessive number of alignment tasks,
which easily surpassed the concurrency limit of the FaaS service. This
prevented the invocation of the index correction functions, resulting in a
deadlock. Additionally, relying on a single instance of Redis proved to be a
bottleneck.

To address the aforementioned issues, our implementation consists only of
asynchronous tasks that can be scheduled independently while treating FaaS as
a work queue. Tasks involving synchronization are divided into two steps,
ensuring that there are no blocking calls during execution. However, due to
the stateless nature of FaaS, intermediate data must be stored in S3 storage
for subsequent functions. Consequently, our pipeline is highly scalable.
Through object storage for intermediate data, we provide temporal decoupling
of tasks. Strict parallelism is no longer required, making the concurrency
guarantees of FaaS systems sufficient. On the negative side, this approach
results in increased data movements, potentially leading to longer execution
times and higher costs.

##  5\. Intra-function parallelism

With FaaS providers enabling the configuration of multiple vCPUs per function,
new possibilities arise for exploiting this capability by incorporating
parallel data processing within each function instance. In our pipeline, the
gem-mapper tool used in the Alignment step allows to leverage multiple threads
for parallel processing. Thus, our focus lies in investigating the possible
performance gains and challenges presented by this approach.

We can discern the following implications of intra-function parallelism.

(i) Utilizing intra-function parallelism reduces function invocations while
maintaining parallelism, and enhances data processing efficiency for multi-CPU
tasks. This, in turn, lowers the total function invocation count, creating
greater room for invoking more functions and processing larger datasets while
staying within the FaaS platform’s concurrency limit.  (ii) Intra-function
parallelism efficiency is bounded to the parallelizable parts of the process
(Amdahl’s Law (Rodgers, 1985)), which requires a thorough workflow analysis to
avoid over-provisioning functions for tasks that won’t benefit from
multiprocessing.  (iii) Assigning more data to each function results in an
increase of transferred data volume per function. This can result in higher
overheads caused by the high latency of object storage.

We evaluated different multi-CPU function configurations to assess the
performance of intra-function parallelism of the Alignment step. We measured
run time and cost across all configurations using a fixed input sample data
volume. This data was evenly distributed across all functions, each configured
with the same memory and CPU.

Figure 2. Multiprocessing performance analysis on AWS Lambda using different
CPU-Memory configurations

Figure 2 presents the results of the experiment regarding the cost-performance
relationship. Configurations with a higher number of invocations and greater
vCPU resources exhibit reduced execution times, albeit at the cost of higher
expenses. However, these configurations require a larger number of function
invocations. For large input datasets, this may surpass the FaaS platform’s
maximum concurrency limit (e.g., 1000 in AWS Lambda), potentially hampering
speedup and overall performance gains. Conversely, configurations featuring
more vCPUs and larger data volumes per function tend to be more cost-effective
but slower.

Figure 3. Average data transfer times for a single function.

For a more comprehensive analysis of this behavior, we quantified the average
run time of storage access and parallel processing across all configurations.
The results are displayed in Figure 3. As the number of parallel functions
decreases, data download time increases due to a larger data partition size
for each function, consuming more execution time. On the other hand,
increasing vCPUs per function leads to decreased execution time, indicating
efficient intra-function parallelism utilization of the genomic process.

Efficiently leveraging multi-CPU functions is linked to the function-executing
process’s parallelism capacity, which varies for different workloads and
cannot be universally applied. In this specific case, a configuration with 4
vCPUs and 7 GB of memory strikes a good balance between performance and cost.

##  6\. Stateful data movements

The reduce phase (Variant Calling) of our pipeline involves the merging of
mpileup files pertaining to the same FASTA chunk, then proceesed by SiNPle
(Ferretti et al., 2019) to provide variant calls, which are then concatenated.
To ensure scalability and due to function memory constraints, we need to
distribute the computation across different functions. Thus, an initial
shuffling phase is needed to determine data partitions for each function.

The complexity of the shuffling lies in the nature of the mpileup format,
which is composed of reads from the different genome positions found in FASTQ
that have been compared to the reference FASTA. This process necessitates
grouping genome positions (indexes) across multiple mpileup files pertaining
to the same FASTA reference chunk. However, not all indexes may exist in every
file, and their quantity can vary. The shuffling algorithm establishes a set
of ranges of indices that will form the partitions. The size of the partitions
is defined by the memory allocated for each serverless function. This size can
be variable depending on the number of functions to invoke capable of
processing all the intermediate mpileup sets. To determine the partition size,
we sampled lines to calculate the minimum and maximum size, so we can
calculate an approximation of the size of the ranges in bytes to fit into the
function’s memory. Small index ranges are grouped together in a function, in
order to keep the number of invocations low and avoid reaching the concurrency
limit. Once the indexes have been distributed, the respective reduce functions
will be invoked. Each function will fetch the ranges from the mpileup file set
assigned to it and then apply the processing using SiNPle.

We leverage S3 SELECT for the shuffling process, as it allows to perform SQL
queries over the contents of the structured text format of the mpileup files.
First, we query each mpileup set to extract the column of indexes of each one
of the files to be able to compute lines and to establish in a simple way the
ranges of indexes for each one of the partitions. Second, each reduce function
queries the file to extract directly the lines that are within the range of
indexes established in the previous phase. This service helps us avoid the
need for additional processing, which would introduce computational complexity
and time overhead. Specifically, this simplifies the map phase when pre-
processing indexes and during the reduce phase when creating partitions. On
the down side, the implementation we propose is not transparent to multiple
cloud providers making it totally dependent on Amazon Web Services.

##  7\. Evaluation

After a comprehensive review of the serverless Variant Calling pipeline
architecture and its implementation details, we will now evaluate two key
aspects:  (i) Performance evaluationwhere we will compare its execution time
to the single-node HPC version, and  (ii) Scalability evaluationwhere we will
use a larger dataset to assess the architecture’s scalability.

Figure 4. Tasks run time and concurrency.

_Experiments setup:_ The HPC configuration consists of a Dell PowerEdge R440
node equipped with two Intel Xeon Gold 5118 CPUs @ 2.30GHz, with a total of 24
cores (48 threads) and 32GB of memory. For the serverless AWS Lambda
execution, we utilized a Docker function runtime configured with $8192$ MB of
memory (which provides 4 vCPUs) and an ephemeral disk size of $10240$ MB. We
used Amazon S3 for the storage of the input and intermediate data, with a
bucket located in the same region (us-east-1) as the serverless functions.

###  7.1. Performance

As a first benchmark, we used the Trypanosome
genome444https://tritrypdb.org/tritrypdb/app/downloads/Current_Release/TbruceiTREU927/fasta/data/
and sequence reads
SRR6052133555https://trace.ncbi.nlm.nih.gov/Traces/?view=run_browser&acc=SRR6052133&display=download
from NCBI.

Table 1. Comparison of execution times between HPC and Serverless for the Variant Calling stages. _Stage_ |  _HPC_ |  _Serverless_
---|---|---
Genome Indexing |  0 min 14.20 s |  0 min 9.81 s
Alignment |  0 min 14.20 s |  0 min 48.10 s
Index correction |  - |  0 min 7.63 s
Generate mpileup |  51 min 15.79 s |  1 min 6.55 s
Index shuffle |  - |  0 min 10.73 s
Variant Calling |  54 min 5.04s s |  0 min 27.82 s
Total |  106 min 8.21 s |  2 min 50.64 s

Table 1 presents execution times for the stages detailed in Section 2. The HPC
version efficiently leverages multithreading for genome indexing and alignment
processes but experiences performance bottlenecks in the mpileup and variant
calling phases as they predominantly consist of sequential code. Conversely,
the serverless setup optimally distributes workloads, resulting in a
significant reduction in overall execution time even for sequential code.

###  7.2. Scalability

Next, we wanted to assess the scalability of our architecture by using a
larger input dataset such as the
hg19666http://hgdownload.cse.ucsc.edu/goldenpath/hg19/bigZips/ FASTA reference
genome with $20$ partitions and the SRR15068323 FASTQ read set with $60$
partitions from the NCBI archive
777https://trace.ncbi.nlm.nih.gov/Traces/?view=run_browser&acc=SRR15068323&display=download.
This represents a realistic benchmark of variant calling in the field of human
genomics. The files have an uncompressed size of 3.0 GB and 5.4 GB,
respectively.

Figure 4 illustrates the pipeline’s tasks execution time and concurrency over
time, using 1200 tasks for the alignment stage and 112 for the reduce stage.
Thanks to the asynchronous code redesign, we can comfortably launch more tasks
than the available default concurrency of 1000 functions for Lambda. In this
case, we launch a total of 1200 parallel functions distributed in two flights
of 1000 and 200 invocations respectively. This implies that we can scale and
launch even more functions if necessary for larger datasets, regardless of the
concurrency limit of the FaaS platform. Furthermore, the FaaS platform adapts
to the highly elastic parallelism throughout the pipeline for optimal cost-
effectiveness, leaving no idle resources under-utilized. The serverless
architecture allowed us to scale this workload, and achieve an execution time
of 7 minutes and 16.57 seconds. For reference, the same workload ran for
$\approx 14$ hours on the HPC setup. The total execution cost was 10.68 USD
for AWS Lambda GB-sec billing, with an additional 0.0644 USD for 32.2GB data
scanned in S3 SELECT.

##  8\. Insights

The main insights we gained from migrating this variant-calling bioinformatics
task to a serverless architecture are:

  1. (1)

Optimizing performance-cost with multi-CPU functions requires balancing intra-
and inter-function parallelism. Achieving a compromise is a non-trivial task,
since allocating more CPUs per function is beneficial for tasks that leverage
multiprocessing, but at the same time it also significantly increases
execution costs.

  2. (2)

Effective data partitioning is vital to exploit the massive parallelism
offered by FaaS. Data indexing can help in reducing the costs of partitioning,
since it prevents writing duplicated partitioned data back to storage.
However, a partition size still has to be explicitly configured, which adds
further complexity.

  3. (3)

An asynchronous and non-blocking architecture is strictly necessary to
guarantee scalability and to avoid the concurrency limits of serverless
platforms. This requires, however, a deep understanding of the pipeline, and
may result in a complex code re-architecting.

  4. (4)

Although object storage offers scalability for intermediate workflow data
storage, it suffers from slow performance and elevated costs. Its affordable
storage cost per GB is completely outweighed by the very expensive blocking
I/O from Lambda runtime billing. Specifically, in the experiment detailed in
Section 7.2, the accumulated I/O time accounts for 18% of the total runtime
during the alignment stage and 34% during the reduce stage, considering the
blocking S3 SELECT calls. Despite object storage being the state-of-practice
for serverless workloads, there still remains a need for a high-performance,
scalable storage service in the Cloud for short-term data storage.

  5. (5)

Serverless data processing services for object storage, like S3 SELECT,
provide a nice synergy with serverless FaaS workloads. Yet, S3 SELECT is
limited to a closed set of SQL expressions over semi-structured files (CSV,
JSON or Parquet). We expect Cloud providers to introduce new services similar
to S3 Select, but aimed to facilitate scientific data management in serverless
architectures.

##  9\. Related Work

The complexity of genomic workloads has led to the exploration of using
serverless computing in several works. In (John et al., 2021), the authors use
SWEEP to execute a serverless variant calling workflow. SWEEP defines
workflows as a directed acyclic graph (DAG) of stateless tasks and manages the
deployment of these tasks in serverless backends (functions or containers).
Their approach focuses on parallelizing and distributing workflows based on
the number of samples — in contrast, our approach is centered on designing a
distributed algorithm to parallelize and scale the analysis of individual
samples, ensuring scalability for large sample sizes.

Other works also explored serverless computing for various genomic algorithms,
like alignment (Niu et al., 2019), SNP genotyping data mining (Crespo-Cepeda
et al., 2019) or RNA visualization (Lee et al., 2019). These works parallelize
the algorithms, but they consist of a single parallel serverless stage, using
object storage for input and output data. In the case of (Hung et al., 2020)
they employed serverless for genome alignment but used virtual machines for
stateful computations. In our approach, we implement a fully serverless
architecture with different parallel processing stages using object storage
for data shuffling and stateful data movements.

##  10\. Conclusion

In this article, we have reviewed the implications of moving a single-node HPC
workload to a serverless architecture. Although we succeeded in the goal of
scalability, it required daunting engineering work to wrap the original code
and make it suitable for FaaS. Limitations such as the stateless nature of
functions or concurrency and resource constraints proved to be challenging. As
a result, we offer the pipeline in open source for future optimization works
and as a highly scalable tool for genomic variant calling.

###### Acknowledgements.

This work has been partially funded by the EU Horizon programme under grant
agreements No. 101086248, No. 101092644, No. 101092646, No. 101093110 and by
the Spanish Ministry of Science and Innovation and State Research Agency
(Agencia Estatal de Investigación) under grant agreement PID2019-106774RB-C22.
Aitor Arjona is a URV Martí Franquès grant fellow.

## References

  * (1)
  * Archive (2022) European Nucleotide Archive. 2022\.  Statistics.  https://www.ebi.ac.uk/ena/browser/about/statistics.
  * Arjona et al. (2023) Aitor Arjona, Arnau Gabriel-Atienza, Sara Lanuza-Orna, Xavier Roca-Canals, Ayman Bourramouss, Tyler K. Chafin, Lucio Marcello, Paolo Ribeca, and Pedro García-López. 2023\.  Scaling a Variant Calling Genomics Pipeline with FaaS. In _Proceedings of the 9th International Workshop on Serverless Computing_ (Bologna, Italy) _(WoSC ’23)_. Association for Computing Machinery, New York, NY, USA, 59–64.  https://doi.org/10.1145/3631295.3631403
  * Barcelona-Pons and García-López (2021) Daniel Barcelona-Pons and Pedro García-López. 2021.  Benchmarking parallelism in FaaS platforms.  _Future Generation Computer Systems_ 124 (2021), 268–284.  https://doi.org/10.1016/j.future.2021.06.005
  * Burgin et al. (2023) J. Burgin, A. Ahamed, C. Cummins, R. Devraj, K. Gueye, D. Gupta, V. Gupta, M. Haseeb, M. Ihsan, E. Ivanov, S. Jayathilaka, V. Balavenkataraman Kadhirvelu, M. Kumar, A. Lathi, R. Leinonen, M. Mansurova, J. McKinnon, C. O’Cathail, J. rio, S. Pesant, N. Rahman, G. Rinck, S. Selvakumar, S. Suman, S. Vijayaraja, Z. Waheed, P. Woollard, D. Yuan, A. Zyoud, T. Burdett, and G. Cochrane. 2023\.  The European Nucleotide Archive in 2022.  _Nucleic Acids Res_ 51, D1 (Jan 2023), D121–D125.
  * Crespo-Cepeda et al. (2019) Rodrigo Crespo-Cepeda, Giuseppe Agapito, Jose Vazquez-Poletti, and Mario Cannataro. 2019.  Challenges and Opportunities of Amazon Serverless Lambda Services in Bioinformatics. 663–668.  https://doi.org/10.1145/3307339.3343462
  * Di Tommaso et al. (2017) Paolo Di Tommaso, Maria Chatzou, Evan W Floden, Pablo Prieto Barja, Emilio Palumbo, and Cedric Notredame. 2017.  Nextflow enables reproducible computational workflows.  _Nature biotechnology_ 35, 4 (2017), 316–319.
  * Ferretti et al. (2019) L. Ferretti, C. Tennakoon, A. Silesian, G. Freimanis, and P. Ribeca. 2019.  SiNPle: Fast and Sensitive Variant Calling for Deep Sequencing Data.  _Genes_ (2019).
  * Hellerstein et al. (2018) Joseph M. Hellerstein, Jose Faleiro, Joseph E. Gonzalez, Johann Schleier-Smith, Vikram Sreekanti, Alexey Tumanov, and Chenggang Wu. 2018\.  Serverless Computing: One Step Forward, Two Steps Back.  arXiv:1812.03651 [cs.DC]
  * Hung et al. (2020) Ling-Hong Hung, Xingzhi Niu, Wes Lloyd, and Ka Yee Yeung. 2020\.  Accessible and interactive RNA sequencing analysis using serverless computing.  _bioRxiv_ (2020).  https://doi.org/10.1101/576199 arXiv:https://www.biorxiv.org/content/early/2020/10/03/576199.full.pdf
  * John et al. (2021) Aji John, Kathleen Muenzen, and Kristiina Ausmees. 2021\.  Evaluation of serverless computing for scalable execution of a joint variant calling workflow.  _PLOS ONE_ 16, 7 (07 2021), 1–12.  https://doi.org/10.1371/journal.pone.0254363
  * Jonas et al. (2017) Eric Jonas, Qifan Pu, Shivaram Venkataraman, Ion Stoica, and Benjamin Recht. 2017.  Occupy the Cloud: Distributed Computing for the 99%. In _Proceedings of the 2017 Symposium on Cloud Computing_ (Santa Clara, California) _(SoCC ’17)_. Association for Computing Machinery, New York, NY, USA, 445–451.  https://doi.org/10.1145/3127479.3128601
  * Lee et al. (2019) Benjamin D Lee, Michael A Timony, and Pablo Ruiz. 2019\.  DNAvisualization.org: a serverless web tool for DNA sequence visualization.  _Nucleic Acids Res_ 47, W1 (July 2019), W20–W25.
  * Marco-Sola et al. (2012) Santiago Marco-Sola, Michael Sammeth, Roderic Guigó, and Paolo Ribeca. 2012.  The GEM mapper: fast, accurate and versatile alignment by filtration.  _Nature methods_ 9, 12 (2012), 1185–1188.
  * Niu et al. (2019) Xingzhi Niu, Dimitar Kumanov, Ling-Hong Hung, Wesley Lloyd, and Ka Yee Yeung. 2019.  Leveraging Serverless Computing to Improve Performance for Sequence Comparison. 683–687.  https://doi.org/10.1145/3307339.3343465
  * Pu et al. (2019) Qifan Pu, Shivaram Venkataraman, and Ion Stoica. 2019\.  Shuffling, Fast and Slow: Scalable Analytics on Serverless Infrastructure. In _16th USENIX Symposium on Networked Systems Design and Implementation (NSDI 19)_. USENIX Association, Boston, MA, 193–206.  https://www.usenix.org/conference/nsdi19/presentation/pu
  * Rodgers (1985) David P. Rodgers. 1985\.  Improvements in Multiprocessor System Design.  _SIGARCH Comput. Archit. News_ 13, 3 (jun 1985), 225–231.  https://doi.org/10.1145/327070.327215
  * Sampe et al. (2021) Josep Sampe, Pedro Garcia-Lopez, Marc Sanchez-Artigas, Gil Vernik, Pol Roca-Llaberia, and Aitor Arjona. 2021.  Toward Multicloud Access Transparency in Serverless Computing.  _IEEE Software_ 38, 1 (2021), 68–74.  https://doi.org/10.1109/MS.2020.3029994
  * Shankar et al. (2018) Vaishaal Shankar, Karl Krauth, Qifan Pu, Eric Jonas, Shivaram Venkataraman, Ion Stoica, Benjamin Recht, and Jonathan Ragan-Kelley. 2018.  numpywren: serverless linear algebra.  arXiv:1810.09679 [cs.DC]
  * Spillner et al. (2018) Josef Spillner, Cristian Mateos, and David A. Monge. 2018\.  FaaSter, Better, Cheaper: The Prospect of Serverless Scientific Computing and HPC. In _High Performance Computing_ , Esteban Mocskos and Sergio Nesmachnow (Eds.). Springer International Publishing, Cham, 154–168.