Cross-Species Read Assignment

Cross-Species Read Assignment: Why Reads in Conserved Genes Are Sorted by Margin, Not by Truth
Zetobit · Bioinformatics Insight Series Model Systems
BIOINFORMATICS INSIGHT SERIES Cross-Species Read Assignment Why Reads in Conserved Genes Are Sorted by Margin, Not by Truth mouse wins human wins tie conserved orthologues Kanna Nandakumar, PhD ZETOBIT
Bioinformatics Insight Series

Cross-Species Read Assignment

Why Reads in Conserved Genes Are Sorted by Margin, Not by Truth

A patient-derived xenograft is two organisms in one tube. Before any expression value, mutation call or signature exists, every read has to be assigned to one of them. Almost all of that sorting is easy, and the part that isn’t is concentrated in exactly the places where the biology is most conserved — which is to say, in the genes most likely to be doing something interesting.

The mechanics are simple enough to state in a sentence. Align each read to the human reference and to the mouse reference, compare the two alignments, and give the read to whichever genome fits it better. Where the fit is decisively better on one side, the answer is a fact about the molecule. Where the two scores are within a point or two of each other, the answer is a fact about the rule.

This is not a marginal technical wrinkle. It is the step that determines what a xenograft transcriptome contains, and the same structure recurs anywhere two genomes share a sample: dual RNA-seq of an infected tissue, co-culture with feeder cells, cell lines grown in mouse-derived matrix.

How much host is actually there

More than most protocols assume, and it varies enormously between samples. Resected PDX tissue can carry up to 70–80% murine DNA through no error in dissection, because mouse stromal and immune cells infiltrate the graft2. In metastatic models the ceiling is higher still: one group reported PDX RNA-seq from lung, liver and brain metastases running as high as 99% mouse reads3.

Across 21 pediatric solid tumour PDXs profiled by three modalities, the median human read fraction was 98% for whole-exome, 93% for whole-genome and 94% for RNA-seq — exome capture excludes some mouse sequence, which flatters WES relative to the others1. The human fractions in WES and WGS from the same models correlated at R = 0.97, which supports a useful conclusion: contamination is set at sampling, not at sequencing1. It is a property of the tumour and how it was taken, and it will differ between your treated and untreated arms if those arms differ in stromal content.

Why the margin is thin exactly where it matters

Human and mouse orthologous gene products are on average about 85% identical at the protein level3, and the nucleotide picture in coding regions is correspondingly close. The consequences at read scale are large. Sequencing five mouse liver samples with a human exome kit, one group found the resulting reads mappable to 49% of the human genome and to 409 cancer genes, and catalogued 1,207,556 mouse-specific alleles that align to the human reference — of which 467,232, or 38.7%, were highly sensitive to contamination and are a recurring source of false somatic mutations in public datasets2.

Half the human genome is close enough to mouse that a mouse read can land on it. What decides the outcome is not whether the read can align to both, but by how much one alignment beats the other.

Transcriptome data has it worse than genomic data, for a reason worth spelling out. Coding sequence is the most conserved part of both genomes, and RNA-seq samples almost nothing else — so where whole-genome data draws most of its reads from divergent intergenic and intronic sequence that separates cleanly, RNA-seq draws its reads disproportionately from the fraction where the two species agree. The consequence shows up directly in the benchmarks: more reads are classified as ambiguous in RNA-seq than in either exome or genome data from the same models1. The assay whose output is most often read as a quantitative measurement is the one with the largest undecidable fraction.

The rule, stated plainly

Disambiguate compares alignment scores between the two BAM files and, in the event of a tie, uses edit distance; reads that still cannot be separated — or that aligned to neither genome — are labelled ambiguous4. XenofilteR works from edit distance between read and reference6. Bamcmp compares mapping scores between the two alignments10. Xenome takes a different route, classifying reads by 25-mer content into human, mouse, both, neither and ambiguous5; xengsort is a faster modern implementation of that idea7.

The xengsort authors make the point that matters about the alignment-based family: in principle these tools all do the same thing, and large differences between them arise from the alignment parameters rather than from the tool itself7. That is borne out empirically — in the 21-PDX benchmark, Disambiguate and bamcmp both produced better classifications when run on STAR alignments than on HISAT2 alignments1. The aligner is not a substrate for the disambiguation step. It is part of it.

