XenofilteR: computational deconvolution of mouse and human reads in tumor xenograft sequence data

Background Mouse xenografts from (patient-derived) tumors (PDX) or tumor cell lines are widely used as models to study various biological and preclinical aspects of cancer. However, analyses of their RNA and DNA profiles are challenging, because they comprise reads not only from the grafted human cancer but also from the murine host. The reads of murine origin result in false positives in mutation analysis of DNA samples and obscure gene expression levels when sequencing RNA. However, currently available algorithms are limited and improvements in accuracy and ease of use are necessary. Results We developed the R-package XenofilteR, which separates mouse from human sequence reads based on the edit-distance between a sequence read and reference genome. To assess the accuracy of XenofilteR, we generated sequence data by in silico mixing of mouse and human DNA sequence data. These analyses revealed that XenofilteR removes > 99.9% of sequence reads of mouse origin while retaining human sequences. This allowed for mutation analysis of xenograft samples with accurate variant allele frequencies, and retrieved all non-synonymous somatic tumor mutations. Conclusions XenofilteR accurately dissects RNA and DNA sequences from mouse and human origin, thereby outperforming currently available tools. XenofilteR is open source and available at https://github.com/PeeperLab/XenofilteR. Electronic supplementary material The online version of this article (10.1186/s12859-018-2353-5) contains supplementary material, which is available to authorized users.


Background
Cancer research heavily relies on model systems such as cell lines. These cell lines have typically been cultured for decades and only partially recapitulate the genetic features of patient tumors [1]. More advanced clinical cancer models are the cell line-derived xenograft and patient-derived xenografts (PDX) [2]. With this approach, either a cancer cell line or a patient tumor sample is injected or transplanted into a host, generally immunodeficient mice. In these xenografts, the complex interactions between the tumor and its microenvironment are at least partially recapitulated, as is the heterogeneity in tumors in the case of PDX [3][4][5][6][7][8]. For these reasons, xenograft models might serve as a better proxy for human tumor samples and have become indispensable for development, validation and optimization of cancer treatment regimens [1,2,9]. Despite its limitations [8,10], the wide applicability of PDX, and more generally of tumor xenografts, is reflected by tens of thousands publications describing numerous biological, mechanistic and preclinical applications [11][12][13][14][15][16].
In spite of this tremendous popularity, sequence analysis of RNA or DNA from tumor xenograft and PDX samples is challenging: the sequence data contain not only DNA and RNA from the grafted human tumor cells but also from the mouse, mostly due to infiltrating stromal cells [17]. When sequencing the combined 'bulk tumor' DNA, sequence reads originating from the mouse result in false positive single nucleotide variants (SNV) when calling mutations [18]. Similar challenges are observed when sequencing RNA: beside false positive SNVs, the gene expression levels are often obscured by reads that derive from mouse cells [19]. Despite mouse-derived sequence reads representing a potential source of bias in sequence analysis of tumor xenografts, the number of tools to solve this important issue is surprisingly limited.
Some solutions have been proposed to bioinformatically remove mouse host sequences from the analysis. The most straightforward method is to map all reads first to the mouse reference genome. Sequence reads failing to map are remapped to the human reference, which is followed by standard downstream analyses [20]. A major disadvantage of this method is that human reads from evolutionary conserved regions will also map to the mouse reference genome. Such reads are inadvertently removed from further analysis, which erodes the read depth and thus sensitivity of variant detection in DNA sequencing. Similarly, it erodes gene expression estimates (counts) when sequencing RNA.
An improved version of this concept, developed for RNA sequence data but also applicable to DNA sequence data, uses a so-called k-mer approach with a mixed mouse/human reference set [19]. This method catalogs for every possible sequence of length k, its presence in the human and mouse reference genome sequences. If a k-mer is unique to one reference, its occurrence in sequencing data indicates the species' origin. Distinction between conserved regions, which are also the most problematic in cross strain filtering, would require long k-mers. However, k-mer elongation rapidly increases computer memory requirements and is therefore less feasible.
Deconvolution based on the alignments of sequence reads to a human and mouse reference genome separately has also been proposed [21,22]. This method utilizes the alignment scores of each sequence read to the mouse and human reference genome to categorize reads as human or mouse. Both methods shows a much better performance as compared to filtering for reads that do not map to the mouse reference genome [19]. However, the number of supported, open-source solutions are limited and improvements in accuracy and ease of use are necessary.
The challenges in the analysis of sequence data from xenografts and the limited availability of tools motivated us to firstly provide a detailed study into the effect of mouse reads on subsequent analyses. Furthermore, we set out to develop a method for accurate filtering on species' origin using a procedure that is easily applicable in bio-informatics pipelines to improve analysis of DNA and RNA sequence data from xenografts.

