SIFTS勘正-平芜编程栈

很久之前的一篇博客，参考：

如何量化1个PDB结构所描述的事件？

里面简单提到了一个SIFTS的计算方法，不走寻常SIFTS数据库，仅仅是从pdb或mmcif本身结构的index，我们biopython解析之后拿到 Biopython 的 res.id 是 (hetflag, auth_seq_id, insertion_code)，但是用7W1M简单测试了一下，我们当时测试的是chain H，也就是CTCF蛋白，res.id[1]符合实际结构建模部分序列index，

但是PDB/mmCIF 规范里没有任何东西保证 auth_seq_id 与 UniProt 坐标一致，

也就是说有可能是7w1m的投稿人正好按按照uniprot给chain H编号，也就是auth_seq_id 和 uniprot_index之间没有必然联系。

比如说，我这里随便举一个例子：PDB 2lyz

如果按照之前biopython解析结构数据，拿到auth index的话，就是1-129

我们到pdb中查看也没有问题

但是实际序列比对之后，也就是用uniprot的序列mapping之后

这一段其实是 P00698 19–147。

我们目前的做法依然是需要依赖SIFTS数据库的，下面是草稿

# ══════════════════════════════════════════════════════════════ # SIFTS PDB ↔ UniProt 逐残基映射 —— 独立模块 # 依赖: 标准库 + httpx（缺失时回退 requests） # ══════════════════════════════════════════════════════════════ """SIFTS PDB <-> UniProt per-residue mapping (standalone, no fragment logic).""" import gzip import random import re import time import xml.etree.ElementTree as ET from pathlib import Path from typing import Any, Dict, List, Optional try: import httpx except ImportError: # pragma: no cover httpx = None try: import requests except ImportError: # pragma: no cover requests = None def http_get_raw(url: str, retries: int = 5, timeout: float = 60.0) -> Any: """GET url, return raw bytes, or "NOT_FOUND" / "BAD_REQUEST".""" headers = {"User-Agent": "IDRInter-SIFTS-test/1.0"} for attempt in range(retries): try: if httpx is not None: res = httpx.get(url, timeout=timeout, follow_redirects=True, headers=headers) elif requests is not None: res = requests.get(url, timeout=timeout, allow_redirects=True, headers=headers) else: raise RuntimeError("install httpx or requests") code = int(res.status_code) if code == 200: return res.content if code in (400, 401, 403): return "BAD_REQUEST" if code in (404, 410): return "NOT_FOUND" if code in (429, 500, 502, 503, 504) and attempt < retries - 1: time.sleep(min(5.0, 0.5 * (2 ** attempt))) continue except Exception: if attempt < retries - 1: time.sleep(min(5.0, 0.5 * (2 ** attempt))) continue return "NOT_FOUND" return "NOT_FOUND" def convert_aa_3to1(code_3: str) -> str: mapping = { 'ASP': 'D', 'GLU': 'E', 'PHE': 'F', 'LYS': 'K', 'ASN': 'N', 'GLN': 'Q', 'ARG': 'R', 'TRP': 'W', 'TYR': 'Y', 'PYL': 'O', 'SEC': 'X', 'ALA': 'A', 'CYS': 'C', 'GLY': 'G', 'HIS': 'H', 'ILE': 'I', 'LEU': 'L', 'MET': 'M', 'PRO': 'P', 'SER': 'S', 'THR': 'T', 'VAL': 'V', } return mapping.get(str(code_3).upper(), str(code_3)[0] if code_3 else '') def _safe_int(val: Any) -> Optional[int]: """'91' -> 91, '91A' -> 91, 'null' -> None.""" s = str(val).strip() if s.lower() in ("null", "none", "nan", ""): return None m = re.match(r"^\d+", s) if not m: return None return int(m.group(0)) def _parse_sifts_xml_to_rows(xml_text: str) -> List[List[str]]: """SIFTS XML -> 7-column rows: [PDB_ID, PDB_chain, PDB_res, PDB_pos, Uniprot_ID, Uniprot_res, Uniprot_pos] """ root = ET.fromstring(xml_text) pdb_id = root.get("dbAccessionId", "") def _local(tag: str) -> str: return tag.rsplit("}", 1)[-1] rows: List[List[str]] = [] for residue in root.iter(): if _local(residue.tag) != "residue": continue pdb_ref = uni_ref = None for child in residue: if _local(child.tag) != "crossRefDb": continue src = child.get("dbSource") if src == "PDB" and pdb_ref is None: pdb_ref = child elif src == "UniProt" and uni_ref is None: uni_ref = child if pdb_ref is None or uni_ref is None: continue rows.append([ pdb_id, pdb_ref.get("dbChainId", ""), convert_aa_3to1(pdb_ref.get("dbResName", "")), pdb_ref.get("dbResNum", ""), uni_ref.get("dbAccessionId", ""), uni_ref.get("dbResName", ""), uni_ref.get("dbResNum", ""), ]) return rows def download_sifts_mapping( pdb_id: str, output_dir: Path, max_retries: int = 5, wait_time: float = 5.0, ) -> Optional[Path]: """Download SIFTS XML, convert to a TSV cache, return its path.""" pdb_lower = pdb_id.lower() output_dir = Path(output_dir) output_dir.mkdir(parents=True, exist_ok=True) out_path = output_dir / f"{pdb_lower}.tsv" if out_path.exists() and out_path.stat().st_size > 0: return out_path url = f"https://www.ebi.ac.uk/pdbe/files/sifts/{pdb_lower}.xml.gz" for _ in range(max_retries): raw = http_get_raw(url) if raw in ("NOT_FOUND", "BAD_REQUEST") or not raw: time.sleep(random.uniform(1, wait_time)) continue try: xml_text = gzip.decompress(raw).decode("utf-8") except Exception: xml_text = raw.decode("utf-8", errors="ignore") rows = _parse_sifts_xml_to_rows(xml_text) if not rows: time.sleep(random.uniform(1, wait_time)) continue with open(out_path, "w") as f: for row in rows: f.write("\t".join(row) + "\n") return out_path return None def parse_sifts_mapping( sifts_path: Path, chain_id: str, uniprot_ids: List[str], ) -> Dict[str, List]: """Per-residue mapping for one auth chain (and the given UniProt accessions). uni_pos: UniProt positions, including residues not observed in the structure. pdb_pos: author residue numbers; None where the residue is missing in the structure. """ sifts_path = Path(sifts_path) if not sifts_path.exists(): return {"uni_pos": [], "pdb_pos": []} uni_set = {str(u).split("-")[0].strip() for u in uniprot_ids if u} chain = str(chain_id).strip() uni_pos: List[int] = [] pdb_pos: List[Optional[int]] = [] try: with open(sifts_path, "r") as f: lines = f.read().splitlines() except Exception: return {"uni_pos": [], "pdb_pos": []} for line in lines: if not line.strip(): continue cols = line.split("\t") if len(cols) < 7: continue if cols[1].strip() != chain: continue if cols[4].split("-")[0].strip() not in uni_set: continue u = _safe_int(cols[6]) if u is None: continue uni_pos.append(u) pdb_pos.append(_safe_int(cols[3])) return {"uni_pos": uni_pos, "pdb_pos": pdb_pos} def list_chain_uniprots(sifts_path: Path) -> Dict[str, List[str]]: """{auth_chain: sorted UniProt accessions} from the SIFTS TSV.""" result: Dict[str, set] = {} with open(sifts_path, "r") as f: for line in f: cols = line.rstrip("\n").split("\t") if len(cols) < 7: continue result.setdefault(cols[1].strip(), set()).add(cols[4].split("-")[0].strip()) return {k: sorted(v) for k, v in result.items()} def load_sifts_mapping( sifts_path: Path, chain_id: str, uniprot_ids: Optional[List[str]] = None, ) -> Dict[str, List]: """parse_sifts_mapping with automatic UniProt lookup for the chain.""" if uniprot_ids is None: uniprot_ids = list_chain_uniprots(sifts_path).get(str(chain_id).strip(), []) return parse_sifts_mapping(sifts_path, chain_id, uniprot_ids)