A third strategy sidesteps the comparison by removing it: build one concatenated human-plus-mouse reference and align to it once, letting the aligner’s own mapping-quality machinery resolve the competition internally. That approach performs well — it was among the top methods in an independent benchmark2 — but it does not abolish the decision. It relocates it into the aligner’s multi-mapping and MAPQ rules, where it is even less visible and where the parameters governing it were tuned for a different problem entirely.

None of this means xenograft data are unreliable, and it is worth being precise about how well the tools do work. Disambiguate’s own evaluation on simulated exome data assigned only 0.06% and 0.01% of reads to the wrong species, with a higher but still small ambiguous fraction — and the authors noted plainly that both rates are higher for RNA-seq than for DNA4. In the 21-model benchmark, every protocol tested removed mouse reads effectively, differences between them were usually within 5% of the human reads, and swapping between five different mouse strain references changed classification by less than 0.5%1. These are good tools solving a hard problem well. The argument here is about the residue they leave, and about where that residue lands.

ONE READ IN A CONSERVED EXON read 100 bp aligned to human reference alignment score 148 · edit distance 2 aligned to mouse reference alignment score 146 · edit distance 3 margin = 2 out of 150 Where the margin is this small, the assignment is a policy — and the policies differ. FOUR CONVENTIONS, FOUR MATRICES Keep with the graft The read counts toward human expression. Inflates conserved genes in proportion to host content. Assign to the host The read is filtered out as contamination. Deflates conserved genes even in a pure sample. Hold out as ambiguous The read goes to a third pile that is reported. Undercounts, visibly — the loss is measurable. Discard silently The read never reaches the count matrix. Undercounts, invisibly — no record in the output. All four are defensible. None is wrong. Applied to the same FASTQ they produce four different expression matrices, and the genes that differ between them are the same genes every time.
Figure 1. Schematic; the scores are illustrative. The four policies shown are all in use across current tools and pipelines — and the choice between them is usually made implicitly, by installing a tool rather than by deciding anything.

What that does to a single gene

The Xenome paper contains the clearest demonstration I know of, and it is worth reading carefully because it runs in the direction people do not expect. In a pure human cDNA dataset — where every read is, by construction, human — the authors compared their k-mer classifier against an alignment-based analysis on the gene MYH3. The alignment-based approach classified 29 reads as human. Xenome classified 2,713. Every read in the gap had been placed in the “both” category by the alignment method, because MYH3 is highly conserved between the two species and the reads aligned equally well to each5.

Two orders of magnitude, in one gene, in a sample containing no mouse at all. Whatever expression value MYH3 carried into the downstream analysis was a property of the sorting method, not of the transcriptome.

That the effect scales with conservation has been measured directly. In a 2018 preprint benchmarking several approaches on simulated xenograft RNA-seq stratified by orthologue similarity, naive mapping to the human genome and Disambiguate each misassigned more than 30% of murine reads to human in genes above 95% sequence similarity; the rate fell to around 12% and then below 10% in less homologous orthologues, while Xenome misclassified fewer than one read per thousand even in the most homologous set9. Those figures come from simulated data in a preprint and should be treated as indicative rather than exact. The gradient is the finding: error rate is a function of conservation.

This is not noise. It is a bias whose magnitude is set by each gene’s evolutionary history, so it lands on the same genes in every sample, in every replicate, in the same direction.

A random error averages out across replicates. This one does not. It reproduces.

The aggregate picture is reassuring and misleading at the same time

Compare expression with and without host-read removal across 21 PDX models and the Spearman correlations run from 0.980 to 0.999, median 0.9971. On that evidence you would conclude the problem is negligible. It is negligible for most genes, and that is precisely why it is dangerous.

The same study identified 821 genes that behaved as regression outliers in at least two samples — genes whose expression was inflated by unfiltered mouse reads. They were significantly enriched in stromal signatures (myogenesis, pancreas beta cells), immune signatures (allograft rejection, inflammatory response), and, more surprisingly, KRAS signalling. Signature scores for those sets were almost uniformly higher in the unsorted samples, and the number of affected genes tracked the contamination level across samples1. Meanwhile the genes most resistant to contamination were the ones you would nominate as cancer-cell-intrinsic: cell division, TP53 signalling, DNA repair, metabolism1.

So the artefact sorts itself into the tumour microenvironment story. A PDX study reporting an immune or stromal expression programme, or a difference in one between arms, is reporting on exactly the gene sets that host reads inflate — and doing so at a global correlation of 0.997.

