Analisis Esensial dengan phyloseq — Script & Penjelasan

Dokumen ini memuat script R lengkap untuk analisis esensial menggunakan phyloseq dan penjelasan langsung di bawah tiap blok kode.

🧪

Skrip ini mengasumsikan Anda sudah menjalankan tahap DADA2 sehingga tersedia: ASV_table.csv, 99_taxonomy_table.csv, dan ASV_sequences.fasta di folder output/dada2_outputs/.

0) Setup

suppressPackageStartupMessages({
  library(phyloseq)
  library(Biostrings)
  library(ggplot2)
  library(vegan)      # PERMANOVA (adonis2), diversity helpers
  library(DECIPHER)   # optional: alignment for tree
  library(phangorn)   # optional: ML tree for UniFrac
})

Paket yang dipanggil:

phyloseq untuk integrasi data mikrobioma.
Biostrings untuk manipulasi sekuens DNA.
ggplot2 untuk visualisasi.
vegan untuk analisis keanekaragaman & PERMANOVA.
DECIPHER dan phangorn opsional bila ingin membuat pohon filogenetik untuk UniFrac.

Path Input & Output

projdir      <- "./../.."
in_dada2_dir <- file.path(projdir, "output/dada2_outputs")
outdir       <- file.path(projdir, "output/phloseq_outputs")
dir.create(outdir, showWarnings = FALSE, recursive = TRUE)

in_dada2_dir: folder tempat hasil DADA2 disimpan.
outdir: semua hasil analisis phyloseq akan disimpan di sini.

otu_csv      <- file.path(in_dada2_dir, "ASV_table.csv")
tax_csv      <- file.path(in_dada2_dir, "99_taxonomy_table.csv")
asv_fasta    <- file.path(in_dada2_dir, "ASV_sequences.fasta")
metadata_csv <- file.path(projdir, "data", "sample_metadata.csv")
stopifnot(file.exists(otu_csv), file.exists(tax_csv), file.exists(asv_fasta))

Input wajib: ASV table, taxonomy table, dan ASV fasta. Metadata CSV opsional.

1) Load Data & Buat Objek phyloseq

otu_df <- read.csv(otu_csv, row.names = 1, check.names = FALSE, stringsAsFactors = FALSE)
rownames(otu_df) <- trimws(rownames(otu_df))
colnames(otu_df) <- trimws(colnames(otu_df))
otu_df[] <- lapply(otu_df, function(x) as.numeric(as.character(x)))

Membaca tabel ASV (baris = sampel, kolom = ASV) dan memastikan semua nilai berupa numerik.

tax_df <- read.csv(tax_csv, check.names = FALSE)
rownames(tax_df) <- tax_df$ASV
tax_df$ASV <- NULL

Membaca tabel taksonomi; setiap baris adalah ASV dengan rank (Phylum, Genus, dll.).

common_asvs <- intersect(colnames(otu_df), rownames(tax_df))
otu_df <- otu_df[, common_asvs, drop=FALSE]
tax_df <- tax_df[common_asvs, , drop=FALSE]

Menyamakan ASV yang muncul di kedua tabel agar konsisten.

asv_seqs <- readDNAStringSet(asv_fasta)
asv_seqs <- asv_seqs[names(asv_seqs) %in% common_asvs]
asv_seqs <- asv_seqs[match(common_asvs, names(asv_seqs))]

Mengambil sekuens DNA untuk setiap ASV (opsional, diperlukan untuk pohon filogenetik).

if (file.exists(metadata_csv)) {
  meta_df <- read.csv(metadata_csv, check.names = FALSE)
  rownames(meta_df) <- meta_df$SampleID
} else {
  meta_df <- data.frame(SampleID = rownames(otu_df), Group = "GroupA", row.names = rownames(otu_df))
}

Membaca metadata sampel. Jika tidak ada, skrip membuat dummy metadata dengan semua sampel = GroupA.

common_samples <- intersect(rownames(otu_df), rownames(meta_df))
otu_df  <- otu_df[common_samples, , drop=FALSE]
meta_df <- meta_df[common_samples, , drop=FALSE]
 
OTU  <- otu_table(as.matrix(otu_df), taxa_are_rows = FALSE)
TAX  <- tax_table(as.matrix(tax_df))
SAMP <- sample_data(meta_df)
ps   <- phyloseq(OTU, TAX, SAMP)
saveRDS(ps, file.path(outdir, "01_phyloseq_input.rds"))

Membangun objek phyloseq (ps) yang menggabungkan data ASV, taksonomi, dan metadata.

2) Eksplorasi Awal (mia/miaViz)

tse <- mia::convertFromPhyloseq(ps)
tse <- mia::transformAssay(tse, method = "relabundance")

