!!!!! Good candidate for the workshop for the clinicians !!!!!

“First, I generate ps_rarefied from ps_filt. Then, for cleaner composition plots, I create ps_abund / ps_abund_rel by filtering taxa for plotting (e.g., keeping taxa with mean relative abundance > 0.1%; resulting in ~95 taxa across 239 samples). Is it correct to compute alpha diversity and beta diversity using ps_abund (the taxa-filtered object)?”

Not really. ps_abund / ps_abund_rel (taxa filtered “for plotting”) is generally not the right object for alpha- or beta-diversity, unless you explicitly want diversity calculated only within that filtered subset. Luckily, I directly using alpha- and beta diversity from qiime2-result and MicrobiotaProcess –> which means we should feed MicrobiotaProcess with ps_filt (not rarefied) and for differential analysis also using ps_filt (not rarefied)

For that sentence, I mean these specific objects:

• For Jaccard / unweighted UniFrac (presence–absence / detection-sensitive):
  ✅ compute beta diversity on ps_rarefied
  = the rarefied-count phyloseq object created from ps_filt (with the full taxon set, i.e., no plotting-only taxon filter).

• For Bray–Curtis (abundance-based):
  ✅ compute beta diversity on a non-rarefied transformed object, typically:

  • ps_rel = relative abundances computed from ps_filt (not rarefied), or
  • a Hellinger-transformed version derived from ps_filt (if that is your chosen workflow).

• For CLR / Aitchison distance:
  ✅ compute beta diversity on a CLR-transformed version of the non-rarefied counts from ps_filt
  (often named ps_clr / otu_clr, created from ps_filt + pseudocount), not on ps_abund.

The “plotting-filtered object” I’m warning against is:

• ps_abund / ps_abund_rel = taxa filtered for visualization (e.g., mean relative abundance > 0.1%).

Why it’s not “correct” for diversity

Filtering taxa by mean relative abundance (e.g., >0.1%) removes many low-abundance / low-prevalence taxa and changes the community profile.

Alpha diversity

• Observed richness will be artificially lower (because taxa were removed).
• Shannon/Simpson can also change because you altered the abundance distribution. ✅ Therefore, compute alpha diversity on ps_rarefied (unfiltered taxa), not on ps_abund.

Beta diversity

This depends on the distance metric, but heavy taxon filtering can distort results:

• Presence/absence distances (e.g., Jaccard, unweighted UniFrac) are very sensitive to removing taxa → filtering can strongly change distances and clustering.
• Abundance-weighted distances (e.g., Bray–Curtis, weighted UniFrac) can also shift after filtering; mild filters (e.g., remove singletons) may be acceptable, but mean rel. abundance >0.1% can be quite aggressive (95 taxa suggests strong pruning). ✅ Best practice: compute beta diversity on ps_rarefied (for Jaccard/unweighted UniFrac) or on non-rarefied transformed data (for Bray/CLR workflows), not on a plotting-filtered object.

Recommended workflow (clean separation of purposes)

For analysis (diversity)

• Alpha diversity: ps_rarefied (derived from ps_filt)
• Beta diversity:
  • Jaccard / unweighted UniFrac: ps_rarefied
  • Bray–Curtis: usually from ps_rel (or Hellinger), often without rarefaction
  • Aitchison/CLR: CLR-transformed (non-rarefied) data (no rarefaction)

For visualization (composition plots only)

• Optional cleaner composition plots: ps_abund / ps_abund_rel derived from ps_rel (taxa filtered to make stacked bars/heatmaps readable)

When would ps_abund be acceptable for beta diversity?

Only if your question is explicitly:

“How do samples differ considering only dominant taxa (≥0.1% mean abundance)?”

Then it’s a valid different analysis, but you should label it clearly as “dominant-taxa-only beta diversity.”

If you tell me which beta distance(s) you use (Bray, Jaccard, weighted/unweighted UniFrac, CLR/Aitchison), I can recommend the best object + transformation for your exact setup.

---

The correct answer is: beta diversity does not universally require rarefaction.

Why this confusion happens

Historically, many microbiome pipelines used rarefied counts for beta diversity to “standardize depth.” That’s why you often see “beta diversity → rarefy” in older workflows. But in modern practice, rarefaction for beta diversity is optional and sometimes not recommended, because it discards data and adds randomness.

