HUMAnN: Depth-normalization with spike-ins

Downstream analysis and statistics
1

Hello everyone,

We’re using spike-ins (cells added to the sample at a known amount) as the basis for an absolute quantification pipeline on shotgun metagenomics data from stool samples, and using HUMAnN to obtain functional pathway quantification. We’ve seen that read depth varies a lot between samples, and the sample size per sample is also different, leading to differences in DNA input between samples. Before finalizing our approach I wanted to check whether this is already addressed by existing tools or documented recommendations.

So here the specific question: read depth normalization at the community (unstratified) level. We know humann_renorm_table can normalize gene family and pathway values for sequencing depth, but this returns compositional values, each sample’s total is forced to a fixed sum, which discards genuine differences in absolute magnitude between samples. Since we want to preserve those absolute differences using the spike-in as our depth-normalization method instead, is there an existing recommended approach within HUMAnN for this: depth-normalizing community-level pathway abundances while keeping absolute magnitude intact?

We were thinking on the following approach: derive a scaling factor (per sample) from two known spike-in species using RPK measured from a separate alignment against the spike-in references (bowtie2 for example). We’d then use the spike-in free reads and run HUMAnN. And finally, apply the scaling factor to gene family/pathway RPKs to get absolute copy numbers, then divide by stool mass (mg) collected.

Is there any other recommended approach, or prior discussion of best practices, we should be aware of?

Thank you very much in advance, and thank you for your work and for sharing it with the community!

Paula

2

Tagging @franzosa in case he has thoughts on this.

I think an example would look something like this in absolute units (copy numbers and cells):

Sample 1:

  • Taxon A: 3 cells, each with 1 gene 1, 9 gene 2
  • Spike-in: 10 cells, each with 1 gene 3 (I think the number of unique genes and copy number in the spike-in doesn’t actually matter, assuming it’s the same between samples which presumably it would be if it’s coming from the same stock.)

Sample 2:

  • Taxon A: 5 cells, each with 2 gene 1, 8 gene 2
  • Spike-in: 10 cells, each with 1 gene 3 (as before since the spike-in is the same in all samples)

What you want to show is that sample 1 has 3 total copies of Taxon A gene 1 and sample 2 has 10 total copies of Taxon A gene 1. What you actually have to work with are (1) knowledge that the spike-in was equal in absolute value across the samples and (2) the stratified relative abundances:

Sample 1:

  • Taxon A gene 1 (RA_A1): 3 / (3 + 18 + 10)
  • Taxon A gene 2 (RA_A2): 18 / (3 + 18 + 10)
  • Spike-in gene 3 (RA_S3): 10 / (3 + 18 + 10)

Sample 2:

  • Taxon A gene 1 (RA_A1): 10 / (10 + 16 + 10)
  • Taxon A gene 2 (RA_A2): 16 / (10 + 16 + 10)
  • Spike-in gene 3 (RA_S3): 10 / (10 + 16 + 10)

(HUMAnN will produce the RA_A1, RA_A2, and RA_S3 for each sample since those are the stratified abundances.)

Therefore, if you take the ratio of each stratified abundance vs. the spike-in relative abundance, you’ll get RA_A1/RA_S3 = (3 / (3 + 18 + 10)) / (10 / (3 + 18 + 10)) = 3/10 in the first sample and RA_A1/RA_S3 = (10 / (10 + 16 + 10)) / (10 / (10 + 16 + 10)) = 10/10 in the second sample.

Even if you don’t know anything about the spike-in gene’s copy number or absolute abundance, you could correctly conclude that the ratio of the stratified gene copy number in the first to second sample is (3/10) / (10/10) = 3/10. Since you can do the same thing for all your samples, this is enough to do differential abundance analysis or most other analyses. It also has the advantage of not being dependent on which spike-in gene you’re taking the ratios relative to, and in fact it might be better to take the ratio relative to the spike-in’s whole abundance to reduce variation sample-to-sample. However, if you knew more about the gene and the spike-in’s abundance (e.g., here the spike-in gene is single copy and we know the absolute abundance spiked in), you could actually take the 3/10 and 10/10 to conclude you had 3 total copies and 10 total copies of the gene in the 2 samples respectively since we had 10 total copies of the spike-in single copy gene in both.

To sum that up, if you have a gene you know the copy number of in your spike-in and you know the number of cells/CFUs you spiked in, you can use this math to get out the total copy number of the rest of the genes for the rest of the species. If you don’t know the copy number of any spike-in genes or you think those values will have too-high variance (because their abundance is close to 0), you can use the relative abundance of the whole spike-in (e.g. from MetaPhlAn) instead in the ratio. The resulting values will be directly comparable across all samples on the absolute scale; they’ll just be scaled by one global unknown value relative to the true total copy numbers.

Let me know if that makes sense,

Will

1 Like
3

Hi Will, thank you for your answer; that fully makes sense, and it’s exactly the kind of conversion we were planning to apply.

What we were thinking, though, is that there may be an issue then with the pathway abundances output. Let me explain; when we run HUMAnN on reads that still include the spike-in, we do get the stratified per-species gene abundances (that we need for the calculation above). However, as I understand, the spike-in genes also make it into the pathway abundance step. That would mean that the spike-in genes also contribute to the pathway completeness, potentially turning a pathway “complete” when the spike-in fills the missing steps, a false positive pathway.

So my questions are:

  1. Is there a HUMAnN parameter that lets us exclude genes from specific species (i.e. the spike-in) before the pathway abundance step, so gene families are still computed with the spike-in present (for scaling step) but pathway completeness is computed only from the real community?

  2. If not, would this workaround make sense to you:

  • remove spike-in reads from the sample entirely, run HUMAnN on the spike-in free reads to get clean pathway abundances for the real community

  • separately, align spike-in reads with bowtie2 to get their own gene abundances to calculate the scaling factor

    Would then, the gene abundances from the separate bowtie2 alignment be comparable and consistent with HUMAnN gene families’ output?

Thanks again!

4

A few thoughts:

  1. HUMAnN produces both the per-taxon stratified and cross-taxon aggregated gene abundances. I think the math I put earlier would hold on the per-taxon stratified set, which should be uninfluenced by the fact that the spike-in genes are also in the HUMAnN run. However, if you’re looking at the cross-taxon aggregated gene abundances, the presence of the spike-in would indeed be an issue.
  2. Off the top of my head I don’t know of a parameter to exclude specific genes from specific species, but the docs or Eric might have more information.
  3. If you’re still trying to get the cross-taxon aggregated gene abundances without the spike-in, running HUMAnN on the spike-in free reads should work. I think there’s an option to save the per-read mapping file from HUMAnN which could help you deplete these reads too (i.e. run HUMAnN on everything, remove the reads that map to your spike-in, then run HUMAnN again), but you’d have to check the docs on this.
  4. To get the scaling factors, unless you have really low abundance of the spike-in or it’s not in a reference database, it’s probably better to estimate the spike-in abundance with something like MetaPhlAn 4 rather than a direct mapping. First, this means you don’t have to worry about things like gene- and genome-size normalization that would otherwise be issues when getting a per-species relative abundance. Second, we’ve found in internal testing that doing a direct bowtie2 map of the reads to the spike-in genome can substantially overstate the abundance of the spike-in since other species’ reads get recruited. This is much less of an issue in a taxonomic abundance method that competitively maps all the reads against a large database rather than a single genome.

Let me know if that makes sense,

Will