ILLUSTRATIVE EXAMPLE - ALL DATA ARE FICTIONAL AND FOR DEMONSTRATION PURPOSES ONLY

1. Introduction

DEMO CLIENT commissioned Sycomia to carry out the bioinformatic analysis of a 16S rRNA metagenomic dataset. Sycomia was responsible for the analysis of the sequencing data provided by the client, from raw FASTQ files to the results presented in this report. 

In this study, faecal samples from 40 subjects were collected in two groups, Control and Treated (n = 20 each), and the V3-V4 hypervariable regions of the 16S rRNA gene were sequenced on an Illumina MiSeq platform. The objectives of the analysis were:

• to assess the quality of the sequencing data and generate a high-confidence profile of amplicon sequence variants (ASVs).
• to describe the taxonomic composition of the faecal microbiota in both groups.
• to determine whether the intervention affects microbial diversity within samples (alpha diversity) and community structure between samples (beta diversity).

2. Material & Methods

2.1. Samples and sequencing data

The analysis was performed on demultiplexed paired-end FASTQ files (Illumina MiSeq, 2 × 300 bp) provided by the client, corresponding to 40 human faecal samples (Control, n = 20; Treated, n = 20). Samples metadata are detailed in Table 1.

2.2. Bioinformatic workflow

Data were processed with nf-core/ampliseq v2.18.0 (Figure 1), a community-curated Nextflow pipeline run in containers, which pins the version of every tool and ensures full reproducibility. The software versions used are listed in Table 2 and are also exported with the results (software_versions.yml and MultiQC report).

2.3. Quality control and read processing

Raw reads were assessed with FastQC and summarised with MultiQC. Primer sequences were removed with Cutadapt and read pairs in which either primer was not found were discarded. Reads were then processed with DADA2: forward and reverse reads were truncated at 270 and 220 bp respectively, based on quality profiles, and reads with more than two expected errors were discarded. Error models were learned from the data, Amplicon Sequence Variants (ASVs) were inferred independently on forward and reverse reads, read pairs were merged and chimeric sequences were removed.

2.4. Taxonomic classification

Amplicon Sequence Variants were classified against the SBDI-GTDB reference database using a naive Bayesian classifier. ASVs assigned to mitochondria or chloroplasts, or unassigned at kingdom level, were removed before downstream analyses. Relative abundance tables were computed at each taxonomic rank, from phylum to species.

2.5. Diversity analysis

2.6. Statistical significance

Unless otherwise stated, tests were two-sided and a p-value or False Discovery Rate (FDR) below 0.05 was considered significant. Alpha diversity indices were compared between groups with the Mann‑Whitney U test. Differences in community composition were tested by PERMANOVA on each distance matrix, with 999 permutations.

2.7. Bioinformatics tools and versions

3. Results

3.1. Sequencing data quality

All 40 samples produced high-quality data, with a median of 101,238 read pairs per sample (range 78,412-131,999). The median proportion of bases with a Phred score ≥ 30 was 93.9% for forward reads and 88.6% for reverse reads. Quality remained high along the full length of forward reads and most of the reverse reads, with the moderate 3′ decline typical of 2 × 300 bp chemistry (Figure 2, Figure 3). Full per-sample metrics (read counts, Q30, GC content, duplication) are provided in Table 3.

GC content profiles were multimodal and highly consistent across samples (Figure 4), a pattern commonly observed in amplicon sequencing data, where reads from different taxa contribute distinct GC compositions. No additional peak was observed that would suggest contamination by host DNA or PhiX. Per-base sequence content, GC content, duplication levels and overrepresented sequences are commonly flagged as unusual by FastQC for amplicon libraries: all reads share the same primer and target region, so many reads are identical copies of the same amplicon sequence variant. 

The high proportion of duplicate reads reported by FastQC is expected for amplicon libraries (Figure 5). Unlike shotgun sequencing, where each read samples a different, essentially random position across a large genome, every read here originates from the same short, fixed target region (the V3-V4 locus, ~460 bp) amplified from a limited pool of template sequences. With most reads drawn from a handful of numerically dominant ASVs, a large fraction are exact copies of one another. Duplication levels of 90% or more are therefore routine for amplicon libraries and do not reflect PCR over-amplification or a loss of library complexity, as they would in a shotgun library.

3.2. Read processing and ASV reference

Primer sequences were then removed with Cutadapt. Read pairs in which either the forward or the reverse primer could not be found were discarded, since their orientation and amplicon boundaries cannot be established reliably, and retaining them would introduce ASVs of variable length at the denoising step. This filter removed only a small fraction of read pairs, consistently across samples (Figure 6), indicating efficient and specific amplification of the V3-V4 target: a high discard rate at this step would point to primer mismatches, off-target amplification or degraded input material.

