Roary Pangenome

# Install Roary via conda/mamba (recommended)
mamba install -c conda-forge -c bioconda roary

# Verify installation
roary --version
# 3.13.0

# Verify dependent tools
which cd-hit blastp mcl bedtools mafft
# /opt/conda/bin/cd-hit
# /opt/conda/bin/blastp
# /opt/conda/bin/mcl
# /opt/conda/bin/bedtools
# /opt/conda/bin/mafft

# Install Python parsing dependencies
pip install pandas matplotlib seaborn biopython dendropy

decisions:
  - id: D1
    param: inputAnnotations
    kind: derived
    source: upstream
    ask: "Which GFF3 files, and were they all produced by the same annotation tool and version?"
    default: "the annotation stage output"

  - id: D2
    param: identityThreshold
    kind: required
    source: user
    ask: "How similar must two proteins be to count as the same gene? Lower values merge diverged orthologs; higher values split them."
    default: "95%"

  - id: D3
    param: coreDefinition
    kind: required
    source: user
    ask: "In what fraction of isolates must a gene appear before it counts as core?"
    default: "99%"

  - id: D4
    param: paralogSplitting
    kind: required
    source: user
    ask: "Should duplicated genes within a genome be split into separate families, or kept together?"
    default: "split"

  - id: D5
    param: coreAlignment
    kind: required
    source: upstream
    ask: "Is a concatenated core-gene alignment needed for a downstream phylogeny?"
    default: "not produced - it is the slowest stage by a wide margin"

  - id: D6
    param: alignmentSpeed
    kind: optional
    source: user
    depends_on: [D5]
    ask: "Align every gene accurately, or use the fast route sufficient for tree building?"
    default: "fast alignment"
    skip_if: "no core alignment requested"

  - id: D7
    param: clusterInflation
    kind: optional
    source: user
    ask: "How tightly should the clustering step group families?"
    default: 1.5

  - id: D8
    param: maxClusters
    kind: optional_conditional
    source: data
    ask: "Is the expected number of gene families above the default cap?"
    default: 50000

  - id: D9
    param: processes
    kind: never_ask
    source: data
    reason: "Affects runtime only, not the gene families"
    default: "min(8, available_cores)"

# Run Roary on all GFF3 files in current directory; emit core gene alignment
roary -e --mafft -p 8 -o pangenome -f roary_out/ *.gff

# Inspect summary statistics
cat roary_out/summary_statistics.txt
# Core genes        (99% <= strains <= 100%) 2823
# Soft core genes   (95% <= strains <  99%)   78
# Shell genes       (15% <= strains <  95%)  1542
# Cloud genes       ( 0% <= strains <  15%)  3104
# Total genes                                 7547

# Create a dedicated conda environment
mamba create -n roary_env -c conda-forge -c bioconda roary mafft fasttree python=3.11 -y
mamba activate roary_env

# Verify Roary
roary --version
# 3.13.0

# Confirm pipeline dependencies are reachable
for tool in cd-hit blastp mcl bedtools mafft FastTree; do
    if command -v $tool >/dev/null; then
        echo "OK: $tool"
    else
        echo "MISSING: $tool"
    fi
done

# Install Python parsing tools
pip install pandas matplotlib seaborn biopython dendropy

