import sys,os import json import string import traceback from optparse import OptionParser import requests import glob import random import urllib3 urllib3.disable_warnings() #import util import subprocess from collections import Counter from datetime import datetime BASE_DIR = "/data/projects/glygen/downloads/rcsb_pdb/temp" STEP_2_DIR = BASE_DIR + "/step_2" STEP_3_DIR = BASE_DIR + "/step_3" # Working folder: only the auth-mismatch logs are written here now STEP_4_DIR = BASE_DIR + "/step_4" # Final datasets go into a dated release folder (MM_DD_YYYY, e.g. 06_11_2026) # inside RELEASES_DIR, and RELEASES_DIR/current is a symlink to the latest one. RELEASES_DIR = "/data/projects/glygen/downloads/rcsb_pdb" FINAL_FILE_SUFFIX = "_glycosites_rcsb_pdb.csv" # _glycosites_rcsb_pdb.csv def lookup_mapping(sifts_map, pdb_id, label_chain, label_pos, auth_chain): """ SIFTS' by_chain_pos index is keyed on ".|" -- SIFTS organizes residues by its own PDBe/label numbering (equivalent to mmCIF label_seq_id), with the auth chain/position attached as attributes on each residue. So the correct exact-match key combines the auth CHAIN with the LABEL position, not an auth position -- there's no reliable auth position available from RCSB to use here (see note in step_2: connect_target's auth_seq_id turned out to report the glycan's own residue number, not the protein's, so it isn't used at all). The auth position for output comes entirely from the matched SIFTS entry itself (entry["auth_resnum"]). Tries, in order: 1. auth_chain (from step_2's entity-level asym_ids/auth_asym_ids crosswalk) + label_pos -- "auth_exact" 2. label_chain + label_pos -- "label_exact" (covers cases where the auth-chain crosswalk was unavailable, or SIFTS' own auth chain happens to equal the label chain anyway) 3. label_pos alone, across every auth chain copy of the entity -- "position_fallback" -- only accepted when every candidate at that position agrees on the UniProt residue, to avoid guessing between genuinely different chain copies of a homo-oligomer. Returns (list_of_entries, match_type) where match_type is one of "auth_exact", "label_exact", "position_fallback", or None. """ if auth_chain: auth_key = "%s.%s|%s" % (pdb_id, auth_chain, label_pos) if auth_key in sifts_map["by_chain_pos"]: return sifts_map["by_chain_pos"][auth_key], "auth_exact" label_key = "%s.%s|%s" % (pdb_id, label_chain, label_pos) if label_key in sifts_map["by_chain_pos"]: return sifts_map["by_chain_pos"][label_key], "label_exact" pos_key = "%s|%s" % (pdb_id, label_pos) candidates = sifts_map["by_pos_only"].get(pos_key) if not candidates: return None, None uniprot_targets = set((c.get("uniprot_ac"), c.get("uniprot_resnum")) for c in candidates) if len(uniprot_targets) == 1: return candidates, "position_fallback" return None, None # genuinely ambiguous between chain copies -- skip # Synthetic-site rules (see synthetic_reason) use this table: # standard amino acids, PDB 3-letter code -> UniProt 1-letter code. # Used to skip glycosylation sites where the PDB residue is NOT the residue UniProt # has at that position (engineered mutations, variants, conflicts): those are real # in that structure but are not true glycosylation sites of the UniProt protein. # PDB residues that are not standard amino acids (e.g. MSE, SEP) are never # skipped, because SIFTS maps modified residues to their parent amino acid. # Run with --keep-mismatches to switch the skip off. THREE_TO_ONE = { "ALA": "A", "ARG": "R", "ASN": "N", "ASP": "D", "CYS": "C", "GLN": "Q", "GLU": "E", "GLY": "G", "HIS": "H", "ILE": "I", "LEU": "L", "LYS": "K", "MET": "M", "PHE": "F", "PRO": "P", "SER": "S", "THR": "T", "TRP": "W", "TYR": "Y", "VAL": "V", } SKIP_MISMATCHES = "--keep-mismatches" not in sys.argv def residue_mismatch(entry): """True if the PDB residue is a standard amino acid that differs from UniProt's.""" pdb_one = THREE_TO_ONE.get((entry.get("pdb_aa") or "").upper()) uni_one = (entry.get("uniprot_aa") or "").upper() return bool(pdb_one and uni_one and pdb_one != uni_one) def synthetic_reason(entry, glycosylation_type): """ Return why this site looks like a synthetic (engineered) glycosylation site, or "" if it looks genuine. residue_mismatch the glycosylated PDB residue is not the UniProt residue engineered_sequon N-linked: the PDB construct has an N-X-S/T sequon that UniProt's sequence does not have, because - UniProt's S/T slot (two on) is not S/T/C (S/T mutated in), or - UniProt has Pro as X (N-P-S/T is not a sequon) and the construct does not (Pro mutated out) sequon_not_in_uniprot N-linked: the S/T slot (or the X) is not mapped to this UniProt entry at all (it sits in a tag, linker or fusion partner), so the sequon does not exist in the real protein The neighbour fields come from step_3 (SIFTS, unobserved residues included). Maps without neighbour fields are never judged on these rules. """ if residue_mismatch(entry): return "residue_mismatch" if glycosylation_type == "N-linked" and entry.get("uniprot_ac") and entry.get("n2_state"): pdb_n2 = THREE_TO_ONE.get((entry.get("n2_pdb_aa") or "").upper(), "") if pdb_n2 in ("S", "T"): if entry.get("n1_state") == "unmapped" or entry.get("n2_state") == "unmapped": return "sequon_not_in_uniprot" uni_n2 = (entry.get("n2_uniprot_aa") or "").upper() if entry.get("n2_state") == "mapped" and uni_n2 and uni_n2 not in ("S", "T", "C"): return "engineered_sequon" pdb_n1 = THREE_TO_ONE.get((entry.get("n1_pdb_aa") or "").upper(), "") uni_n1 = (entry.get("n1_uniprot_aa") or "").upper() if entry.get("n1_state") == "mapped" and uni_n1 == "P" and pdb_n1 and pdb_n1 != "P": return "engineered_sequon" return "" FINAL_HEADER = [ "uniprotkb_ac", "glycosylation_site_uniprotkb", "glycosylation_type", "amino_acid_chain", "glycosylation_site_pdb", "amino_acid", "saccharide", "pdb_glycan_entity_id", "pdb_extended_id", "xref_id", "xref_key", "src_xref_key", "src_xref_id", "amino_acid_chain_auth", "glycosylation_site_auth", "carb_name", "host_organism", "host_org_cell_line", "host_org_cell", "host_org_organ", "notes", ] # Flagged in the "notes" column but KEPT (the other reasons are dropped). NOTE_ONLY = { "sequon_not_in_uniprot": "Possibly an engineered N-sequon site", } def update_current_symlink(release_name): """ Point RELEASES_DIR/current at the release folder just written (the equivalent of `rm current; ln -s $version current`). The link target is relative (just the folder name), so it keeps working if the whole directory is moved. The new link is created under a temporary name and renamed over the old one, so `current` is never missing or half-updated while something is reading it. If `current` exists but is a real file/folder rather than a symlink, it is left alone. """ link = os.path.join(RELEASES_DIR, "current") if os.path.lexists(link) and not os.path.islink(link): print("WARNING: %s exists and is not a symlink -- not touching it, 'current' NOT updated." % link) return False tmp_link = "%s.tmp_%d" % (link, os.getpid()) if os.path.lexists(tmp_link): os.remove(tmp_link) os.symlink(release_name, tmp_link) os.replace(tmp_link, link) print("Updated symlink: %s -> %s" % (link, release_name)) return True def extract_data(tax_name, release_dir): """Returns True if the final dataset for this species was written.""" map_file = "%s/%s.sifts_map.json" % (STEP_3_DIR, tax_name) if not os.path.exists(map_file): print("No SIFTS map found for %s (missing %s), skipping." % (tax_name, map_file)) return False with open(map_file) as FR: sifts_map = json.load(FR) in_file = "%s/%s.csv" % (STEP_2_DIR, tax_name) out_file = "%s/%s%s" % (release_dir, tax_name, FINAL_FILE_SUFFIX) mismatch_file = "%s/%s.auth_mismatches.csv" % (STEP_4_DIR, tax_name) fallback_count = 0 mismatch_count = 0 checked_count = 0 aa_skipped_rows = 0 aa_skipped_sites = set() aa_file = "%s/%s.residue_mismatches.csv" % (STEP_4_DIR, tax_name) AW = open(aa_file, "w") AW.write("\"%s\"\n" % ("\",\"".join( ["src_xref_id", "amino_acid_chain", "glycosylation_site_pdb", "pdb_residue", "uniprotkb_ac", "glycosylation_site_uniprotkb", "uniprot_residue", "reason", "pdb_residue_plus1", "uniprot_residue_plus1", "plus1_state", "pdb_residue_plus2", "uniprot_residue_plus2", "plus2_state"]))) # Drop diagnostics (counted per UNIQUE site, not per output row) sifts_pdbs = set(k.split("|")[0] for k in sifts_map["by_pos_only"]) all_sites = set() kept_sites = set() dropped = {} # site -> reason dropped_file = "%s/%s.dropped_sites.csv" % (STEP_4_DIR, tax_name) FW = open(out_file, "w") FW.write("\"%s\"\n" % ("\",\"".join(FINAL_HEADER))) MW = open(mismatch_file, "w") MW.write("\"%s\"\n" % ("\",\"".join( ["src_xref_id", "amino_acid_chain", "glycosylation_site_pdb", "glycosylation_site_auth_rcsb", "glycosylation_site_auth_sifts", "match_type"]))) with open(in_file, "r") as FR: header_line = FR.readline().strip() step2_header = header_line[1:-1].split("\",\"") for line in FR: line = line.strip() if not line: continue row = line[1:-1].split("\",\"") row_dict = dict(zip(step2_header, row)) pdb_id = row_dict["src_xref_id"].strip().lower() label_chain = row_dict["amino_acid_chain"] label_pos = row_dict["glycosylation_site_pdb"] auth_chain = row_dict["amino_acid_chain_auth"] site_id = (pdb_id, label_chain, auth_chain, label_pos) all_sites.add(site_id) row_written = False entries, match_type = lookup_mapping(sifts_map, pdb_id, label_chain, label_pos, auth_chain) if not entries: if site_id not in kept_sites: dropped[site_id] = ("pdb_not_in_sifts" if pdb_id not in sifts_pdbs else "residue_not_found_or_ambiguous") continue if match_type == "position_fallback": fallback_count += 1 for entry in entries: uniprotkb_ac = entry.get("uniprot_ac") glycosylation_site_uniprotkb = entry.get("uniprot_resnum") glycosylation_site_auth = entry.get("auth_resnum") # SIFTS' own matched auth chain is authoritative; fall back to # step_2's RCSB-crosswalk value only if SIFTS didn't supply one. amino_acid_chain_auth = entry.get("auth_chain_id") or row_dict["amino_acid_chain_auth"] synth = synthetic_reason(entry, row_dict["glycosylation_type"]) note = NOTE_ONLY.get(synth, "") keep_with_note = bool(note) if row_dict.get("chemically_synthetic"): syn_note = ("protein entity is chemically synthesised (RCSB entity.src_method = syn); " "the glycosylation may be synthetic") note = syn_note if not note else note + "; " + syn_note if synth and (SKIP_MISMATCHES or keep_with_note): drop_it = SKIP_MISMATCHES and not keep_with_note if drop_it: aa_skipped_rows += 1 if site_id not in kept_sites: dropped.setdefault(site_id, synth) site_key = (uniprotkb_ac, glycosylation_site_uniprotkb, row_dict["glycosylation_type"]) if site_key not in aa_skipped_sites: aa_skipped_sites.add(site_key) AW.write("\"%s\"\n" % ("\",\"".join([ row_dict["src_xref_id"], label_chain, label_pos, entry.get("pdb_aa") or "", str(uniprotkb_ac), str(glycosylation_site_uniprotkb), entry.get("uniprot_aa") or "", synth, entry.get("n1_pdb_aa") or "", entry.get("n1_uniprot_aa") or "", entry.get("n1_state") or "", entry.get("n2_pdb_aa") or "", entry.get("n2_uniprot_aa") or "", entry.get("n2_state") or ""]))) if drop_it: continue if uniprotkb_ac is None or glycosylation_site_uniprotkb is None: # residue exists in the structure but isn't covered by a # UniProt segment in SIFTS -- nothing to report for this row if site_id not in kept_sites: dropped.setdefault(site_id, "no_uniprot_mapping") continue kept_sites.add(site_id) dropped.pop(site_id, None) # Cross-check: step_2 also carries an RCSB-sourced auth position # (connect_target.auth_seq_id, from the polymer-side struct_conn -- # unlike the branched-side block, this one is expected to be the # protein's real auth number). SIFTS remains the value actually # written to the output; this just flags any disagreement between # the two independent sources for later review, rather than # silently trusting either one. rcsb_auth = row_dict.get("glycosylation_site_auth", "") sifts_auth = "" if glycosylation_site_auth is None else str(glycosylation_site_auth) if rcsb_auth and sifts_auth: checked_count += 1 if rcsb_auth != sifts_auth: mismatch_count += 1 MW.write("\"%s\"\n" % ("\",\"".join([ row_dict["src_xref_id"], label_chain, label_pos, rcsb_auth, sifts_auth, match_type, ]))) newrow = [ str(uniprotkb_ac), str(glycosylation_site_uniprotkb), row_dict["glycosylation_type"], row_dict["amino_acid_chain"], row_dict["glycosylation_site_pdb"], row_dict["amino_acid"], row_dict["saccharide"], row_dict["pdb_glycan_entity_id"], row_dict["pdb_extended_id"], row_dict["xref_id"], row_dict["xref_key"], row_dict["src_xref_key"], row_dict["src_xref_id"], amino_acid_chain_auth, sifts_auth, row_dict["carb_name"], row_dict["host_organism"], row_dict["host_org_cell_line"], row_dict["host_org_cell"], row_dict["host_org_organ"], note, ] FW.write("\"%s\"\n" % ("\",\"".join(newrow))) FW.close() MW.close() AW.close() if aa_skipped_rows: print(" %s: skipped %d rows (%d distinct sites) that look synthetic (mutated residue or engineered sequon) -- see %s" % ( tax_name, aa_skipped_rows, len(aa_skipped_sites), aa_file)) else: os.remove(aa_file) # --- drop report ------------------------------------------------------- reasons = Counter(dropped.values()) print(" %s: %d unique step_2 sites -> %d kept, %d dropped %s" % ( tax_name, len(all_sites), len(kept_sites), len(dropped), dict(reasons))) if dropped: with open(dropped_file, "w") as DW: DW.write("\"%s\"\n" % ("\",\"".join( ["src_xref_id", "amino_acid_chain", "amino_acid_chain_auth", "glycosylation_site_pdb", "reason"]))) for (pid, lc, ac, lp), why in sorted(dropped.items()): DW.write("\"%s\"\n" % ("\",\"".join([pid, lc, ac, lp, why]))) if not mismatch_count: # no point keeping an empty mismatch-log file around os.remove(mismatch_file) if fallback_count: print(" %s: %d rows needed the position-only SIFTS fallback (no auth_chain/label_chain exact match)" % (tax_name, fallback_count)) if mismatch_count: print(" %s: %d/%d auth-position cross-checks disagreed between RCSB and SIFTS -- see %s" % ( tax_name, mismatch_count, checked_count, mismatch_file)) print(" %s: wrote %s" % (tax_name, out_file)) return True def main(): # Dated release folder, named once for the whole run (MM_DD_YYYY) release_name = datetime.now().strftime("%m_%d_%Y") release_dir = os.path.join(RELEASES_DIR, release_name) os.makedirs(release_dir, exist_ok=True) os.makedirs(STEP_4_DIR, exist_ok=True) print("Writing final datasets to %s" % release_dir) seen = {} written = 0 file_list = glob.glob("%s/*.sifts_map.json" % STEP_3_DIR) for in_file in file_list: tax_name = os.path.basename(in_file).replace(".sifts_map.json", "") if tax_name not in seen: if extract_data(tax_name, release_dir): written += 1 seen[tax_name] = True subprocess.getoutput("chmod -R 775 %s/ %s/" % (release_dir, STEP_4_DIR)) # Only move `current` if this run actually produced datasets, so a failed # run can never leave it pointing at an empty folder. if written: update_current_symlink(release_name) else: print("WARNING: no datasets were written -- 'current' symlink NOT updated.") return if __name__ == '__main__': main() print("Script finished successfully.", flush=True)