Which analyses care, and which don’t

Effect of unremoved host reads, by analysis
AnalysisEffectObserved
Expression, genome-wide Minimal Spearman 0.980–0.999 between filtered and unfiltered1
Expression, stromal and immune gene sets Systematically inflated 821 outlier genes enriched in immune and stromal signatures1
Somatic mutation calling Severely inflated Mean 1,127 exonic mutations unfiltered vs 38 filtered1; 419 true SNVs became 107,826 in an in silico 3:1 mix6
Mutational signatures Distorted SBS46, an artefact signature, at exposure 0 filtered vs 0.18 unfiltered1
Gene fusions Minimal No significant difference in fusion counts1
Copy number Minimal Profiles nearly identical in all but one of 21 models1

Most figures are from a single 21-model pediatric solid tumour cohort and one in silico mixing experiment; magnitudes will differ with tumour type, host content and library chemistry. The pattern — aggregate and structural measures robust, per-gene and per-variant measures fragile — is the part that generalises.

The mutation result has a useful corollary. Because the artefacts arise from a shared source, they recur: pairwise overlap of artefact calls between genetically unrelated PDX models averaged 12.6% by Jaccard index, against 0.63% for genuine mutations1. Recurrence across unrelated models is therefore a usable filter, and blacklists of human-genome-aligned mouse alleles exist for the purpose2.

It is also worth stating what the filtering tools cost. In simulated benchmarks from a PDX consortium, all five tools tested exceeded 99% precision, but recall varied: XenofilteR recovered 96.60% of true human reads in whole-exome data and 89.63% in RNA-seq, while BBSplit achieved the highest precision without a recall penalty8. A recall of 89.63% means one in ten genuine human reads discarded — and, per the conservation gradient above, not one in ten from each gene.

The same problem, opposite default

Dual RNA-seq of an infected host faces an identical decision with a different convention. The standard pipeline aligns reads in parallel against host and pathogen references, quantifies the reads that map equally well to both — the field calls them cross-mappings — and discards them, using only unequivocally assigned reads for quantification11. Discarding is the sane default when the two organisms are a mammal and a bacterium, because genuine ambiguity is rare and the discarded pile is small.

But the rate is not a constant. As reviewers of plant–parasite dual RNA-seq put it, cross-mapping arises from insufficient divergence between the interacting organisms, and the proportion of cross-mapped reads depends on how far apart they are evolutionarily12. Human and bacterium, human and mouse, and plant and parasitic plant are three very different regimes running through pipelines that look the same in a methods section. Two adjacent fields have settled on opposite defaults for the same ambiguous read — xenograft tools mostly keep or reassign it, dual RNA-seq discards it — and neither default is derived from the divergence of the particular pair of organisms in the experiment.

The fix the field already described

There is an obvious better answer, and it was written down in 2012. The Xenome authors noted that where reads are ambiguous, the unambiguous reads belonging to the same homologous gene could be used to inform the classification, and pointed at expectation-maximisation methods as the basis for such an extension — while being clear that they had not implemented it5.

That is the same machinery transcript quantifiers already use for multi-mapping reads within a single genome: rather than forcing each ambiguous read to one location, distribute it probabilistically according to the abundance implied by the reads that could be assigned. Cross-species assignment is the multi-mapping problem with a species label attached, and the field solved the within-genome version more than a decade ago. The cross-species version is still largely resolved by a hard comparison of two numbers, one read at a time, with no memory of what the rest of the gene looked like.

What to do about it

  1. Treat host fraction as a per-sample covariate. It varies from a few percent to nearly everything, it is set at sampling, and it will correlate with treatment arm whenever treatment changes stromal content. Report it with every sample and check whether it differs between groups before interpreting any expression difference.
  2. Report the ambiguous pile, not just the two species. The count of reads that could not be separated is the direct measure of how much of the result is policy. RNA-seq generates more of them than DNA does1.
  3. Record the aligner and its parameters as part of the disambiguation method. “Filtered with Disambiguate” is not a method description; the same tool on STAR and on HISAT2 gives different answers17.
  4. Never compare expression across samples processed with different tools or aligners. Tool identity behaves as a batch variable, and it acts on a specific, reproducible set of genes rather than diffusely.
  5. Run a second tool as a sensitivity analysis and report the genes that move. This is cheap, and the list of movers is informative in itself — it is the conserved fraction of your result.
  6. Apply extra scrutiny to immune, stromal and microenvironment findings in xenograft data. These are the gene sets host reads inflate. A finding there needs the host fraction, the tool, and ideally an orthogonal confirmation before it is a finding.
  7. Use cross-model recurrence and mouse-allele blacklists for mutation artefacts, and expect artefact calls to sit at low variant allele fraction1.
  8. In single-cell data, classify cells rather than reads. Per-cell host read fraction is strongly bimodal — in one PDX sample only 0.6% of cells fell in the ambiguous middle — which makes cell-level assignment far better posed than read-level assignment1.
  9. In dual RNA-seq, state the cross-mapping policy and justify it for your species pair, rather than inheriting a default calibrated on a much more divergent pair.