# Verify GFF3 files include the embedded FASTA section
mkdir -p roary_input/
cp annotations/*/sample*.gff roary_input/

for GFF in roary_input/*.gff; do
    if grep -q "^##FASTA" "$GFF"; then
        echo "OK: $GFF"
    else
        echo "MISSING ##FASTA in: $GFF"
    fi
done

# Roary uses GFF filename (without .gff) as the sample column name.
# Rename or symlink to ensure unique, descriptive names.
cd roary_input/
ls *.gff | head
# strain_A.gff  strain_B.gff  strain_C.gff

# Sanity-check sample size and unique names
from pathlib import Path

gff_files = sorted(Path("roary_input/").glob("*.gff"))
print(f"GFF files: {len(gff_files)}")

names = [g.stem for g in gff_files]
duplicates = [n for n in names if names.count(n) > 1]
assert not duplicates, f"Duplicate sample names: {set(duplicates)}"

# Show first 5
for g in gff_files[:5]:
    size_mb = g.stat().st_size / 1e6
    print(f"  {g.name}: {size_mb:.1f} MB")

# Run Roary on the prepared GFF directory
# -e: extract core gene alignment per gene, then concatenate
# --mafft: use MAFFT for alignment (faster than default PRANK)
# -p 8: 8 parallel processes
# -i 95: min BLASTP percentage identity (default 95)
# -cd 99: % isolates a gene must be in to be "core" (default 99)
roary -e --mafft \
    -p 8 \
    -i 95 \
    -cd 99 \
    -o pangenome \
    -f roary_out/ \
    roary_input/*.gff

# Expected runtime:
#   ~50 genomes:  10–20 min
#   ~500 genomes: 4–8 hours

ls roary_out/
# accessory_binary_genes.fa            gene_presence_absence.csv
# accessory_binary_genes.fa.newick     gene_presence_absence.Rtab
# accessory.header.embl                number_of_conserved_genes.Rtab
# accessory.tab                        number_of_genes_in_pan_genome.Rtab
# blast_identity_frequency.Rtab        number_of_new_genes.Rtab
# clustered_proteins                   number_of_unique_genes.Rtab
# core_accessory.header.embl           pan_genome_reference.fa
# core_accessory.tab                   summary_statistics.txt
# core_gene_alignment.aln

import pandas as pd
import numpy as np

df = pd.read_csv("roary_out/gene_presence_absence.csv", low_memory=False)
print(f"Gene families: {len(df)}")
print(f"Columns: {len(df.columns)}")
# Standard metadata columns end with 'Inference', then sample columns follow
metadata_cols = list(df.columns[:14])
sample_cols = list(df.columns[14:])
print(f"Samples: {len(sample_cols)}")

# Build a binary presence/absence matrix (1 = locus tag present, 0 = empty)
binary = df[sample_cols].notna().astype(int)
binary.index = df["Gene"]
print(f"Binary matrix shape: {binary.shape}")

# Pan-genome partition by frequency
n_samples = len(sample_cols)
freq = binary.sum(axis=1) / n_samples
core      = freq[freq >= 0.99]
soft_core = freq[(freq >= 0.95) & (freq < 0.99)]
shell     = freq[(freq >= 0.15) & (freq < 0.95)]
cloud     = freq[freq < 0.15]

print(f"\nPan-genome partition (n={n_samples} genomes):")
print(f"  Core      (>=99 %): {len(core):>5}")
print(f"  Soft-core (95–99 %): {len(soft_core):>5}")
print(f"  Shell    (15–95 %): {len(shell):>5}")
print(f"  Cloud      (<15 %): {len(cloud):>5}")
print(f"  Total            : {len(freq):>5}")

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

df = pd.read_csv("roary_out/gene_presence_absence.csv", low_memory=False)
sample_cols = list(df.columns[14:])
n = len(sample_cols)
binary = df[sample_cols].notna().astype(int)
freq = binary.sum(axis=1)

fig, axes = plt.subplots(1, 2, figsize=(11, 4))

# Histogram of gene frequencies
axes[0].hist(freq, bins=range(1, n + 2), color="#2980B9", edgecolor="white")
axes[0].axvline(0.99 * n, color="red", linestyle="--", label="99 % core cutoff")
axes[0].axvline(0.15 * n, color="orange", linestyle="--", label="15 % shell cutoff")
axes[0].set_xlabel("Number of genomes containing the gene")
axes[0].set_ylabel("Number of gene families")
axes[0].set_title("Pan-genome frequency distribution")
axes[0].legend()

# Pie chart of pan-genome partition
freqp = freq / n
core_n      = (freqp >= 0.99).sum()
soft_n      = ((freqp >= 0.95) & (freqp < 0.99)).sum()
shell_n     = ((freqp >= 0.15) & (freqp < 0.95)).sum()
cloud_n     = (freqp < 0.15).sum()
axes[1].pie([core_n, soft_n, shell_n, cloud_n],
            labels=[f"Core ({core_n})", f"Soft-core ({soft_n})",
                    f"Shell ({shell_n})", f"Cloud ({cloud_n})"],
            colors=["#27AE60", "#3498DB", "#F39C12", "#95A5A6"],
            autopct="%1.1f%%", startangle=90)
axes[1].set_title(f"Pan-genome partition (n={n} genomes)")

plt.tight_layout()
plt.savefig("pangenome_distribution.png", dpi=150, bbox_inches="tight")
print(f"Saved pangenome_distribution.png  (total families: {len(freq)})")

# FastTree expects an aligned multi-FASTA; Roary produces this with `-e`
FastTree -nt -gtr -nosupport \
    < roary_out/core_gene_alignment.aln \
    > roary_out/core_gene_tree.nwk

echo "Tree built: roary_out/core_gene_tree.nwk"
# (Optional) IQ-TREE for support values:
# iqtree -s core_gene_alignment.aln -m GTR+G -bb 1000 -nt 8

# Visualize the FastTree result alongside the gene presence/absence matrix
import dendropy
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np

tree = dendropy.Tree.get(path="roary_out/core_gene_tree.nwk", schema="newick")
leaves = [leaf.taxon.label for leaf in tree.leaf_node_iter()]
print(f"Tree leaves: {len(leaves)}")
print(f"Tree length: {tree.length():.4f} substitutions/site")

# Print first 10 leaves with their root-to-tip distances
for leaf in list(tree.leaf_node_iter())[:10]:
    dist = leaf.distance_from_root()
    print(f"  {leaf.taxon.label:<25}  d={dist:.4f}")

# Use the bundled python helper from the Roary install
# (download from: https://github.com/sanger-pathogens/Roary/blob/master/contrib/roary_plots/roary_plots.py)
python roary_plots.py \
    roary_out/core_gene_tree.nwk \
    roary_out/gene_presence_absence.csv \
    --format png \
    --labels

ls *.png
# pangenome_frequency.png
# pangenome_matrix.png
# pangenome_pie.png

import pandas as pd

df = pd.read_csv("roary_out/gene_presence_absence.csv", low_memory=False)
sample_cols = list(df.columns[14:])
n = len(sample_cols)
binary = df[sample_cols].notna().astype(int)
binary.index = df["Gene"]

freq = binary.sum(axis=1)
cloud_mask = freq < (0.15 * n)
shell_mask = (freq >= 0.15 * n) & (freq < 0.95 * n)

per_genome = pd.DataFrame({
    "genes_total":   binary.sum(axis=0),
    "cloud_genes":   binary.loc[cloud_mask].sum(axis=0),
    "shell_genes":   binary.loc[shell_mask].sum(axis=0),
})
per_genome = per_genome.sort_values("cloud_genes", ascending=False)
print(per_genome.head(10).to_string())

per_genome.to_csv("per_genome_accessory_counts.csv")
print("\nSaved per_genome_accessory_counts.csv")