Metabarcoding Community Analysis
Faiz Fadhlirobby
2026-07-23
1 Introduction
Metabarcoding is a molecular technique that uses targeted DNA sequencing to identify multiple species simultaneously within a complex environmental or bulk sample. By extracting genetic material from sources like soil or water and comparing it to reference databases, it allows researchers to rapidly profile hidden biodiversity without capturing individual organisms.
Microbial communities are essential because they maintain the fundamental life support systems of our planet, driving nutrient cycles, supporting human health, and sustaining agriculture. They do not live in isolation but form complex networks that perform vital ecological and biological functions that single organisms cannot achieve alone.
Within the human body, one of the most critical of these networks is
the gut microbiome, which is profoundly shaped by long term nutritional
habits. However, the microbial response to food is rarely universal. A
diet that sustains a healthy microbial balance in one population may
cause severe dysbiosis in another, as an individual’s geographical
origin, cultural lifestyle, and established microbial baseline heavily
dictate their biological suitability to different foods. To explore this
dynamic interaction between origin and dietary tolerance, this analysis
utilizes metabarcoding data from the dietswap study
(O’Keefe et al.). By examining the gut communities of rural Africans and
African Americans before and after exchanging their traditional
high-fiber and Western diets, we can mathematically quantify how deeply
our population-level characteristics and immediate nutritional shifts
interact to alter human microbial ecology.
2 Import Library
Here, the researchers utilized data from R packages:
microbiome and phyloseq for microbiome
analysis, tidyverse for data tidying, and finally
vegan to calculate various diversity indices.
library(microbiome)
library(phyloseq)
library(tidyverse)
library(vegan)Once all the libraries have been imported, you simply need to call
the dietswap data.
data("dietswap")3 Data Explanation
The dietswap data is contained within the
microbiome package and has likely already been formatted as
a phyloseq object. This data structure generally integrates
three core elements: a matrix of microbial taxa abundance per sample
(otu_table), variables characterizing the subjects (sample_data), and a
taxonomic reference classifying each taxon into broader biological
hierarchical levels (tax_table).
dietswap## phyloseq-class experiment-level object
## otu_table() OTU Table: [ 130 taxa and 222 samples ]
## sample_data() Sample Data: [ 222 samples by 8 sample variables ]
## tax_table() Taxonomy Table: [ 130 taxa by 3 taxonomic ranks ]The dietswap dataset comprises 130 distinct taxonomic
entries at specific levels, based on 222 samples that vary in both their
characteristics and ecological attributes.
We have a primary tax_otu table containing 130 rows of
taxa and 222 columns of samples, specific to the treatment types. Each
value represents the abundance of a given taxon within a sample.
otu_table(dietswap)[1:5,1:10]## OTU Table: [5 taxa and 10 samples]
## taxa are rows
## Sample-1 Sample-2 Sample-3 Sample-4 Sample-5
## Actinomycetaceae 0 1 0 1 0
## Aerococcus 0 0 0 0 0
## Aeromonas 0 0 0 0 0
## Akkermansia 18 97 67 256 21
## Alcaligenes faecalis et rel. 1 2 3 2 2
## Sample-6 Sample-7 Sample-8 Sample-9 Sample-10
## Actinomycetaceae 0 0 0 0 0
## Aerococcus 0 0 0 0 0
## Aeromonas 0 0 0 0 0
## Akkermansia 16 26 30 19 125
## Alcaligenes faecalis et rel. 2 2 2 2 7Metadata contains detailed information about the samples. We can determine what kind of research can be developed based on this metadata.
head(meta(dietswap))## subject sex nationality group sample timepoint
## Sample-1 byn male AAM DI Sample-1 4
## Sample-2 nms male AFR HE Sample-2 2
## Sample-3 olt male AFR HE Sample-3 2
## Sample-4 pku female AFR HE Sample-4 2
## Sample-5 qjy female AFR HE Sample-5 2
## Sample-6 riv female AFR HE Sample-6 2
## timepoint.within.group bmi_group
## Sample-1 1 obese
## Sample-2 1 lean
## Sample-3 1 overweight
## Sample-4 1 obese
## Sample-5 1 overweight
## Sample-6 1 obesemetadata <- meta(dietswap)
metadata <- metadata %>%
mutate(timepoint = as.factor(timepoint),
timepoint.within.group = as.factor(timepoint.within.group))
summary(metadata)## subject sex nationality group sample timepoint
## azh : 6 female:102 AAM:123 DI:72 Length:222 1:38
## azl : 6 male :120 AFR: 99 ED:75 Class :character 2:37
## byn : 6 HE:75 Mode :character 3:38
## cxj : 6 4:37
## dwc : 6 5:35
## eve : 6 6:37
## (Other):186
## timepoint.within.group bmi_group
## 1:112 lean :56
## 2:110 overweight:76
## obese :90
##
##
##
## The following are the definitions for each column:
subject: Anonymized alphabetic identifier representing the individual participant.
sex: The biological sex of the participant.
nationality: The geographic and cultural cohort of the participant, designated as either native rural South Africans (AFR) or African Americans (AAM).
group: The specific clinical phase during which the specimen was collected. Categories include Home Environment (HE) for baseline dietary habits, Endoscopy (ED) for clinical pre-colonoscopy baselines, and Dietary Intervention (DI) for the active diet-swap phase.
sample: The unique alphanumeric tag assigned to each biological specimen.
timepoint: The chronological sequence of the sample collection across the entire duration of the study, scaled from 1 to 6.
timepoint.within.group: The relative sub-sampling sequence (coded as 1 or 2) denoting early or late collection within a specific clinical phase.
bmi_group: A categorical classification of the participant’s Body Mass Index, stratified into lean, overweight, or obese.
4 Relative Abundance Analysis
5 Converting to Relative Abundance
In metabarcoding, each fecal or gut sample typically yields a different total number of DNA reads (sequencing depth).
For example:
A sample from Subject 1 might yield a total of 10,000 bacterial reads.
A sample from Subject 2 might yield a total of 50,000 bacterial reads.
If you compare the absolute numbers directly, the comparison becomes biased and unfair (Subject 2 would appear to have more bacteria for any given taxon simply because the machine sequenced more DNA from that sample).
By converting the data to relative abundance, you standardize all samples to the same scale (where the sample total equals 1 or 100%).
DS.rel <- transform_sample_counts(
dietswap,
function(x) x / sum(x)
)5.1 Aggregation to the Phylum Level
The taxa that have the highest hierarchy in this data are at the phylum level. To make it easier to see abundance, we will see abundance from the phylum level
DS.phylum <- tax_glom(
DS.rel,
taxrank = "Phylum"
)
top10 <- names(sort(
taxa_sums(DS.phylum),
decreasing = TRUE
))[1:10]
DS.top10 <- prune_taxa(
top10,
DS.phylum
)
p <- plot_bar(DS.top10, fill = "Phylum")
p +
geom_bar(aes(color = Phylum, fill = Phylum), stat = "identity", position = "stack") +
facet_wrap(~ group, scales = "free_x") +
theme_minimal() +
labs(
x = "Samples",
y = "Relative Abundance"
) +
theme(
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
legend.position = "bottom",
strip.text = element_text(face = "bold", size = 12),
panel.grid.major.x = element_blank(),
panel.grid.minor = element_blank()
)
Based on the plot above, across the DI, ED, and HE diets, only the phyla Firmicutes and Bacteroidetes constitute the majority of the gut microbiota; they are consistently abundant in all diets.
6 Visualisasi Relative Abundance by Group
This plot shows the average abundance across all samples within each diet group.
df <- psmelt(DS.top10)
df.mean <- df %>%
group_by(
group,
Phylum
) %>%
summarise(
MeanAbundance = mean(Abundance),
.groups = "drop"
)
ggplot(df.mean, aes(x = group, y = MeanAbundance, fill = Phylum)) +
geom_bar(stat = "identity") +theme_bw() +
theme(axis.text.x =element_text(angle=45,hjust = 1))
While the previous abundance plot showed no apparent differences in abundance across groups, these group-specific plots reveal that the DI and HE diets resulted in high levels of Firmicutes, whereas the ED diet yielded a greater abundance of Bacteroidetes. Bacterial populations other than these two were minimal and constituted only a very small minority.
7 Alpha Diversity Analysis
When considering diversity indices, three popular ones stand out:
Observed: The number of ASVs found in each habitat.
Shannon :The Shannon index relates to the level of species diversity within a community. A high Shannon value indicates a complex community, taking into account both the number of ASVs and their evenness.
Simpson: The Simpson index measures dominance. A community dominated by a single taxon will yield a different value compared to a community with even distribution.
plot_richness(dietswap,
measures = c("Observed", "Shannon", "Simpson"),
x = "group",
color = "nationality") +
# Adds a boxplot underneath the points to show medians and quartiles
geom_boxplot(alpha = 0.7, outlier.shape = NA) +
# A clean theme
theme_bw() +
# Professional labels
labs(x = "Study Phase",
color = "Nationality")
From the plot, In Observed Diversity, the rural Africans (AFR, blue) have a slightly higher median number of unique taxa compared to the African Americans (AAM, red) across all phases. However, the boxes overlap significantly, meaning the raw number of different bacteria in their guts is relatively similar. In Shannon and Simpson, we see a massive, structural difference. The AAM group (red) has notably higher Shannon and Simpson scores than the AFR group (blue) in every single study phase. This is driven by the heavy dominance of specific fiber degrading taxa (like Prevotella) in the AFR cohort, which lowers their overall evenness score despite high richness.
8 Beta Diversity Analysis
Beta diversity analysis is a measuring the degree of difference or similarity in species composition between two or more environments or samples. One of the most frequently used metrics is Bray-Curtis.
9 Principal Coordinate Analysis (PCoA)
Since beta diversity utilizes a distance matrix, PCoA is used to help interpret beta diversity in two dimensions.
ordination <- ordinate(dietswap,method = "PCoA",distance = "bray")
plot_ordination(dietswap, ordination, color = "group", shape = "nationality") +
theme_bw() +
geom_point(size = 3) 
Based on the plot above, PC1 suggests that nationality has a significant impact on similarity. This can be observed in the positioning of the circles (AAM) in the upper corner, while the triangles (AFR) are located on the left.
10 Non-metric MultiDimensional Scaling (NMDS)
As an alternative to PCoA which entails various assumptions NMDS offers a non-metric approach by ranking data from most similar to most different.
ordu.nmds <- ordinate(dietswap,method = "NMDS",distance = "bray")plot_ordination(dietswap, ordu.nmds, color = "group", shape = "nationality") +
theme_bw() +
geom_point(size = 3, alpha = 0.8)
In the PCoA, there is extensive overlap across both groups and nationalities; in contrast, NMDS effectively distinguishes by nationality, with only 2 to 5 samples deviating from their original nationality.
11 PERMANOVA
Unlike classical MANOVA, PERMANOVA utilizes distance matrices (such as Bray-Curtis or Euclidean distance) and permutation algorithms to assess significance. Its key advantages include being free from distributional assumptions, making it highly suitable for skewed, zero-inflated, or qualitative data.
meta <- data.frame(sample_data(dietswap))
bc.dist <- phyloseq::distance(
dietswap,
method = "bray"
)
adonis2(bc.dist ~ group, data = meta)## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = bc.dist ~ group, data = meta)
## Df SumOfSqs R2 F Pr(>F)
## Model 2 1.036 0.02629 2.9563 0.006 **
## Residual 219 38.386 0.97371
## Total 221 39.422 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1From this test, it was proven (0.06 P-Value) that the type of diet affects the bacterial ecosystem in the intestine. But it only applies to 2.6% of popularity and that is very small. So we need to try to combine the type of diet and the type of race
adonis2(bc.dist ~ nationality * group, data = meta)## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = bc.dist ~ nationality * group, data = meta)
## Df SumOfSqs R2 F Pr(>F)
## Model 5 8.478 0.21505 11.835 0.001 ***
## Residual 216 30.944 0.78495
## Total 221 39.422 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1A significant shift occurs when integrating diet type with ethnicity, with R² values ranging from 2% to 21%. This indicates that 21% of the variation in gut bacterial populations is heavily dependent on the combination of diet and ethnicity. Consequently, dietary approaches should be tailored to specific ethnic groups to achieve optimal results.
12 Dominant Taxa
First, we must select the 20 most dominant taxa at the genus level.
# 1. Agglomerate the relative abundance data to the Genus rank
DS.genus <- tax_glom(DS.rel, taxrank = "Genus")
# 2. Identify the top 20 Genera
top20 <- names(sort(taxa_sums(DS.genus), decreasing = TRUE))[1:20]
# 3. Prune the phyloseq object to only keep those top 20
DS.top20 <- prune_taxa(top20, DS.genus)
# 4. CRITICAL STEP: Melt the new Genus-level object into 'df'
df <- psmelt(DS.top20)To answer which taxa drive the differences between habitats, we can use a heatmap like the one below.
# 1. Update the summarization to include nationality
df.mean <- df %>%
group_by(nationality, group, Genus) %>%
summarise(MeanAbundance = mean(Abundance), .groups = "drop") %>%
# Create a combined column for the x-axis (e.g., "AAM_HE", "AFR_DI")
mutate(Cohort = paste(nationality, group, sep = "_")) %>%
# Reorder the Genus factor based on overall abundance rather than alphabetically
mutate(Genus = fct_reorder(Genus, MeanAbundance, .fun = sum))
# 2. Plot the improved heatmap
ggplot(df.mean, aes(x = Cohort, y = Genus, fill = MeanAbundance)) +
# Add white borders between the tiles for a cleaner look
geom_tile(color = "white", linewidth = 0.5) +
# Using pseudo_log_trans is safer than log10 if you have exact 0s in your data
scale_fill_viridis_c(trans = scales::pseudo_log_trans(base = 10),
option = "magma", # 'magma' often provides better contrast for heatmaps
name = "Mean Rel.\nAbundance") +
theme_minimal() +
labs(
x = "Cohort (Nationality_Phase)",
y = NULL) +
theme(
# Make x-axis labels bold and readable
axis.text.x = element_text(angle = 45, hjust = 1, face = "bold"),
# Scientifically correct formatting: Italicize Genus names
axis.text.y = element_text(face = "italic"),
# Remove background grid lines since the tiles provide the structure
panel.grid = element_blank()
)
This demonstrates that rural Africans naturally harbor massive colonies of Prevotella in their guts. (Biological note: Prevotella are fiber-degrading bacteria; they thrive on a traditional African diet rich in complex carbohydrates and vegetables.)
Observe what happens when this African group is switched to a Western-style diet (high in meat and fat, low in fiber). The box labeled AFR_DI immediately shifts to pink/dark purple. This means that, in a matter of weeks, their Prevotella populations starved and plummeted drastically!
In contrast, the African-American (AAM) group shows dark purple boxes across all phases (AAM_HE, AAM_ED, AAM_DI). Their guts accustomed to a Western diet from the start do not harbor significant Prevotella populations.
Thank you for reading and Happy Learning! I hope this analysis provides a solid foundation for your future microbiome research. Keep exploring, and best of luck with your data journey!