import sys, os import json import glob import csv import subprocess # --------------------------------------------------------------------------- # Configuration # --------------------------------------------------------------------------- BASE_DIR = "/data/projects/glygen/downloads/rcsb_pdb/temp" IN_DIR = BASE_DIR + "/step_1" OUT_DIR = BASE_DIR + "/step_2" # role -> glycosylation_type GLYCO_TYPE_MAP = { "N-Glycosylation": "N-linked", "O-Glycosylation": "O-linked", "C-Mannosylation": "C-linked", "S-Glycosylation": "S-linked", } HEADER = [ "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", "chemically_synthetic", ] # --------------------------------------------------------------------------- # Helpers # --------------------------------------------------------------------------- def map_glyco_type(role): """O-Glycosylation -> O-linked, C-Mannosylation -> C-linked, etc.""" if role is None: return "" return GLYCO_TYPE_MAP.get(role, role.replace("Glycosylation", "linked")) def format_amino_acid(comp_id): """ASN -> Asn (capitalise first letter only).""" if comp_id is None: return "" return comp_id.capitalize() def pick_saccharide(reference_identifiers): """GlyTouCan first, else GlyCosmos, else GlyGen accession.""" buckets = {"GlyTouCan": [], "GlyCosmos": [], "GlyGen": []} if reference_identifiers: for o in reference_identifiers: name = o.get("resource_name") acc = o.get("resource_accession") if name in buckets and acc: buckets[name].append(acc) for src in ["GlyTouCan", "GlyCosmos", "GlyGen"]: if buckets[src]: return buckets[src][0] return "" def match_key(role, chain, seq_id, comp_id): """Raw-value key shared between step-2 sites and step-4 glycan instances.""" return (str(role), str(chain), str(seq_id), str(comp_id)) def extract_host_organism_info(polymer_entity): """ Given one entry from obj["polymer_entities"], return the four host-organism output columns for that entity. rcsb_entity_host_organism and entity_src_gen are both lists -- an entity can have more than one source record (e.g. a fusion construct expressed from more than one gene/host). Rules: - data missing (key absent / empty list) -> "" (blank) - a single source record with a null scientific name -> host_organism = "null" - more than one source record (either list) -> host_organism = "fusion constructs", the other three columns left blank """ result = {"host_organism": "", "host_org_cell_line": "", "host_org_cell": "", "host_org_organ": ""} host_org_list = polymer_entity.get("rcsb_entity_host_organism") or [] src_gen_list = polymer_entity.get("entity_src_gen") or [] if len(host_org_list) > 1 or len(src_gen_list) > 1: result["host_organism"] = "fusion constructs" return result if host_org_list: ncbi_name = host_org_list[0].get("ncbi_scientific_name") result["host_organism"] = "null" if ncbi_name is None else ncbi_name # else: no host organism data at all -- leave "" (missing) if src_gen_list: rec = src_gen_list[0] result["host_org_cell_line"] = rec.get("pdbx_host_org_cell_line") or "" result["host_org_cell"] = rec.get("pdbx_host_org_cell") or "" result["host_org_organ"] = rec.get("pdbx_host_org_organ") or "" return result # --------------------------------------------------------------------------- # Core extraction # --------------------------------------------------------------------------- def build_glycan_map(obj): """Map raw site key -> (pdb_glycan_entity_id, saccharide) from branched_entities.""" glycan_map = {} for be in (obj.get("branched_entities") or []): ent_id = be.get("rcsb_id", "") cont = be.get("rcsb_branched_entity_container_identifiers") or {} saccharide = pick_saccharide(cont.get("reference_identifiers")) for inst in (be.get("branched_entity_instances") or []): for conn in (inst.get("rcsb_branched_struct_conn") or []): role = conn.get("role") if role is None: # skip sugar-sugar internal bonds continue cp = conn.get("connect_partner") or {} key = match_key(role, cp.get("label_asym_id"), cp.get("label_seq_id"), cp.get("label_comp_id")) glycan_map[key] = (ent_id, saccharide) return glycan_map SYN_FLAGGED = {"sites": 0, "entries": set()} def is_chemically_synthetic(pe): """ True if RCSB records this polymer entity as chemically synthesised (flagged, not dropped). Reads rcsb_polymer_entity.src_method (or entity.src_method) from the step 1 GraphQL query; a value starting with "syn" counts. If the field is absent nothing is treated as synthetic. """ for block in ("rcsb_polymer_entity", "entity"): val = ((pe.get(block) or {}).get("src_method") or "").strip().lower() if val: return val.startswith("syn") return False def build_sites(obj): """ One site per polymer struct_conn -- every glycosylation connection gets a row here, whether or not a matching branched glycan entity is found later in build_glycan_map (that lookup is a separate, optional join in extract_data, defaulting to blank saccharide/pdb_glycan_entity_id rather than dropping the site). """ sites = [] for pe in (obj.get("polymer_entities") or []): synthetic = is_chemically_synthetic(pe) cont = pe.get("rcsb_polymer_entity_container_identifiers") or {} asym = cont.get("asym_ids") or [] auth = cont.get("auth_asym_ids") or [] chain_to_auth = {a: (auth[i] if i < len(auth) else a) for i, a in enumerate(asym)} host_info = extract_host_organism_info(pe) for inst in (pe.get("polymer_entity_instances") or []): for conn in (inst.get("rcsb_polymer_struct_conn") or []): role = conn.get("role") # keep only glycosylation; skip disulfide bridges, metal coordination, etc. if role not in GLYCO_TYPE_MAP: continue cp = conn.get("connect_partner") or {} ct = conn.get("connect_target") or {} chain = ct.get("label_asym_id") seq_pdb = ct.get("label_seq_id") aa_comp = ct.get("label_comp_id") sites.append({ "key": match_key(role, chain, seq_pdb, aa_comp), "glycosylation_type": map_glyco_type(role), "amino_acid_chain": chain or "", "glycosylation_site_pdb": "" if seq_pdb is None else str(seq_pdb), "amino_acid": format_amino_acid(aa_comp), "carb_name": cp.get("label_comp_id") or "", "amino_acid_chain_auth": chain_to_auth.get(chain, chain or ""), "glycosylation_site_auth": ct.get("auth_seq_id") or "", "host_organism": host_info["host_organism"], "host_org_cell_line": host_info["host_org_cell_line"], "host_org_cell": host_info["host_org_cell"], "host_org_organ": host_info["host_org_organ"], "chemically_synthetic": "yes" if synthetic else "", }) if synthetic: SYN_FLAGGED["sites"] += 1 SYN_FLAGGED["entries"].add((obj.get("entry") or {}).get("id") or "?") return sites def build_xref_list(obj, src_xref_id): """Original xref logic: PubMed -> DOI per citation, fallback to PDB id.""" xref_id_list = [] if "citation" in obj and obj["citation"] is not None: for o in obj["citation"]: pubmed_flag = False if o.get("pdbx_database_id_PubMed") is not None: xref_id_list.append({"xref_id": o["pdbx_database_id_PubMed"], "xref_key": "protein_xref_pubmed"}) pubmed_flag = True if not pubmed_flag and o.get("pdbx_database_id_DOI") is not None: xref_id_list.append({"xref_id": o["pdbx_database_id_DOI"], "xref_key": "protein_xref_doi"}) if not xref_id_list: xref_id_list.append({"xref_id": src_xref_id, "xref_key": "protein_xref_pdb"}) return xref_id_list def extract_data(tax_name): out_file = os.path.join(OUT_DIR, "%s.csv" % tax_name) with open(out_file, "w", newline="") as fh: writer = csv.writer(fh, quoting=csv.QUOTE_ALL) writer.writerow(HEADER) for in_file in glob.glob(os.path.join(IN_DIR, "%s.*.json" % tax_name)): doc = json.load(open(in_file)) for obj in doc["data"]["entries"]: src_xref_key = "protein_xref_pdb4glycans" src_xref_id = obj["entry"]["id"] # pdb_extended_id from database_2 pdb_extended_id_list = [] for o in (obj.get("database_2") or []): acc = o.get("pdbx_database_accession") if acc: pdb_extended_id_list.append(acc) pdb_extended_id_list = list(set(pdb_extended_id_list)) or [""] xref_id_list = build_xref_list(obj, src_xref_id) sites = build_sites(obj) glycan_map = build_glycan_map(obj) for pdb_extended_id in pdb_extended_id_list: for xo in xref_id_list: xref_key, xref_id = xo["xref_key"], str(xo["xref_id"]) for s in sites: ent_id, saccharide = glycan_map.get(s["key"], ("", "")) row = { "glycosylation_type": s["glycosylation_type"], "amino_acid_chain": s["amino_acid_chain"], "glycosylation_site_pdb": s["glycosylation_site_pdb"], "amino_acid": s["amino_acid"], "saccharide": saccharide, "pdb_glycan_entity_id": ent_id, "pdb_extended_id": pdb_extended_id, "xref_id": xref_id, "xref_key": xref_key, "src_xref_key": src_xref_key, "src_xref_id": src_xref_id, "amino_acid_chain_auth": s["amino_acid_chain_auth"], "glycosylation_site_auth": s["glycosylation_site_auth"], "carb_name": s["carb_name"], "host_organism": s["host_organism"], "host_org_cell_line": s["host_org_cell_line"], "host_org_cell": s["host_org_cell"], "host_org_organ": s["host_org_organ"], "chemically_synthetic": s["chemically_synthetic"], } writer.writerow([row[c] for c in HEADER]) subprocess.getoutput("chmod -R 775 %s" % OUT_DIR) return def main(): os.makedirs(OUT_DIR, exist_ok=True) seen = {} for in_file in glob.glob(os.path.join(IN_DIR, "*.json")): tax_name = in_file.split("/")[-1].split(".")[0] if tax_name not in seen: extract_data(tax_name) seen[tax_name] = True return if __name__ == '__main__': main() if SYN_FLAGGED["sites"]: print("Flagged %d glycosylation connections in %d entries whose protein is chemically synthetic " "(entity.src_method = syn): %s" % (SYN_FLAGGED["sites"], len(SYN_FLAGGED["entries"]), ", ".join(sorted(SYN_FLAGGED["entries"])[:20])), flush=True) print("Script finished successfully.", flush=True)