Finding the nitrogen-fixation (nif) operon with Pynteny¶

Motivation: GitHub issue — "This could be great for identifying operons. Do you have any examples of using this tool to identify common operons? For example, nif operons or others involved in nitrogen fixation."

Biological nitrogen fixation — reducing atmospheric N₂ to ammonia — is carried out by nitrogenase, encoded by the compact, strongly conserved nifH–nifD–nifK operon. That co-located, co-oriented gene cluster is exactly the signal Pynteny is built to detect.

In this notebook we will:

screen six genomes (three diazotrophs, three non-fixers) for the nif operon;
learn to tune the three knobs of a synteny search — strand, gene spacing, and order;
meet three real-world complications along the way: paralogous HMMs, an operon on a megaplasmid, and an operon split by an 11-kb excision element.

Everything runs from the small inputs committed next to this notebook — the six nif HMMs and their metadata under data/ — so you do not need the full 432 MB PGAP database.

Setup¶

In [1]:

Copied!





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

# Resolve paths relative to this notebook
BASE     = Path.cwd()
DATA     = BASE / "data"
HMM_DIR  = DATA / "hmms"
HMM_META = DATA / "nif_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)

# Keep the notebook output tidy: 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)
print("HMMs committed for this case study:")
print("\n".join(p.name for p in sorted(HMM_DIR.glob('*.HMM'))))
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

# Resolve paths relative to this notebook
BASE     = Path.cwd()
DATA     = BASE / "data"
HMM_DIR  = DATA / "hmms"
HMM_META = DATA / "nif_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)

# Keep the notebook output tidy: 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)
print("HMMs committed for this case study:")
print("\n".join(p.name for p in sorted(HMM_DIR.glob('*.HMM'))))

HMMs committed for this case study:
TIGR01282.1.HMM
TIGR01283.1.HMM
TIGR01285.1.HMM
TIGR01286.1.HMM
TIGR01287.1.HMM
TIGR01290.1.HMM

1. The biology: the nitrogenase gene panel¶

We use six HMMs from the NCBI PGAP/TIGRFAM collection — the three structural genes plus three FeMo-cofactor biosynthesis genes. The metadata file maps each HMM accession to a gene symbol, which lets us write the synteny structure in readable nifH/nifD/nifK terms (via --gene_ids) instead of raw accessions.

In [2]:

Copied!

meta = pd.read_csv(HMM_META, sep="\t")
meta
meta = pd.read_csv(HMM_META, sep="\t")
meta

Out[2]:

	#ncbi_accession	label	product_name	gene_symbol	ec_numbers
0	TIGR01287.1	nifH	nitrogenase iron protein	nifH	1.18.6.1
1	TIGR01282.1	nifD	nitrogenase molybdenum-iron protein alpha chain	nifD	1.18.6.1
2	TIGR01286.1	nifK	nitrogenase molybdenum-iron protein subunit beta	nifK	1.18.6.1
3	TIGR01283.1	nifE	nitrogenase iron-molybdenum cofactor biosynthesis protei...	nifE	NaN
4	TIGR01285.1	nifN	nitrogenase iron-molybdenum cofactor biosynthesis protei...	nifN	NaN
5	TIGR01290.1	nifB	nitrogenase cofactor biosynthesis protein NifB	nifB	NaN

⚠️ Built-in gotcha — paralogy. nifE/nifN are ancient duplications of the structural genes nifD/nifK, so their HMMs cross-hit: one peptide is often matched by several models. Throughout this notebook we pass best_hmm_wins=True, which keeps each peptide's single highest-scoring HMM and stops the cross-hits from silently breaking the synteny filter.

2. The genome panel¶

Three diazotrophs from three different phyla, and three non-fixers — including one deliberately counter-intuitive negative.

In [3]:

Copied!

genomes = (pd.read_csv(BASE / "genomes.tsv", sep="\t")
             .rename(columns=lambda c: c.lstrip("#")))
genomes[["name", "role", "lineage", "lifestyle", "nif_expectation"]]
genomes = (pd.read_csv(BASE / "genomes.tsv", sep="\t")
             .rename(columns=lambda c: c.lstrip("#")))
genomes[["name", "role", "lineage", "lifestyle", "nif_expectation"]]

Out[3]:

	name	role	lineage	lifestyle	nif_expectation
