Finding the SusC–SusD polysaccharide-utilization pair with Pynteny¶
Members of the phylum Bacteroidota harvest complex polysaccharides using
polysaccharide-utilization loci (PULs). At the heart of every PUL sits a tandem pair of
outer-membrane proteins — susC (a TonB-dependent transporter) and susD (a
substrate-binding lipoprotein) — that together form the "pedal-bin" import machine.
susC is a great illustration of why sequence similarity alone is a weak annotation signal:
it belongs to the large SusC/RagA/OmpW family of TonB-dependent receptors found across many
phyla, so a lone susC-like hit tells you little. The susC–susD tandem, on the other
hand, is a crisp syntenic signature of a genuine PUL — exactly the kind of genomic-context
signal Pynteny is built to exploit.
In this notebook we show, across three Bacteroidota and three non-Bacteroidota genomes, that:
susCalone is promiscuous — it lights up across phyla;- the
susC–susDtandem is specific to Bacteroidota; and - relaxing the strand constraint reveals the full PUL repertoire — dozens of loci in a gut Bacteroides.
The two HMMs and their metadata are committed under data/, so you do not need the full
432 MB PGAP database.
Setup¶
import warnings
warnings.filterwarnings("ignore", message=".*IProgress.*") # cosmetic tqdm-in-notebook notice
from pathlib import Path
import urllib.request, gzip, shutil, logging
import pandas as pd
from pynteny.api import Build, Search
BASE = Path.cwd()
DATA = BASE / "data"
HMM_DIR = DATA / "hmms"
HMM_META = DATA / "sus_hmm_meta.tsv"
GENOMES = DATA / "genomes" # downloaded below (git-ignored)
PEPTIDES = DATA / "peptides" # built below (git-ignored)
RESULTS = BASE / "results"
for d in (GENOMES, PEPTIDES, RESULTS):
d.mkdir(parents=True, exist_ok=True)
# Negatives legitimately match nothing, which Pynteny logs at ERROR level; silence <= ERROR so
# expected "no hits" cases don't print scary red lines.
logging.disable(logging.ERROR)
pd.set_option("display.max_colwidth", 60)
# susC/RagA is a large, promiscuous family: use a permissive reporting threshold and let the
# synteny filter (not the bit-score) provide the specificity.
HMMSEARCH_ARGS = "-E 1e-10"
print("HMMs committed for this case study:",
", ".join(p.name for p in sorted(HMM_DIR.glob('*.HMM'))))
HMMs committed for this case study: NF033071.0.HMM, TIGR04056.1.HMM
1. The two genes¶
Two HMMs from the NCBI PGAP collection. We give each a clean gene_symbol so the synteny
structure reads in plain susC/susD terms (via --gene_ids).
pd.read_csv(HMM_META, sep="\t")
| #ncbi_accession | label | product_name | gene_symbol | ec_numbers | |
|---|---|---|---|---|---|
| 0 | TIGR04056.1 | susC | SusC/RagA family TonB-linked outer membrane protein | susC | NaN |
| 1 | NF033071.0 | susD | starch-binding outer membrane lipoprotein SusD | susD | NaN |
2. The genome panel¶
Three Bacteroidota (expected to carry PULs) and three non-Bacteroidota controls. Two of
the controls (Pseudomonas, E. coli) deliberately carry susC-like TonB-dependent receptors
but no susD partner — the crux of the specificity story.
genomes = (pd.read_csv(BASE / "genomes.tsv", sep="\t")
.rename(columns=lambda c: c.lstrip("#")))
genomes[["name", "role", "lineage", "lifestyle", "sus_expectation"]]
| name | role | lineage | lifestyle | sus_expectation | |
|---|---|---|---|---|---|
| 0 | Bacteroides_thetaiotaomicron_VPI5482 | positive | Bacteroidota | human gut symbiont | dozens of susC-susD PULs (model Sus organism) |
| 1 | Bacteroides_fragilis_NCTC9343 | positive | Bacteroidota | human gut symbiont | many susC-susD PULs |
| 2 | Flavobacterium_johnsoniae_UW101 | positive | Bacteroidota | soil/freshwater glider | many susC-susD PULs (environmental Bacteroidota) |
| 3 | Pseudomonas_aeruginosa_PAO1 | negative | Pseudomonadota (Gamma) | soil/opportunist | susC-like TonB receptors present, but NO susD tandem |
| 4 | Escherichia_coli_MG1655 | negative | Pseudomonadota (Gamma) | enteric | a few susC-like TonB receptors, no susD tandem |
| 5 | Bacillus_subtilis_168 | negative | Bacillota (Firmicutes) | soil saprophyte | no SusC/SusD system |
3. Build labelled peptide databases¶
pynteny build translates every ORF (Prodigal) and writes a FASTA whose record labels encode
each gene's contig, position and strand — the positional information the synteny filter
reads. Already-built genomes are skipped, so re-running is cheap.
NCBI = "https://ftp.ncbi.nlm.nih.gov/genomes/all"
def fetch_and_build(name, refseq_path):
fna, faa = GENOMES / f"{name}.fna", PEPTIDES / f"{name}.faa"
if faa.exists():
return faa
if not fna.exists():
base = refseq_path.rsplit("/", 1)[-1]
gz = fna.with_suffix(".fna.gz")
urllib.request.urlretrieve(f"{NCBI}/{refseq_path}/{base}_genomic.fna.gz", gz)
with gzip.open(gz, "rb") as fi, open(fna, "wb") as fo:
shutil.copyfileobj(fi, fo)
gz.unlink()
Build(data=fna, outfile=faa).run()
return faa
for _, row in genomes.iterrows():
fetch_and_build(row["name"], row["refseq_path"])
n = sum(1 for line in open(PEPTIDES / f"{row['name']}.faa") if line.startswith(">"))
print(f"{row['name']:<40} {n:>5} ORFs")
Bacteroides_thetaiotaomicron_VPI5482 4964 ORFs Bacteroides_fragilis_NCTC9343 4351 ORFs Flavobacterium_johnsoniae_UW101 5173 ORFs Pseudomonas_aeruginosa_PAO1 5715 ORFs Escherichia_coli_MG1655 4319 ORFs Bacillus_subtilis_168 4242 ORFs
4. Search A — susC alone (no genomic context)¶
>susC
A bare search for the SusC/RagA receptor family.
DISPLAY_COLS = ["contig", "gene_number", "locus", "strand", "gene_symbol"]
def synteny_search(genome, struct, *, unordered=False, subdir="tmp"):
# Run a Pynteny synteny search; return the hits DataFrame (empty if no match).
outdir = RESULTS / subdir / genome
outdir.mkdir(parents=True, exist_ok=True)
try:
hits = Search(
data=PEPTIDES / f"{genome}.faa", synteny_struc=struct, gene_ids=True,
unordered=unordered, hmm_dir=HMM_DIR, hmm_meta=HMM_META,
hmmsearch_args=HMMSEARCH_ARGS, outdir=outdir,
).run()
return hits.hits
except SystemExit: # Pynteny exits when an HMM matches nothing -> a negative
return pd.DataFrame()
hitsA = {g: synteny_search(g, ">susC", subdir="A_susC_alone") for g in genomes["name"]}
pd.DataFrame({"role": genomes.set_index("name")["role"],
"susC_alone_hits": {g: len(df) for g, df in hitsA.items()}})
| role | susC_alone_hits | |
|---|---|---|
| Bacteroides_thetaiotaomicron_VPI5482 | positive | 69 |
| Bacteroides_fragilis_NCTC9343 | positive | 46 |
| Flavobacterium_johnsoniae_UW101 | positive | 32 |
| Pseudomonas_aeruginosa_PAO1 | negative | 6 |
| Escherichia_coli_MG1655 | negative | 0 |
| Bacillus_subtilis_168 | negative | 0 |
susC is promiscuous: it scores dozens of hits in the Bacteroidota, but it also lights
up in Pseudomonas aeruginosa — whose genome encodes many SusC-like TonB-dependent receptors
that have nothing to do with polysaccharide-utilization loci. On these hits alone you could
not confidently annotate a real susC. That is the annotation ambiguity synteny will resolve.
5. Search B — the strict susC–susD tandem¶
>susC 0 >susD
> requires the (+) strand; 0 requires the two genes to be immediately adjacent — the
canonical PUL core.
hitsB = {g: synteny_search(g, ">susC 0 >susD", subdir="B_tandem_strict")
for g in genomes["name"]}
pd.DataFrame({"role": genomes.set_index("name")["role"],
"tandem_member_hits": {g: len(df) for g, df in hitsB.items()},
"= susC-susD pairs": {g: len(df) // 2 for g, df in hitsB.items()}})
| role | tandem_member_hits | = susC-susD pairs | |
|---|---|---|---|
| Bacteroides_thetaiotaomicron_VPI5482 | positive | 48 | 24 |
| Bacteroides_fragilis_NCTC9343 | positive | 34 | 17 |
| Flavobacterium_johnsoniae_UW101 | positive | 24 | 12 |
| Pseudomonas_aeruginosa_PAO1 | negative | 0 | 0 |
| Escherichia_coli_MG1655 | negative | 0 | 0 |
| Bacillus_subtilis_168 | negative | 0 | 0 |
The tandem is specific. Every Bacteroidota genome yields many adjacent susC–susD
pairs; Pseudomonas — which had 6 lone-susC hits — drops to zero, because none of its
TonB receptors sits next to a susD. Synteny converted an ambiguous family hit into a confident
PUL call. Here are the first few pairs in Bacteroides thetaiotaomicron, alternating
susC,susD along the contig:
hitsB["Bacteroides_thetaiotaomicron_VPI5482"].sort_values("gene_number")[DISPLAY_COLS].head(8)
| contig | gene_number | locus | strand | gene_symbol | |
|---|---|---|---|---|---|
| 0 | NC_004663.1 | 32 | (22906, 26169) | pos | susC |
| 24 | NC_004663.1 | 33 | (26182, 28011) | pos | susD |
| 1 | NC_004663.1 | 144 | (139982, 143212) | pos | susC |
| 25 | NC_004663.1 | 145 | (143223, 144962) | pos | susD |
| 2 | NC_004663.1 | 197 | (200124, 203564) | pos | susC |
| 26 | NC_004663.1 | 198 | (203588, 204979) | pos | susD |
| 3 | NC_004663.1 | 213 | (222704, 225646) | pos | susC |
| 27 | NC_004663.1 | 214 | (225665, 227200) | pos | susD |
6. Search C — the tandem on either strand: the full PUL repertoire¶
susC 0 susD # no strand symbols, unordered=True
The strict (+)-strand search only catches the roughly half of PULs encoded on the forward
strand. Dropping the strand and order constraints recovers them all — revealing how many
SusC–SusD loci each genome really carries.
hitsC = {g: synteny_search(g, "susC 0 susD", unordered=True, subdir="C_tandem_any_strand")
for g in genomes["name"]}
# Count PULs as the number of susC loci recovered (each PUL contributes one susC).
def n_susC(df):
return int((df["gene_symbol"] == "susC").sum()) if len(df) else 0
pd.DataFrame({
"role": genomes.set_index("name")["role"],
"(+)-strand PULs": {g: len(df) // 2 for g, df in hitsB.items()},
"PULs (any strand)": {g: n_susC(df) for g, df in hitsC.items()},
}).reindex(genomes["name"])
| role | (+)-strand PULs | PULs (any strand) | |
|---|---|---|---|
| name | |||
| Bacteroides_thetaiotaomicron_VPI5482 | positive | 24 | 49 |
| Bacteroides_fragilis_NCTC9343 | positive | 17 | 37 |
| Flavobacterium_johnsoniae_UW101 | positive | 12 | 25 |
| Pseudomonas_aeruginosa_PAO1 | negative | 0 | 0 |
| Escherichia_coli_MG1655 | negative | 0 | 0 |
| Bacillus_subtilis_168 | negative | 0 | 0 |
Roughly double the pairs appear once both strands count: ~48 PULs in B. thetaiotaomicron, ~36 in B. fragilis, ~25 in Flavobacterium. The gut Bacteroides — which must process the huge diversity of dietary and host glycans — carry the largest PUL arsenals, exactly as the biology predicts. Every non-Bacteroidota control stays firmly at zero.
7. Summary¶
summary = pd.DataFrame({
"role": genomes.set_index("name")["role"],
"susC_alone": {g: len(df) for g, df in hitsA.items()},
"tandem_strict (members)": {g: len(df) for g, df in hitsB.items()},
"PULs (any strand)": {g: n_susC(df) for g, df in hitsC.items()},
}).reindex(genomes["name"])
summary
| role | susC_alone | tandem_strict (members) | PULs (any strand) | |
|---|---|---|---|---|
| name | ||||
| Bacteroides_thetaiotaomicron_VPI5482 | positive | 69 | 48 | 49 |
| Bacteroides_fragilis_NCTC9343 | positive | 46 | 34 | 37 |
| Flavobacterium_johnsoniae_UW101 | positive | 32 | 24 | 25 |
| Pseudomonas_aeruginosa_PAO1 | negative | 6 | 0 | 0 |
| Escherichia_coli_MG1655 | negative | 0 | 0 | 0 |
| Bacillus_subtilis_168 | negative | 0 | 0 | 0 |
Takeaways¶
- Single-gene hits can mislead.
susCbelongs to a widespread TonB-receptor family (SusC/RagA/OmpW); a lone hit — e.g. the six in Pseudomonas — does not mean "polysaccharide utilization locus". - Synteny supplies the missing specificity. The adjacent
susC–susDtandem is found in the Bacteroidota and nowhere in the controls, turning an ambiguous family match into a confident functional call. - Tune strand to the question.
>susC 0 >susDcharacterises the (+)-strand loci precisely; dropping the strand symbols (susC 0 susD,unordered) inventories the whole PUL repertoire — which scales with lifestyle (gut Bacteroides ≫ environmental Flavobacterium).
Reproducible CLI equivalent: run_case_study.sh ·
Data provenance & full write-up: README.md. HMMs from NCBI PGAP
(susC = TIGR04056.1, susD = NF033071.0); genomes from NCBI RefSeq (accessions in
genomes.tsv). This case study generalises the marine-metagenome Sus example from the Pynteny
docs to a curated, reproducible genome panel.