✕  Close

1 Introduction

1.1 Background

  • This project uses a longitudinal 16S rRNA amplicon sequencing public data to characterize how the murine gut microbiome changes over time and whether these trajectories differ between male and female mice.
  • Samples were collected repeatedly from the same animals across early and late time windows, enabling within-mouse comparisons while leveraging biological replication across individuals.
  • We will generate amplicon sequence variants (ASVs), assign taxonomy, and quantify both within-sample diversity (alpha diversity) and between-sample community shifts (beta diversity).
  • The primary analysis contrasts early (Day0–25) and late (Day124–175) timepoints, with secondary evaluations of sex effects and sex-by-time interactions to assess whether community maturation patterns are consistent or diverge by sex.

1.2 Aims

  1. Characterize longitudinal community change
    • Quantify early-to-late shifts in microbiome structure using beta diversity and ordination
    • Summarize within-mouse stability versus between-mouse variability across time
  2. Test early vs late differences with replication
    • Compare alpha diversity and community composition between early (Day0–25) and late (Day124–175) time windows
    • Evaluate how robust the early/late signal is across mice (biological replicates)
  3. Assess sex effects on microbiome trajectories
    • Estimate the main effect of sex on community structure and diversity
    • Test for sex-by-time-window interaction (do males and females shift differently from early to late?)

1.3 Materials

Group Male Female Total
Early (Day0–25) 102 91 193
Late (Day124–175) 69 68 137
Total 171 159 330
  • Data source
  • Sequencing format
    • Paired-end Illumina MiSeq FASTQ files generated from 16S rRNA amplicons
    • Filenames encode mouse ID, sex (F/M), and day (e.g., M6D141, F3D8)

1.4 Sample name and color key

  • Nomelclature: [ M (male) | F (female) ][ mouse number ]D[ day ]
  • Examples:
    • M6D141: Male mouse 6 at Day 141
    • F3D8: Female mouse 3 at Day 8
  • Color keys:
    • Sex, N=2: ◼ M (male) and ◼ F (female)
    • Group, N=2: ◼ E (early) and ◼ L (late)

2 Environment

  • OS: macOS 26.2
  • Platform: aarch64-apple-darwin20
  • Software: R, FastQC, MultiQC, cutadapt, VSEARCH
  • R pakcages: dada2 (v1.36.0), Biostrings (v2.76.0), ShortRead (v1.66.0), data.table (v1.18.0), ggplot2 (v4.0.1), patchword (v1.3.2), RColorBrewer (v1.1-3), plotly (v4.11.0)

3 Preprocessing

  • The following workflow was adapted from the DADA2 tutorial and modified as needed to generate this quality analysis report.

3.1 FASTQ quality assessment

  • The ridge plot below shows per-base quality scores across read positions. Each x-axis position includes two ridges: forward reads (left) and reverse reads (right)
  • Interpretation: Forward reads have higher overall quality. Reverse read quality drops sharply around 157 bp, whereas forward reads remain consistently above a Phred score of 30

3.2 Adator trimming

  • Optional primer/adapter trimming using cutadapt
$ cutadapt \
  -g FORWARD_PRIMER_SEQ \
  -G REVERSE_PRIMER_SEQ \
  --discard-untrimmed \
  -o sample_R1.trimmed.fastq.gz \
  -p sample_R2.trimmed.fastq.gz \
  sample_R1.fastq.gz sample_R2.fastq.gz

3.3 Preparation in the R environment

  • Set the working directory to import FASTQ files and save all outputs generated in this report
# Do not run
base_dir <- "WORKING_DIRECTORY" # path to your project folder
setwd(base_dir)

3.4 Load R packages

  • Required packages are listed in Section 2
library(dada2) # ASV inference and DADA2 pipeline functions
library(Biostrings) # DNA/RNA sequence handling (FASTA I/O, string ops)
library(ShortRead) # FASTQ input and quality assessment utilities
library(data.table) # fast, memory-efficient tabular data manipulation
library(ggplot2) # static plotting with the grammar of graphics
library(patchwork) # combine multiple ggplot figures into layouts
library(RColorBrewer) # color palettes (e.g., for plots/heatmaps)
library(plotly) # interactive versions of plots (hover/zoom)
library(phyloseq) # microbiome data structures + analysis/visualization