我们按照上面的方法测试一下，现在真实的数据应该如何

CACHE_DIR = Path("/data2/test_sifts_cache") PDBS = ["7w1m", "4ins", "2lyz"] def to_ranges(positions): """连续位置压成区间: [223,224,225,270] -> [(223,225),(270,270)]""" out = [] for p in sorted(positions): if out and p == out[-1][1] + 1: out[-1][1] = p else: out.append([p, p]) return [tuple(r) for r in out] for pdb in PDBS: tsv = download_sifts_mapping(pdb, CACHE_DIR) if tsv is None: print(f"{pdb}: SIFTS 下载失败") continue print("=" * 80) print(f"{pdb} -> {tsv}") for chain, unis in sorted(list_chain_uniprots(tsv).items()): m = load_sifts_mapping(tsv, chain, unis) if not m["uni_pos"]: continue observed = [u for u, p in zip(m["uni_pos"], m["pdb_pos"]) if p is not None] missing = sum(p is None for p in m["pdb_pos"]) print(f" chain {chain} uniprot={','.join(unis)}") print(f" mapped={len(m['uni_pos'])} observed={len(observed)} missing={missing}") print(f" covered_uni={to_ranges(observed)}") print()

输出如下

================================================================================ 7w1m -> /data2/test_sifts_cache/7w1m.tsv chain A uniprot=Q14683 mapped=1233 observed=399 missing=834 covered_uni=[(2, 201), (1027, 1225)] chain B uniprot=Q9UQE7 mapped=1217 observed=528 missing=689 covered_uni=[(1, 256), (915, 1060), (1092, 1217)] chain C uniprot=O60216 mapped=631 observed=248 missing=383 covered_uni=[(10, 92), (154, 171), (321, 394), (558, 630)] chain D uniprot=Q8WVM7 mapped=1258 observed=967 missing=291 covered_uni=[(86, 1052)] chain E uniprot=Q6KC79 mapped=1467 observed=1226 missing=241 covered_uni=[(1193, 1216), (1231, 1280), (1293, 1357), (1380, 1475), (1484, 1505), (1524, 1629), (1646, 1690), (1708, 1729), (1746, 1987), (1998, 2372), (2389, 2471), (2533, 2628)] chain H uniprot=P49711 mapped=727 observed=322 missing=405 covered_uni=[(223, 231), (265, 577)] ================================================================================ 4ins -> /data2/test_sifts_cache/4ins.tsv chain A uniprot=P01315 mapped=21 observed=21 missing=0 covered_uni=[(88, 108)] chain B uniprot=P01315 mapped=30 observed=30 missing=0 covered_uni=[(25, 54)] chain C uniprot=P01315 mapped=21 observed=21 missing=0 covered_uni=[(88, 108)] chain D uniprot=P01315 mapped=30 observed=30 missing=0 covered_uni=[(25, 54)] ================================================================================ 2lyz -> /data2/test_sifts_cache/2lyz.tsv chain A uniprot=P00698 mapped=129 observed=129 missing=0 covered_uni=[(19, 147)]