0	Azotobacter_vinelandii_DJ	positive	Gammaproteobacteria	free-living aerobe	canonical nifHDK + nifENB cluster
1	Sinorhizobium_meliloti_1021	positive	Alphaproteobacteria	legume symbiont	nifHDKE on the pSymA megaplasmid
2	Nostoc_PCC7120	positive	Cyanobacteria	heterocyst-forming	nifHDK interrupted by the 11-kb nifD excision element
3	Escherichia_coli_MG1655	negative	Gammaproteobacteria	enteric	no nitrogenase
4	Bacillus_subtilis_168	negative	Bacillota (Firmicutes)	soil saprophyte	no nitrogenase
5	Klebsiella_pneumoniae_342	negative	Gammaproteobacteria	endophyte/opportunist	this strain lacks nif (N-fixation is strain-specific)

3. Build labelled peptide databases¶

pynteny build translates every ORF (with Prodigal) and writes a FASTA whose record labels encode each gene's contig, position and strand — the positional information the synteny filter later reads. We download each genome from NCBI RefSeq and build it (skipping any already built, so re-running the notebook is cheap).

In [4]:

Copied!





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']:<34} {n:>5} ORFs")
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']:<34} {n:>5} ORFs")

Azotobacter_vinelandii_DJ           4856 ORFs
Sinorhizobium_meliloti_1021         6336 ORFs
Nostoc_PCC7120                      6199 ORFs
Escherichia_coli_MG1655             4319 ORFs
Bacillus_subtilis_168               4242 ORFs
Klebsiella_pneumoniae_342           4668 ORFs

4. Search A — the strict structural operon¶

>nifH 0 >nifD 0 >nifK

> requires the gene on the (+) strand; 0 means the genes must be immediately adjacent (zero genes between). This is the most specific possible query for the canonical nitrogenase operon — it demands exact collinearity.

In [5]:

Copied!





DISPLAY_COLS = ["contig", "gene_number", "strand", "gene_symbol", "locus", "product"]

def synteny_search(genome, struct, *, unordered=False, best_hmm_wins=True, 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, best_hmm_wins=best_hmm_wins,
            hmm_dir=HMM_DIR, hmm_meta=HMM_META, outdir=outdir,
        ).run()
        return hits.hits
    except SystemExit:           # Pynteny exits when an HMM matches nothing -> a negative
        return pd.DataFrame()

STRICT = ">nifH 0 >nifD 0 >nifK"
hitsA = {g: synteny_search(g, STRICT, subdir="A_strict_nifHDK") for g in genomes["name"]}

pd.DataFrame({"role": genomes.set_index("name")["role"],
              "A_strict_nifHDK_hits": {g: len(df) for g, df in hitsA.items()}})
DISPLAY_COLS = ["contig", "gene_number", "strand", "gene_symbol", "locus", "product"]

def synteny_search(genome, struct, *, unordered=False, best_hmm_wins=True, 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, best_hmm_wins=best_hmm_wins,
            hmm_dir=HMM_DIR, hmm_meta=HMM_META, outdir=outdir,
        ).run()
        return hits.hits
    except SystemExit:           # Pynteny exits when an HMM matches nothing -> a negative
        return pd.DataFrame()

STRICT = ">nifH 0 >nifD 0 >nifK"
hitsA = {g: synteny_search(g, STRICT, subdir="A_strict_nifHDK") for g in genomes["name"]}

pd.DataFrame({"role": genomes.set_index("name")["role"],
              "A_strict_nifHDK_hits": {g: len(df) for g, df in hitsA.items()}})

Out[5]:

	role	A_strict_nifHDK_hits
Azotobacter_vinelandii_DJ	positive	3
Sinorhizobium_meliloti_1021	positive	3
Nostoc_PCC7120	positive	0
Escherichia_coli_MG1655	negative	0
Bacillus_subtilis_168	negative	0
Klebsiella_pneumoniae_342	negative	0

Two of the three diazotrophs light up; every non-fixer stays dark. Here is the operon in Azotobacter vinelandii — three genes, same strand, perfectly adjacent:

In [6]:

Copied!

hitsA["Azotobacter_vinelandii_DJ"][DISPLAY_COLS]
hitsA["Azotobacter_vinelandii_DJ"][DISPLAY_COLS]

Out[6]:

	contig	gene_number	strand	gene_symbol	locus	product
0	NC_012560.1	129	pos	nifH	(136759, 137631)	nifH
1	NC_012560.1	130	pos	nifD	(137758, 139236)	nifD
2	NC_012560.1	131	pos	nifK	(139337, 140908)	nifK