Read losses were then distributed across the successive DADA2 steps (Table 4, Figure 7). Quality filtering and denoising each removed a small fraction of reads, while read pair merging accounted for the largest single loss, as expected for a ~460 bp amplicon sequenced 2 × 300 bp, where the overlap between forward and reverse reads is limited. Chimeras and off-target sequences were rare: mitochondrial, chloroplast and unassigned sequences accounted for a median of 0.22%, 0.04% and 0.16% of reads, respectively. Overall, a median of 82.9% of raw read pairs was retained, a high yield for V3-V4 data, and a total of 353 ASVs were obtained.

Rarefaction curves assess whether sequencing depth was sufficient to describe each community. Each sample is progressively subsampled to increasing numbers of reads, and a diversity index is recomputed at every depth. The resulting curve shows how much of the community is captured as more reads are considered. A curve that is still rising at the maximum depth indicates that taxa remain undetected and that the sample is under-sequenced, whereas a curve that flattens shows that additional reads no longer reveal new information. Curves are shown here as the mean across samples of each group, with error bars giving the standard deviation.

Rarefaction was performed for two complementary diversity indices (Figure 8). Observed features counts the number of distinct ASVs detected and is the most demanding criterion, since rare taxa are only sampled at high depth. The Shannon index combines the number of ASVs and the evenness of their abundances, and therefore stabilises earlier, as it is dominated by the more abundant taxa. The Shannon index plateaued from approximately 10,000 reads onward, and observed features from approximately 25,000 reads, well below the rarefaction depth of 66,000 reads, in both groups. Sequencing depth was therefore sufficient to capture both the richness and the structure of each community, and the diversity estimates that follow are not limited by sequencing effort.

3.3. Taxonomic composition

The 353 ASVs were assigned to 8 phyla (Figure 9). The community was dominated by Firmicutes (mean 45.8% in Control and 47.9% in Treated) and Bacteroidota (44.1% and 33.9%), followed by Actinobacteriota (4.2% and 10.2%), Verrucomicrobiota and Proteobacteria. At genus level, the most abundant taxa were Bacteroides, Faecalibacterium, Prevotella, Bifidobacterium, Roseburia (Figure 10 and Figure 11). Profiles showed marked inter-individual variability, including a subset of samples enriched in Prevotella, a pattern frequently described in human gut microbiome studies.

3.4. Alpha diversity

Alpha diversity was assessed with four indices capturing distinct properties of the communities: richness alone (observed features), richness weighted by phylogenetic breadth (Faith's PD), and richness combined with, or reduced to, the evenness of abundances (Shannon index and Pielou's evenness). Richness and phylogenetic diversity were comparable between groups (observed features, median 158 vs 154, p = 0.551; Faith's PD, 12.29 vs 12.27, p = 0.766), whereas both the Shannon index (5.91 vs 5.68, p = 0.005) and Pielou's evenness (0.806 vs 0.773, p = 0.020) were significantly lower in the Treated group (Figure 12, Table 5). Since evenness decreased while richness remained unchanged, this pattern reflects a redistribution of abundances among the taxa already present, rather than a loss of taxa from the community.

3.5. Beta diversity

Beta diversity was assessed with four distance metrics differing in whether they weight taxon abundance and whether they account for phylogeny: Jaccard (presence/absence), Bray-Curtis (abundance), and their phylogenetic counterparts, unweighted and weighted UniFrac. Group separation was strongest with the abundance-based metrics, Bray-Curtis (R² = 0.091, p = 0.001) and weighted UniFrac (R² = 0.191, p = 0.002). With the presence/absence metrics, the effect was marginal or absent: Jaccard reached borderline significance with a three- to six-fold smaller effect size (R² = 0.030, p = 0.046), and unweighted UniFrac was not significant (R² = 0.031, p = 0.116; Figure 13, Table 6). This gradient indicates that the intervention shifts the relative abundances of taxa already shared across samples rather than changing which taxa are present, consistent with the loss of evenness observed above. The higher R² obtained with weighted UniFrac further suggests that the taxa driving this shift are phylogenetically clustered rather than dispersed across the tree.

4. Discussion & Conclusion

The sequencing data were of high quality, with high base quality, a high yield of usable reads and very low levels of chimeric and off-target sequences. Sequencing depth was sufficient to capture the diversity of every sample. All 40 samples met the sequencing-depth criterion and were therefore retained for downstream diversity analyses.

The faecal microbiota showed a composition typical of the human gut, dominated by Firmicutes and Bacteroidota. The intervention did not modify the number of taxa present, but reduced the evenness of the community and shifted its structure. Because the group effect was detected with abundance-based metrics (Bray-Curtis, weighted UniFrac) but not, or only weakly, with presence/absence metrics (Jaccard, unweighted UniFrac), the intervention appears to act mainly on the relative proportions of taxa already present rather than by introducing or eliminating taxa. The stronger signal observed with weighted UniFrac indicates that differences in abundance-weighted phylogenetic community structure contribute to the observed separation between groups.

Identifying the taxa responsible for this shift, and its potential functional consequences, requires differential abundance testing and functional prediction, which are included in the Advanced analysis.

5. Deliverables