我们一个一个来验证：

首先是7w1m

* Q14683：chain A，covered_uni=[(2, 201), (1027, 1225)]

正好未建模区域为202-1026

* Q9UQE7：chain B，[(1, 256), (915, 1060), (1092, 1217)]

* chain C，covered_uni=[(10, 92), (154, 171), (321, 394), (558, 630)]

* chain D，covered_uni=[(86, 1052)]

* chain E，covered_uni=[(1193, 1216), (1231, 1280), (1293, 1357), (1380, 1475), (1484, 1505), (1524, 1629), (1646, 1690), (1708, 1729), (1746, 1987), (1998, 2372), (2389, 2471), (2533, 2628)]

* chain H，covered_uni=[(223, 231), (265, 577)]

其余的例子也就不举了，

import pandas as pd pdb, chain = "7w1m", "H" tsv = download_sifts_mapping(pdb, CACHE_DIR) print("UniProt:", list_chain_uniprots(tsv)[chain]) m = load_sifts_mapping(tsv, chain) pairs = [(u, p) for u, p in zip(m["uni_pos"], m["pdb_pos"]) if p is not None] print("mapped residues :", len(m["uni_pos"])) print("observed residues:", len(pairs)) print("missing residues :", sum(p is None for p in m["pdb_pos"])) print("covered_uni :", to_ranges([u for u, _ in pairs])) print("author_pos_range :", (pairs[0][1], pairs[-1][1])) pd.DataFrame({"uniprot_pos": m["uni_pos"], "pdb_pos": m["pdb_pos"]})

我们这里再回到biopython解析的逻辑，我们之前用的就是res.id[1]

def to_ranges(positions): """连续位置压成区间: [223,224,225,270] -> [(223,225),(270,270)]""" out = [] for p in sorted(positions): if out and p == out[-1][1] + 1: out[-1][1] = p else: out.append([p, p]) return [tuple(r) for r in out] model = structure[0] # Get the first model for chain in model: print(f"Chain ID: {chain.id}") list_res = [] for residue in chain: list_res.append(residue.id[1]) list_range = to_ranges(list_res) print(f" Residue range: {list_range}") # resname and residue number

我们可以看到，7w1m这个例子大部分是准的，但是多了一些东西

这些都是非氨基酸，都是结构文件里带进来的其他小分子之类，

我们公平的使用下述biopython代码，如果是res.id[1]逻辑的话

from Bio.PDB import PDBParser from Bio.PDB.Polypeptide import is_aa parser = PDBParser(QUIET=True) # QUIET=True to suppress warnings about missing atoms structure = parser.get_structure("example", "/data2/7W1M.pdb") def to_ranges(positions): """连续位置压成区间: [223,224,225,270] -> [(223,225),(270,270)]""" out = [] for p in sorted(positions): if out and p == out[-1][1] + 1: out[-1][1] = p else: out.append([p, p]) return [tuple(r) for r in out] model = structure[0] # Get the first model for chain in model: aa = [r for r in chain if is_aa(r, standard=True)] print(f"Chain ID: {chain.id}") print(f" Residue range: {to_ranges(r.id[1] for r in aa)}")

这个就是当时偶然测试没有问题的例子：

我们用其他例子、程序检验一下

* 2lyz

这里问题和前面说的一样，我们如果拿这段序列去做uniprot注释，以为是1-129，结果就会把1-129的uniprot序列特征注释到实际上是19-147的区域上，就会有偏差，就有问题了

* 4ins

我们这个例子也要实际pdb比对一下：1D PFV: 4INS

再到uniprot中验证，确实不是开头的1-21

其他chain我们保险起见也一一查看一下

所以其实很明确了，

auth_* 就是指作者（投稿人）自己搞的一套编号、chain id等等，

这个一般在解析之后和结构数据库里的chain id、residue id都不一定对应上，

biopython 这里就是一个坑，