Sinorhizobium meliloti matches identically — but look at the contig column: the operon lives on NC_003037.1, the pSymA megaplasmid, not the chromosome. Pynteny does not care which replicon a syntenic block sits on.

In [7]:

Copied!

hitsA["Sinorhizobium_meliloti_1021"][DISPLAY_COLS]
hitsA["Sinorhizobium_meliloti_1021"][DISPLAY_COLS]

Out[7]:

	contig	gene_number	strand	gene_symbol	locus	product
0	NC_003037.1	489	pos	nifH	(453556, 454449)	nifH
1	NC_003037.1	490	pos	nifD	(454549, 456051)	nifD
2	NC_003037.1	491	pos	nifK	(456142, 457683)	nifK

Nostoc scored 0 here. That is not a bug — its operon is on the (−) strand and is interrupted, so it cannot match a strict (+)-strand, zero-gap query. To find it we relax the search.

5. Search B — the robust, order- & strand-agnostic panel¶

nifB 80 nifH 80 nifD 80 nifK 80 nifE 80 nifN

No strand symbols (strand-agnostic) and unordered=True (order-agnostic) ask only: "are these six genes clustered within ~80 ORFs of one another, in any arrangement?" This is the recommended screen for "does this genome fix nitrogen?".

In [8]:

Copied!





PANEL = "nifB 80 nifH 80 nifD 80 nifK 80 nifE 80 nifN"
hitsB = {g: synteny_search(g, PANEL, unordered=True, subdir="B_panel_unordered")
         for g in genomes["name"]}

pd.DataFrame({"role": genomes.set_index("name")["role"],
              "B_panel_unordered_hits": {g: len(df) for g, df in hitsB.items()}})
PANEL = "nifB 80 nifH 80 nifD 80 nifK 80 nifE 80 nifN"
hitsB = {g: synteny_search(g, PANEL, unordered=True, subdir="B_panel_unordered")
         for g in genomes["name"]}

pd.DataFrame({"role": genomes.set_index("name")["role"],
              "B_panel_unordered_hits": {g: len(df) for g, df in hitsB.items()}})

Out[8]:

	role	B_panel_unordered_hits
Azotobacter_vinelandii_DJ	positive	12
Sinorhizobium_meliloti_1021	positive	6
Nostoc_PCC7120	positive	7
Escherichia_coli_MG1655	negative	0
Bacillus_subtilis_168	negative	0
Klebsiella_pneumoniae_342	negative	0

