MIToS-generated mmCIF templates fail in ColabFold/AlphaFold due to missing PDBx categories

[diegozea]

Summary

When generating a template mmCIF from a PDB using MIToS (e.g. selecting a single chain, then writing with MMCIFFile / write_file), the resulting mmCIF can be syntactically valid but still fails in ColabFold/AlphaFold template mode with:

ValueError: mmCIF file ... is missing required field _chem_comp.id

Root cause

AlphaFold’s mmCIF parser (alphafold/data/mmcif_parsing.py) expects several PDBx/mmCIF categories that are commonly present in wwPDB mmCIF files, notably:

• _chem_comp.id and _chem_comp.type (used to classify monomers as “peptide”)
• _entity_poly_seq.* (polymer sequence)
• _struct_asym.* (mapping between entity IDs and chain IDs)
• plus header fields like _entry.id and _exptl.method

MIToS' current mmCIF writing path (via BioStructures.MMCIFDict(residues; ...) and writemmcif) primarily emits _atom_site.* information. In addition, when starting from PDB input, residue identifiers often have empty PDBe numbering, which may lead to _atom_site.label_seq_id being written as "." rather than an integer. AlphaFold’s parser later casts label/auth sequence IDs to integers, so this can be another failure point after _chem_comp is addressed.

Expected behavior

It would be very helpful if MIToS documentation (or helper utilities) could:

1. warn that “atom_site-only” mmCIFs may not be sufficient for AlphaFold/ColabFold templating, and/or
2. provide a small utility to “upgrade” a minimal mmCIF to include the minimal additional PDBx categories required by AlphaFold.

Repro (typical)

• Start from a PDB
• Use MIToS to select a chain / rename residues
• Save as mmCIF using MMCIFFile + write_file
• Use that mmCIF as a ColabFold template → error above

[diegozea]

Workaround

Patch MIToS mmCIF into an AlphaFold-friendly template mmCIF ( @ZulieD )

This workaround reads the MIToS-generated mmCIF (typically missing _chem_comp, _entity_poly_seq, _struct_asym, etc.), then:

• assigns a sequential PDBe residue number per chain so _atom_site.label_seq_id is numeric
• merges label=true and label=false mmCIF dict outputs to include both label_* and auth_* atom_site identifiers
• adds minimal AlphaFold-required PDBx categories: _entry, _exptl, _struct_asym, _entity_poly_seq, _chem_comp
• uses MIToS.Utils.THREE2ONE to decide if a residue name is “peptide-like” (otherwise mark as non-polymer), with an option to force all monomers to peptide

using MIToS.PDB
using BioStructures
using MIToS.Utils: THREE2ONE

# Parse numeric part of PDB residue number like "34A" -> 34
_auth_num(resnum::AbstractString) = something(
    tryparse(Int, replace(strip(resnum), r"[A-Za-z]" => "")),
    typemax(Int),
)

# Parse insertion code like "34A" -> 'A', "34" -> ' '
_ins_char(resnum::AbstractString) = begin
    m = match(r"[A-Za-z]$", strip(resnum))
    m === nothing ? ' ' : m.match[1]
end

# Treat residue as peptide-like if it appears in MIToS' 3-letter mapping table
_is_peptide_resname(resname::AbstractString) =
    haskey(THREE2ONE, uppercase(strip(resname)))

# Best-effort element guess if missing
function _guess_element(atomname::AbstractString)
    s = uppercase(strip(atomname))
    isempty(s) && return "?"
    if s[1] in ('C','N','O','S','P','H')
        return string(s[1])
    end
    if length(s) ≥ 2
        two = s[1:2]
        if two in ("ZN","FE","MG","NA","CL","BR","MN","CO","NI","CU","CD","HG","PB","SR","CS","BA")
            return two
        end
    end
    return string(s[1])
end

"""
    patch_mmcif_for_alphafold(infile, outfile;
        entry_id=nothing,
        exptl_method="computational model",
        force_peptide_types=false,
        fill_missing_elements=true)

Patch a MIToS-generated mmCIF (often atom_site-only) into a minimal PDBx/mmCIF
that AlphaFold/ColabFold can use as a template.

Returns `outfile`.
"""
function patch_mmcif_for_alphafold(infile::AbstractString, outfile::AbstractString;
    entry_id::Union{Nothing,String}=nothing,
    exptl_method::AbstractString="computational model",
    force_peptide_types::Bool=false,
    fill_missing_elements::Bool=true,
)
    # 1) Read residues from mmCIF
    residues = read_file(infile, MMCIFFile; label=true)

    # Sort residues for stable numbering/sequence tables
    sort!(residues, by = r -> (
        tryparse(Int, r.id.model) === nothing ? 1 : parse(Int, r.id.model),
        r.id.chain, _auth_num(r.id.number), _ins_char(r.id.number)
    ))

    # 2) Assign PDBe_number sequentially per chain -> numeric label_seq_id
    chain_counter = Dict{String,Int}()
    for r in residues
        c = r.id.chain
        n = get!(chain_counter, c, 0) + 1
        chain_counter[c] = n

        r.id = PDBResidueIdentifier(
            string(n),      # PDBe_number (drives label_seq_id)
            r.id.number,    # PDB number (auth_seq_id)
            r.id.name, r.id.group, r.id.model, r.id.chain
        )

        if fill_missing_elements
            r.atoms = [
                isempty(a.element) ?
                    PDBAtom(a.coordinates, a.atom, _guess_element(a.atom),
                           a.occupancy, a.B, a.alt_id, a.charge) :
                    a
                for a in r.atoms
            ]
        end
    end

    # 3) Build mmCIFDict twice (label/auth) and merge so both sets exist
    mm_label = BioStructures.MMCIFDict(residues; label=true)
    mm_auth  = BioStructures.MMCIFDict(residues; label=false)
    for (k, v) in mm_auth
        if !haskey(mm_label, k)
            mm_label[k] = v
        end
    end
    mm = mm_label

    # 4) Add minimal non-atom_site categories AlphaFold expects
    eid = entry_id === nothing ? splitext(basename(infile))[1] : entry_id
    mm["_entry.id"] = [eid]
    mm["_exptl.method"] = [String(exptl_method)]

    chains = unique([r.id.chain for r in residues])
    chain2entity = Dict{String,String}(c => string(i) for (i, c) in enumerate(chains))
    mm["_struct_asym.id"] = chains
    mm["_struct_asym.entity_id"] = [chain2entity[c] for c in chains]

    # entity_poly_seq
    ent = String[]; num = String[]; mon = String[]
    for c in chains
        seqnum = 0
        for r in residues
            r.id.chain == c || continue
            seqnum += 1
            push!(ent, chain2entity[c])
            push!(num, string(seqnum))
            push!(mon, r.id.name)
        end
    end
    mm["_entity_poly_seq.entity_id"] = ent
    mm["_entity_poly_seq.num"] = num
    mm["_entity_poly_seq.mon_id"] = mon

    # chem_comp
    monomers = unique(mon)
    sort!(monomers)
    mm["_chem_comp.id"] = monomers
    mm["_chem_comp.type"] = [
        (force_peptide_types || _is_peptide_resname(m)) ?
            "L-peptide linking" : "non-polymer"
        for m in monomers
    ]
    mm["_chem_comp.name"] = monomers  # optional

    # 5) Write patched CIF
    open(outfile, "w") do io
        BioStructures.writemmcif(io, mm)
    end

    return outfile
end

Usage example:

patch_mmcif_for_alphafold("input.cif", "patched.cif"; force_peptide_types=true)

(Setting force_peptide_types=true is useful if residue names are nonstandard but still peptide-like.)