The shape of the error

The recurring structure in this series is a measurement that is real and an inference attached to it that is not. This case adds a twist: the ambiguity here is genuine. When a read from a conserved exon aligns equally well to both genomes, no amount of care recovers which nucleus it came from, because the information required to answer was never in the molecule. The read is 100 bases long and those 100 bases exist in both species.

So the tool cannot be blamed for guessing. What can be examined is what happens to the guess afterwards: it enters a count matrix as an integer, indistinguishable from a read that was assigned with a margin of forty. The matrix has one column per sample and one row per gene, and nowhere to record that some of its entries were decided by two points of alignment score.

A xenograft expression value is a claim about a human tumour. It is produced by a procedure that also had to decide, thousands of times, which of two genomes a read belonged to — and the decisions it found hardest were the ones about the genes the two species share.

References

  1. Bhandari M, He F, Rogojina A, et al. Benchmarking mouse contamination removing protocols in patient-derived xenografts genomic profiling. npj Precision Oncology 2025;9:113. nature.com/articles/s41698-025-00902-z
  2. Jo SY, Kim E, Kim S. Impact of mouse contamination in genomic profiling of patient-derived models and best practice for robust analysis. Genome Biology 2019;20:231. genomebiology.biomedcentral.com — s13059-019-1849-2
  3. Chromatin conformation capture (Hi-C) sequencing of patient-derived xenografts: analysis guidelines. bioRxiv 2020 (preprint; not peer reviewed). Cited here for its summary of published host-content and orthologue-identity figures. biorxiv.org — 2020.10.17.343814
  4. Ahdesmäki MJ, Gray SR, Johnson JH, Lai Z. Disambiguate: an open-source application for disambiguating two species in next generation sequencing data from grafted samples. F1000Research 2016;5:2741. f1000research.com/articles/5-2741
  5. Conway T, Wazny J, Bromage A, et al. Xenome — a tool for classifying reads from xenograft samples. Bioinformatics 2012;28(12):i172–i178. academic.oup.com/bioinformatics/28/12/i172
  6. Kluin RJC, Kemper K, Kuilman T, et al. XenofilteR: computational deconvolution of mouse and human reads in tumor xenograft sequence data. BMC Bioinformatics 2018;19:366. bmcbioinformatics.biomedcentral.com — s12859-018-2353-5
  7. Zentgraf J, Rahmann S. Fast lightweight accurate xenograft sorting. Algorithms for Molecular Biology 2021;16:2. link.springer.com — s13015-021-00181-w
  8. Systematic establishment of robustness and standards in patient-derived xenograft experiments and analysis. Cancer Research 2020;80(11):2286–2297. aacrjournals.org — Cancer Res 80(11):2286
  9. SMAP: exploiting high-throughput sequencing data of patient-derived xenografts. bioRxiv 2018 (preprint; not peer reviewed). biorxiv.org — 440008
  10. Khandelwal G, Girotti MR, Smowton C, et al. Next-generation sequencing analysis and algorithms for PDX and CDX models. Molecular Cancer Research 2017;15(8):1012–1016. pubmed.ncbi.nlm.nih.gov/28442585
  11. Westermann AJ, Barquist L, Vogel J. Resolving host–pathogen interactions by dual RNA-seq. PLoS Pathogens 2017;13(2):e1006033. pubmed.ncbi.nlm.nih.gov/28207848
  12. Refining dual RNA-seq mapping: sequential and combined approaches in host–parasitic plant dynamics. Frontiers in Plant Science 2024;15:1483717. frontiersin.org — fpls.2024.1483717
Zetobit, LLC · Bioinformatics consulting · Lexington, KY · zetobit.com
Next
Next

cfDNA Fragmentomics