3.5 Import FASTQ files

  • Enumerate R1 files in the FASTQ directory, derive unique sample IDs from filenames, and (if paired-end) validate that each sample has a matching R2 file
  • List FASTQ files, derive sample IDs, and validate paired-end matches
fastq_dir <- file.path(base_dir, "FASTQ", "raw")
r1_suffix <- "_R1.fastq.gz"
r2_suffix <- "_R2.fastq.gz"

r1_files <- sort(list.files(fastq_dir, pattern = paste0(r1_suffix, "$"), full.names = TRUE))
stop_if_not(length(r1_files) > 0, paste0("No R1 files matching '*", r1_suffix, "' found in: ", fastq_dir))

sample_ids <- vapply(r1_files, derive_sample_id_from_r1, character(1))
stop_if_not(!anyDuplicated(sample_ids), "Duplicate sample IDs derived from R1 filenames.")

names(r1_files) <- sample_ids

if (isTRUE(is_paired_end)) {
        r2_files <- file.path(fastq_dir, paste0(sample_ids, r2_suffix))
        missing_r2 <- !file.exists(r2_files)
        stop_if_not(
                !any(missing_r2),
                paste0("Missing R2 files for: ", paste(sample_ids[missing_r2], collapse = ", "))
        )
        names(r2_files) <- sample_ids
} else {
        r2_files <- character(0)
}

3.6 Read quality

3.6.1 Raw reads

  • Per-base quality profiles from the original demultiplexed FASTQs, before any trimming or filtering.
if (isTRUE(is_paired_end)) {
        qp_raw <- filter_qp_inputs(unname(r1_files), unname(r2_files))

        p_raw_R1 <- plotQualityProfile(qp_raw$R1, aggregate = TRUE) +
                labs(subtitle = "Read 1")

        p_raw_R2 <- plotQualityProfile(qp_raw$R2, aggregate = TRUE) +
                labs(subtitle = "Read 2")

        print(p_raw_R1 + p_raw_R2)
} else {
        p_raw_R1 <- plotQualityProfile(unname(r1_files), aggregate = TRUE)
        print(p_raw_R1)
}

3.6.2 After cutadapt

  • Quality profiles after cutadapt processing, showing read quality after primer/adaptor removal and any configured trimming.
if (isTRUE(do_cutadapt)) {
        cutadapt_dir <- base::file.path(base_dir, "FASTQ", "cutadapt")
        ca_r1 <- sort(list.files(cutadapt_dir, pattern = paste0(r1_suffix, "$"), full.names = TRUE))
        ca_ids <- vapply(ca_r1, derive_sample_id_from_r1, character(1))

        if (isTRUE(is_paired_end)) {
                ca_r2 <- file.path(cutadapt_dir, paste0(ca_ids, r2_suffix))
                stop_if_not(all(file.exists(ca_r2)), "Missing some cutadapt R2 files.")

                qp_afterqc <- filter_qp_inputs(ca_r1, ca_r2)

                p_afterqc_R1 <- plotQualityProfile(qp_afterqc$R1, aggregate = TRUE) +
                        labs(subtitle = "After trimming (cutadapt) — Read 1")

                p_afterqc_R2 <- plotQualityProfile(qp_afterqc$R2, aggregate = TRUE) +
                        labs(subtitle = "After trimming (cutadapt) — Read 2")

                print(p_afterqc_R1 + p_afterqc_R2)

                filt_in_r1 <- ca_r1
                filt_in_r2 <- ca_r2
        } else {
                qp_afterqc <- list(R1 = ca_r1, R2 = NULL)

                p_afterqc_R1 <- plotQualityProfile(qp_afterqc$R1, aggregate = TRUE) +
                        labs(subtitle = "After trimming (cutadapt) — Read 1")

                print(p_afterqc_R1)

                filt_in_r1 <- ca_r1
                filt_in_r2 <- NULL
        }
} else {
        filt_in_r1 <- unname(r1_files)
        filt_in_r2 <- if (isTRUE(is_paired_end)) unname(r2_files) else NULL
}