When rarefaction is a good idea for beta diversity

Use ps_rarefied when your beta metric is very sensitive to library size, especially:

• Presence/absence distances: Jaccard, unweighted UniFrac These are strongly influenced by detection of rare taxa, which depends on depth. Rarefaction is a common way to make them comparable.
• If your field/reviewer expectations specifically require rarefied beta diversity for certain metrics.

When rarefaction is not necessary (often better avoided)

Prefer non-rarefied approaches when you use:

• Bray–Curtis on relative abundance (or Hellinger-transformed data) Rarefaction can reduce power and add noise; relative abundance already removes the “sum” effect (though compositionality remains).
• Aitchison distance / PCA via CLR transform (Euclidean on CLR) This is a standard compositional approach and typically does not use rarefaction.
• DESeq2 VST / rlog for ordination-like PCA (count-model-based normalization) Again: no rarefaction.

Important note about your original code

Doing rarefaction → compositional (“relative abundance”) is usually redundant for beta diversity/ordination, because after rarefaction every sample has the same total, and compositional transformation just divides by that constant.

Corrected “rule of thumb”

• Alpha diversity: ps_rarefied ✅ (common)
• Beta diversity:
  • Unweighted UniFrac / Jaccard: ps_rarefied ✅ (often recommended)
  • Bray–Curtis / ordination on abundances: ps_rel or Hellinger ✅ (rarefaction optional)
  • Aitchison (CLR): CLR-transformed (non-rarefied) ✅ (no rarefaction)

If you tell me exactly which beta-diversity distance(s) you compute (Bray–Curtis? weighted/unweighted UniFrac? Jaccard? Aitchison/CLR?), I can recommend the cleanest “use this object + this transformation” setup.

---

Doing rarefaction and then compositional normalization (relative abundance) is often redundant.

Short answer: Don’t need both — rarefaction + compositional normalization is redundant in most cases. Doing rarefaction and then compositional normalization (relative abundance) is often redundant.

# RAREFACTION
set.seed(9242)  # This will help in reproducing the filtering and nomalisation.
ps.ng.tax <- rarefy_even_depth(ps.ng.tax, sample.size = 6389)
total <- 6389

# NORMALIZE number of reads in each sample using median sequencing depth.
total = median(sample_sums(ps.ng.tax))
#> total
#[1] 42369
standf = function(x, t=total) round(t * (x / sum(x)))
ps.ng.tax = transform_sample_counts(ps.ng.tax, standf)
ps.ng.tax_rel <- microbiome::transform(ps.ng.tax, "compositional")

saveRDS(ps.ng.tax, "./ps.ng.tax.rds")

What happens mathematically

1) Rarefaction forces every sample to the same library size (e.g., 6,389 reads). 2) Compositional normalization divides each sample by its total so it sums to 1.

After rarefaction, every sample has the same total, so compositional normalization is essentially: relative_abundance = rarefied_counts / constant_depth. That step is valid, but typically unnecessary.

Practical implication

• If your goal is relative abundance–based visualization or many beta-diversity analyses, you can usually skip rarefaction and compute relative abundances (or other transformations) on the full data.
• If your goal is alpha diversity, rarefaction is commonly used, but then you typically keep that rarefied object for alpha diversity only (not as a general-purpose normalized dataset).

So: rarefaction alone can be fine for specific tasks (especially alpha diversity), but rarefaction + compositional normalization together is usually not needed.

---

Mapping new ps_* objects to the old ps.ng.tax* objects (and overwrite notes)