All three diazotrophs are now recovered (and A. vinelandii's two nitrogenase copies give 12 hits), while every non-fixer is still rejected. Here is the Nostoc neighbourhood, sorted along the contig:

In [9]:

Copied!

nostoc = hitsB["Nostoc_PCC7120"].sort_values("gene_number")
nostoc[DISPLAY_COLS]
nostoc = hitsB["Nostoc_PCC7120"].sort_values("gene_number")
nostoc[DISPLAY_COLS]

Out[9]:

	contig	gene_number	strand	gene_symbol	locus	product
2	NC_003272.1	1457	neg	nifB	(1694529, 1694942)	nifB
4	NC_003272.1	1458	neg	nifN	(1695055, 1696389)	nifN
6	NC_003272.1	1459	neg	nifE	(1696389, 1697831)	nifE
8	NC_003272.1	1461	neg	nifK	(1698743, 1700281)	nifK
10	NC_003272.1	1475	neg	nifD	(1711821, 1713263)	nifD
0	NC_003272.1	1476	neg	nifH	(1713396, 1714283)	nifH
3	NC_003272.1	1539	neg	nifB	(1776670, 1778097)	nifB

The Nostoc twist — a 55-year-old discovery, visible in the coordinates¶

Reading the (−) strand operon 5'→3' (high gene number to low): nifH → nifD → big gap → nifK → nifE → nifN → nifB. nifD and nifK are not adjacent — they are separated by ~13 ORFs. That gap is the famous 11-kb nifD excision element: in vegetative cells it interrupts the operon, and the recombinase XisA excises it only during heterocyst differentiation, reconstituting a functional nifHDK. The genome was sequenced from vegetative DNA, so the interruption is right there in the gene coordinates.

In [10]:

Copied!





def span(locus):
    # Return (start, end) from a locus value that may be a tuple or a "(start, end)" string.
    if isinstance(locus, str):
        a, b = locus.strip("() ").split(",")
        return int(a), int(b)
    return int(locus[0]), int(locus[1])

g = nostoc.set_index("gene_symbol")
nifK_end   = span(g.loc["nifK", "locus"])[1]
nifD_start = span(g.loc["nifD", "locus"])[0]
print(f"nifK ends at   {nifK_end:>9,} bp")
print(f"nifD starts at {nifD_start:>9,} bp")
print(f"intergenic gap = {nifD_start - nifK_end:,} bp  (~the 11-kb nifD excision element)")
def span(locus):
    # Return (start, end) from a locus value that may be a tuple or a "(start, end)" string.
    if isinstance(locus, str):
        a, b = locus.strip("() ").split(",")
        return int(a), int(b)
    return int(locus[0]), int(locus[1])

g = nostoc.set_index("gene_symbol")
nifK_end   = span(g.loc["nifK", "locus"])[1]
nifD_start = span(g.loc["nifD", "locus"])[0]
print(f"nifK ends at   {nifK_end:>9,} bp")
print(f"nifD starts at {nifD_start:>9,} bp")
print(f"intergenic gap = {nifD_start - nifK_end:,} bp  (~the 11-kb nifD excision element)")

nifK ends at   1,700,281 bp
nifD starts at 1,711,821 bp
intergenic gap = 11,540 bp  (~the 11-kb nifD excision element)

6. Search C — recovering Nostoc as an ordered operon¶

Once we tell Pynteny the truth — (−) strand (<) and a nifK–nifD gap wide enough to span the element — the operon reappears as a clean collinear block:

<nifK 15 <nifD 0 <nifH

In [11]:

Copied!

nostoc_C = synteny_search("Nostoc_PCC7120", "<nifK 15 <nifD 0 <nifH",
                          subdir="C_nostoc_strand_tuned")
nostoc_C[DISPLAY_COLS]
nostoc_C = synteny_search("Nostoc_PCC7120", "<nifK 15 <nifD 0 <nifH",
                          subdir="C_nostoc_strand_tuned")
nostoc_C[DISPLAY_COLS]

Out[11]:

	contig	gene_number	strand	gene_symbol	locus	product
1	NC_003272.1	1461	neg	nifK	(1698743, 1700281)	nifK
2	NC_003272.1	1475	neg	nifD	(1711821, 1713263)	nifD
0	NC_003272.1	1476	neg	nifH	(1713396, 1714283)	nifH

7. Summary¶

In [12]:

Copied!





summary = pd.DataFrame({
    "role":              genomes.set_index("name")["role"],
    "A_strict_nifHDK":   {g: len(df) for g, df in hitsA.items()},
    "B_panel_unordered": {g: len(df) for g, df in hitsB.items()},
}).reindex(genomes["name"])
summary
summary = pd.DataFrame({
    "role":              genomes.set_index("name")["role"],
    "A_strict_nifHDK":   {g: len(df) for g, df in hitsA.items()},
    "B_panel_unordered": {g: len(df) for g, df in hitsB.items()},
}).reindex(genomes["name"])
summary

Out[12]:

	role	A_strict_nifHDK	B_panel_unordered
name
Azotobacter_vinelandii_DJ	positive	3	12
Sinorhizobium_meliloti_1021	positive	3	6
Nostoc_PCC7120	positive	0	7
Escherichia_coli_MG1655	negative	0	0
Bacillus_subtilis_168	negative	0	0
Klebsiella_pneumoniae_342	negative	0	0

Every diazotroph detected, every non-fixer rejected.

Takeaways¶

Synteny beats independent gene hits. Requiring nifH–nifD–nifK co-located is far more specific than three separate HMM searches — the difference between "has a nitrogenase-like gene" and "has a nitrogenase operon".
Tune three knobs to the biology: strand (> < / none), max gene spacing (the integers), and order (unordered). Start permissive to screen, tighten to characterise.
Mind paralogues. With cross-hitting models like nifD/nifE and nifK/nifN, use best_hmm_wins=True.
Presence is strain-specific, not genus-level. Klebsiella pneumoniae 342 carries no nif even though the genus is a textbook nitrogen fixer — only a genome-level search settles it.

Reproducible CLI equivalent: run_case_study.sh · Data provenance & full write-up: README.md. HMMs from NCBI PGAP; genomes from NCBI RefSeq (accessions in genomes.tsv). The Nostoc nifD element: Golden, Robinson & Haselkorn (1985) Nature 314, 419–423.