\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[ngerman,polish,english]{babel}



\usepackage{float}








\begin{document}

\title{GeneMiner: a tool for extracting phylogenetic markers from
next-generation sequencing data}



\author[1]{Pulin Xie}%
\author[1]{Yongling Guo}%
\author[2]{Yue Teng}%
\author[3]{Wenbin Zhou}%
\author[1]{Yan Yu}%
\affil[1]{Sichuan University}%
\affil[2]{Beijing Institute of Microbiology and Epidemiology}%
\affil[3]{The University of North Carolina at Chapel Hill}%


\vspace{-1em}



  
  \date{June 14, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
The advancement of next-generation sequencing (NGS) technologies has
been revolutionary for the field of evolutionary biology. This
technology has led to an abundance of available genomes and
transcriptomes for researchers to mine. Specifically, researchers can
mine for various types of molecular markers that are vital for
phylogenetic, evolutionary, and ecological studies. Numerous tools have
been developed to extract these molecular markers from NGS data.
However, due to an insufficient number of well-annotated reference
genomes for non-model organisms, it remains challenging to obtain these
markers accurately and efficiently. Here, we present GeneMiner, an
improved and expanded version of our previous tool, Easy353. GeneMiner
combines the reference-guided de Bruijn graph assembly with seed
self-discovery and greedy extension. Additionally, it includes a
verification step using a parameter-bootstrap method to reduce the
pitfalls associated with using a relatively distant reference. Our
results using both experimental and simulation data showed GeneMiner can
accurately acquires phylogenetic molecular markers for plants using
transcriptomic, genomic, and other NGS data. GeneMiner is designed to be
user-friendly, fast, and memory efficient. Further, it is compatible
with Linux, Windows, and macOS. All source codes are publicly available
on GitHub for easy accessibility and transparency
(https://github.com/yyscu/GeneMiner).%
\end{abstract}%



\sloppy


\textbf{Introduction}

DNA markers from the nuclear and plastid genomes have been widely
applied for phylogenetic, evolutionary, and ecological studies in the
last few decades (Palmer et al., 1985; Palmer \& Thompson, 1982; Palmer
\& Zamir, 1982). Researchers have identified numerous single-to-low-copy
nuclear genes (Li et al., 2008; Li et al., 2017; Small et al., 2004; Wu
et al., 2006) to estimate the phylogeny for seed plants from
transcriptomic and genomic data. Building phylogenetic frameworks from
different types of next-generation sequencing (NGS) data and large-scale
molecular markers is becoming commonplace and fundamental for other
applied biological studies (Wang et al., 2020; Wen et al., 2020).
Molecular markers can be obtained from a variety of sequencing
resources, including RNA-seq (Wang et al., 2009), Hyb-Seq (Weitemier et
al., 2014), and shallow whole-genome sequencing (genome skimming)
(Straub et al., 2012). For example, Allen et al. (2017) obtained
single-copy orthologs from whole genome sequencing. Zhang et al. (2019)
extracted single-copy orthologs and ultra conserved elements from genome
skimming. Further, Liu et al. (2021) captured single-copy nuclear genes,
organellar genomes, and nuclear ribosomal DNA from deep genome skimming
data. aTRAM (Allen et al., 2015, 2018) exploits a BLAST-based iterative
search-and-assemble approach to extract specific genes from NGS data
(Allen et al., 2017). HybPiper (Johnson et al., 2016) and HybPhyloMaker
(F\selectlanguage{ngerman}ér \& Schmickl, 2018) can filter reads by mapping to the reference
using BWA/bowtie2 and subsequently assemble those reads into contigs. We
have also designed software for mining NGS data in the past, our
previous tool, Easy353 (Zhang et al., 2022), enables researchers to mine
Angiosperms353 (Johnson et al., 2019) genes from transcriptome and
enriched genome based on the reference-guided de Bruijn graph. However,
the aforementioned workflows and tools still present the following
challenges for retrieving phylogenetic markers: (1) Some markers have
lower/higher coverage than others, leading to unever read coverage; (2)
The stability and accuracy of assembly results depend on the chosen
reference sequence; (3) The putative paralogs in the assembly results
can lead to misestimation of branch lengths; (4) They require
high-performance computing servers and advanced bioinformatics skills.

In the realm of phylogenetics, molecular markers are characteristically
succinct, with their genomic arrangement frequently deemed
inconsequential. Single-to-low-copy orthologous genes, a small subset of
genomic data, are often used for phylogenetic studies at the genus level
or higher taxonomic levels. It is unnecessary to assemble complete and
sophisticated genome sequencing data to obtain these genes. With these
considerations, we introduce GeneMiner: a pipeline designed for the
extraction of phylogenetic markers from short reads NGS datasets. This
pipeline employs our proprietary reference-guided de Bruijn graph
construction algorithm. Our algorithm deliberately circumvents the need
for independent assembly tools such as SPAdes or Velvet. Compared to
other available tools, GeneMiner can captures gene fragments from both
transcriptome and genome skimming raw sequencing data more quickly,
accurately, and comprehensively. Importantly, GeneMiner can achieve all
of this on personal computers. Compared to our previous tool, Easy353,
GeneMiner contains several innovative features. These include: (1) no
restriction on the type of molecular markers, supporting the direct use
of sequences in GenBank format as reference sequences; (2) a
verification method to evaluate the accuracy of recovered target genes;
(3) an optimized weighted node model to accommodate distantly related
reference sequences; (4) a collection of new methods such as
re-filtering, re-assembly, and soft boundary to improve assembly
capability. Additionally, GeneMiner boasts excellent cross-platform
compatibility, supporting Windows, Mac, and Linux operating systems,
provides a user-friendly GUI interface for Windows and Mac users (Figure
1-A), and has distinct computational parameters that improve accuracy
over other tools in this category (Figure 1-B).

\textbf{Materials and methods}

\textbf{Workflow of GeneMiner}

GeneMiner is a Python-based tool used to mine phylogenetic markers from
NGS data. It consists of three major steps (Figure 2 A-C). First, it
collects reads that reliably map to the reference of closely related
taxa using a k-mer filtering method. Next, it uses a hybrid assembly
algorithm, a reference-guided de Bruijn graph assembly with seed
self-discovery and a greedy extension method, to recover the target
genes. Finally, it employs a mutation model and repeated resampling to
statistically evaluate the impact of the reference on the assembly and
confirm the assembly results.

\textbf{Reads filtering}

GeneMiner stores reads with high similarity to the reference of closely
related taxa. The more closely related the input sequences are to the
reference sequence, the better the results of the program. The process
includes two stages: building and retrieving the reference hash table
(Figure 2-A). To create the reference hash table, we split the reference
sequences of multiple genes are split into a k-mers dictionary, which is
simply a dictionary of length-k strings of nucleotide. A record of the
position information, the occurrence number, and the gene name of all
k-mers are added in the hash table. The k-mer
size,\emph{k\textsubscript{f}} , is set by the user. To use the
reference hash table, raw reads from NGS data (from Illumina, Roche-454,
ABI, or other sequencing platforms) in FASTQ format are used as input.
Similar to the reference sequences, each raw read is also split into
k-mers where the same \emph{k\textsubscript{f}} is k-mer sized as
previosuly designated by the user. If the k-mer of a read matches a
k-mer in the reference hash table, the whole read will be kept and
assigned to the target gene's filtered dataset. In this process,
GeneMiner can concurrently process multiple genes and samples, leading
to an accelerated analysis speed.

\textbf{Assembly}

Similar to Easy353, GeneMiner uses a reference-guided de Brujin graph
assembly with seed self-discovey and a greedy extension method (Figure
2-B). Using this design, a k-mer is initially chosen as a seed, acting
as the starting point for the assembly and it is self-discovered from
the datasets. After identifying the seed, a greedy extension technique
is applied, progressively appending the most probable next k-mer at each
stage. This is completed through utilization of information from the de
Bruijn graph and k-mer count. The reference sequence is utilized
throughout the program to ensure the assembled sequence resembles the
reference sequence.

The GeneMiner method is updated from our Easy353 program in the
following two ways: seed selection and the method for extending the seed
into contigs. For choosing the seed, GeneMiner splits both of the
reference sequences and filtered reads into k-mers
using\emph{k\textsubscript{a}} (different from the
\emph{k\textsubscript{f}}used in the filtering) and builds k-mer sets as
\emph{R} and \emph{T} , respectively. The intersection of the two sets
becomes the potential seed set. Seeds are weighted, and those with
higher weights are prioritized. In Easy353, only the number of reads
were considered (\emph{W\textsubscript{seed}} = \emph{C
\textsubscript{read}} ). However, in GeneMiner,
\emph{W\textsubscript{seed}} is jointly determined by both reference
sequences and reads. The weight is calculated as\(W_{\text{seed}}=C_{\text{ref}}\ast\operatorname{}{{(C}_{\text{read}}+}1)\),
where \emph{C\textsubscript{read}} represents the counts of the k-mer in
the set \emph{T} , and \emph{C\textsubscript{ref}} represents the counts
of the k-mer in the set \emph{R.} For extending the seed into contigs,
each k-mer is taken as a node and assigned a weighted score according to
its frequency and position in the reference sequences. In Easy353, the
node-weighted score only considered the position of a k-mer in the
reference sequences, and the weight formula used was\(W_{\text{node\ }}={C_{\text{read}}}^{(1-P_{R})}\).
However, in GeneMiner, the node-weighted score is calculated
as\(\text{\ \ }W_{\text{node\ }}={C_{\text{read}}\ast e}^{(1-|P_{R}-\ P_{T}|)}\), where \emph{C\textsubscript{read}} represents the
occurrences of the k-mer in the set \emph{T} , \emph{P\textsubscript{R}}
represents the average percentage distance of the k-mer in the set
\emph{R} , and\emph{P\textsubscript{T}} represents the percentage
distance of the k-mer in target gene. If a k-mer appears at multiple
positions in a reference sequence, the \emph{W\textsubscript{node}} is
downgraded to\emph{C\textsubscript{read}} to avoid the introducing of
insertion errors or duplication issues. In the case of a given seed, any
secondary branches branching off from the main path, such as ''bubbles''
in the graph, will be disregarded during sequence reconstruction as long
as they have a lower score.

\textbf{Verifying results}

One valuable step, which is often overlooked, when designing tools for
NGS data, is verification of the results. GeneMiner employs a novel
verification method based on mutation models and repeated resampling to
assess the influence of the reference and also confirm the stability of
assembly results statistically (Figure 2-C). The mutation model is
generated by calculating both the variation rate and the indel rate in
the reference sequences. Based on the initially generated contig
(\emph{Seq\textsubscript{1}} ), mutations and resamples are applied
using the mutation model to obtain \emph{ref\textsubscript{1}}
,\emph{ref\textsubscript{2}} \ldots{}. \emph{ref\textsubscript{n}} . All
simulated \emph{ref\textsubscript{1}} , \emph{ref\textsubscript{2}}
\ldots{}\emph{ref\textsubscript{n}} are implemented as the new reference
for the whole assembly process and get each new recovered target gene is
noted as \emph{Seq\textsubscript{2}} . A bootstrap score is then used to
assess the stability of the assembly results, calculated
as\(\frac{N}{L}\ast 100\), where \emph{N} represents the matched base number
in the pairwise alignment \emph{Seq\textsubscript{1}}
and\emph{Seq\textsubscript{2}} , and \emph{L} represents the alignment
length. A recovered target gene with a high \emph{bootstrap
score}indicates a stable assembly and a probably well-chosen reference
data. Conversely, a recovered target gene with a low \emph{bootstrap
score}requires further consideration for reference selection and is
probably not able to be used for phylogenetic studies.

\textbf{Parameter and performance testing}

To evaluate the overall performance of GeneMiner, we conducted two
rounds of evaluation, Test I and Test II. Initially, we input short-read
NGS data (sequence data, 50-300bp) and the sequences of target genes
from closely related taxa (reference sequences) into GeneMiner. We then
compared the results from GeneMiner with the known results, referred to
as the gold standard.

For Test I, we used two public RNA-seq data sets from \emph{Arabidopsis
thaliana} and \emph{Oryza sativa} as the sequencing files (can be found
in the Data Accessibility section). The reference sequences were derived
from Angiosperms353 of \emph{Brassicaceae} (53,379 sequences)
and\emph{Poaceae} (224,049 sequences) , excluding \emph{A. thaliana}
and\emph{O. sativa} . Angiosperms353 is a set of targeted sequencing
probes designed to reconstruct the Tree of Angiosperms. It was developed
based on 353 putatively single-copy protein-coding genes from the One
Thousand Plant Transcriptomes Initiative project (Johnson et al., 2019).
The gold standard for this test were the Angiosperms353 genes that are
specific to these two species. Thus our gold standard for this test
includes with 349 genes from \emph{A. thaliana} and 347 genes from
\emph{O. sativa.}Within each collection of reference sequences, the
variability between the gold standard gene and its closest reference
sequence ranges from 0-15.4\%, and 0-18.3\%, respectively. Both the
reference and the gold standard were downloaded from the Kew Tree of
Life Explorer at (\emph{https://treeoflife.kew.org).} Performance of
aTRAM, Easy353, GeneMiner and HybPiper were compared, evaluating speed,
memory efficiency, and disk space usage using similar parameters in Test
I. Based on the difference between the recovered target genes and the
gold standard, we categorized the results into four different levels as
follows: \emph{level1} (100\% identity and and 90-100\% coverage,
high-quality), \emph{level2} (100\% identity and 0-90\% coverage,
medium-quality), \emph{level3} (only with indels, medium-quality),
and\emph{level4} (with substitutions, low-quality).

While Test I used real data, Test II utilized simulation datasets. We
simulated different input data in order to examine the impact of data
sources and program parameters on the results, including sequencing
depth, reference and parameter-bootstrap. We constructed simulation
datasets ranging from 1 to 100x based on the transcriptome of \emph{A.
thaliana} from Phytozome version 12 (release TAIR10) (Goodstein et al.,
2012) using NGSNGS version 0.9 (Henriksen et al., 2023). We used two two
different sets of reference sequences, one set being real reference
sequences, the same as in Test I with Brassicaceae, and the other set
being simulated reference sequences. The simulated reference sequences
were generated using custom scripts based on the gold standard to
construct 31 groups of reference sequences with variances ranging from 0
to 30\%. Each group contained 347 genes, with each gene having 50
sequences. The gold standard was also constructed from the transcriptome
of \emph{A. thaliana} , including a total of 349 angiosperm genes. The
scripts used to construct the datasets and the completed simulated
reference sequences can be found in the Supplementary Materials.

In Test II, we first used data of varying depths and real reference
sequences, with and without a Weight-node model to test the impact of
using data of different depths on the GeneMiner results. In the test
that did not use the Weight-node model, the weights of the nodes were
set equal to the counts of k-mers generated by the filtered reads. We
then conducted a grid test using data with depths ranging from 1 to 100x
and simulated reference sequences with variances ranging from 0\% to
30\%. This round of evaluation looked at the impact of the distance
between reference sequences and target sequences on the results. Based
on the evaluation results, we assessed the capability of
parameter-bootstrap to validate results using data with a depth of 50x
and 0-30\% variance. In Test II, we used the same four levels as
outlined for Test I to categorize the assembly results. All codes,
reference sequences, and test data involved in Test I and Test II are
available on the following GitHub page:
https://github.com/yyscu/GeneMiner-Test.

\textbf{Results}

\textbf{Performance comparison of different tools}

In Test I, we evaluated the operational performance and computational
results of GeneMiner with those of aTRAM, Easy353, and HybPiper (Table
1). Figure 1-B displays the average values for the speed, memory usage,
and disk space usage required for each tool when analyzing two
species:\emph{A. thaliana} and \emph{O. sativa.} GeneMiner exhibited the
fastest speed, being on average 1.8, 4.52, and 42.62 times faster than
Easy353, HybPiper, and aTRAM, respectively. GeneMiner also exhibited
disk usage under 1 GB and higher memory usage than both HybPiper and
aTRAM, ranging between (6.91-7.78 GB). For the assembly results (Figure
3-A), GeneMiner acquired the highest number of bases (382,013 and
417,976), provided the highest accuracy (average 99.17\%) and sequence
coverage (average 99.92\%), and the lowest indel rate (average 0.64\%)
and substitution rate (average 0.20\%). In terms of the number of genes
acquired, for both species, aside from aTRAM which performed poorly on
\emph{A. thaliana} , securing only 303 genes, all other methods managed
to capture more than 90\% of the genes. Among them, GeneMiner had the
highest recovery rate, recovering 349 and 346 genes from \emph{A.
thaliana} and\emph{O. sativa} , respectively. The number of high-quality
sequences, rated as level 1 (303 and 313), was significantly higher than
the other software packages, and the number of low-quality sequences,
rated as level 4 (20 and 18), was significantly lower compared to the
other tools.

\textbf{Features of GeneMiner}

The results for assessing the impact of sequencing depth on assembly
performance are stored in Table S1 (Weighted-node) and Figure 3-B. As
the sequencing depth increased from 1 to 100x, the number of
high-quality results (\emph{level1} ) increased from 24 to 297 while the
sum of high-quality and medium-quality (\emph{levels 1\&2} ) results
increased from 270 to 330. The overall accuracy rate also improved
significantly with increasing sequencing depth, rising from 82.84\% to
99.47\%. The assembly performance was observed to be adequate when the
sequencing depth exceeded 20x, with statistics gradually trending
towards convergence as it surpassed 50x (Figure 3-B). Notably, any
increase in sequencing deepth above 50x only had a marginal effect on
improving the assembly. Compared to using counts of k-mers, the
employment of the Weighted-node model increases the number of
high-quality results. Additionally, in instances where the depth is
above 10x, both the indel rate and substitution rate are notably reduced
(Figure 3-B).

The results of the grid test with different variances of reference
sequences from 0 to 30\% and depth variations from 1 to 100x are stored
in Table S2. The heatmap of the results for high-quality genes obtained
(\emph{level 1} ) can be seen in Figure 3-D, and the heatmaps
for\emph{levels 2} through \emph{4,} indel rate, substitution rate, and
recovered genes are in Figure S1. The results in Figure 3-D show that
when the sequencing depth is below 10x, or the variance in reference
sequences is above 18\%, the number of high-quality genes obtained is
less than 50\% of the total number of genes, and the number of
successfully recovered genes is also greatly reduced. These data
demonstrate that the sequencing depth and the variance in the reference
sequences are the most important factors that affect gene recovery.

The results of the evaluation of the bootstrap parameter are shown in
Table S3. For each set of results, we conducted a detailed breakdown
based on bootstrap scores, dividing them into five parts: {[}0, 50),
{[}50,90), {[}90, 99), {[}99,100), and 100. We then compared the changes
in the number of genes of different qualities, substitution rate, and
other parameters after removing genes with different scores. The results
indicate that a higher bootstrap score (above 90) contains the majority
of high-quality (\emph{level 1} ) and medium-quality (\emph{level 2 \&
3} ) genes and excludes some low-quality (\emph{level 4} ) genes. This
trend is particularly evident when the variance in reference sequences
is greater than 11\%. However, if users select a bootstrap score that is
too strict, for example, removing all sequences with a bootstrap score
that is not equal to 100, many high-quality sequences will be lost. Of
note, when the variance in reference sequences is less than or equal to
15\%, removing genes with a bootstrap score lower than 90 hardly
sacrifices high-quality results and can effectively eliminate
low-quality results that are caused by highly variable reference
sequences. We visually compared the results of removing those with a
bootstrap score lower than 90 with the untreated results (Figure 3-C).
These results reveal that when the variance in reference sequences is
above 11\%, filtering through the bootstrap score can effectively reduce
base variations caused by incorrect assembly.

\textbf{Discussion}

\textbf{Pre-assembly filtering}

One key factor that makes GeneMiner faster in comparison to other
software is its use of k-mer-based reads filtering (see Methods
section), as other tools employ sequence alignment-based methods, such
as BLAST, for read filtering. Although the k-mer based filtering method
is quick, it has several drawbacks. Firstly, when the value of k is
small, different sequences may generate the same k-mer, leading to a
higher false-positive rate. Simultaneously, the k-mer method may
encounter difficulties when dealing with genes containing a large amount
of complex structural variation and repetitive regions. To address the
consequence of a small k, we designed the re-filtering method to
iteratively (automatically) filter the labeled filtered datasets using
the gradient rising \emph{k\textsubscript{f}} (the length of the k-mer
in the reads filtering step). This re-filtering method prevents the
slowdown of the subsequent assembly caused by excessively large,
filtered datasets. Additionally, a larger \emph{k\textsubscript{f}} is
recommended for repeat regions. Regarding the second issue, we suggest
that GeneMiner should only be used for genes that do not have complex
structural variation and repetitive regions. Typically, genes used as
molecular phylogenetic markers have simpler structures and do not
contain repetitive regions.

\textbf{Assembly in GeneMiner}

Genes used for phylogenetic studies are typically shorter than the
average gene, allowing us to employ more computational processes during
the assembly to enhance the accuracy of the assembly. Compared to our
previously developed program Easy353, GeneMiner has implemented the
enhancements and improvements outlined below.

\textbf{(1) Seed selection}

As the initial step of sequence extension, seed selection is crucial in
determining the accuracy and completeness of an assembly. To simplify
this process, GeneMiner quickly scans through sequencing data to
automatically identify and select appropriate seeds without any manual
work. This saves and effort while ensuring optimal assembly. We assume
that k-mers, which occur at high frequencies in both the reference and
sequencing data, are conserved regions with a higher probability of
occurrence in the target genes. To achieve an unbiased selection of
candidate seeds, we apply a weighted seed model that accounts for the
k-mer counts in both the reference and sequencing data. This model
assigns a weighted score to each candidate seed, with a stronger
preference toward conserved regions in the reference, to avoid the risk
of high-frequency false positives or repeat regions in the sequencing
data. If the resulting assembly from a particular seed candidate is
unsatisfactory according to the assembly length and completeness,
GeneMiner will select a new candidate seed to optimize assembly
performance.

(2) \textbf{Weighted node model}

De Bruijn graph is the foundational methodology for almost all
short-read genome and transcriptome assembly tools (Bao et al., 2014;
Cameron et al., 2017; Chang et al., 2015; Li et al., 2017; Pandey et
al., 2017). In the field of genomics, each node in a de Bruijn graph
stands for a k-mer. These nodes are connected through directed edges
when their (k-1) long suffixes match another node's (k-1) long prefix.
The k-mers are most often derived from unassembled DNA sequencing reads.
The key concepts of GeneMiner include utilizing de Bruijn graphs to
establish connections between k-mers. We employ a weighted node model
that combines information from both reference sequences and the input
reads to guide seed selection and node connection and use depth-first
search and stacks to enable efficient seed greedy extension and
backtracking. The weighted node model encompasses both seed selection
and node-to-node connection choices, with a distinct emphasis on
assigning weighted scores. By taking into account both the reads and
reference sequences information, the model importantly balances the
impact of sequencing errors and reference bias on assembly. As a result,
the weighted node model helps make optimal decisions for large numbers
of non-unique node-to-node connections or seed selections.

\textbf{(3) Re-assembly and soft boundary}

During the subsequent assembly or re-assembly processes, the accuracy of
the bases gradually diminishes as the extension length along both sides
of the recovered target gene increases. However, this decline is not
sudden but rather a gradual decrease within a certain buffer length. To
address this, we introduced a compensation mechanism called \emph{soft
boundary. Soft boundary} is a vague boundary between the regions with
low and high read coverage; this buffer ensures that the accuracy of
phylogenetic information remains relatively stable with only subtle
changes in its edge. To support maximum performance of the program, we
suggest setting the \emph{soft boundary} to half the read length.

\textbf{Parameter-bootstrap verification}

Being able to accurately determine phylogeny depends on the accurate and
efficient acquisition of target genes. At this time, selecting a
suitable reference from closely related taxa for non-model species still
remains challenging for researchers. To address this issue, we developed
a parameter-bootstrap verification solution in GeneMiner. This method
generates a new set of simulated reference data by employing a mutation
model and repeated re-sampling. GeneMiner can then include the newly set
of simulated reference data as input to assemble the target gene again.
We observe target genes with \emph{bootstrap scores} below 90 tend to
have an unstable assembly, and therefore may not be reliable reference
choices. While target genes with bootstrap scores above 90 can also
exhibit an unstable assembly, their overall indel and substitution rates
are comparatively lower. Our statistical analysis led us to select
target genes with bootstrap scores of 90 or higher, which allowed us to
mitigate the potential impact of unstable assembly. This indicates that
the parameter-bootstrap verification method is effective in evaluating
assembly results and can guide reference sequence selection. GeneMiner
offers improved sensitivity to false substitutions compared to
insertions and deletions, which excludes accounting for insertions and
deletions in reference data. Users should exercise caution when
inserting and removing large portions of the assembly results, although
these large insertions and deletions have minimal impact on the
phylogenetic results after being corrected. As an added measure to
ensure accuracy when using GeneMiner, we suggest that users employ a
software such as trimAL (Capella-Gutierrez et al., 2009) to edit the
sorting file before constructing a phylogenetic tree.

\textbf{The impact of sequencing depth}

Although GeneMiner can extract many phylogenetic markers from different
types of NGS data, we cannot guarantee 100\% accuracy in identifying
target genes and eliminating incorrect assembly. Obtaining the best
results from GeneMiner is dependent on sufficient sequencing depth of
the input data and closely related reference sequences. In Test II, when
the depth is less than 20x, the number of high-quality assembly results
significantly decreases. In practice, if users are utilizing
transcriptome data, although there is no consistent depth on account of
significant expression differences among target genes, the coverage is
often much more than 20x due to the relatively small size of
transcripts. Thus, utilizing transcriptome data for gene mining should
generally yield satisfactory results.

If users are employing shallow genome sequencing data, the minimum depth
that GeneMiner needs will depend on the type of genes to be mined. For
single-to-low-copy nuclear genes, at least 10x depth is needed, while
coverage of 20x or more is recommended. For mining chloroplast,
mitochondrial, or high-copy-number ribosomal genes (e.g., commonly used
ITS genes), the minimum sequencing depth can be as low as 1x. This is
because organelles like chloroplasts and mitochondria are far more
abundant in cells than nuclei, providing higher coverage during
sequencing. However, we recommend using deeper sequencing data whenever
it is possible. Furthermore, although GeneMiner can process long-read
sequencing data without encountering errors, we do not recommend its
utilization, as the pipeline is unable to harness the full benefits of
longer reads.

\textbf{The impact of reference sequences}

In Test II, when the variance in the reference sequences 18\%, the
number of high-quality assembly results significantly decreases. In our
test datasets, each gene included 50 reference sequences with the same
mutation rate, which mildly increased the software's tolerance to
variance. In practice, one does not know the actual mutation rate
between the reference sequence and the target sequence in advance,
although it can be estimated through the pairwise comparison results for
the reference sequence. For the best performance of GeneMiner, we
suggest that the average mutation rate in the pairwise comparison of
reference sequences should not exceed 15\%.

Most often, users may have only very limited reference sequences
available to them . If this is the case, we suggest that users use all
the genes of species from the same family as reference sequences,
without giving special consideration to the mutation rate of those
sequences. Particularly, if there is only a single or few very close
reference sequences available, the tolerance of GeneMiner to reference
sequences will significantly decrease. This is because the k-mer based
pre-filtering method is less effective at providing correct reads based
on limited reference sequences. Additionally, the assembly process might
fail to acquire the correct seed and allocate enough weight to the
correct results. In such cases, if providing more reference sequences is
not possible, users should try to reduce the values
of\emph{k\textsubscript{f}} and \emph{k\textsubscript{a}} in the program
(corresponding to ''k1'' and ''k2'', respectively, in GeneMiner). This
adjustment should result in increased output and longer contigs, but
this might also increase the assembly errors in the results.

In addition, paralogous genes can also impact the results. This is
because the paralogs can lead to significant changes in the branch
length and even the topology of the phylogenetic tree. The most
prevalent method of identification of orthologs predominantly relies on
sequence similarity and sequence length, such as reciprocal best hits
(RBH), where orthologs are assumed if two genes from different genomes
find each other as the best hit in the other genome. For GeneMiner, we
use reference-guided assembly and constraint of sequence length to
obtain putative orthologs. If there are multiple very similar paralogous
genes in both the reference sequence and the target species, the target
gene may be incorrectly assembled from different paralogous genes.
Therefore, we recommend that researchers should clean up the paralogous
genes in the reference sequence in advance when providing it, to obtain
more accurate results. Similarly, GeneMiner is not designed or robust
enough to handle polyploids. Researchers should aim to use sequencing
data from diploid species whenever possible. Additionally, GeneMiner can
export contigs with degenerate bases, which can be further improved
using our PPD script (Zhou et al., 2022). This feature allows us to
accurately identify putative paralogs that may be overlooked by other
tools. PPD identifies paralogs by identifying shared heterozygosity at a
locus between individuals and heterozygosity at a locus within
individuals.

\textbf{Computational Performance}

Finally, we compared the average performance of four tools in terms of
speed, memory efficiency, disk space usage, gene recovery, and accuracy
(Figure 3-A). GeneMiner demonstrated a significant speed advantage in
Test I, outperforming Easy353 and HybPiper, while aTRAM lagged due to
its reliance on BLAST tasks. GeneMiner also stood out for its efficient
memory usage, which users can easily manage by selecting different
subsets of reference sequences. The principal reason for memory
consumption when using GeneMiner is construction of the reference hash
table. The memory usage of this program depends on the variety of genes
and the similarity among the same genes from different species. When
there is high similarity between reference sequences of the same gene,
variations in the number of reference genes will not significantly
elevate memory usage. In Test I, Poaceae and Brassicaceae comprised
224,049 sequences from 347 species and 53,379 sequences from 349
species, respectively. Despite a higher quantity of reference genes in
Poaceae compared to Brassicaceae, the peak memory usage for Poaceae was
lower. This is mainly attributed to the lower variability of the 353
genes within the Poaceae family compared to Brassicaceae. In practice,
if the reference sequences only come from a closely related genus,
GeneMiner only requires approximately 0.33 GB of memory usage per
sample. Regarding disk space usage, GeneMiner was the second-best
performer in Test I, utilizing an average of 0.30 GB, slightly higher
than Easy353 (0.225 GB). In conclusion, GeneMiner's speed, memory, and
disk usage are highly conducive for small servers and even personal
computers, facilitating the ease of gene retrieval tasks. GeneMiner's
user-friendly computational attributes make it a desirable tool for
various users working with gene acquisition tasks.

\textbf{Conclusion}

In summary, GeneMiner is a convenient and efficient tool that utilizes a
reference-guided de Bruijn graph with seed self-discovery and a greedy
extension for phylogenetic marker extraction from non-model organisms
using NGS data. GeneMiner offers a highly optimized approach for
validating assembly results through a range of optimization techniques.
With outstanding speed and optimal resource utilization, GeneMiner
supports multi-core parallel computing and allows for limited memory
consumption. It can easily be used on both high-performance computing
clusters and standard personal computers. GeneMiner is open-source, easy
to install, user-friendly, and compatible with all widely used operating
systems (i.e., Linux, Windows, and macOS). Overall, this software fills
a need in the field for easy-to-use tools for mining NGS data from
non-model organisms.

\textbf{Data Accessibility}

The two public datasets for \emph{A. thaliana} (SRA: SRR18391637)
and\emph{O. sativa} (SRA: SRR9663069) are available for download from
the European Nucleotide Archive (http://www.ebi.ac.uk). All codes,
reference sequences, and test data involved in Test I and Test II can be
found at the following site: https://github.com/yyscu/GeneMiner-Test.
All source codes of GeneMiner are publicly available on GitHub at
https://github.com/yyscu/GeneMiner.

\textbf{Acknowledgements}

This work was supported by the National Natural Science Foundation of
China (Grant No. 32071666, 32271552) and the Science \& Technology
Fundamental Resources Investigation Program (Grant No. 2022FY101000).

\textbf{Figure 1.} A. The interface of GeneMiner-GUI. B. Performance
comparison between GeneMiner and aTRAM, Easy353, and HybPiper. The bar
graph represents the average values for the speed, memory usage, and
disk space usage used by these four tools when analyzing two
species,\emph{A. thaliana} and \emph{O. sativa} .

\textbf{Figure 2.} The workflow of GeneMiner. A. Filtering reads. B.
Assembly. C. Verification.

\textbf{Figure 3.} A. A comparison of the results obtained by aTRAM,
Easy353, GeneMiner, and HybPiper for Angiosperms353 genes in \emph{A.
thaliana} (left, light-colored bars) and \emph{O. sativa} (right,
dark-colored bars). B. Comparison of mining results for Angiosperms353
genes in simulated sequencing data at various depths with (right,
dark-colored bars) and without (Depth only, left, light-colored bars)
the application of the Weighted-node model. C. Comparative outcomes of
mining Angiosperms353 genes at 50x coverage using simulated sequencing
data with different degrees of variability, with (right, dark-colored
bars) and without (left, light-colored bars) parameter-bootstrap
verification. D. Heatmap depicting the quantity of high-quality data
(\emph{level 1} ) obtained by GeneMiner using simulated sequencing data
and reference sequences with different degrees of variability at diverse
gradient depths.

\textbf{Figure S1.} The heatmaps for levels 2 through 4, indel rate,
substitution rate, and recovered genes.

\textbf{References}

Allen, J. M., Boyd, B., Nguyen, N., Vachaspati, P., Warnow, T., Huang,
D. I., Grady, P. G. S., Bell, K. C., Cronk, Q. C. B., Mugisha, L.,
Pittendrigh, B. R., Soledad Leonardi, M., Reed, D. L., \& Johnson, K. P.
(2017). Phylogenomics from Whole Genome Sequences Using
aTRAM.\emph{Systematic Biology} , syw105.
https://doi.org/10.1093/sysbio/syw105

Allen, J. M., Huang, D. I., Cronk, Q. C., \& Johnson, K. P. (2015).
aTRAM - automated target restricted assembly method: A fast method for
assembling loci across divergent taxa from next-generation sequencing
data. \emph{BMC Bioinformatics} , \emph{16} (1), 98.
https://doi.org/10.1186/s12859-015-0515-2

Allen, J. M., LaFrance, R., Folk, R. A., Johnson, K. P., \& Guralnick,
R. P. (2018). aTRAM 2.0: An Improved, Flexible Locus Assembler for NGS
Data. \emph{Evolutionary Bioinformatics} , \emph{14} , 117693431877454.
https://doi.org/10.1177/1176934318774546

Bao, E., Jiang, T., \& Girke, T. (2014). AlignGraph: Algorithm for
secondary de novo genome assembly guided by closely related
references.\emph{Bioinformatics} , \emph{30} (12), i319--i328.
https://doi.org/10.1093/bioinformatics/btu291

Cameron, D. L., Schröder, J., Penington, J. S., Do, H., Molania, R.,
Dobrovic, A., Speed, T. P., \& Papenfuss, A. T. (2017). GRIDSS:
Sensitive and specific genomic rearrangement detection using positional
de Bruijn graph assembly. \emph{Genome Research} , \emph{27} (12),
2050--2060. https://doi.org/10.1101/gr.222109.117

Capella-Gutierrez, S., Silla-Martinez, J. M., \& Gabaldon, T. (2009).
trimAl: A tool for automated alignment trimming in large-scale
phylogenetic analyses. \emph{Bioinformatics} , \emph{25} (15),
1972--1973. https://doi.org/10.1093/bioinformatics/btp348

Chang, Z., Li, G., Liu, J., Zhang, Y., Ashby, C., Liu, D., Cramer, C.
L., \& Huang, X. (2015). Bridger: A new framework for de novo
transcriptome assembly using RNA-seq data. \emph{Genome Biology}
,\emph{16} (1), 30. https://doi.org/10.1186/s13059-015-0596-2

Chikhi, R., \& Rizk, G. (2013). \emph{Space-efficient and exact de
Bruijn graph representation based on a Bloom filter} . 9.
https://doi.org/10.1186/1748-7188-8-22

Cock, P. J. A., Antao, T., Chang, J. T., Chapman, B. A., Cox, C. J.,
Dalke, A., Friedberg, I., Hamelryck, T., Kauff, F., Wilczynski, B., \&
de Hoon, M. J. L. (2009). Biopython: Freely available Python tools for
computational molecular biology and bioinformatics.\emph{Bioinformatics}
, \emph{25} (11), 1422--1423.
https://doi.org/10.1093/bioinformatics/btp163

Fér, T., \& Schmickl, R. E. (2018). HybPhyloMaker: Target Enrichment
Data Analysis From Raw Reads to Species Trees. \emph{Evolutionary
Bioinformatics} , \emph{14} , 1176934317742613.
https://doi.org/10.1177/1176934317742613

Goodstein, D. M., Shu, S., Howson, R., Neupane, R., Hayes, R. D., Fazo,
J., \ldots{} \& Rokhsar, D. S. (2012). Phytozome: a comparative platform
for green plant genomics. \emph{Nucleic acids research} , 40(D1),
D1178-D1186.

Johnson, M. G., Gardner, E. M., Liu, Y., Medina, R., Goffinet, B., Shaw,
A. J., Zerega, N. J. C., \& Wickett, N. J. (2016). HybPiper: Extracting
coding sequence and introns for phylogenetics from high-throughput
sequencing reads using target enrichment. \emph{Applications in Plant
Sciences} , \emph{4} (7), 1600016. https://doi.org/10.3732/apps.1600016

Johnson, M. G., Pokorny, L., Dodsworth, S., Botigué, L. R., Cowan, R.
S., Devault, A., Eiserhardt, W. L., Epitawalage, N., Forest, F., Kim, J.
T., Leebens-Mack, J. H., Leitch, I. J., Maurin, O., Soltis, D. E.,
Soltis, P. S., Wong, G. K., Baker, W. J., \& Wickett, N. J. (2019). A
Universal Probe Set for Targeted Sequencing of 353 Nuclear Genes from
Any Flowering Plant Designed Using k-Medoids Clustering.\emph{Systematic
Biology} , \emph{68} (4), 594--606.
https://doi.org/10.1093/sysbio/syy086

Li, D., Huang, Y., Leung, C.-M., Luo, R., Ting, H.-F., \& Lam, T.-W.
(2017). MegaGTA: A sensitive and accurate metagenomic gene-targeted
assembler using iterative de Bruijn graphs. \emph{BMC Bioinformatics}
,\emph{18} (S12), 408. https://doi.org/10.1186/s12859-017-1825-3

Li, M., Wunder, J., Bissoli, G., Scarponi, E., Gazzani, S., Barbaro, E.,
Saedler, H., \& Varotto, C. (2008). Development of COS genes as
universally amplifiable markers for phylogenetic reconstructions of
closely related plant species. \emph{Cladistics} , \emph{24} (5),
727--745. https://doi.org/10.1111/j.1096-0031.2008.00207.x

Li, Z., De La Torre, A. R., Sterck, L., Cánovas, F. M., Avila, C.,
Merino, I., Cabezas, J. A., Cervera, M. T., Ingvarsson, P. K., \& Van de
Peer, Y. (2017). Single-Copy Genes as Molecular Markers for Phylogenomic
Studies in Seed Plants. \emph{Genome Biology and Evolution} ,
\emph{9}(5), 1130--1147. https://doi.org/10.1093/gbe/evx070

Liu, B., Ma, Z., Ren, C., Hodel, R. G. J., Sun, M., Liu, X., Liu, G.,
Hong, D., Zimmer, E. A., \& Wen, J. (2021). Capturing single\selectlanguage{english}-copy
nuclear genes, organellar genomes, and nuclear ribosomal DNA from deep
genome skimming data for plant phylogenetics: A case study in
Vitaceae.\emph{Journal of Systematics and Evolution} , \emph{59} (5),
1124--1138. https://doi.org/10.1111/jse.12806

Palmer, J. D., Jorgensen, R. A., \& Thompson, W. F. (1985). CHLOROPLAST
DNA VARIATION AND EVOLUTION IN PISUM: PATTERNS OF CHANGE AND
PHYLOGENETIC ANALYSIS. \emph{Genetics} , \emph{109} (1), 195--213.
https://doi.org/10.1093/genetics/109.1.195

Palmer, J. D., \& Thompson, W. F. (1982). Chloroplast DNA rearrangements
are more frequent when a large inverted repeat sequence is
lost.\emph{Cell} , \emph{29} (2), 537--550.
https://doi.org/10.1016/0092-8674 (82)90170-2

Palmer, J. D., \& Zamir, D. (1982). Chloroplast DNA evolution and
phylogenetic relationships in \emph{Lycopersicon} . \emph{Proceedings of
the National Academy of Sciences} , \emph{79} (16), 5006--5010.
https://doi.org/10.1073/pnas.79.16.5006

Pandey, P., Bender, M. A., Johnson, R., \& Patro, R. (2017). deBGR: An
efficient and near-exact representation of the weighted de Bruijn
graph.\emph{Bioinformatics} , \emph{33} (14), i133--i141.
https://doi.org/10.1093/bioinformatics/btx261

Pevzner, P. A., Tang, H., \& Waterman, M. S. (2001). An Eulerian path
approach to DNA fragment assembly. \emph{Proceedings of the National
Academy of Sciences} , \emph{98} (17), 9748--9753.
https://doi.org/10.1073/pnas.171285098

Henriksen, R. A., Zhao L., Korneliussen T.S. (2023). NGSNGS:
next-generation simulator for next-generation sequencing
data.\emph{Bioinformatics} , 39(1), btad041.

Schulz, M. H., Zerbino, D. R., Vingron, M., \& Birney, E. (2012). Oases:
Robust de novo RNA-seq assembly across the dynamic range of expression
levels. \emph{Bioinformatics} , \emph{28} (8), 1086--1092.
https://doi.org/10.1093/bioinformatics/bts094

Simpson, J. T., Wong, K., Jackman, S. D., Schein, J. E., Jones, S. J.
M., \& Birol, \selectlanguage{polish}İ. \selectlanguage{english}(2009). ABySS: A parallel assembler for short read
sequence data. \emph{Genome Research} , \emph{19} (6), 1117--1123.
https://doi.org/10.1101/gr.089532.108

Small, R. L., Cronn, R. C., \& Wendel, J. F. (2004). Use of nuclear
genes for phylogeny reconstruction in plants. \emph{Australian
Systematic Botany} , \emph{17} (2), 145. https://doi.org/10.1071/SB03015

Sohn, J., \& Nam, J.-W. (2016). The present and future of \emph{de
novo}whole-genome assembly. \emph{Briefings in Bioinformatics} , bbw096.
https://doi.org/10.1093/bib/bbw096

Straub, S. C. K., Parks, M., Weitemier, K., Fishbein, M., Cronn, R. C.,
\& Liston, A. (2012). Navigating the tip of the genomic iceberg:
Next-generation sequencing for plant systematics. \emph{American Journal
of Botany} , \emph{99} (2), 349--364.
https://doi.org/10.3732/ajb.1100335

Wang, X., Liang, D., Jin, W., Tang, M., Shalayiwu, Liu, S., \& Zhang, P.
(2020). Out of Tibet: Genomic Perspectives on the Evolutionary History
of Extant Pikas. \emph{Molecular Biology and Evolution} , \emph{37} (6),
1577--1592. https://doi.org/10.1093/molbev/msaa026

Wang, Z., Gerstein, M., \& Snyder, M. (2009). RNA-Seq: A revolutionary
tool for transcriptomics. \emph{Nature Reviews Genetics} , \emph{10}
(1), 57--63. https://doi.org/10.1038/nrg2484

Weitemier, K., Straub, S. C. K., Cronn, R. C., Fishbein, M., Schmickl,
R., McDonnell, A., \& Liston, A. (2014). Hyb-Seq: Combining Target
Enrichment and Genome Skimming for Plant
Phylogenomics.\emph{Applications in Plant Sciences} , \emph{2} (9),
1400042. https://doi.org/10.3732/apps.1400042

Wen, J., Xie, D.-F., Price, M., Ren, T., Deng, Y.-Q., Gui, L.-J., Guo,
X.-L., \& He, X.-J. (2021). Backbone phylogeny and evolution of
Apioideae (Apiaceae): New insights from phylogenomic analyses of
plastome data. \emph{Molecular Phylogenetics and Evolution} , \emph{161}
, 107183. https://doi.org/10.1016/j.ympev.2021.107183

Wen, J., Yu, Y., Xie, D.-F., Peng, C., Liu, Q., Zhou, S.-D., \& He,
X.-J. (2020). A transcriptome-based study on the phylogeny and evolution
of the taxonomically controversial subfamily Apioideae
(Apiaceae).\emph{Annals of Botany} , \emph{125} (6), 937--953.
https://doi.org/10.1093/aob/mcaa011

Wu, F., Mueller, L. A., Crouzillat, D., P\selectlanguage{ngerman}étiard, V., \& Tanksley, S. D.
(2006). Combining Bioinformatics and Phylogenetics to Identify Large
Sets of Single-Copy Orthologous Genes (COSII) for Comparative,
Evolutionary and Systematic Studies: A Test Case in the Euasterid Plant
Clade. \emph{Genetics} , \emph{174} (3), 1407--1420.
https://doi.org/10.1534/genetics.106.062455

Zerbino, D. R., \& Birney, E. (2008). Velvet: Algorithms for de novo
short read assembly using de Bruijn graphs. \emph{Genome Research}
,\emph{18} (5), 821--829. https://doi.org/10.1101/gr.074492.107

Zhang, F., Ding, Y., Zhu, C.-D., Zhou, X., Orr, M. C., Scheu, S., \&
Luan, Y.-X. (2019). Phylogenomics from low-coverage whole-genome
sequencing. \emph{Methods in Ecology and Evolution} , \emph{10} (4),
507--517. https://doi.org/10.1111/2041-210X.13145

Zhang, Z., Xie, P., Guo, Y., Zhou, W., Liu, E., \& Yu, Y. (2022).
Easy353: A Tool to Get Angiosperms353 Genes for Phylogenomic
Research.\emph{Molecular Biology and Evolution} , \emph{39} (12),
msac261. https://doi.org/10.1093/molbev/msac261

Zhou, W., Soghigian, J., \& Xiang, Q.-Y. (Jenny). (2022). A New Pipeline
for Removing Paralogs in Target Enrichment Data. \emph{Systematic
Biology} , \emph{71} (2), 410--425.
https://doi.org/10.1093/sysbio/syab044

\textbf{Hosted file}

\verb`Table 1.xlsx`  available at \url{https://authorea.com/users/607221/articles/683413-geneminer-a-tool-for-extracting-phylogenetic-markers-from-next-generation-sequencing-data}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-1/Figure-1}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-2/Figure-2}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Figure-3/Figure-3}

\end{center}
\end{figure}

\selectlanguage{english}
\FloatBarrier
\end{document}