ps_* object	What it contains	Old R-script equivalent (ps.ng.tax*)	Overwrite notes for old object(s)	Taxonomic composition	Alpha diversity	Beta diversity	DESeq2 differential abundance
ps_raw	Raw imported phyloseq object (integer counts; as imported)	(no direct equivalent)	—	❌	❌	❌	⚠️ only if you also apply sample filtering first
ps_base	ps_raw + taxonomy + sample metadata aligned (still raw counts)	(closest to) ps.ng.tax before overwrite	Old ps.ng.tax BEFORE overwrite: integer count table (absolute abundance).	❌	❌	❌	⚠️ only if you also apply sample filtering first
ps_pruned	Optional subset of ps_base (e.g., remove unwanted samples by ID/pattern); still raw counts	(subset of) ps.ng.tax before overwrite	Old ps.ng.tax BEFORE overwrite: same as above, but without any sample subsetting unless you did it.	❌	❌	❌	⚠️ only if you also apply low-depth filtering
ps_filt	Filtered samples (low-depth samples removed) + taxa with nonzero totals; absolute counts	ps.ng.tax before overwrite (but with explicit low-depth filtering)	Old ps.ng.tax BEFORE overwrite: raw integer counts; new ps_filt is the “cleaned” version after sample-depth QC.	✅ as a starting point (but plot on ps_rel)	✅ as input to rarefaction	✅ as input to rarefaction	✅ as input to ps_deseq
ps_rel	Relative abundance (compositional) computed from ps_filt	ps.ng.tax_rel (conceptually)	In the old script, ps.ng.tax_rel is “relative abundance of ps.ng.tax”, but its meaning depends on whether it was computed before or after ps.ng.tax was overwritten.	✅ primary	❌	❌	❌
ps_abund	Absolute counts after “plotting taxa filter” (e.g., mean rel. abundance > 0.1%), derived from ps_filt via ps_rel	ps.ng.tax_abund	Old ps.ng.tax_abund was created by filtering taxa using mean relative abundance (> 0.001) and then pruning counts.	✅ (if you want cleaner plots)	❌	❌	❌ (not recommended)
ps_abund_rel	Relative abundance computed from ps_abund (filtered taxa set)	ps.ng.tax_abund_rel	Old ps.ng.tax_abund_rel was relative abundance of the filtered-taxa object.	✅ (clean composition plots)	❌	❌	❌
ps_rarefied	Rarefied counts from ps_filt (even depth)	ps.ng.tax after overwrite	Old ps.ng.tax was overwritten 1× by rarefy_even_depth(ps.ng.tax, sample.size = 41764) → after this, ps.ng.tax no longer meant raw counts; it meant rarefied counts.	❌	✅ primary	✅ primary	❌
ps_deseq	Non-rarefied integer counts from ps_filt + optional count-based taxon prefilter (e.g., total ≥ 10)	(no direct equivalent)	Old ps.ng.tax_abund is not a good DESeq2 analogue because it used a mean-relative-abundance filter; ps_deseq uses count-based prefiltering (optional) and keeps integer counts.	❌	❌	❌	✅ primary

Overwrite summary (old script):

• ps.ng.tax was overwritten 1 time:
  • Before overwrite: ps.ng.tax = absolute abundance (raw integer counts).
  • After overwrite: ps.ng.tax = rarefied counts, produced by rarefy_even_depth(ps.ng.tax, sample.size = 41764).

ps_* object	What it contains	Taxonomic composition	Alpha diversity	Beta diversity	DESeq2 differential abundance
ps_raw	Raw imported phyloseq object (integer counts; as imported)	❌ (not recommended directly)	❌	❌	⚠️ only if you also apply sample filtering first
ps_base	ps_raw + taxonomy + sample metadata aligned (still raw counts)	❌	❌	❌	⚠️ only if you also apply sample filtering first
ps_pruned	Optional subset of ps_base (e.g., remove unwanted samples by ID/pattern); still raw counts	❌	❌	❌	⚠️ only if you also apply low-depth filtering
ps_filt	Filtered samples (low-depth samples removed) + taxa with nonzero totals; absolute counts	✅ as a starting point (but plot on ps_rel)	✅ as the input to rarefaction	✅ as the input to rarefaction	✅ as the input to ps_deseq
ps_rel	Relative abundance (compositional) computed from ps_filt	✅ primary	❌	❌	❌
ps_abund	Absolute counts after “plotting taxa filter” (e.g., mean rel. abundance > 0.1%), derived from ps_filt via ps_rel	✅ (if you want cleaner plots)	❌	❌	❌ (not recommended)
ps_abund_rel	Relative abundance computed from ps_abund (filtered taxa set)	✅ (clean composition plots)	❌	❌	❌
ps_rarefied	Rarefied counts from ps_filt (even depth)	❌	✅ primary	✅ primary	❌
ps_deseq	Non-rarefied integer counts from ps_filt + optional count-based taxon prefilter (e.g., total ≥ 10)	❌	❌	❌	✅ primary