Mengubah phyloseq ke TreeSummarizedExperiment (TSE) dan menghitung relative abundance.

p_genus <- miaViz::plotAbundance(tse, abund_values = "relabundance", rank = "Genus") +
  theme_minimal(11) +
  labs(title = "Genus composition (relative)", x = "Sample", y = "Relative abundance") +
  scale_y_continuous(labels = scales::percent_format(1), limits = c(0,1))
ggsave(file.path(outdir, "02_abundance_genus.png"), p_genus, dpi = 300)

Membuat barplot komposisi genus dan menyimpan ke PNG. Sama untuk level Phylum.

Genus composition

Phylum composition

3) Alpha Diversity

alpha_df <- data.frame(
  SampleID = sample_names(ps),
  Observed = phyloseq::estimate_richness(ps, measures = "Observed")[,1],
  Shannon  = phyloseq::estimate_richness(ps, measures = "Shannon")[,1],
  Simpson  = phyloseq::estimate_richness(ps, measures = "Simpson")[,1],
  Group    = sample_data(ps)$Group
)
write.csv(alpha_df, file.path(outdir, "02_alpha_diversity.csv"), row.names = FALSE)

Menghitung 3 metrik keanekaragaman alfa: Observed, Shannon, dan Simpson.
Disimpan ke CSV.

p_alpha <- function(metric) {
  ggplot(alpha_df, aes(x = Group, y = .data[[metric]])) +
    geom_boxplot(outlier.shape = NA) +
    geom_jitter(width = 0.15, alpha = 0.6) +
    labs(x = "Group", y = metric, title = paste(metric, "diversity")) +
    theme_bw(base_size = 12)
}
ggsave(file.path(outdir, "02a_alpha_observed.png"), p_alpha("Observed"))

Membuat plot boxplot + jitter untuk setiap metrik alpha diversity.

4) Pohon Filogenetik (Opsional, untuk UniFrac)

do_tree <- TRUE
if (do_tree) {
  aln   <- DECIPHER::AlignSeqs(asv_seqs, processors = 0)
  phdat <- phangorn::phyDat(as(aln, "matrix"), type="DNA")
  dm    <- phangorn::dist.ml(phdat)
  tree0 <- phangorn::NJ(dm)
  fit   <- phangorn::pml(tree0, data = phdat)
  fit   <- phangorn::optim.pml(fit, model="GTR", optInv=TRUE, optGamma=TRUE, rearrangement="stochastic")
  tree  <- fit$tree
  ps    <- merge_phyloseq(ps, phy_tree(tree))
  saveRDS(ps, file.path(outdir, "03_phyloseq_with_tree.rds"))
}

Jika do_tree = TRUE, skrip akan:

Align sekuens ASV.
Hitung jarak evolusi (dist.ml).
Buat pohon Neighbor-Joining, optimasi Maximum Likelihood (model GTR).
Simpan ps baru dengan pohon.

5) Beta Diversity + PERMANOVA

ps_rel <- transform_sample_counts(ps, function(x) x / sum(x))
dist_bray <- phyloseq::distance(ps_rel, method = "bray")
ord_bray  <- ordinate(ps_rel, method = "PCoA", distance = dist_bray)

Bray–Curtis: jarak berbasis kelimpahan relatif antar sampel.
PCoA: ordination untuk visualisasi jarak.

p_bray <- plot_ordination(ps_rel, ord_bray, color="Group") +
  geom_point(size=3, alpha=0.85) +
  theme_bw(base_size = 12) +
  labs(title = "PCoA (Bray-Curtis)")
ggsave(file.path(outdir, "04a_pcoa_bray.png"), p_bray, width=6.5, height=5, dpi=300)

Menyimpan gambar PCoA Bray-Curtis.

meta_for_adonis <- as(sample_data(ps_rel), "data.frame")
adon_bray <- adonis2(as.matrix(dist_bray) ~ Group, data = meta_for_adonis, permutations = 999)
sink(file.path(outdir, "04b_permanova_bray.txt")); print(adon_bray); sink()

Melakukan PERMANOVA (adonis2) untuk menguji apakah komposisi komunitas berbeda signifikan antar grup.

Ringkasan Output

02_abundance_genus.png & 02_abundance_phylum.png: barplot komposisi taksonomi.
02a_alpha_observed.png, 02b_alpha_shannon.png, 02c_alpha_simpson.png: plot alpha diversity.
03_phyloseq_with_tree.rds: phyloseq dengan pohon filogenetik.
04a_pcoa_bray.png: plot PCoA Bray-Curtis.
04b_permanova_bray.txt: hasil PERMANOVA.

📊

Output ini bisa langsung dipakai untuk laporan, publikasi, atau workshop sebagai analisis mikrobioma standar.

3 Amplicon Sequence Variants (DADA2)About Us