Mustafa Abdallah

Ashraf Mahgoub

Saurabh Bagchi

Somali Chaterji

Athena: Automated Tuning of Genomic Error Correction Algorithms using Language Models

The performance of most error-correction algorithms that operate on genomics reads is dependent on the proper choice of its configuration parameters, such as the value of k in k -mer based techniques. In this work, we target the problem of finding the best values of these configuration parameters to optimize error correction and indirectly improve genome assembly. We perform this in an adaptive manner, adapted to different datasets and to EC tools, due to the observation that different configuration parameters are optimal for different datasets, i.e., from different platforms and species and vary with the EC algorithm being applied. We use language modeling techniques from the Natural Language Processing (NLP) domain in our algorithmic suite, Athena, to automatically tune the performance-sensitive configuration parameters. Through the use of N-Gram and Recurrent Neural Network (RNN) language modeling, we validate the intuition that the Error Correction (EC) performance can be computed quantitatively and efficiently using the "perplexity" metric, repurposed from NLP. After training the language model, we show that the perplexity metric calculated for runtime data has a strong negative correlation with the correction of the erroneous NGS reads. Therefore, we use the perplexity metric to guide a hill climbing-based search, converging toward the best k-value. Our approach is suitable for both de novo and comparative sequencing (resequencing), eliminating the need for a reference genome to serve as the ground truth. True Calculate Perplexity metric (P) on 𝑆 [𝑘 ], 𝑆 [𝑘 δ] and 𝑆 [𝑘 + δ] EC Tool (e.g., Lighter, Blue,…etc.) Genomic Reads from Sequencer (Uncorrected) 𝐷 Smaller set (S') of genomic reads (Uncorrected)

Introduction

Rapid advances in next-generation sequencing (NGS) technologies, with the resulting drops in sequencing costs, offer unprecedented opportunities to characterize genomes across the tree of life. While NGS techniques allow for rapid parallel sequencing, they are more error-prone than Sanger reads and generate different error profiles, e.g., substitutions, insertions, and deletions. The genome analysis workflow needs to be carefully orchestrated so errors in reads are not magnified downstream. Consequently, multiple error-correction (EC) techniques have been developed for improved performance for applications ranging from de novo variant calling to differential expression, iterative k-mer selection for improved genome assembly [START_REF] Mahadik | Scalable genomic assembly through parallel de bruijn graph construction for multiple k-mers[END_REF] and for long-read error correction [START_REF] Szalay | De novo sequencing and variant calling with nanopores using poreseq[END_REF][START_REF] Sameith | Iterative error correction of long sequencing reads maximizes accuracy and improves contig assembly[END_REF].

Importantly, the values of the performance-sensitive configuration parameters are dependent not only on the dataset but also on the specific EC tool (Table 2, Table 5 in Appendix). The performance of many EC algorithms is highly dependent on the proper choice of configuration parameters, e.g., k-value (length of the substring) in k -spectrum-based techniques. Selecting different Figure 1: Overview of Athena's workflow. First, we train the language model using the entire set of uncorrected reads for the specific dataset. Second, we perform error correction on a subsample from the uncorrected reads using an EC tool (e.g., Lighter or Blue) and a range of k-values. Third, we compute perplexity of each corrected sample, corrected with a specific EC tool, and decide on the best k-value for the next iteration, i.e., the one corresponding to the lowest perplexity metric because EC quality is anti-correlated with the perplexity metric. This process continues until the termination criteria are met. Finally, the complete set of reads is corrected with the best k-value found and then used for evaluation.

k-values has a trade-off such that small values increase the overlap probability between reads, however, an unsuitably small value degrades EC performance because it does not allow the algorithm to discern correct k-mers from erroneous ones. In contrast, unsuitably high k-values decrease the overlap probability and hurts EC performance because most k-mers will now appear unique. The k-mers that appear above a certain threshold frequency, and are therefore expected to be legitimate, are solid k-mers, the others are called insolid or untrusted k-mers. In k-spectrum-based methods, the goal is to convert insolid k-mers to solid ones with a minimum number of edit operations. Thus, an adaptive method for finding the best k-value and other parameters is needed for improved EC performance, and in turn, genome assembly.

Many existing EC solutions (e.g., Reptile [START_REF] Yang | Reptile: representative tiling for short read error correction[END_REF], Quake [START_REF] Kelley | Quake: quality-aware detection and correction of sequencing errors[END_REF], Lighter [START_REF] Song | Lighter: fast and memory-efficient sequencing error correction without counting[END_REF], Blue [START_REF] Greenfield | Blue: correcting sequencing errors using consensus and context[END_REF]) require users to specify k. The best value is usually found by exploration over k-values [START_REF] Kao | Echo: a reference-free short-read error correction algorithm[END_REF] and evaluating performance metrics, e.g., EC Gain, Alignment Rate, Accuracy, Recall, and Precision. However, a reference genome is typically needed to serve as ground truth for such evaluations, making this approach infeasible for de novo sequencing tasks.To the best of our knowledge, existing tools leave the best parameter choice to the end user [START_REF] Peng | Idba-a practical iterative de bruijn graph de novo assembler[END_REF][START_REF] Mahadik | Scalable genomic assembly through parallel de bruijn graph construction for multiple k-mers[END_REF], although Kmergenie provides intuitive abundance histograms or heuristics to guide k-value selection, albeit for genome assembly. Further, Kmergenie does not provide the optimal k-value for different tools, only accounting for the dataset for computing the abundance histograms. We find that this is an incorrect simplification, e.g., the optimal k-values for Blue and Lighter in our evaluation vary (Tables 2,5, and 6). Within assembly also, it is well known that different k-values are optimal for different assembly pipelines [START_REF] Limasset | Toward perfect reads: self-correction of short reads via mapping on de bruijn graphs[END_REF]. Thus, the need for a data-driven method to select the optimal k-value is crucial. For example, we found that k=25 gives the best EC gain performance for D5 using Lighter. However, if the same k-value was used for D3, the EC gain is only 65% compared to the gain of 95.34% when using the optimal k-value (Table 5 in Appendix). The read on the left is an erroneous read selected from dataset #3 (D3), while the read on the right is the same read, after correction with Lighter. When using language modeling to compute the perplexity for both reads, we notice that the read on the right has a lower perplexity value (15.2), relative to the erroneous read (77.72), as the sequence of k-mers after correction has a higher probability of occurrence. Also notice that the probability of a sequence of k-mers depends on both their frequencies and their relative order in the read, which allows the perplexity metric to capture how likely it is to observe this k-mer given the neighboring k-mers in a given read.

Our solution, Athena1 finds the best value of the configuration parameters for correcting errors in genome sequencing, such as the value of k in k-mer based methods. Further, Athena does not require access to a reference genome to perform its function of determining the optimal parameter configuration2 . Athena, like other EC tools, leverages the fact that NGS reads have the property of reasonably high coverage, 30X-150X coverage depth is commonplace. From this, it follows that the likelihood of correct overlaps for a given portion of the genome will outnumber the likelihood of erroneous ones [START_REF] Yang | Reptile: representative tiling for short read error correction[END_REF]. Athena captures the underlying semantics of NGS reads from different datasets such that the language model (LM) can learn the complexities of a language and can also predict future sequences. This is integral to traditional NLP tasks such as speech recognition, machine translation, or text summarization. LM's goal is to learn a probability distribution over sequences of symbols and words, in this case specialized to the dataset and to the EC tool of choice for optimal performance. Therefore, it can be used to identify sequences that do not fit the semantics of the language and thus contribute to a high value of the "Perplexity metric". This personalization is conducive to our problem because the LM then learns the patterns pertaining to the specific genome.

In our context, we use LM to estimate the probabilistic likelihood that some observed sequence is solid or insolid, in the context of the specific genome. We train an LM using the original dataset (with uncorrected reads), and then use this trained model to compute the Perplexity metric when the EC tool is run with a specific configuration, i.e., for specific parameters. We show empirically that the perplexity metric is inversely correlated to the quality of error correction. Crucially, the perplexity metric evaluation does not require the computationally expensive alignment to a reference genome. Through a stochastic optimization method, we evaluate the search space to pick the best configuration parameter value to guide EC-k in k-mer-based methods and the Genome Length (GL) in the RACER EC tool. Moreover, most EC tools are evaluated based on their direct ability to reduce error rates rather than to improve genome assembly. Although the assembly can benefit substantially from EC tools, this can also sometime negatively affect assembly quality due to conversion of benign errors into damaging errors [START_REF] Heydari | Evaluation of the impact of illumina error correction tools on de novo genome assembly[END_REF]. In our evaluation, we see that the EC improvement due to Athena also leads to higher quality assembly (Table 3). In summary, this paper makes the following contributions.

1. Athena develops a Language Model-based suite to efficiently determine the best configuration parameter value for any existing Error Correction tool.

2. We introduce a likelihood-based metric, the Perplexity metric repurposed from NLP, to evaluate EC quality that does not need a gold standard genome. This is the first such use of this metric in the bioinformatics domain. 3. We compare and contrast two LM variants of Athena, N-gram and RNN-based. Through this, we show that N-Gram modeling can be faster to train while char-RNN provides similar accuracy to N-gram, albeit with significantly lower memory footprint, conducive to multitenant analysis pipelines e.g., MG-RAST for metagenomics processing [START_REF] Meyer | Mg-rast version 4lessons learned from a decade of low-budget ultra-highthroughput metagenome analysis[END_REF]. 4. We apply Athena to 3 EC tools: Lighter [START_REF] Song | Lighter: fast and memory-efficient sequencing error correction without counting[END_REF], Blue [START_REF] Greenfield | Blue: correcting sequencing errors using consensus and context[END_REF], and RACER [START_REF] Ilie | Racer: Rapid and accurate correction of errors in reads[END_REF] on 5 real datasets, with varied error rates and read lengths. We show that Athena was successful in finding either the best parameters (k for Lighter and Blue, and GenomeLength for RACER) or parameters that perform within 0.27% in overall alignment rate to the best values using exhaustive search against a reference genome. We couple EC, with best parameter selection by Athena, to the Velvet assembler and find that it improves assembly quality (NG50) by a Geometric Mean of 10.3X across the 5 datasets.

2 Our Solution: Athena

Application of Language Models

We use two different LM variants in our Athena suite. We describe them next. N-Gram Language Models: We train an N-Gram model [START_REF] Brown | Class-based n-gram models of natural language[END_REF], which is word-based, from the input set of reads before correction. This N-Gram model needs word-based segmentation of the input read as a pre-processing phase. Then, we use this trained LM to evaluate EC performance.

RNN Language Models:

The second technique is RNN-based LM [START_REF] Kombrink | Recurrent neural network based language modeling in meeting recognition[END_REF], using different RNN architectures, e.g., standard RNNs, LSTMs, and GRUs. These models can be trained either as word-based models or character-based models. We train our RNN variant as character-based models to avoid having to make the decision about how to segment the genomic string, as we have to do for the N-Gram model. More details about theses LM techniques is given in the Appendix, Sec 7.1. Contrasting our 2 LM variants: Although training N-Gram LMs is much faster relative to RNNbased models (3-5 minutes for N-Gram vs. 10-11 hours for RNN-based across the 5 datasets), they still have the requirement of splitting a read into words of specific length. Further, RNN-based models have much lower memory footprint and storage requirements relative to N-Gram. This is because N-Gram models need to store conditional probabilities in large tables with an average size of 0.6-1.1 GB across the 5 datasets. In contrast, RNNs only need to store the network architecture and weights with an average size of 3-5 MB across the 5 datasets. Also, we used an LSTM variant of Athena and found that it took 3 times longer to train and test but gave insignificant improvement in perplexity over RNN.

Intuition for the Use of the Perplexity metric

Our design uses the Perplexity metric (PP) to provide an accurate, and importantly quick, estimation of the EC quality with the current configuration parameter value(s). The Perplexity metric is based on the LM trained on the entire original (uncorrected) dataset. The Perplexity metric is then calculated on the dataset of the corrected reads (entire dataset for N-Gram and 1% for RNN) to measure the EC performance. It measures how well LM can predict the next element in an input stream. Suppose the input stream is H and the next element is e. Then, the Perplexity metric is inversely proportional to the probability of seeing "e" in the stream, given the history H for the learned model. This method works because there is a high negative correlation of the Perplexity metric with both EC metrics-Alignment Rate and EC Gain. Given this anti-correlation, we can rely on the Perplexity metric as an evaluation function, and apply a simple search technique (e.g., hill climbing) to find the best k-value for a given dataset. In this description, for simplicity of exposition, we use the k-value in k-mer based techniques as an example of Athena-tuned configuration parameter. However, Athena can tune any other relevant configuration parameter in EC algorithms and we experimentally show the behavior with another parameter-Genome Length-in the RACER tool. Figure 2 shows an example how Perplexity can evaluate the likelihood of a sequence of k-mers using their frequencies and contextual dependencies. In this example, we notice that the corrected read set (i.e., on the right) has a considerably lower Perplexity value (15.2), relative to the erroneous set (77.72). Thus, our intuition that the Perplexity metric reflects the correctness of the read dataset is validated through the observed negative relationship. For an N-Gram language model, perplexity of a sentence is the inverse probability of the test set, normalized by the number of words. It is clear from (1) that minimizing perplexity is the same as maximizing the probability of the observed set of m words from W 1 to W m . For RNN, perplexity is measured as the exponential of the mean of the cross-entropy loss (CE) as shown in [START_REF] Szalay | De novo sequencing and variant calling with nanopores using poreseq[END_REF], where ŷ is the predicted next character-the output of the RNN-and |V | is the vocabulary size used during training.

P P (W ) = m 1 P (W1, W2, ...., Wm) ≈ m 1 m i=1 P (Wi|Wi-1, ..., Wi-n) (1) CE(y, ŷ) = - |V | i=1 yi log( ŷi). (2)
P erplexity = fn(LM, DC , k) (3) kopt = argmin k i fn(LM, DC , ki) (4)

Search through the Parameter Space

Our objective is to find the best k-value that will minimize the Perplexity of the corrected dataset. We assume that Perplexity can be modeled by the function in [START_REF] Sameith | Iterative error correction of long sequencing reads maximizes accuracy and improves contig assembly[END_REF]. Here, LM: trained language model, D C : corrected genomic dataset, and k: the configuration parameter we wish to tune. f n is a discrete function as k-values are discrete, and therefore, its derivative is not computable. Thus, a gradient-based optimization technique will not work. Hence, we use a simple hill-climbing technique to find the value of k that gives the minimum value of f n , for the given LM and D C (4).

The following pseudo-code describes the steps used for finding the best k-value for a given dataset. We start with Algorithm 1, which invokes Algorithm 2 multiple times, each time with a different starting value. We begin by training an LM on the original uncorrected read set (D 0 ). Second, we assume that the best value of k lies in a range from A to B (initially set to either the tool's recommended range, or between 1 and L, where L is the read size). We apply an existing EC algorithm (Lighter, Blue, or RACER in our evaluation) with different initial values (k 0 ,..., k m ) ∈ (A, B), to avoid getting stuck in local minima, going through multiple iterations for a given initial value. We evaluate the Perplexity for the corrected dataset with current values of k-k i , and its neighbors, (k i -δ and k i + δ). In each iteration, we apply hill-climbing search to identify the next best value of k i for the following iteration. The algorithm terminates whenever the Perplexity relative to k i is less than the perplexities of both its neighbors or the maximum number of (userdefined) iterations is reached. However, all shown results are with respect to only one initial value (i.e., m = 1 in k1, k2

Time and Space Complexity

Because we apply hill climbing search to find the best value of k, the worst-case time complexity of the proposed algorithm is L, the upper bound of the range of k-values to search for. For the space complexity, Athena only needs to save the Perplexity values of previously investigated k-values, which is also linear in terms of L.

Evaluation with Real Datasets

In this section, we evaluate Athena variants separately by correcting errors in 5 real datasets and evaluating the quality of the resultant assembly. We implement the N-Gram model using the SRILM toolkit [START_REF] Stolcke | Srilm -an extensible language modeling toolkit[END_REF]. SRILM is an open source toolkit that supports building statistical N-Gram LMs, in addition to evaluating the likelihood of new sequences using the Perplexity metric. For the RNN LM implementation, we build on the TensorFlow platform [START_REF] Abadi | Tensorflow: A system for large-scale machine learning[END_REF]. Specifically, we utilized Character-level RNN models, char-RNN for short, which takes a chunk of text and models the probability distribution of the next character in the sequence, given a sequence of previous characters. After correction, we run the Bowtie2 aligner [START_REF] Langmead | Fast gapped-read alignment with bowtie 2[END_REF] and measure the Alignment Rate and the Error Correction Gain (EC Gain). A higher value for either metric implies superior error correction. We do a sweep through a range of k-values and measure the alignment rate to determine if the Athena-generated k-value is optimal or its distance from optimality. For interpreting the execution time results, our experiments were performed on Dell Precision T3500 Workstation, with 8 CPU cores, each running at 3.2GHZ, 12GB RAM, and Ubuntu 16.04 Operating System. We use 3 EC tools, in pipeline mode with Athena, namely, Lighter, Blue, and RACER. Blue uses a k-mer consensus to target different kinds of errors, e.g., substitution, deletion and insertion errors, as well as uncalled bases (N). This improves the performance of both alignment and assembly [START_REF] Greenfield | Blue: correcting sequencing errors using consensus and context[END_REF]. In contrast, Lighter is much faster as it uses only a sample of k-mers to perform correction. Third, RACER uses a different configuration parameter distinct from the k-value, specifically Genome Length (GL), and we are able to tune GL as well. Our ability to tune any of these EC algorithm's parameters is in line with our vision and ongoing work to design extensible blocks of software to expedite algorithmic development in bioinformatics [START_REF] Mahadik | Sarvavid: a domain specific language for developing scalable computational genomics applications[END_REF]. Incidentally, we started using another popular EC tool, Reptile, but it only allowed for a smaller range of k-values, beyond which it ran into out-of-memory errors. Hence, to demonstrate results with the full range of k values, we restricted ourselves to Lighter, Blue, and RACER.

Our datasets are Illumina short reads [Table 1], used in multiple prior studies (e.g., [START_REF] Yang | Reptile: representative tiling for short read error correction[END_REF][START_REF] Schrder | Shrec: a short-read error correction method[END_REF]). For these, there exist ground-truth reference genomes, which we use to evaluate the EC quality. The five datasets have different read lengths (from 36bp to 100bp) and different error rates (from < 3% to 43%).

Optimal Parameter Selection

The results of using Athena with Lighter, Blue, and RACER tools are shown in Table 2 for each dataset. The value of k found through exhaustive testing (or, GL for RACER) along with the EC quality is given first. Then it is noted when the Athena-selected configuration matches this theoretical best. We find that in 25 of 30 cases (83.3%) the configuration found by Athena matches the theoretical best. In the remaining cases, our chosen parameters are within 0.27% in overall alignment rate to the best values. We see that the optimal parameter value almost always corresponds to the lowest perplexity scores computed using Athena's language modeling (Table 5, Table 6). Further, the anti-correlation between perplexity and alignment rate holds for both optimal and non-optimal k-values. This shows that our hypothesis is valid across a range of kvalues. We notice that the feasible range of k-values in Blue is [START_REF] Schrder | Shrec: a short-read error correction method[END_REF][START_REF] Van Aken | Automatic database management system tuning through large-scale machine learning[END_REF], distinct from Lighter's. Another interesting observation is that the optimal k-values are different across the two different EC tools, Lighter and Blue, for the same dataset, as observed before [START_REF] Song | Lighter: fast and memory-efficient sequencing error correction without counting[END_REF]). Athena can be applied to a different configuration parameter, GL for the RACER tool, in line with our design as a generalpurpose tuning tool.

N-Gram Language Model Results

We start by training an N-Gram language model from the original dataset. We divide each read into smaller segments (words) of length L s (set to 7 by default). A careful reader may be concerned that selecting the best value of L s is a difficult problem in itself, of the same order of difficulty as our original problem. Fortunately, this is not the case and L s selection turns out to be a much simpler task. From domain knowledge, we know that if we use L s = 4 or less, the frequencies of the different words will be similar as this increases the overlap probability between the generated words, thus reducing the model's discriminatory power, while a large value will mean that the model's memory footprint increases. We find that for a standard desktop-class machine with 32GB of memory, L s = 8 is the maximum that can be accommodated. Further, we find that the model performance is not very sensitive in the range (5-7), so we end up using L s = 7. The same argument holds for selecting a history of words, and we use a tri-gram model (history of 3 words) for all our experiments. Second, we compare the perplexity metric for datasets corrected with different k-values and compare the perplexity metric (without a reference genome) to the alignment rate (using a reference genome). We always report the average perplexity, which is just the total perplexity averaged across all words. Our results show a high negative correlation (≤ -0.946) between the two metrics on the 5 datasets, as shown in Table 3. To reiterate, the benefit of using the perplexity metric is that it can be computed without the ground truth and even where a reference genome is available, it is more efficient than aligning to the reference and then computing the alignment rate.

Char-RNN Language Model Results

For training our RNN, we used the "tensorflow-char-rnn" library [START_REF] Abadi | Tensorflow: A system for large-scale machine learning[END_REF]. After parameter tuning, we use the following for our experiments: 2 hidden layers with 300 neurons per layer, output layer with size 4 (i.e.,, corresponding to the four base pairs), mini-batch size 200, and learning rate 2e -3 respectively. For each of the 5 datasets, we used 90% for training and 10% for validation, with no overlap, for coming up with the optimal RNN architecture.

For our char-RNN results, we find that the perplexity metric has a strong negative relation to the overall alignment rate [Table 3], with the absolute value of the correlation always greater than 0.93. Here we have to sample the uncorrected set for calculating the perplexity measure because using an RNN to perform the calculation is expensive. This is because it involves, for predicting each character, doing 4 feed-forward passes (corresponding to the one-hot encodings for A, T, G, or C), each through 600 neurons. Empirically, for a test sample size of 50K, this translates to approximately 30 minutes on a desktop-class machine. In the experiments with the real datasets, we use 50K samples with uniform sampling, and in synthetic experiments, we use 100K samples (i.e., only 1% of the dataset). More importantly though, the strong quantitative relationship between perplexity and the EC quality is maintained throughout.

Comparison with Self-tuning EC Tool

Here, we compare Athena with the EC tool, Fiona [START_REF] Schulz | Fiona: a parallel and automatic strategy for read error correction[END_REF], which estimates its parameters internally. The purpose of this comparison is to show that Athena can tune k-mer-based approaches (RACER specifically for this experiment) to achieve comparable performance to suffix array-based approaches (e.g., Fiona), reducing the gap between the two approaches. [START_REF] Yang | A survey of error-correction methods for nextgeneration sequencing[END_REF] and [START_REF] Molnar | Correcting illumina data[END_REF] show a similar comparison between different EC approaches concluding that the automatic selection of configuration parameters, based on the datasets, is crucial for EC performance. However, they do not perform such parameter tuning automatically. Table 3 presents the overall alignment rate for our 5 evaluation datasets, calculated after doing correction by Fiona. We notice that RACER, when tuned with Athena, outperforms automatic tuning by Fiona in 3 of the 5 datasets, while they are equal in one dataset. Finally, Fiona is only better on D 5 by 0.05%, which may be deemed insignificant.

Impact on Assembly Quality

Here we show the impact on genome assembly quality of using an EC tool tuned with Athena. We use Velvet [START_REF] Zerbino | Velvet: algorithms for de novo short read assembly using de bruijn graphs[END_REF] to perform the assembly and QUAST [START_REF] Gurevich | quality assessment tool for genome assemblies[END_REF] to evaluate the assembly quality. We compare the NG50 before and after correction done by RACER using the best GL found by Athena. The results (Table 3) show a significant improvement on NG50 by 2.26X, 46X, 4X, 7.27X, and 15X. For D1, the improvement is the lowest, since D1 contains the least fraction of errors to start with. These improvements are consistent with what was reported in [START_REF] Heydari | Evaluation of the impact of illumina error correction tools on de novo genome assembly[END_REF][START_REF] Greenfield | Blue: correcting sequencing errors using consensus and context[END_REF], which measured improvement due to the use of error correction tools with manually tuned configuration parameters. Consider that in our problem statement, we are trying to search through a space of configuration parameters in order to optimize a metric (EC Gain or Alignment Rate). The search space can be large and since the cost of searching shows up as a runtime delay, it is important to reduce the time that it takes to evaluate that metric of each search point. In today's state-of-the-art, to find the best value of a configuration parameter [START_REF] Chikhi | Informed and automated k-mer size selection for genome assembly[END_REF][START_REF] Mahadik | Scalable genomic assembly through parallel de bruijn graph construction for multiple k-mers[END_REF], the method would be to pick a single k, run the EC tool with that value, then perform alignment (with one of several available tools such as Bowtie2), and finally compute the metric for that value. In contrast, with Athena, to explore one point in the search space, we run the EC algorithm with the k-value, and then compute the Perplexity metric, which does not involve the time consuming alignment step. Here, we evaluate the relative time spent in exploring one point in the search space using Athena vis-à-vis the current state-ofthe-art. The result is shown in Table 4. For this comparison, the alignment is done by Bowtie2 [START_REF] Langmead | Fast gapped-read alignment with bowtie 2[END_REF]. We find that using the baseline approach, each step in the search takes respectively 6.2X, 4.8X, and 4.7X, 3.6X, and 5.82X longer for the 5 datasets. Further, while we use the hill-climbing technique to search through the space, today's baseline methods use exhaustive search, such as in Lighter [START_REF] Song | Lighter: fast and memory-efficient sequencing error correction without counting[END_REF] and thus the end-to-end runtime advantage of Athena will be magnified.

Search time improvement with Athena

Related Work

Error Correction Approaches: EC tools can be mainly divided into three categories: kspectrum based, suffix tree/array-based, and multiple sequence alignment-based (MSA) methods. Each tool takes one or more configuration parameters. While we have experimented with Athena applied to the first kind, with some engineering effort, it can be applied to tuning tools that belong to the other two categories. Language Modeling in Genomics: In the genomics domain, LM was used in [START_REF] Ganapathiraju | Comparative n-gram analysis of whole-genome sequences[END_REF] to find the characteristics of organisms in which N-Gram analysis was applied to 44 different bacterial and archaeal genomes and to the human genome. In subsequent work, they used N-Gram-based LM for extracting patterns from whole genome sequences. Others [START_REF] Coin | Enhanced protein domain discovery by using language modeling techniques from speech recognition[END_REF] have used language modeling to enhance domain recognition in protein sequences. For example, [START_REF] King | ngloc: an n-gram-based bayesian method for estimating the subcellular proteomes of eukaryotes[END_REF] has used N-Gram analysis specifically to create a Bayesian classifier to predict the localization of a protein sequence over 10 distinct eukaryotic organisms. RNNs can be thought of as a generalization of Hidden Markov Models (HMMs) and HMMs have been applied in several studies that seek to annotate epigenetic data. For example, [START_REF] Song | Spectacle: fast chromatin state annotation using spectral learning[END_REF] presents a fast method using spectral learning with HMMs for annotating chromatin states in the human genome. Thus, we are seeing a steady rise in the use of ML techniques traditionally used in NLP being used to make sense of genomics data. Automatic Parameter Tuning: The field of computer systems has had several successful solutions for automatic configuration tuning of complex software systems. Our own work [START_REF] Mahgoub | Rafiki: a middleware for parameter tuning of nosql datastores for dynamic metagenomics workloads[END_REF] plus others [START_REF] Van Aken | Automatic database management system tuning through large-scale machine learning[END_REF] have shown how to do this for distributed databases, while other works have done this for distributed computing frameworks like Hadoop [START_REF] Bei | RFHOC: A Random-Forest Approach to Auto-Tuning Hadoop's Configuration[END_REF][START_REF] Li | Mronline: Mapreduce online performance tuning[END_REF] or cloud configurations [START_REF] Alipourfard | Cherrypick: Adaptively unearthing the best cloud configurations for big data analytics[END_REF]. We take inspiration from them but our constraints and requirements are different (such as, avoiding reliance on ground truth corrected sequences).

Discussion

The space to search for finding the optimal configuration is non-convex in general. Therefore, it is possible to get stuck in a local minima, and hence, we use multiple random initializations. However, in our evaluation, we find that a single initialization suffices. Some EC tools have a number of performance-sensitive configuration parameters with interdependencies. There, systems such as Rafiki [START_REF] Mahgoub | Rafiki: a middleware for parameter tuning of nosql datastores for dynamic metagenomics workloads[END_REF] can encode the dependencies, while relying on Athena's LM to compute the corresponding performance metric, converging toward optimal parameters. We can use Athena to optimize the k-value in DBG-based assemblers as well, though there will be different optimization passes since the optimal values are likely to be different for the error correction and assembly stages.

Conclusion

The performance of most EC tools for NGS reads is highly dependent on the proper choice of its configuration parameters, e.g., k-value selection in k -mer based techniques as shown in Table 2. Using our Athena suite, we target the problem of automatically tuning these parameters using language modeling techniques from the NLP domain without the need for a ground truth genome. Through N-Gram and char-RNN language modeling, we validate the intuition that the EC performance can be computed quantitatively using the "perplexity" metric, which then guides a hill climbing-based search toward the best k-value. We evaluate Athena with 5 different real datasets, plus, with synthetically injected errors. We find that the predictive performance of the perplexity metric is maintained under all scenarios, with absolute correlation scores higher than 0.93. Further, using the perplexity metric, Athena can search for and arrive at the best k-value, or within 0.27% of the assembly quality obtained using brute force.

Error Correction and Evaluation

The majority of error correction tools share the following intuition: high-fidelity sequences (or, solid sequences) can be used to correct errors in low-fidelity sequences (or, in-solid sequences). However, they vary significantly in the way they differentiate between solid and in-solid sequences. For example, [START_REF] Yang | Reptile: representative tiling for short read error correction[END_REF] corrects genomic reads containing insolid k-mers using a minimum number of edit operations such that these reads contain only solid k-mers after correction. The evaluation of de novo sequencing techniques rely on likelihood-based metrics such as ALE and CGAL, without relying on the availability of a reference genome. On the other hand, comparative sequencing or resequencing, such as to study structural variations among two genomes, do have reference genomes available.

Language Modeling

To increase the accuracy of detecting words in speech recognition, language modeling techniques have been used to see which word combinations have higher likelihood of occurrence than others, thus improving context-based semantics. Thus, language modeling is being used in many applications such as speech recognition, text retrieval, and many NLP applications. The main task of these statistical models is to capture historical information and predict the future sequences based on that information. N-Gram-Based Language Modeling. This type of modeling is word-based. The main task that N-Gram based models [START_REF] Brown | Class-based n-gram models of natural language[END_REF] have been used for is to estimate the likelihood of observing a word W i , given the set of previous words W 0 , • • • W i-1 , estimated using the following equation:

P (W 0 , W 1 , ..., W m ) = m i=1 P (W i |W i-1 , ..., W 1 ) ≈ m i=1 P (W i |W i-1 , ..., W i-n ) (5)
Where n represents the number of history words the model uses to predict the next word. Obviously, a higher n results in better prediction, at the cost of higher training time resulting from a more complex model.

Char-RNN-Based Language Modeling. Recurrent neural network (RNN) is a very popular class of neural networks for dealing with sequential data, frequently encountered in the NLP domain. The power of RNN is that each neuron or unit can use its internal state memory to save information from the previous input and use that state, together with the current input, to determine what the next output should be. Character-level RNN models, char-RNN for short, operate by taking a chunk of text and modeling the probability distribution of the next character in the sequence, given a sequence of previous characters. This then allows it to generate new text, one character at a time [START_REF] Graves | Generating sequences with recurrent neural networks[END_REF]. As shown in Fig. 3, RNNs consist of three main layers: Input Layer, Hidden Layer, and Output Layer. First, Input Layer takes x t vector, which is input at a time step t, usually a one-hot encoding vector of the t th word or character of the input sentence. Second, Hidden Layer consists of the hidden state at the same time step s t , which represents the memory of this network. It is calculated as a non-linear function f (e.g., tanh ) of the previous hidden state s t-1 and the input at current time step x t with the following relation:

s t = f (U x t + W s t-1 ). ( 6
)
where, W is a matrix that consists of hidden weights of this hidden layer. Finally, Output Layer consists of a vector o t , which represents the output at step t and contains prediction probabilities for the next character in the sentence. Formally, its length equals the size of the vocabulary and is calculated using a softmax function. Backpropagation was used to train the RNN to update weights and minimize the error between the observed and the estimated next word. For Deep RNN architectures, there are multiple parameters that affect the performance of the model. The two main parameters are: Number of Hidden Layers and Number of Neurons per Layer. For Char-RNN language modeling, vocabulary would include the four nucleotide bases as characters A, C, G, and T. Each input is a one-hot encoding vector for the four nucleotides. Each output vector at each time step also has the same dimension. Perplexity of the Language Model. Perplexity is a measurement of how well a probability distribution predicts a sample. In NLP, perplexity is one of the most effective ways of evaluating the goodness of fit of a language model since a language model is a probability distribution over entire sentences of text [START_REF] Azzopardi | Investigating the relationship between language model perplexity and ir precision-recall measures[END_REF]. For example, 5 per word perplexity of a model translates to the model being as confused on test data as if it had to select uniformly and independently from 5 possibilities for each word. Thus, a lower perplexity indicates that language model is better at making predictions.

Detailed Results

In this section, we show a detailed version of the results presented in Table 2.

Evaluation with Synthetically Injected Errors

Here we experiment on data sets where we synthetically inject errors of three kinds -insertion, deletion, and substitution. The real data sets used in our evaluation in Section 3 belonged to Illumina and therefore had primarily substitution errors (about 99% of all errors). Hence our synthetic injections are meant to uncover if the relationship between perplexity and error rate holds for the other error types. We start by randomly collecting 100K short reads from the reference genome (i.e.,, without any error) for each of the two organisms used in the real data sets-E. coli (D1, D2) and Acinetobacter (D3). Afterward, we inject errors of each type as follows: is likely because D3 had a higher natural ambient rate of error and hence the additional injection does not cause as much increase in the perplexity metric.

Figure 2 :

 2 Figure 2: An example showing how the perplexity metric encodes errors in genomic reads.The read on the left is an erroneous read selected from dataset #3 (D3), while the read on the right is the same read, after correction with Lighter. When using language modeling to compute the perplexity for both reads, we notice that the read on the right has a lower perplexity value (15.2), relative to the erroneous read (77.72), as the sequence of k-mers after correction has a higher probability of occurrence. Also notice that the probability of a sequence of k-mers depends on both their frequencies and their relative order in the read, which allows the perplexity metric to capture how likely it is to observe this k-mer given the neighboring k-mers in a given read.

Figure 3 :

 3 Figure 3: Structure of a Recurrent Neural Network consisting of a chain of repeating modules of neural networks. The RNN predicts the next character at each time step (t + 1), knowing the history of characters in previous time steps.

Figure 4 :

 4 Figure 4: N-Gram (Figures A and B) and RNN (Figures C and D) Perplexity metric for different types of synthetic errors: Indels and Substitution errors, and a mixture of the three for E. coli str reference genome (FiguresA and C) and Acinetobacter sp. reference genome (FiguresB and D). We compare two versions of such errors: high and low error rates.

  Take k as the argmin of K i values picked at step 5 and do a complete correction on the entire dataset D 0 and return Dc. Evaluate Perplexity of k i and its neighbors (k i -δ, k i + δ).2-If fn(LM, S , k i ) < fn(LM, S , k i -δ)&&fn(LM, S , k i ) < fn(LM, S , k i + δ), return k i as k * 3-Else, set k i+1 = value of neighbor with lower fn value.

	Algorithm 1 Correct Set of Reads Input: Dataset: D 0 , Read Length: L, Maximum number of iterations: IterM ax Output: Corrected Dataset: Dc 1-Train a Language Model (LM) using D 0 . 2-Select random sample S from D 0 . 3-Pick k 0 ,k 1 ,...,km: initial values of k in the range of (1, L). 4-for i in the range (0, m) 5-Call Find Optimal K i (S , L, LM, IterMax, k i ) 6-Algorithm 2 Find Optimal k 4-if IterMax > 0 Do: 5-decrement IterMax and Call Find Optimal K (S', L, LM, IterMax , k i+1 ) 6-else
	7-return best k i so far

, • • • , km). Input: Data set: S , Read Length: L, Language Model: LM, MaxNumOfIterations: IterMax, Initial value of k: k i Output: Bestk : k * 1-

Table 1 :

 1 Datasets' description with coverage, number of reads, read lengths, and the reference genome.

	Dataset Coverage #Reads Read Length	Genome Type	Reference Genome Accession Number
	D1	80X	20.8M	36 bp	E. coli str. K-12 substr	NC 000913	SRR001665
	D2	71X	7.1M	47 bp	E. coli str. K-12 substr	NC 000913	SRR022918
	D3	173X	18.1M	36 bp	Acinetobacter sp. ADP1	NC 005966	SRR006332
	D4	62X	3.5M	75 bp	B. subtilis	NC 000964.3	DRR000852
	D5	166X	7.1M	100 bp	L. interrogans C sp. ADP1	NC 005823.1	SRR397962

Table 2 :

 2 Comparison of Lighter, Blue, and RACER using 5 datasets. This is for finding the best k-value (GL for RACER) using Athena variants vs. exhaustive search. We find either the optimal value or within 0.27% (over Alignment Rate) and within 8.5% (EC Gain) of the theoretical best (in the worst case), consistent with the reported results by Lighter (Figure5in[START_REF] Song | Lighter: fast and memory-efficient sequencing error correction without counting[END_REF]). We also notice that for RACER, GL found by Athena is within 3% of the reference GL (except for the RNN model with D5, which still achieves very close performance for both Alignment Rate and EC Gain.)

	Dataset	Exhaustive Search	With Athena (RNN)	With Athena (N-gram)
					Lighter
		Selected	Alignment	EC	Selected	Alignment	EC	Selected	Alignment	EC
		k	Rate (%)	Gain (%)	k	Rate (%)	Gain (%)	k	Rate(%)	Gain (%)
	D1	k=17	98.95%	96.3%	Same as Exhaustive Search	Same as Exhaustive Search
	D2	k=15	61.42%	73.8%	k=17	61.15%	80.1%		Same as RNN
	D3	k=15	80.44%	86.78 %	k=17	80.39%	95.34%	Same as Exhaustive Search
	D4	k=17	93.95%	89.87%	Same as Exhaustive Search	Same as Exhaustive Search
	D5	k=17	92.15%	81.7%	k = 25	92.09%	83.8%	Same as Exhaustive Search
					Blue
	D1	k=20	99.53%	99%	k=25	99.29%	98.6%	Same as Exhaustive Search
	D2	k=20	57.44%	4.61%	Same as Exhaustive Search	Same as Exhaustive Search
	D3	k=20	84.17%	99.2%	Same as Exhaustive Search	Same as Exhaustive Search
	D4	k=20	95.31%	98.5%	Same as Exhaustive Search	Same as Exhaustive Search
	D5	k=20	92.33%	88.9%	Same as Exhaustive Search	Same as Exhaustive Search
					RACER
	D1	GL=4.7M	99.26%	84.8%	Same as Exhaustive Search	Same as Exhaustive Search
	D2	GL=4.7M	81.15%	92.9%	Same as Exhaustive Search	Same as Exhaustive Search
	D3	GL=3.7M	84.11%	88.27%	Same as Exhaustive Search	Same as Exhaustive Search
	D4	GL=4.2M	95.33%	97%	Same as Exhaustive Search	Same as Exhaustive Search
	D5	GL=4.2M	92.29%	81.63%	GL=20M	92.28%	80.5%	Same as Exhaustive Search

Table 3 :

 3 Comparison of Overall Alignment Rate between Fiona's and RACER's (with and without Athena's tuning). Columns 5 & 6 demonstrate the strong anti-correlation values between Perplexity and Alignment Rate. The last two columns show the assembly quality (in terms of NG50) before and after correction by Racer tuned with Athena. Improvements in NG50 are shown between parentheses

	Dataset	Correlation (N-Gram)	Correlation (RNN)	Fiona + Bowtie2 (Alignment Rate)	RACER w/o Athena + Bowtie2 (Alignment Rate)	RACER w/ Athena + Bowtie2 (Alignment Rate)	NG50 of Velvet w/o EC	NG50 of Velvet w/ (Racer+Athena)
	D1	-0.977	-0.938	99.25%	85.01%	99.26%	3019	6827 (2.26X)
	D2	-0.981	-0.969	73.75%	58.66%	81.15%	47	2164 (46X)
	D3	-0.982	-0.968	83.12%	80.79%	84.11%	1042	4164 (4X)
	D4	-0.946	-0.930	95.33%	93.86%	95.33%	118	858 (7.27X)
	D5	-0.970	-0.962	92.34%	90.91%	92.29%	186	2799 (15X)

Table 4 :

 4 Search time comparison for estimating the perplexity metric with Athena (N-gram) for a point in search space vs. estimating overall alignment rate with Bowtie2

				Dataset
	D1		D2	D3	D4		D5
	Athena Bowtie2 Athena Bowtie2 Athena Bowtie2 Athena Bowtie2 Athena Bowtie2
	1m 38s 10m 5s	49s	3m 53s 1m 39s 7m 50s	52s	3m 8s	1m 40s 9m 42s

Just as Athena is the Greek Goddess of wisdom and a fierce warrior, we wish our technique to unearth the genomic codes underlying disease in a fearless war against maladies.

In our evaluation, we use Bowtie2 for alignment and measure alignment rate as a metric to evaluate Athena. However, alignment is not needed for Athena to work.

1. Deletion: We select the index of injection (position from which to delete sequence), which follows a uniform distribution U(0,l -d), where l is the length of the read and d is the length to delete.

2. Insertion: Similar to deletion errors, the index of insertion follows a uniform distribution U(0,l -I), where l is the length of the read and I is the length of inserted segment. Moreover, each base pair in the inserted segment follows a uniform distribution over the four base pairs (A, C, G, or T).

3. Substitution: For this type of synthetic errors, we select k positions of substitutions following a uniform distribution U(0,l), where l denotes the length of the read. Second, we substitute each base pair of the k positions with another base pair following a uniform distribution over the four base pairs (A, C, G, or T). Note that there is one-in-four chance that the substitution results in the same base pair.

We test the proposed approach against two levels of synthetic errors: low error rate (Uniform in [START_REF] Mahadik | Scalable genomic assembly through parallel de bruijn graph construction for multiple k-mers[END_REF][START_REF] Kelley | Quake: quality-aware detection and correction of sequencing errors[END_REF] bp in error per read), and high error rate (Uniform in [START_REF] Song | Lighter: fast and memory-efficient sequencing error correction without counting[END_REF][START_REF] Limasset | Toward perfect reads: self-correction of short reads via mapping on de bruijn graphs[END_REF] bp in error per read). In all cases, the language model is trained using the original data set, mentioned in Table 1, with the natural ambient rate of error. Figure 4 shows the results for the perplexity metric in the data set after error injection (i.e.,, without doing correction), for the three error types and an equi-probable mixture of the three. The results show that for all types of errors, the direct quantitative relation between error rate and perplexity holds-the perplexity values are higher for high error rate injection. One observation is that the values of perplexity for insertion and substitution errors are higher than for deletion errors. For example, insertion and substitution perplexities are 3X & 4.7X the deletion perplexity for high error rate data sets, and 2X & 1.5X for low error rate data sets. This is expected because in deletion errors, the segment of the read after the index of deletion is a correct segment and hence is expected by the language model. On the other hand, for insertion and substitutions, the added synthetic base pairs will most likely construct wrong segments. Such segments will have very low counts in the language model and therefore produce higher overall perplexities. Another observation is that the perplexities for the injection into data set D1 are higher than for D3. This


Orna Almogi

Lena Dankin
email: lenadankin@tau.ac.il

Nachum Dershowitz

Lior Wolf

A Hackathon for Classical Tibetan

Keywords: Tibetan, Buddhist studies, hackathon, stemming, segmentation, intertextual alignment, text classification

We describe the course of a hackathon dedicated to the development of linguistic tools for Tibetan Buddhist studies. Over a period of five days, a group of seventeen scholars, scientists, and students developed and compared algorithms for intertextual alignment and text classification, along with some basic language tools, including a stemmer and word segmenter.

I INTRODUCTION

In February 2016, a group of four Tibetologists (from the University of Hamburg), one digital humanities scholar (from Europe), and twelve computer scientists (from Israel and Europe) got together in Kibbutz Lotan in the Arava region of Israel with the stated goal of developing algorithmic methods for advancing Tibetan Buddhist textual studies. Participants were either recruited by the organizers or responded to an announcement on several mailing lists. See Figure 1. Most of the computer scientists had background in machine learning, and a few of them also had experience with natural language processing (NLP) research, but without any prior experience with Tibetan texts. The computer scientist organizers were quite familiar with programming workshops and contests and thought that the challenges presented by Tibetan texts would pose an ideal opportunity to explore the hackathon format. The hackathon is a short and intense event where computer scientists collaborate to develop software. For that purpose, it was essential to recruit as many software developers as possible. Some of the recruited students had participated in other hackathons. The plan was to have a focused hacking event, with specific goals to work towards, goals that had been provided by the Tibetan scholars.

The six-hour drive down from Tel Aviv (including a stop to admire desert flora) afforded an opportunity for everyone to get to know each other. The back seats of the van were piled high with computer equipment and the kibbutz was to provide the necessary fast internet connection. The isolation of the kibbutz created an intense working environment and encouraged long hours; the stark natural beauty of the location contributed to a shared sense of tranquility of purpose.

Several hackathons have been conducted for the purpose of the development of tools for digital humanities, before and after ours. In June 2015, the eTRAP team organized a hackathon for text reuses (http://www.etrap.eu/tutorials/2015-goettingen/). Twenty-three participants from fifteen different institutes worked on the detection of textual reuses across data in different languages and from different genres, using the TRACER tool [START_REF] Büchler | Informationstechnische Aspekte des Historical Text Re-use[END_REF][START_REF] Büchler | Scaling historical text re-use[END_REF]. More recently, in May 2017, another hackathon took place in Helsinki, which brought together historians, linguists, psychologists, and computer scientists to work on four different tasks, including the analysis of the written media with the political elite (https://www.helsinki.fi/en/helsinki-centre-for-digital-humanities/helsinki-digital-humanitieshackathon-2017-dhh17). In November 2017, The National Library of Israel hosted a 24-hour hackathon dedicated to the goal of developing state-of-the-art tools and applications for national cultural treasures using a iiif server for their large collection of images (https://hackathon.nli.org.il/index-en.html).

The two main tasks that confronted our group that week (February 14-18) were (1) to develop algorithms for finding intertextual parallels that are only approximately the same, and (2) to experiment with algorithmic classification methods for identifying authorship and style. In both cases, the concern was centered on language issues specific to Tibetan.

After a quick lesson in Tibetan, the Buddhist canon, and modern Tibetan encoding conventions for the benefit of the less knowledgeable, the group split into four loose teams, devoted to the following goals: (A) dataset preparation; (B) language tool development; (C) intertextual alignment; and (D) text classification. We describe each of these efforts in turn in the sections that follow. Each team consisted of a few computer scientists, chosen based on the individual background, experience, and interests, plus a Tibetan scholar who provided annotated data sets and analyzed results. Twice a day we held synchronization round-ups, where each team briefed everyone about their progress, discussed their next steps, and raised problems they stumbled across.

II PRELIMINARIES

Tibetan is a monosyllabic language (Tibetan morphemes normally consist of one syllable) belonging to the Tibeto-Burman branch of the Sino-Tibetan family. The language is ergative, with a plethora of (usually monosyllabic) grammatical particles, which are often omitted. Occasionally, the same syllable can be written using one of several orthographic variations, for example, sogs and stsogs. In the case of verbs, the syllable has various inflectional forms that are often homophones, a fact that can result in variants in reading due to scribal errors or lack of standardization. An example of such inflectional forms is sgrub, bsgrubs, bsgrub, sgrubs (present, past, future and imperative, respectively), all of which are homophones. The intransitive form of the verb offers even more inflectional forms that yield homophones with their transitive counterpart, 'grub and grub (present/future and past, respectively). See [START_REF] Beyer | The Classical Tibetan Language[END_REF]] for details about the language.

The Tibetan Buddhist canon consists of two parts: the Kangyur (bKa' 'gyur), which commonly comprises 108 volumes containing what is believed by tradition to be the Word of the Buddha, texts that were mostly translated directly from the Sanskrit original (with some from other languages and others indirectly via Chinese); and the Tengyur (bsTan 'gyur), commonly comprising about 210 volumes consisting of canonical commentaries, treatises, and various kinds of manuals that were written in the seventh to thirteenth centuries and likewise mostly translated from Sanskrit, with some works from other languages and a few originally written in Tibetan. Overall, this corpus contains 77 million occurrences (tokens) of 81,000 different syllable types. The average transcribed syllable length is 3.5 and the average number of syllables in a single document is 18201.

III HACKATHON TASKS

A. Dataset preparation

A prerequisite for the main goals of the hackathon was data with which to work, that is, texts to compare and classify. For this, we took Tibetan Buddhist texts obtained from various sources. These included the Tibetan Buddhist canon in digital form (we used a modified form of the ACIP files of the Kangyur and Tengyur provided by Paul Hackett of Columbia University) and several sets of autochthonous Tibetan Buddhist texts of various authors (compiled by Eric Werner of Universität Hamburg).

In addition, it was necessary to prepare test suites with manually prepared "gold standard" answers, so that the performance of algorithms for finding parallel passages and for classifying texts could be measured. The passages were selected from various sources, particularly from (a) two doxographical texts ('grub mtha'), the gZhung lugs rnam 'byed by Phywa pa Chos kni sengge (1109-1169) and the 'Grub mtha' mdzod by Klong chen pa Dri med 'od zer (1308-1364), the latter including "borrowed" passages from the former [START_REF] Werner | Phywa-pa Chos-kyi-seng-ge's (1109-1169) depiction of Mahāyāna philosophy: A critical edition and annotated translation of the chapters on Yogācāra and Mādhyamaka philosophy from the gZhung lugs rnam 'byed, a doxography of the twelfth century[END_REF], and (b) Rong zom Chos kyi bzang po's (11th c.) collected writings, which features numerous cases of parallel passages. These Tibetan works were provided in textual form, transcribed according to the Wylie convention [START_REF] Turrell | A standard system of Tibetan transcription[END_REF]. In this system, Tibetan is transliterated into Latin characters without diacritics; thus various Tibetan letters are represented by two or three Latin consonants. The decision to work with transliterated texts was made partly because they were the ones available at the time, but also because the computer scientists didn't understand Tibetan script, so this transliteration made it possible for them to progress quickly without the need to acquire a new alphabet. The texts had to be "cleaned" by removing sigla and by standardizing punctuation.

B. Language tools

Since syllables having the same base form may take many different surface forms, stemming is a crucial stage in almost every text-processing task one would like to perform in Tibetan, as for many other languages. So, to support present and future analysis of Tibetan texts, developing a stemmer was one of the first orders of business.

Usually, in Indo-European and Semitic languages, stemming is performed on the word level. However, in Tibetan, in which multisyllabic words are not separated by spaces or other marks, a syllable-based stemming mechanism is required even in order to segment the text into lexical items. Stemming is not the same as (grammatical) lemmatization, and the stemming process can result in a stem that is not itself a lexical entry in a dictionary. Moreover, unlike Indo-European languages, stemming of Tibetan is mostly relevant to verbs and verbal nouns (which are common in the language). Despite being inaccurate in some cases, stemming (for Tibetan, as for other languages) can improve tasks such as word segmentation and the detection of intertextual parallels [START_REF] Enjamin Klein | Finding inexact quotations within a Tibetan uddhist corpus[END_REF]. Even for Tibetan words consisting of more than one syllable, stemming each "substantial" syllable (i.e. excluding grammatical particles) makes sense since all the inflections are embedded at the syllable level. For instance, the words brtag dbyad (analysis) and brtags dpyad (analyzed) are stemmed to rtog dpyod (to analyze, analysis).

The stemmer we developed is a rule-based application that works in the following manner: first, the syllable is divided into a sequence of Tibetan letters. This stage is required because the Wylie transliteration scheme represents some Tibetan letters by more than one character (e.g. zh, tsh). There is, fortunately, no ambiguity in the process of segmentation into Tibetan letters. By design, the transliteration ensures that whenever a sequence of two or three characters represents a single letter, it cannot also be interpreted in context as a sequence of distinct Tibetan letters.

For the analysis of the Tibetan syllable we used an octuple (8-component) scheme: Each Tibetan syllable should contain one core letter and one vowel. Other positions (subscript, superscript, coda, prescript, postscript, and appended particle) are not obligatory. Each position contains a single letter, except for that of the appended particle, which can be any of six syllables. The "stem" of a syllable is defined by us as consisting of the core letter or stacked letter (which, in turn, consists of the core letter and a superscript or a subscript, or both), the vowel (syllabic contractions contain two vowels at most), and the coda (if extant). Syllables can be considered stemmically identical if these are consistent, despite additions or omissions of a prescript and/or a postscript. The final stage of the stemming is normalization, since there are groups of Tibetan letters that can be replaced one with another without changing the basic meaning of the syllable (in inflectional forms). Since the goal is to group all syllables that are ultimately stemmically identical into one and the same stem, we normalized all tuples according to an elaborate set of rules.

The stemmer, as described, extracts the information encoded in each Wylie transliterated syllable and makes it explicit. An important task, given two syllables, is to evaluate their stemmic similarity. Some substitutions can be considered silent or synonymous; others change the meaning completely; and there is a continuous spectrum in between. Metric learning algorithms were used to assess the relative importance of each substitution.

Another important language task is word segmentation, that is, grouping syllables into words (lexical units). Since no spaces or special characters are used to mark word boundaries, the reader has to rely on language models to detect the word boundaries. As opposed to the stemming task, we had recourse to an annotated corpus for the segmentation task, that is, a word-segmented corpus, with which it was possible to train a supervised model. The training data that was used, consisting of 37,000 sentences, was obtained from the Tibetan in Digital Communication project (http://larkpie.net/tibetancorpus).

The approach taken at the hackathon was based on a flavor of recurrent neural networks (RNNs) called "long short-term memory" (LSTM) [START_REF] Hochreiter | Long short-term memory[END_REF]. LSTMs have been used in the past for word segmentation of Chinese text [START_REF] Chen | Long short-term memory neural networks for Chinese word segmentation[END_REF]. The tuple representation of syllables was used for this purpose; see details in [Almogi et al., 2016]. Several LSTM setups were compared; the best configuration yielded an F1 score of 0.95. In addition, a more traditional algorithm, the conditional random field (CRF), was applied to the data, yielding a lower F1 score of 0.89. This technique was previously applied on Tibetan script in [START_REF] Liu | Tibetan word segmentation as syllable tagging using conditional random field[END_REF].

It bears noting that our efforts to train a word2vec model [START_REF] Mikolov | Efficient estimation of word representations in vector space[END_REF] to represent Tibetan syllables did not result in a solid representation, in the sense that pairs of vectors with high (cosine) similarity did not usually represent synonyms. For that reason, the vector representation that was developed for the stemmer was also essential for the word segmentation task.

Both the stemmer and word segmenter have been made publicly available and can be accessed from http://www.cs.tau.ac.il/~nachumd/Tools.html. Additional details may be found in [Almogi et al., 2016].

C. Intertextual alignment

The primary goal of the hackathon was to develop and compare tools for finding parallel passages between Tibetan texts that are the result of either acknowledged citations (with or without attributions) or borrowing (i.e. with no acknowledgement whatsoever). Generally, for determining the history of composition or relative chronology of a text, passages need not match precisely. That is, in addition to the fact that orthographical differences or omission/addition of grammatical particles are of no great significance, it is often the case that cited or borrowed passages are not necessarily reproduced verbatim, but are often slightly paraphrased or shortened, or both. For determining the identity of persons involved in the composition of the text and its transmission-that is, the author, translator, scribe, or editorthe precision of the match is of greater significance, and even variation in orthography or omission/addition of grammatical particles may be relevant. In this regard, however, textual scholars take into consideration that texts were often copied and edited and that through these processes changes could have been introduced into the text, either deliberately-particularly in terms of standardization of orthography and verb inflection, employment of particles, and even substitutions of terminology in cases of archaism-or unintentionally.

Broadly speaking, there are two cases of interest: (a) an approximate alignment of what could be considered to be exactly the same text, that is, an alignment that allows variants that are considered accidental or non-substantial (that is, variations regarding omission/addition or different forms of the same grammatical particles, orthography, inflectional forms in the case of verbs, archaism vs. standardization, and the like), and (b) an approximate alignment of passages that contained the same text but in modified form of some sort, that is, an alignment that allows substantial variants in addition to the non-substantial ones (omission/addition of a substantial syllable, replacement of a substantial syllable by a completely different one, omission/addition of a string of syllables, occurrence of the same syllables in a different order, and the like). To address the problem of substantial variants that could occur also when a (more or less) exact citation or borrowing was intended, that is, such that have been intentionally introduced by either the author himself or by the scribes and editors during the process of transmission, or such that have been unintentionally crept in during the processes of composition and copying, a limited number of substantial variants must be admitted as well.

Three algorithms competed with one another on this task during the hackathon.

1. One algorithm was TRACER [START_REF] Büchler | Informationstechnische Aspekte des Historical Text Re-use[END_REF][START_REF] Büchler | Scaling historical text re-use[END_REF], based on the "bag of words" representation method. TRACER is a general text reuse detection algorithm with a seven-level architecture. Each step is configurable and can be optimized to specific text reuse tasks and corpora. The steps are preprocessing, featuring, selection, scoring, and post-processing. This approach is called feature-based linking, where only text-reuse units with shared features are compared, as opposed to the comparison of the full text of passages, all against all. All passages are compared by comparing the words they contain, ignoring word order.

2. Another method was based on Agents for Actors (AfA) [START_REF] Küster | Agents for Actors: A Digital Humanities framework for distributed microservices for text linking and visualization[END_REF], a "digital humanities framework for distributed microservices for text analysis". AfA was originally developed for the purpose of identifying allusions to Shakespearean passages in transcriptions of dialogues in films (hence "actors" in its name). This algorithm compares passages both on the letter and the word level, and therefore catches variations at the orthographic and formulation levels, respectively. While its primary use is to identify references and allusions in texts, in the hackathon, the algorithm was tested to see how well it can also serve to identify parallel passages for very different types of texts in an unrelated language.

3. The third approach was based on an adaptation of the method of [Barsky et al., 2008], designed for matching DNA subsequences, to our problem, as described in [START_REF] Enjamin Klein | Finding inexact quotations within a Tibetan uddhist corpus[END_REF]. This algorithm looks for "all against all approximate matches" (within some given threshold of difference between passages) by rephrasing the problem as finding maximal paths in a matching graph. That method was modified during the hackathon to work with syllable stems as the basic building block, rather than the individual character level used before. This change improved both the run time and the quality of the results. Since, on average, a syllable has 4 characters, the speedup was two orders of magnitude. As for the results, p@10 ("precision at ten", the fraction of the top ten results that are of relevance) increased from 0.67 to 1, and p@20 increased from 0.37 to 0.63. The improvement were due to the fact that with character-wise alignment syllables can share many letters but have no semantic similarity; see [START_REF] Labenski | Intertextuality in Tibetan texts[END_REF][START_REF] Labenski | Finding Inter-textual Relations in Historical Texts[END_REF]].

An infrastructure subteam, in addition to keeping everything up and running, parallelized the implementation of the third algorithm to run on a Sparc cluster of computers, located at Tel Aviv University. This is necessary for the ultimate goal, considering the large size of the corpus. The idea is simple: divide the texts into overlapping chunks; then run the original algorithm on all chunks in parallel; finally, piece all the results together.

All three algorithms were tested on a test set that was designed during the hackathon. The two doxological texts mentioned above and known to contain many shared passages were chosen, and 24 pairs of parallel passages were manually annotated. Out of the 24 pairs, the TRACER algorithm retrieved 13 pairs, the AFA algorithm retrieved 12 pairs, and the APBT algorithm retrieved 16.

By finding cited or borrowed passages within the corpora of Indo-Tibetan (i.e. translated) and Tibetan (i.e. autochthonous) Buddhist literature, several research questions can be better addressed:

• determining the history of composition of individual texts;

• determining relative chronology of groups of texts;

• determining the intellectual scholarly milieu in which the texts emerged; and

• determining the intellectual history behind the texts (viz. terminology and concepts).

After identifying parallel passages, one can assess the frequencies of letter/syllable/word replacements in the aligned passages of selected texts or text groups. This can serve to help answer further research questions like: determining editorial policies and processes, such as standardization of orthography, standardization of employment of grammatical particles (i.e. according to the so-called sandhi rules); and identifying processes of "revisions" of translated texts.

D. Text classification

The second major task that was addressed at the hackathon was the question of author profiling. While the question as to what extent the issue of authorship can be addressed in the case of translated texts is yet to be looked into carefully, some general research questions related to authorship fall under the purview of machine classification. These include the following: (a) distinguishing between translated texts and autochthonous texts; (b) identifying the period in which a text was composed, viz. Old Tibetan (7-11th c.), Classical Tibetan I (11-14th c.), or Classical Tibetan II (15-20th c.); (c) determining whether a translated canonical work belongs to the early period of translation (snga 'gyur) or the new period (phyi 'gyur); (d) in the case of autochthonous literature, differentiating between the so-called "revealed" texts (texts that are portrayed as having been transmitted supernaturally) versus "composed" texts; and (e) identifying an author's intellectual milieu (e.g. affiliation with a particular school of thought).

A series of experiments were performed on scriptures and treatises, early and late, translated and autochthonous texts. We tried several methods, including bag-of-word features and a perceptron classifier with stochastic gradient descent with features similar to [START_REF] Volansky | n the features of translationese[END_REF], mainly: mean syllable length; mean sentence length; frequency of verbal prefixes and function words; frequency of foreign (Sanskrit) words; and type-to-token ratio. For authorship detection, we first used an automatic word segmenter and then used n-gram frequency and bag-of-words as features. Such a method was shown to be useful in [Koppel et al., 2008]. We didn't advance further in this task, due to a shortage of time.

Both parts of the canon were employed as training data to determine features that are peculiar for the Kangyur, the corpus containing scriptures, on the one hand, and the Tengyur, the corpus containing treatises, commentaries, manuals and the like, on the other. Numerous autochthonous texts, including the entire collected writings of Rong zom Chos kyi bzang po, the entire collected writings of Shākya mchog ldan (1428-1507), several works by Sa kya paṇḍi ta Kun dga' rgyal mtshan (1182-1251), and several texts by Tsong kha pa Blo bzang grags pa (1357-1419) were tested against the translated canonical texts in order to determine features of translated versus autochthonous works. In addition, selected individual texts were tested. For example, Sa skya paṇdi ta's Tshad ma rigs gter was compared with Dharmakīrti's (7th c.) Pramāṇavarttika in Tibetan translation, which enabled a comparison of autochthonous versus translated work on similar topics. The Mañjuśrīnāmasaṅgīti commentary ascribed to Rong zom pa (and at the same time included in the Tengyur as an Indian work in Tibetan translation) was compared with the canon in its entirety, as was the Tengyur alone with other works by Rong zom pa and additional autochthonous works, which provided a comparison of works whose origin has been considered doubtful with translated and autochthonous literature. The classification results are undergoing analysis by the Tibetan scholars.

IV CONCLUSION

The intense hackathon format proved to be quite exhilarating. Towards evening, each group reported on the day's accomplishments and vicissitudes. No single task was actually brought to completion on site, but the saplings were planted, and the ideas and prototype tools have continued to grow and develop in the ensuing weeks.

Based on our experience, we would recommend such a hackathon format for other welldefined interdisciplinary efforts in the computational humanities. It pays to come wellprepared to the event with clear goals and clean test data. And it is crucial to allocate resources for bringing the products and results of the hackathon to a stable and useful state after the event.

As a matter of fact, the authors held a second hackathon one year later (February 2017) on a kibbutz in the Galilee, again for the development of tools for Tibetan Buddhist texts, but this time concentrating on manuscripts and computer-vision aspects.

Figure 1 .

 1 Figure 1. Poster announcement of the hackathon.

Journal of Data Mining and Digital Humanities http://jdmdh.episciences.org

http://jdmdh.episciences.org

Acknowledgements

We thank the staff at Kibbutz Lotan and all the hackathon participants (listed below). This research was supported in part by a grant (#I-145-101.3-2013) from the German-Israeli Foundation for Scientific Research and Development, and by the Khyentse Center for Tibetan uddhist Textual Scholarship, niversit t Hamburg, thanks to a grant by the Khyentse Foundation. N.D.'s research benefitted from a fellowship at the Paris Institute for Advanced Studies (France), with the financial support of the French state, managed by the French National Research Agency's "Investissements d'avenir" program (ANR-11-LABX-0027-01 Labex RFIEA+). Hackathon participants: rna Almogi, Kfir ar, Marco üchler, Lena Dankin, Nachum Dershowitz, Daniel Hershcovich, Yair Hoffman, Marc W. Küster, Daniel Labenski, Peter Naftaliev, Dimitri Pauls, Elad Shaked, Nadav Steiner, Lior Uzan, Dorji Wangchuk, Eric Werner, and Lior Wolf. Participating institutions: Tel Aviv University (School of Computer Science); Universität Hamburg (Khyentse Center for Tibetan Buddhist Textual Scholarship, Department for Indian and Tibetan Studies); Georg-August-Universität Göttingen (Göttingen Centre for Digital Humanities).