Implementation
XenofilteR is an easy-to-use R-package for deconvolution of mouse and human sequence reads form xenograft sequence data. XenofilteR takes a file with 2 bam files (e.g. BWA [23], TopHat [24], STAR [25]) for each sample as input: reads aligned to the human reference and reads aligned to the mouse reference genome (Fig. 1). XenofilteR does not require a specific order of the sequence reads for the input BAM files. Default output of XenofilteR is a new bam file with the sequence reads classified as human. Optionally, a second bam file can be generated with the sequence reads classified as mouse.

Filtering
Sequence reads that only map to a single reference genome are classified to that specific organism. For reads that map to both the human and mouse reference genome the edit distance is calculated by summing soft clips, insertions (both derived from the CIGAR string) and the number of mismatches (bam tag: 'NM') ( Fig. 1). For paired-end sequencing, the edit distance of the forward and reverse read is averaged. Sequence reads with an equal edit distance to mouse as well as human are not in either bam file as these cannot be assigned. Assignment of reads (or read pairs) to either human or mouse is based on the edit distance, with reads having a lower edit distance for the reference genome of a species being classified as originating from that species.
Although sequence reads generally map to one specific location on the genome, some reads can be mapped reasonably well to multiple places on the genome, these mappings are called secondary alignments. In XenofilteR, the edit distance is calculated on the primary alignments only. All secondary alignments are either kept in the filtered output or removed depending on the classification based on the primary alignment. Classification can further be fine-tuned by setting a maximum number for the edit distance (default = 4) and a penalty for unmapped reads in case of paired-end sequencing (default = 8).

