library(qiime2R)
library(ANCOMBC)
library(phyloseq)
library(ggplot2)
library(tidyverse)

# Step 2.1: Create phyloseq object
SVs <- read_qza("{feature_table}")$data
taxonomy <- read_qza("{taxonomy_qza}")$data %>% parse_taxonomy()
metadata <- read_q2metadata("{metadata_tsv}")
tree <- read_qza("{tree_qza}")$data

ps <- phyloseq(otu_table(SVs, taxa_are_rows=TRUE),
               tax_table(as.matrix(taxonomy)),
               sample_data(metadata),
               phy_tree(tree))

# Step 2.2: Rarefaction
ps_rare <- rarefy_even_depth(ps, rngseed=42)

# Step 2.3: Alpha diversity
alpha <- estimate_richness(ps_rare, measures=c("Observed","Shannon","Simpson"))
alpha$sample <- rownames(alpha)
alpha <- merge(alpha, metadata, by.x="sample", by.y="sample_name")

ggplot(alpha, aes(x=condition, y=Shannon, fill=condition)) +
  geom_boxplot() + geom_jitter(width=0.2) +
  theme_bw(base_size=14) + labs(title="Alpha Diversity (Shannon)")
ggsave("05_alpha_diversity.pdf", width=6, height=5)

# Step 2.4: Beta diversity (PCoA)
ord <- ordinate(ps_rare, method="PCoA", distance="bray")
p <- plot_ordination(ps_rare, ord, color="condition", shape="season") +
  geom_point(size=4) + theme_bw(base_size=14) +
  labs(title="Bray-Curtis PCoA")
ggsave("05_pcoa.pdf", p, width=7, height=6)

# Step 2.5: ANCOM-BC differential abundance
abc <- ancombc(phyloseq=ps, formula="condition", group="condition",
               p_adj_method="holm", alpha=0.05)
write.csv(abc$res, "05_ancombc_results.csv")

#!/bin/bash
# Step 3.1: Run PICRUSt2
qiime picrust2 full-pipeline \
  --i-seq {rep_seqs} \
  --i-table {feature_table} \
  --output-dir 06_picrust2 \
  --p-placement-tool sepp \
  --p-threads 8 \
  --p-hsp-method mp \
  --p-max-nsti 2

# Step 3.2: Add descriptions
add_descriptions.py -i 06_picrust2/ec_metagenome_out/pred_metagenome_contrib.tsv -m EC -o 06_picrust2/ec_named.tsv
add_descriptions.py -i 06_picrust2/ko_metagenome_out/pred_metagenome_contrib.tsv -m KO -o 06_picrust2/ko_named.tsv

# Step 3.3: Pathway inference
pathway_pipeline.py -i 06_picrust2/ko_metagenome_out/pred_metagenome_unstrat.tsv.gz -o 06_picrust2/pathways_out -p 8

library(phyloseq)
library(vegan)
library(ggplot2)
library(tidyverse)

# Step 4.1: Load data
ps <- readRDS("{microbiome_ps}")
host <- read.csv("{host_data}")

# Step 4.2: Procrustes analysis
micro_dist <- phyloseq::distance(ps, method="bray")
micro_pcoa <- cmdscale(micro_dist, k=2)
host_pca <- prcomp(host[,-1], scale=TRUE)$x[,1:2]

proc <- protest(micro_pcoa, host_pca, permutations=999)
print(paste("Procrustes m12:", round(proc$t0, 4), "p-value:", proc$signif))

# Step 4.3: Plot Procrustes
proc_df <- data.frame(
  micro_X=micro_pcoa[,1], micro_Y=micro_pcoa[,2],
  host_X=host_pca[,1], host_Y=host_pca[,2]
)

ggplot(proc_df) +
  geom_point(aes(x=micro_X, y=micro_Y), color="#006334", size=3, alpha=0.6) +
  geom_point(aes(x=host_X, y=host_Y), color="#ff6f00", size=3, alpha=0.6) +
  geom_segment(aes(x=micro_X, y=micro_Y, xend=host_X, yend=host_Y),
               color="gray60", alpha=0.4) +
  theme_bw(base_size=14) +
  labs(title=paste("Procrustes: m12 =", round(proc$t0, 3)),
       x="PC1", y="PC2")
ggsave("07_procrustes.pdf", width=8, height=6)

# Step 4.4: Correlation network
corr_results <- data.frame()
for(mg in colnames(host[,-1])) {
  for(tax in taxa_names(ps)[1:50]) {
    cor_val <- cor(host[[mg]], as.vector(otu_table(ps)[tax,]), method="spearman")
    if(abs(cor_val) > 0.5) {
      corr_results <- rbind(corr_results, data.frame(metabolite=mg, taxon=tax, rho=cor_val))
    }
  }
}
write.csv(corr_results, "07_host_microbiome_correlations.csv")