16S Amplicon — Human Body-Site Microbiome (SciFlow)

Classify the mock 16S rRNA community to genus level with the Ribosomal Database Project naive-Bayes classifier (bootstrap confidence >= 0.8). Directly answers 'what microbes are in this sample'.

Verbatim command / script not recorded in this run's stage record.

• input_data_summary.htmlreport · 171 KB
• mock_16s_community.rdpRDP · 42 KB
• moving_pictures_16s.rdpRDP · 5.1 MB
• rdp_report.htmlreport · 9 KB
• summary.jsondata · 2 KB

Pivot the RDP fixrank output into a taxa x samples relative-abundance matrix (mock_abundance.csv) and carry env.csv (sample biome metadata) through for between-group comparison.

# Pivot RDP fixrank output -> taxa x samples abundance matrix (INTEGER COUNTS)
# Output named exactly abundance.csv + env.csv (community_analysis.R hardcodes these)
# FIX(C1 2026-05-31): was dividing by sample_tot to emit relative proportions
# (0.0-1.0), causing detect_data_type() to classify as plant_community and
# skip rarefaction + collapse NMDS to stress=0. Now emits raw integer read counts.
import os, csv, glob, shutil
from collections import defaultdict

INPUT_DIR = os.environ.get("INPUT_DIR", "/work/input")
OUTPUT_DIR = os.environ.get("OUTPUT_DIR", "/work/output")
os.makedirs(OUTPUT_DIR, exist_ok=True)
CONF = 0.8

def sample_of(seqid):
    return seqid.split("_", 1)[0]

def parse_fixrank(path, counts):
    with open(path) as fh:
        for line in fh:
            line = line.rstrip("\n")
            if not line:
                continue
            parts = line.split("\t")
            if len(parts) < 5:
                continue
            sample = sample_of(parts[0])
            best = None
            i = 2
            while i + 2 < len(parts):
                taxon = parts[i].strip().strip('"')
                try:
                    conf = float(parts[i + 2])
                except ValueError:
                    i += 3; continue
                if conf >= CONF and taxon and taxon.lower() != "root":
                    best = taxon
                i += 3
            if best:
                counts[(sample, best)] += 1

counts = defaultdict(int)
rdp_files = glob.glob(os.path.join(INPUT_DIR, "*.rdp"))
print("RDP files found:", rdp_files)
for p in rdp_files:
    parse_fixrank(p, counts)

samples = sorted({s for (s, _t) in counts})
taxa = sorted({t for (_s, t) in counts})
print(f"Parsed {len(taxa)} taxa across {len(samples)} samples")

total_reads = sum(counts.values())
print(f"Total read assignments: {total_reads} (integer counts, rarefaction-ready)")

# Write taxa-as-rows, samples-as-columns (community_analysis.R auto-transposes via
# env-overlap detection or dimension heuristic to get samples-as-rows for analysis)
out = os.path.join(OUTPUT_DIR, "abundance.csv")
with open(out, "w", newline="") as fh:
    w = csv.writer(fh)
    w.writerow(["taxon"] + samples)
    for t in taxa:
        row = [t] + [counts.get((s, t), 0) for s in samples]
        w.writerow(row)
print("Wrote", out)

src_env = os.path.join(INPUT_DIR, "env.csv")
if os.path.exists(src_env):
    shutil.copy2(src_env, os.path.join(OUTPUT_DIR, "env.csv"))
    print("Carried env.csv through")
else:
    print("No env.csv in input dir")

• abundance.csvtable · 48 KB
• env.csvtable · 904 B

Alpha diversity (Shannon/Simpson/observed), beta diversity (Bray-Curtis), and ordination (PCoA/NMDS) across the 6 samples, tested between gut and soil biomes. Produces the interactive community report.

Verbatim command / script not recorded in this run's stage record.

• alpha_diversity.csvtable · 3 KB
• alpha_diversity.htmlreport · 3.6 MB
• alpha_diversity_plot.htmlreport · 3.6 MB
• alpha_diversity_subject.htmlreport · 3.6 MB
• alpha_kruskal_wallis.csvtable · 235 B
• betadisper.htmlreport · 3.6 MB
• ca.htmlreport · 3.6 MB
• cca.htmlreport · 3.7 MB
• community_analysis.RR · 20 KB
• core_microbiome.htmlreport · 395 KB
• dbrda.htmlreport · 3.7 MB
• dca.htmlreport · 3.6 MB
• heatmap.htmlreport · 3.7 MB
• hierarchical_clustering.htmlreport · 3.6 MB
• hill_numbers.htmlreport · 3.6 MB
• indicator_species.htmlreport · 446 KB
• indicator_species_chart.htmlreport · 3.6 MB
• nestedness.htmlreport · 3.6 MB
• nmds.htmlreport · 3.6 MB
• not_applicable.htmlreport · 2 KB
• pairwise_permanova.htmlreport · 382 KB
• pcoa.htmlreport · 3.6 MB
• pcoa_plot.htmlreport · 3.6 MB
• permanova_heatmap_body_site.htmlreport · 3.6 MB

+ 19 more output file(s) — see Raw artifacts.