Fast lightweight accurate xenograft sorting

Abstract

Motivation: With an increasing number of patient-derived xenograft (PDX) models being created and subsequently sequenced to study tumor heterogeneity and to guide therapy decisions, there is a similarly increasing need for methods to separate reads originating from the graft (human) tumor and reads originating from the host species' (mouse) surrounding tissue. Two kinds of methods are in use: On the one hand, alignment-based tools require that reads are mapped and aligned (by an external mapper/aligner) to the host and graft genomes separately first; the tool itself then processes the resulting alignments and quality metrics (typically BAM files) to assign each read or read pair. On the other hand, alignment-free tools work directly on the raw read data (typically FASTQ files). Recent studies compare different approaches and tools, with varying results.

Results: We show that alignment-free methods for xenograft sorting are superior concerning CPU time usage and equivalent in accuracy. We improve upon the state of the art sorting by presenting a fast lightweight approach based on three-way bucketed quotiented Cuckoo hashing. Our hash table requires memory comparable to an FM index typically used for read alignment and less than other alignment-free approaches. It allows extremely fast lookups and uses less CPU time than other alignment-free methods and alignment-based methods at similar accuracy. Several engineering steps (e.g., shortcuts for unsuccessful lookups, software prefetching) improve the performance even further.

Availability: Our software xengsort is available under the MIT license at http://gitlab.com/genomeinformatics/xengsort . It is written in numba-compiled Python and comes with sample Snakemake workflows for hash table construction and dataset processing.

Keywords: Alignment-free method; Cuckoo hashing; Xenograft sorting; k-mer.

Fig. 1

Illustration of (3,4) Cuckoo hashing with 3 hash functions and buckets of size 4. Left: Each key-value pair can be stored at one of up to 12 locations in 3 buckets. For key x, the bucket given by $f_1(x)$ is full, so bucket $f_2(x)$ is attempted, where a free slot is found. Right: If all hb slots are full, key x is placed into one of these slots at random (blue), and the currently present key-value pair is evicted and re-inserted into an alternative slot

Fig. 2

A k-mer is partitioned into its $\ell$-prefix, a middle base and its $\ell$-suffix. Efficient local re-sorting of k-mers according to common $\ell$-prefix and $\ell$-suffix yields groups of k-mers that differ only in their middle base

Fig. 3

Decision rule tree for classifying a DNA fragment from k-mer statistics $(h,h^{\prime },g,g^{\prime },b,x;n)$, meaning number of k-mers of type “host” (h), “weak host” ($h^{\prime }$), “graft” (g), “weak graft” ($g^{\prime }$), “both” (b), and number of k-mers not present in the key-value store (x), respectively; n is the total number of (valid) k-mers in the fragment. We also use weighted scores $S_{\text{host}}:=h+\lfloor h^{\prime }/2\rfloor$ and $S_{\text{graft}}:=g+\lfloor g^{\prime }/2\rfloor$ and thresholds $T_{\text{host}}:=\lfloor n/4\rfloor ,T_{\text{graft}}:=\lfloor n/4\rfloor$ and $T_{\text{both}}:=\lfloor n/4\rfloor$. A fragment is thus classified as “host”, “graft”, “both”, “neither”, or “ambiguous”. Category “ambiguous” is chosen if no other rule applies and no “else” rule is present in a node

Fig. 4

Classification results of different tools (XenofilteR, xenome, xengsort, and partially xengsort with “quick” option) on several datasets: a GIAB human matepair dataset (XenofilteR did not run on this dataset); b Chicken genome; c Human lymphocytic leukemia RNA-seq data; d Patient-derived xenograft (PDX) RNA-seq data. e CPU times on the PDX RNA-seq dataset with different tools and different xengsort parameters (see text)

Fig. 5

Effect of different prefetching levels on the running time of the adversarial chicken genome dataset: no prefetching (0, default), prefetching the second bucket (1), or prefetching both second and third bucket (2). Times are averages over 4 repeated runs

Fig. 6

Effect of shortcut bits for unsuccessful k-mer searches on the running time of the adversarial chicken genome dataset: no extra bits (0, default), one extra bit (1) or two extra bits (2); see “Performance engineering” section. Times are averages over 6 repeated runs