Parallel implementation and computational time
XenofilteR uses functionality from GenomicAlignments and Rsamtools [26] for reading and manipulating bam files. Parallel analysis is implemented in XenofilteR package using BiocParallel. As XenofilteR Fig. 1 Overview of XenofilteR workflow. Sequence reads (fastq) from PDX are mapped with the appropriate aligner (e.g. BWA, Tophat, STAR) to both a human and mouse reference genome. Sequence reads that only map to a single reference genome are classified to that specific organism. For seqeunce reads that map to both the human and mouse reference genome the edit distance is calculated which is defined by the number of base pairs different between the sequence read and the reference genome. Next, XenofilteR classifies the sequence reads as 'human' or 'mouse' based on the edit distance only evaluates the sequences that map to both reference genomes and requires only little information from the bam files, we were able to minimize the CPU time and memory needed for analysis. Xenofil-teR can be run on a desktop computer in single sample mode and in parallel on computer servers. Examples of code to run XenofilteR and further documentation is available at (https://github.com/PeeperLab/XenofilteR).

Results
Mouse sequence reads map to specific regions on the human genome In xenograft models, human tumors are grown in a murine host. Sequence data of these tumor xenografts commonly contain reads that originate from the host. To investigate which genes and exons are likely to be affected by mouse reads, we mapped whole genome DNA sequence data (WGS) of three mouse strains (NOD/ ShiLtJ, BALB/cJ and C57BL/6NJ) [27,28] to a human reference [29]. On average, 0.3% of mouse reads mapped to the human reference genome, of which 18-20% overlapped with an exon of a protein-coding gene. A high correlation was observed in the number of reads mapped to exons between different mouse strains (R 2 = 0.98, Fig. 2a). Mouse reads mapped to specific regions of the genome with~2000 (out of 200.000 exons in total) exons exceeding 100 reads, including exons from known cancer driver genes [30] (Fig. 2b, Additional file 1: Table S1). Mapping of BALB/cJ WGS data to the human reference revealed that 13% of exons have at least a single mouse read mapped, affecting 43% of genes in total (Fig. 2b). For example, out of the ten exons of BCL9, four exons had over 100 mapped reads mapped, the remaining six had only a few reads or none at all (Fig. 2d). Similar results were observed for other cancer-related genes such as PTEN (Fig. 2c).
Also, RNA sequence data of the same three mouse strains (NOD/ShiLtJ, BALB/cJ, C57BL/6NJ) [27,28] were mapped to the human reference genome. As the sequence similarity between mouse and human is highest for the coding regions, the number of RNA sequencing reads that map to the human reference is much higher (4-8% of reads) compared to WGS. The read count per gene from the RNA sequence data correlated (R 2 = 0.52) with read count per gene in the WGS (Fig. 2e), indicating that the same exonic regions are affected with WGS and RNAseq.
Although mouse RNA sequencing and WGS data clearly showed that mouse reads can map to the human reference genome, both methods were performed on the complete RNA and DNA pools of the sample. Whole exome sequencing (WES) on the other hand, includes an enrichment step using baits designed to target exons on the human reference genome. To test the affinity of mouse sequence reads to the human baits, we sequenced eight mouse DNA samples enriched with a human exome kit (Illumina, SureSelect Human Exon Kit 50 Mb capture set, Agilent, G3362). On average, 29.2 million reads were sequenced per sample of which~11% could be mapped to the human reference genome. Furthermore, 85-86% of mapped reads did so to an exon. These findings were highly reproducible, with a high correlation in exon read count between samples (R 2 = 0.93-0.98) but also with the results from WGS (BALB/cJ, R 2 =~0.47), albeit a higher average read count was observed per exon with WES (Fig. 2f ). To summarize, mouse sequence reads map to specific regions on the human genome, an issue that we have observed for RNA sequencing, WGS and WES.

Sequence reads of mouse origin affect downstream analysis of xenografts
In recognition that mouse reads can map to the human reference genome, we set out to determine the effect that these reads have on analyses of eight PDX samples [14]. For each sample, the percentage of mouse stroma was estimated by two pathologists and averaged (Additional file 2: Table S2). Mutation analysis on the WES data of the PDX samples revealed an extremely high number of single nucleotide variants (SNVs), especially in the samples with a high percentage of mouse stroma. A direct comparison of PDX samples containing a high number of mouse sequence reads, mapped to the human reference, revealed that many of the SNVs in the samples overlap with SNV that originate from mouse, for example in one of the exons of PTEN (Fig. 3a).
Genome-wide mutation analysis on the mouse WGS data mapped to the human reference identified 101,068 SNVs (19.5% exonic). Intersection of this list with the lists of SNVs detected in the PDX samples suggested that many SNVs detected in PDX samples are derived from reads that originate from mouse cells. In the PDX sample M005.X1 (~25% mouse stroma), 73,705 SNVs were detected, of which 67,194 overlapped with the 101,068 SNVs from mouse reads mapped to the human reference. The PDX sample M029.X1 (~1% mouse stroma) had a much lower total number of SNVs, only 460 detected SNVs in the PDX samples overlap with the mouse SNVs (Fig. 3b). In conclusion, sequence reads that originate from mouse have a large effect on mutation calling on samples derived from PDX.
The edit distance can be used to classify sequence reads Accurate assignment of reads to either mouse or human is pivotal to assure high quality downstream analyses. Currently available tools generally use the mapping of reads to a combined reference genome or to both genomes as a classification strategy [18,19]. However, due to the sequence similarity between mouse and human, the mapping itself might not provide the optimal distinction between reads of human and mouse origin.
A striking distinction between the alignments to mouse and human reference was the difference in 'edit distance': the number of base pairs in a given mapped  (Fig. 3c). Thus, the edit distance of a sequence read can be used to filter mouse from human sequence reads. Based on these observations, we developed an algorithm, called XenofilteR, which calculates the edit distance for each read that maps to both the human and mouse reference genomes (Fig. 1). The edit distance is calculated by summing soft clips, insertions (both derived from the CIGAR string) and the number of mismatches (bam tag: 'NM'). The reference genome to which a specific sequence read has the lowest edit distance is considered as the species of origin for that read. By differentiating each sequence read in the original input bam files, XenofilteR generates a new bam file, which contains the sequence reads classified as human only. Conversely, XenofilteR can also output the bam file with all reads classified as mouse. XenofilteR is programmed in R and publically available from GitHub (https://github.com/PeeperLab/XenofilteR).

XenofilteR accurately filters mouse reads from human reads from in silico-mixed datasets
To validate this computational method and compare the results to other available methods, we generated fastq files containing both mouse [27,28] and human [29] WGS reads. We generated paired-end and single-end fastq files of different sequence length and multiple percentages of mouse cells ( Fig. 4a and Additional file 3: Table S3). These files were generated for two mouse strains (BALB/cJ, C57BL/6NJ; a full description on how the files were generated is available in the methods section). The combined fastq files were mapped to both human and mouse references (C57BL/6NJ). We applied five tools to the generated data: XenofilteR Strict filtering (filtering of all reads that map to mouse), bamcmp [21], BBsplit [22], Xenome [19] and XenofilteR (all with default settings). Since the origin of each read was known, we could calculate the accuracy of each of the three methods. Because the C57BL/6NJ mouse strain is identical to the mm10 reference genome the most accurate classification was reached with this mouse strain compared to BALB/cJ (Additional file 3: Table S3).
Results from the dataset with mixed human and BALB/ cJ reads strain shows that for all tools true and false positive classification of reads as human depend on both sequence length and on whether sequencing was paired-end, but not on the initial percentage of mouse reads in the mixture (Fig. 4b and Additional file 3: Table S3). Although the Strict filtering method showed the least misclassified mouse reads (0.01%), it was accompanied by a severe decrease in the number of correctly assigned human reads (Fig. 4b). By contrast, both XenofilteR and Xenome correctly identified almost all mouse reads with, respectively, less than 0.02 and 0.04% of mouse sequence reads remaining after filtering. Bamcmp retained the highest number of human reads but also kept a high percentage of mouse sequence reads, especially for the paired-end sequence runs (> 0.20%). Similar results were observed for BBsplit, except that a high number of mouse sequence reads were kept both with single-end and paired-end sequencing ( Fig. 4b and Additional file 3: Table S3).
In addition to the WGS of in silico mixed samples, we also determined the effect of filtering on the detection of somatic variants in a cancer sample. For this purpose, we mixed in silico WES sequence reads of a patient tumor sample M032 [8] with those obtained from mouse WES, with both sequence libraries generated using the same human exome enrichment kit, in a 3:1 ratio. This sample was processed in parallel with Bamcmp, Xenome and XenofilteR. Due to the high number of erroneously filtered sequence reads the performance of both the Strict Filtering method and BBsplit was not further investigated. All three methods were run with default settings followed by mutation calling (GATK). In the original tumor sample, 419 somatic SNVs were detected; in the mixed sample, without exclusion of mouse reads, a total of 107,826 SNVs were observed, comparable to the number of SNVs in PDX sample M005.X1. Filtering with Bamcmp, Xenome or XenofilteR resulted in 547, 449 and 438 SNVs, respectively. The 438 SNVs remaining after XenofilteR filtering included all 419 SNVs identified in the original samples, with almost identical VAFs (Fig. 4c), and an additional 15 false positive SNVs (Fig. 4d). This is an improvement over Bamcmp and Xenome, which both produced more false positives, 128 and 30 respectively (Fig. 4d). In addition, for two SNVs, the VAF was lower after filtering compared to the original tumor (Fig. 4c). Thus, when filtering samples with in silico-mixed mouse and human sequence reads, XenofilteR improves on Bamcmp and Xenome both regarding total number of filtered sequence reads and in retaining mutations of human origin.

XenofilteR accurately filters mouse reads from human reads in PDX samples
In addition to in silico-mixed samples, we tested Xeno-filteR on PDX samples and compared the results to those obtained with the best performing method on the in silico data, Xenome. Patient tumor, normal and PDX were analyzed by WES for three breast cancer samples. Mutations were called on these samples after XenofilteR or Xenome filtering (Fig. 5a). For each SNV identified in the filtered PDX, we traced whether it was either also found in the patient tumor ( Fig. 5a; black), in the matched normal or SNP database ( Fig. 5a; green) or not found in either blood, SNP database or tumor sample ( Fig. 5a; red). This last group represents either false positives or a difference between PDX and tumor (e.g. due to tumor heterogeneity or alternate sequence depth between patient tumor and PDX). However, similar to mutation calling in the in silico-mixed sample, the VAF was much lower for several mutations identified with Xenome compared to XenofilteR. This was reflected not only by the VAF but also by the read counts, on which A C D B Fig. 5 Performance of XenofilteR and Xenome on PDX samples. a: Mutation calling on exome sequence data of a breast cancer PDX sample. The variant allele frequency (VAF) was plotted after filtering with XenofilteR (x-axis) and Xenome (y-axis). Plotted in black are mutations also detected in the patient sample, in green known SNPs and in red SNVs detected in the PDX only. b: Read count of each SNV used to calculate the VAF from A for Xenome and XenofilteR. c: Mutation calling on targeted sequencing of melanoma samples. In green all known SNPs are indicated, in black the remaining SNVs. d: Validation of the SNP rs7121 (GNAS) by Sanger sequencing with human-specific primers the VAF was based (Fig. 4b): they were fewer after filtering with Xenome compared to XenofilteR in almost all cases. This suggests that Xenome might filter too stringently, which results in multiple SNPs and SNVs in the patient tumor to receive a VAF estimate below the true value.
In addition to the three PDX breast samples, we tested ten melanoma PDX samples for which targeted sequencing (using a 360-cancer gene panel) was performed [14]. We again calculated the VAF and indicated known SNPs (Fig. 5c and Additional file 4: Figure S1). Since only PDX were sequenced, no estimate exists for true somatic or germline mutations. Strikingly, and similar to the breast cancer analysis, the VAF of multiple SNVs and SNPs were lower after filtering with Xenome, compared to Xenofil-teR. Again, this suggests that XenofilteR filters are more sensitive, which contributes to its performance.
To further validate these findings, we selected two SNPs with discordant VAFs between XenofilteR and Xenome after filtering. We developed human-specific primers to perform Sanger sequencing on both SNPs. SNP rs7121, located in the gene GNAS, harbored a C > T change, in M041.X1and M046.X1, but not in M043R.X1, in concordance with the WES data. Also, the expected VAF of 50% was observed in the Sanger sequencing in M046.X1 and the VAF of~25% was reflected in the lower peak for T in M041.X1 (Fig. 5d). SNP rs2071313, located in the gene MEN1, showed a G > T change in M041.X1 and M046.X1. Sanger sequencing revealed the SNP in M041.X1 as heterogeneous corresponding to the VAF after filtering with XenofilteR (Additional file 5: Figure S2A). In addition to the lower VAF, the number of sequence reads was much lower after filtering with Xenome, indicative of XenofilteR better representing the real VAF (Additional file 5: Figure S2B). Altogether, we conclude that XenofilteR outperforms Xenome for the analysis of mutation data of mixed human/mouse origin as illustrated by both in silico mixed data and subsequent corroboration in PDX samples from breast cancer and melanoma patients.

XenofilteR allows for filtering of RNA sequencing data
The effect of mouse sequence reads on downstream analysis of PDX samples is not limited to DNA sequencing but affects RNA sequencing also. The method used by XenofilteR, for which classification is based on the edit distance of a read, can also be applied to RNA sequencing data, as the same values to calculate the edit distance are available in the BAM files (CIGAR and the tag NM). To validate whether indeed, filtering of RNA sequence PDX data can be accurately done, we applied XenofilteR on a set of seven PDX samples for which matched patient samples were available [14].
The effect of XenofilteR on the read counts in RNA sequence data was tested using two different samples, one with a high percentage of mouse cells (M005.X1, pathologist estimate was 25% of mouse cells) and one with a low percentage of mouse cells (M019.X1, 1% mouse cells). As expected, the largest difference between filtered and unfiltered read count was observed for sample M005.X1 (Fig. 6a).
Next, we compared the top differentially changed genes (FDR < 0.001) between filtered and unfiltered samples and generated a heat map and cluster analysis including the original patient samples (Fig. 6b). As expected, samples with the highest percentage of mouse cells also showed the highest expression of the selected genes. Most importantly, after filtering with XenofilteR the expression of the selected genes better reflected the expression of the genes in the patient samples.
We also tested XenofilteR on a large data set of 95 melanoma PDX RNA profiles. Although XenofilteR was initially developed to remove infiltrating mouse reads from PDX samples, we investigated whether the we could also use the method to select for mouse reads. For this purpose, XenofilteR was run on this large PDX cohort to remove the reads of human origin, leaving the mouse reads. As expected, considerable variation was observed with regards to the number of sequence reads classified as mouse, with a range from 408,145 to 20,725,475 sequence reads, with on average 6.1% of the total sequence reads classified as mouse (range: 1-35%). Cluster analysis based on the mouse read counts of the top 250 most variable genes showed separation in three clusters with clear expression patterns in specific samples for clusters 1 and 2 (Fig.  6c). Gene Ontology (GO) analysis of cluster 1 (blue) showed that this cluster was highly enriched for genes involved in fat cells and metabolic processes, suggesting the presence of mouse fat cells in this sample (Fig. 6c). We performed the same analysis for the genes in cluster 2 (orange) and found clear enrichment for genes related to muscle cells (Fig. 6c). Both cell types likely represent the predominant components of the murine microenvironment associated with subcutaneous tumor xenografts. Indeed, pathological examination of the H&E stainings confirmed that both fat and muscle cells are abundantly present in these samples (Fig. 6e). We concluded from these data that XenofilteR can be applied to RNA sequencing data as well. Furthermore, we show that gene expression profiles can be generated of exclusively the murine compartment in PDX samples, despite the fact that murine sequence reads represented only a minor fraction of the total number of sequenced reads. Furthermore, based on the murine-specific gene expression profiles, we can identify the predominant cell types surrounding or infiltrating the PDX in the host.

Discussion
High similarity between mouse and human genetics complicates the downstream analysis of both RNA and DNA profiles from tumor xenografts, including PDX models. The key problem is that a substantial percentage of sequence reads of mouse origin can be mapped with high quality to the human reference genome. These murine sequence reads map to conserved regions between mouse and human DNA, mainly exonic regions of the human genome, many of which correspond to cancer-related genes, e.g. PTEN and BCL9. Mutation analysis on PDX WES data revealed thousands of identified SNVs that were false positives, because they arose from mouse sequence reads. Importantly, we show that enrichment steps with human-specific baits also enrich for mouse DNA, likely due to the high DNA sequence similarity between mouse and human in the coding region of the genome, consistent with previous analyses [27][28][29].
Our in-depth analyses of sequence data from PDX samples demonstrate that the edit distance (the number of base pairs in a given mapped read discording with the reference genome) of a sequence read can distinguish the origin of the sequence reads. Sequence reads, aligned to a reference genome by most standard tools, provide the metrics to calculate the edit distance. On the basis of this premise, we developed XenofilteR, a method to classify sequence reads based on edit distances of sequence reads (Additional file 6: Figure S3). Because the edit distance can be calculated from mapped sequence files from the most frequently used sequence mappers, XenofilteR is a highly versatile method that can easily be implemented in current sequence analysis pipelines. Mapping of the sequence data can be done with the trusted and validated mapper of choice (e.g. BWA [23] or Tophat2 [24]), since XenofilteR requires only a CIGAR string and NM-tag, standard values present in the BAM-format [31]. XenofilteR has been implemented in the open-source programming language R. The R-package provides access to the method for investigators without programming knowledge by easy installation [26], extensive documentation and example data. In addition, we provide a Perl implementation to facilitate alternative integration for more advanced programmers of the XenofilteR algorithm in existing pipelines. We have used an earlier version of this algorithm already for the analysis of PDX models by targeted sequencing, both WES and a small cancer gene panel, and it proved to be an important tool for the identification of mutations causing resistance to targeted therapy [14,15].
The validation results of XenofilteR using in silico mixed samples and PDX samples revealed accurate filtering, allowing for reliable VAF and accurate mutation calling. For example, all mutations detected in patient sample M032 (n = 419) were also detected after mixing of the sample with mouse reads and subsequent filtering with XenofilteR. XenofilteR appeared to filter more sensitively in multiple genomic regions than the Xenome algorithm, identifying more mutations with a more accurate VAF. Compared to Bamcmp, XenofilteR filters more thorough, thereby removing more sequence reads of mouse origin without losing somatic mutations. Although no real mutations were missed in the in silico validation, 15 false-positive mutations were identified after filtering with XenofilteR. Interestingly, most of the sequence reads underlying the false positive SNVs and erroneously classified as 'human' were not mapped to the mouse reference genome, despite being of mouse origin. The importance of a good reference genome is also reflected by the in silico admixed WGS sample where classification of the mouse strain C57BL/6NJ yielded fewer mouse reads after filtering compared to the BALB/cJ strain, this because the reference genome is based on the C57BL/6NJ. This exemplifies the importance of good and accurate reference genomes for any filtering to work and suggests a reference genome based on the mouse host will yield the best results.
In addition to removing mouse reads from PDX data, we demonstrate here the feasibility of analyzing the mouse stromal environment associated with (PDX) human tumor xenografts. We show that expression profiles typical for two major mouse cell types (fat cells and muscle cells) can be recovered from the PDX RNA sequence data when isolating for the mouse sequence reads. The high variability in the number of mouse-derived sequence reads per sample will pose a challenge when comparing the gene expression between samples. However, the possibility to analyze mouse gene expression in tumor xenografts will likely help understand interactions between tumor cells and stromal cells.

Conclusions
In conclusion, we present here the XenofilteR algorithm as a solution for the problem of intermingled murine host and human cells in tumor xenografts. XenofilteR can be applied to both DNA and RNA sequencing and uses the edit distance, providing a straightforward and fast implementation that outperforms currently available methods.

PDX samples
The collection and use of human tissue was approved by the Medical Ethical Review Board of the Antoni van Leeuwenhoek. Animal experiments were approved by the animal experimental committee of the institute and performed according to Dutch law. PDX tumor fragments of~5 mm 3 were used for subcutaneous transplantation into NOD.Cg-Prkdcscid Il2rgtm1Wjl/ SzJ (NSG) mice, which was performed under anesthesia. Before reaching the maximum allowed tumor size, mice were sacrificed, tumors were removed, and tumor pieces were (1) fixed in formalin and embedded in paraffin; (2) snap-frozen and stored at − 80 C for further analyses; (3) cryopreserved in 10% fetal calf serum (FCS) in DMSO and stored at − 80 C for additional passages; or (4) re-transplanted into a new set of NSG mice [14].

Next generation sequencing
PDX tumor material was extracted as previously described [14]. Briefly, melanoma PDX samples from 6 patients and 6 matched blood samples were analyzed with WES [14]; 10 melanoma PDX samples were analyzed with targeted sequencing [14]; 6 melanoma PDX samples and matching patient tumor samples were analyzed with RNAseq 13 and 3 breast cancer PDX samples with matching blood normal as well as matching patient tumor were also analyzed with WES [14].
Whole exome sequencing (WES) was performed as described in Kemper et al. [14]. Exome enrichment was performed using the Agilent SureSelect Human Exon Kit 50 Mb capture set (Agilent, G3362). Paired-end 75 reads of targeted-enrichment libraries were sequenced on the HiSeq 2000. Reads were mapped to the Ensembl human reference (hg19) by bwa 0-7.5 with default settings [23]. BAM files were processed using Picard [

Merged mouse with human sequence data
For validation of XenofilteR and comparison of Xenofil-teR with other tools we generated WGS and WES files with artificially mixed mouse and human reads.
WGS data was generated by merging data from WGS from mouse (BALB_cJ And C57 downloaded from: www.sanger.ac.uk/resources/mouse/genomes/) [27,28] and WGS from the 1000Genomes project (NA12878, source: www.internationalgenome.org/data-portal/) [29]. Paired-end sequencing with 100 bp was available for mouse WGS as well as for human WGS. To ensure all reads were mappable, only mapped read pairs were extracted from bam files and converted to fastq using Picard SamTo-Fastq. Read pairs containing adapters or uncertain base calls (N) in either read were filtered.
Using these filtered reads, we generated single-end and paired-end fastq files with 50 bp, 80 bp or 100 bp reads. Mouse and human reads were randomly chosen and mixed in 6 different ratios with 100%, 94%, 88%, 75%, 50% and 0% human reads. All fastq files contained 70 Million reads (SE) or read pairs (PE). The resulting 72 sets were mapped in parallel to GRCh38 and mm10 before analysis with Xeno-filteR, or directly provided to the Xenome pipeline.
The identical strategy was used for the WES data using sequence reads of melanoma tumor (M032) [14] sequence reads from mouse DNA, enriched with a human enrichment kit, in a 4:1 human/mouse ratio.

RNA sequencing
Three sets of RNA sequencing data were used: 3 mouse strains (NOD/ShiLtJ, BALB/cJ, C57BL/6NJ) [27,28], 7 PDX samples for which matched PDX samples were available [14], and 95 PDX samples without matched patient samples. The mouse strains downloaded from The PDX samples with matched patient samples were processed as previously described [14].
For the second set, all 95 PDX samples were uniquely barcoded and pooled into a single stranded library and sequenced. An average of 56 million unique read pairs were sequenced per sample in a range between 33 and 78 million. All samples were mapped with Tophat2 v2.1.0 with the following parameters: --library-type fr-firststrand -g 1 -p 8 -G ENSEMBL_Annotation_v82.gtf. All samples were mapped to human reference GRCh38 (ENSEMBL v82) and mouse reference genome GRCm38 (ENSEMBL v82). The resulting alignments to human and mouse references were provided to XenofilteR (version 1.4, with default settings), to select for reads of mouse origin. Bam files were name sorted with picard tools followed by counting reads with HTseq-count (HTSeq-0.6.1p1) with settings: -m intersection-nonempty -a 10 -i gene_id -s reverse -f bam. Count data generated with HTseq-count was analyzed with DESeq2 [32].

Mutation calling
Variants were called by GATK 2.7-4 using the 'Unified-Genotyper' with default settings, except for "-minIndel-Frac", which was set to 10%. Annotation of the vcf files was performed with annovar (http://annovar.openbioinformatics.org). All variants detected in the germ-line (blood) samples with a Variant Allele Frequency (VAF) over 2% were excluded from further analysis. Variants were further filtered: minimum VAF of 0.1 in at least one of the samples; a minimum of 20× coverage in at least one of the samples; variant positions, listed as a single nucleotide polymorphism (SNP) in the 1000 Genome project, were excluded, except when also present in COSMIC [33]; Variant position were kept only if annotated as exonic by RefSeq (Release 45) and only if the change was non-synonymous.

XenofilteR
All data shown in this paper were analyzed with XenofilteR version 1.4 as available through Github, http://github.com/ PeeperLab/XenofilteR/releases/tree/V1.4. All samples described here were processed with the default settings.