The amount of sequence information in public repositories is growing at a rapid rate. Although these data are likely to contain clinically important information that has not yet been uncovered, our ability to effectively mine these repositories is limited. Here we introduce Sequence Bloom Trees (SBTs), a method for querying thousands of short-read sequencing experiments by sequence, 162 times faster than existing approaches. The approach searches large data archives for all experiments that involve a given sequence. We use SBTs to search 2,652 human blood, breast and brain RNA-seq experiments for all 214,293 known transcripts in under 4 days using less than 239 MB of RAM and a single CPU. Searching sequence archives at this scale and in this time frame is currently not possible using existing tools. The National Institutes of Health (NIH) Sequence Read Archive (SRA) 1 contains ~3 petabases of sequence information that can be used to answer biological questions that single experiments do not have the power to address. However, searching the entirety of such a database for a sequence has not been possible in reasonable computational time. Some progress has been made toward enabling sequence searches on large databases. The NIH SRA provides a sequence search functionality 2 ; however, the search is restricted to a limited number of experiments. Existing full-text indexing data structures such as Burrows-Wheeler transform 3 , FM-index 4 or others 5-7 are currently unable to mine data of this scale. Word-based indices 8, 9 , such as those used by internet search engines, are not appropriate for editdistance-based biological sequence searches. The sequence-specific solution caBLAST and its variants 10-12 require an index of known genomes, genes or proteins, and so cannot search for novel sequences. Further, none of these existing approaches are able to match a query sequence q that spans many short reads. Here, we use an indexing data structure, Sequence Bloom Tree (SBT), to identify all experiments in a database that contain a given query sequence q. A query is an arbitrary sequence, such as a transcript. The SBT index is independent of eventual queries, so the approach is not limited to searching for known sequences, and the index can be efficiently built and stored in limited additional space. It also does not require retaining the original sequence files and can be distributed separately from the data. SBTs are dynamic, allowing insertions and deletions of new experiments. A coarse-grained version of an SBT can be downloaded and subsequently refined as more specific results are needed. They can be searched using low memory for the existence of arbitrary query sequences. We show that SBTs can search large collections of RNA-seq experiments for a given transcript orders of magnitude faster than existing approaches. RESULTS SBTs create a hierarchy of compressed bloom filters 13, 14 , which efficiently store a set of items. Each bloom filter contains the set of k-mers (length-k subsequences) present within a subset of the sequencing experiments. SBTs are binary trees in which the sequencing experiments are associated with leaves, and each node v of the SBT contains a bloom filter that contains the set of k-mers present in any read in any experiment in the subtree rooted at v (Supplementary Fig. 1) . We reduced the space usage by using bloom filters that are compressed by the RRR 15 compression scheme (Online Methods). Hierarchies of bloom filters have been used for data management on distributed systems 16 . However, they have not previously been applied to sequence search, and we find that this allows us to tune the bloom filter error rate much higher than in other contexts (Theorem 2, Online Methods), vastly reducing the space requirements. Bloom filters have also been used for storing implicit de Bruijn graphs 17,18 , and one view of SBTs is as a generalization of this to multiple graphs. We used SBTs to search RNA-seq experiments for expressed isoforms. We built an SBT on 2,652 RNA-seq experiments in the SRA for human blood, breast and brain tissues ( Supplementary Table 1 ). The entire SBT required only 200 GB (2.3% of the size of the original sequencing data) ( Supplementary Table 2 ). For these data, construction of the tree took ≈2.5 min per file (Supplementary Table 3) . These experiments could be searched for a single transcript query in, on average, 20 min (Fig. 1) , using less than 239 MB of RAM with a single thread (Online Methods). We estimate the comparable search time using SRA-BLAST 2 or mapping by STAR 19 to be 2.2 d and 921 d, respectively (Online Methods), though SRA-BLAST and STAR return alignments whereas SBT does not. However, even a very fast aligner such as STAR cannot identify query-containing experiments as fast as SBT. We also tested batches of 100 queries and found SBT was an estimated 4,056 times faster than a batched version of the mapping 3 0 1 approach (Supplementary Fig. 2) . These queries were performed over varying sensitivity threshold θ (the minimum fraction of query k-mers that must exist in order to return a 'hit') as well as the transcripts per million (TPM) threshold used to select the query set ( Supplementary  Figs. 3 and 4) . For approximately half of the queries, the upper levels of the SBT hierarchy provided substantial benefit, particularly on queries that were not expressed in any experiment ( Supplementary  Fig. 5 and Supplementary Table 4 ). SBTs can speed up the use of algorithms, such as STAR or SRA-BLAST, by first ruling out experiments in which the query sequences are not present. This allows the subsequent processing time to scale with the size of the number of hits rather than the size of the database. We first used SBTs to filter the full dataset consisting of 2,652 human blood, breast and brain RNA-seq experiments. We then compared the performance of STAR or SRA-BLAST on the filtered dataset with the time to process the unfiltered dataset with these algorithms. Using SBTs to first filter the data reduced the overall query time of STAR or SRA-BLAST by a factor of ≈3 (Supplementary Fig. 6) . To analyze the accuracy of the SBT filter, we compared the experiments returned by SBT with those in which the query sequence was estimated to be expressed using Sailfish 20 . Because it is impractical to use existing tools to estimate expression over the entire set of experiments, we queried the entire tree, but estimated accuracy on a set of 100 random files on which we ran Sailfish (Fig. 2) . Three collections of representative queries were constructed using Sailfish, denoted by High, Medium and Low, which included transcripts of length >1,000 nt that were likely to be expressed at a higher, medium or low level in at least one experiment contained in the set of 100 experiments on which Sailfish was run. The High set was chosen to be 100 random transcripts with an estimated abundance of >1,000 TPM in at least one experiment. The Medium and Low query sets were similarly chosen randomly from among transcripts with >500 and >100 TPM, respectively. These Sailfish estimates were taken as the ground truth of expression for the query transcripts. Both false positives and false negatives can arise from a mismatch between SBT's definition of present (coverage of k-mers over a sufficient fraction of the query) and Sailfish's definition of expressed (as estimated by read mapping and an expectation-maximization inference). These two definitions are related, but not perfectly aligned, resulting in some disagreement that is quantified by the false-positive rates (FPR) and false-negative rates of Figure 2 . The observed false negatives are primarily driven by a few outlier queries for which the SBT reports no results but their expression is above the TPM threshold as estimated by Sailfish. This is supported by the fact that the average true-positive rate at θ = 0.7 for queries that return at least one file was 96-100%, and the median true-positive rate across all queries was 100% for all but the strictest θ (Fig. 2) . We used SBT to search all blood, brain and breast SRA sequencing runs for the expression of all 214,293 known human transcripts and used these results to identify tissue-specific transcripts ( Supplementary Table 5 and Supplementary Fig. 7) . This search took 3.3 d using a single thread (Supplementary Fig. 8 ). There are presently no search or alignment tools that can solve this scale of sequence search problem in a reasonable time frame, but we estimate an equivalent search using Sailfish would take 92 d. The speed and computational efficiency of SBTs will enable both individual laboratories and sequencing centers to support largescale sequence searches, not just for RNA-seq data, but for genomic and metagenomic collections as well. Researchers could search for conditions from among thousands that are likely to express a given novel isoform or use SBTs to identify metagenomic samples that are likely to contain a particular strain of bacteria. Fast search of this type will be essential to make good use of the ever-growing collection of available sequencing data. Currently, it is difficult to access all the relevant data relating to a particular research question from available sequencing experiments. Individual hospitals, sequencing centers, research consortia and research groups are collecting data at a rapid pace, and face the same difficulty of not being able to test computational hypotheses quickly or to find the relevant conditions for further study. SBTs enable the efficient mining of these data and could be used to uncover biological insights that can be revealed only through the analysis of multiple data sets from different sources. Furthermore, SBTs do not require prior knowledge about sequences of interest, making it possible to identify, for example, the expression of unknown isoforms or long noncoding RNAs. This algorithm makes it practical to search large sequencing repositories and may open up new uses for these rich collections of data. ( 1 5 -t h r e a d ) Time (min) Figure 1 Estimated running times of search tools for one transcript. The SBT per-query time was recorded using a maximum of a single filter in active memory and one thread. The other bars show the estimated time to achieve the same query results using SRA-BLAST and STAR.
The currently most cited line of work on DNA k-mer indices for approximate search are the Sequence Bloom Trees (SBTs) first proposed by Solomon and Kingsford REF .
Fast search of thousands of short-read sequencing experiments
Nature Biotechnology
Nature Biotechnology
Biology Medicine