protein.py 11 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322
  1. # Copyright 2021 DeepMind Technologies Limited
  2. #
  3. # Licensed under the Apache License, Version 2.0 (the "License");
  4. # you may not use this file except in compliance with the License.
  5. # You may obtain a copy of the License at
  6. #
  7. # http://www.apache.org/licenses/LICENSE-2.0
  8. #
  9. # Unless required by applicable law or agreed to in writing, software
  10. # distributed under the License is distributed on an "AS IS" BASIS,
  11. # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  12. # See the License for the specific language governing permissions and
  13. # limitations under the License.
  14. """Protein data type."""
  15. import dataclasses
  16. import io
  17. from typing import Any, Mapping, Optional
  18. import numpy as np
  19. from Bio.PDB import PDBParser
  20. from modelscope.models.science.unifold.data import residue_constants
  21. FeatureDict = Mapping[str, np.ndarray]
  22. ModelOutput = Mapping[str, Any] # Is a nested dict.
  23. # Complete sequence of chain IDs supported by the PDB format.
  24. PDB_CHAIN_IDS = 'ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789'
  25. PDB_MAX_CHAINS = len(PDB_CHAIN_IDS) # := 62.
  26. @dataclasses.dataclass(frozen=True)
  27. class Protein:
  28. """Protein structure representation."""
  29. # Cartesian coordinates of atoms in angstroms. The atom types correspond to
  30. # residue_constants.atom_types, i.e. the first three are N, CA, CB.
  31. atom_positions: np.ndarray # [num_res, num_atom_type, 3]
  32. # Amino-acid type for each residue represented as an integer between 0 and
  33. # 20, where 20 is 'X'.
  34. aatype: np.ndarray # [num_res]
  35. # Binary float mask to indicate presence of a particular atom. 1.0 if an atom
  36. # is present and 0.0 if not. This should be used for loss masking.
  37. atom_mask: np.ndarray # [num_res, num_atom_type]
  38. # Residue index as used in PDB. It is not necessarily continuous or 0-indexed.
  39. residue_index: np.ndarray # [num_res]
  40. # 0-indexed number corresponding to the chain in the protein that this residue
  41. # belongs to.
  42. chain_index: np.ndarray # [num_res]
  43. # B-factors, or temperature factors, of each residue (in sq. angstroms units),
  44. # representing the displacement of the residue from its ground truth mean
  45. # value.
  46. b_factors: np.ndarray # [num_res, num_atom_type]
  47. def __post_init__(self):
  48. if len(np.unique(self.chain_index)) > PDB_MAX_CHAINS:
  49. raise ValueError(
  50. f'Cannot build an instance with more than {PDB_MAX_CHAINS} chains '
  51. 'because these cannot be written to PDB format.')
  52. def from_pdb_string(pdb_str: str, chain_id: Optional[str] = None) -> Protein:
  53. """Takes a PDB string and constructs a Protein object.
  54. WARNING: All non-standard residue types will be converted into UNK. All
  55. non-standard atoms will be ignored.
  56. Args:
  57. pdb_str: The contents of the pdb file
  58. chain_id: If chain_id is specified (e.g. A), then only that chain
  59. is parsed. Otherwise all chains are parsed.
  60. Returns:
  61. A new `Protein` parsed from the pdb contents.
  62. """
  63. pdb_fh = io.StringIO(pdb_str)
  64. parser = PDBParser(QUIET=True)
  65. structure = parser.get_structure('none', pdb_fh)
  66. models = list(structure.get_models())
  67. if len(models) != 1:
  68. raise ValueError(
  69. f'Only single model PDBs are supported. Found {len(models)} models.'
  70. )
  71. model = models[0]
  72. atom_positions = []
  73. aatype = []
  74. atom_mask = []
  75. residue_index = []
  76. chain_ids = []
  77. b_factors = []
  78. for chain in model:
  79. if chain_id is not None and chain.id != chain_id:
  80. continue
  81. for res in chain:
  82. if res.id[2] != ' ':
  83. raise ValueError(
  84. f'PDB contains an insertion code at chain {chain.id} and residue '
  85. f'index {res.id[1]}. These are not supported.')
  86. res_shortname = residue_constants.restype_3to1.get(
  87. res.resname, 'X')
  88. restype_idx = residue_constants.restype_order.get(
  89. res_shortname, residue_constants.restype_num)
  90. pos = np.zeros((residue_constants.atom_type_num, 3))
  91. mask = np.zeros((residue_constants.atom_type_num, ))
  92. res_b_factors = np.zeros((residue_constants.atom_type_num, ))
  93. for atom in res:
  94. if atom.name not in residue_constants.atom_types:
  95. continue
  96. pos[residue_constants.atom_order[atom.name]] = atom.coord
  97. mask[residue_constants.atom_order[atom.name]] = 1.0
  98. res_b_factors[residue_constants.atom_order[
  99. atom.name]] = atom.bfactor
  100. if np.sum(mask) < 0.5:
  101. # If no known atom positions are reported for the residue then skip it.
  102. continue
  103. aatype.append(restype_idx)
  104. atom_positions.append(pos)
  105. atom_mask.append(mask)
  106. residue_index.append(res.id[1])
  107. chain_ids.append(chain.id)
  108. b_factors.append(res_b_factors)
  109. # Chain IDs are usually characters so map these to ints.
  110. unique_chain_ids = np.unique(chain_ids)
  111. chain_id_mapping = {cid: n for n, cid in enumerate(unique_chain_ids)}
  112. chain_index = np.array([chain_id_mapping[cid] for cid in chain_ids])
  113. return Protein(
  114. atom_positions=np.array(atom_positions),
  115. atom_mask=np.array(atom_mask),
  116. aatype=np.array(aatype),
  117. residue_index=np.array(residue_index),
  118. chain_index=chain_index,
  119. b_factors=np.array(b_factors),
  120. )
  121. def _chain_end(atom_index, end_resname, chain_name, residue_index) -> str:
  122. chain_end = 'TER'
  123. return (f'{chain_end:<6}{atom_index:>5} {end_resname:>3} '
  124. f'{chain_name:>1}{residue_index:>4}')
  125. def to_pdb(prot: Protein) -> str:
  126. """Converts a `Protein` instance to a PDB string.
  127. Args:
  128. prot: The protein to convert to PDB.
  129. Returns:
  130. PDB string.
  131. """
  132. restypes = residue_constants.restypes + ['X']
  133. # res_1to3 = lambda r: residue_constants.restype_1to3.get(restypes[r], 'UNK')
  134. def res_1to3(r):
  135. return residue_constants.restype_1to3.get(restypes[r], 'UNK')
  136. atom_types = residue_constants.atom_types
  137. pdb_lines = []
  138. atom_mask = prot.atom_mask
  139. aatype = prot.aatype
  140. atom_positions = prot.atom_positions
  141. residue_index = prot.residue_index.astype(np.int32)
  142. chain_index = prot.chain_index.astype(np.int32)
  143. b_factors = prot.b_factors
  144. if np.any(aatype > residue_constants.restype_num):
  145. raise ValueError('Invalid aatypes.')
  146. # Construct a mapping from chain integer indices to chain ID strings.
  147. chain_ids = {}
  148. for i in np.unique(chain_index): # np.unique gives sorted output.
  149. if i >= PDB_MAX_CHAINS:
  150. raise ValueError(
  151. f'The PDB format supports at most {PDB_MAX_CHAINS} chains.')
  152. chain_ids[i] = PDB_CHAIN_IDS[i]
  153. pdb_lines.append('MODEL 1')
  154. atom_index = 1
  155. last_chain_index = chain_index[0]
  156. # Add all atom sites.
  157. for i in range(aatype.shape[0]):
  158. # Close the previous chain if in a multichain PDB.
  159. if last_chain_index != chain_index[i]:
  160. pdb_lines.append(
  161. _chain_end(
  162. atom_index,
  163. res_1to3(aatype[i - 1]),
  164. chain_ids[chain_index[i - 1]],
  165. residue_index[i - 1],
  166. ))
  167. last_chain_index = chain_index[i]
  168. atom_index += 1 # Atom index increases at the TER symbol.
  169. res_name_3 = res_1to3(aatype[i])
  170. for atom_name, pos, mask, b_factor in zip(atom_types,
  171. atom_positions[i],
  172. atom_mask[i], b_factors[i]):
  173. if mask < 0.5:
  174. continue
  175. record_type = 'ATOM'
  176. name = atom_name if len(atom_name) == 4 else f' {atom_name}'
  177. alt_loc = ''
  178. insertion_code = ''
  179. occupancy = 1.00
  180. element = atom_name[
  181. 0] # Protein supports only C, N, O, S, this works.
  182. charge = ''
  183. # PDB is a columnar format, every space matters here!
  184. atom_line = (
  185. f'{record_type:<6}{atom_index:>5} {name:<4}{alt_loc:>1}'
  186. f'{res_name_3:>3} {chain_ids[chain_index[i]]:>1}'
  187. f'{residue_index[i]:>4}{insertion_code:>1} '
  188. f'{pos[0]:>8.3f}{pos[1]:>8.3f}{pos[2]:>8.3f}'
  189. f'{occupancy:>6.2f}{b_factor:>6.2f} '
  190. f'{element:>2}{charge:>2}')
  191. pdb_lines.append(atom_line)
  192. atom_index += 1
  193. # Close the final chain.
  194. pdb_lines.append(
  195. _chain_end(
  196. atom_index,
  197. res_1to3(aatype[-1]),
  198. chain_ids[chain_index[-1]],
  199. residue_index[-1],
  200. ))
  201. pdb_lines.append('ENDMDL')
  202. pdb_lines.append('END')
  203. # Pad all lines to 80 characters.
  204. pdb_lines = [line.ljust(80) for line in pdb_lines]
  205. return '\n'.join(pdb_lines) + '\n' # Add terminating newline.
  206. def ideal_atom_mask(prot: Protein) -> np.ndarray:
  207. """Computes an ideal atom mask.
  208. `Protein.atom_mask` typically is defined according to the atoms that are
  209. reported in the PDB. This function computes a mask according to heavy atoms
  210. that should be present in the given sequence of amino acids.
  211. Args:
  212. prot: `Protein` whose fields are `numpy.ndarray` objects.
  213. Returns:
  214. An ideal atom mask.
  215. """
  216. return residue_constants.STANDARD_ATOM_MASK[prot.aatype]
  217. def from_prediction(features: FeatureDict,
  218. result: ModelOutput,
  219. b_factors: Optional[np.ndarray] = None) -> Protein:
  220. """Assembles a protein from a prediction.
  221. Args:
  222. features: Dictionary holding model inputs.
  223. fold_output: Dictionary holding model outputs.
  224. b_factors: (Optional) B-factors to use for the protein.
  225. Returns:
  226. A protein instance.
  227. """
  228. if 'asym_id' in features:
  229. chain_index = features['asym_id'] - 1
  230. else:
  231. chain_index = np.zeros_like((features['aatype']))
  232. if b_factors is None:
  233. b_factors = np.zeros_like(result['final_atom_mask'])
  234. return Protein(
  235. aatype=features['aatype'],
  236. atom_positions=result['final_atom_positions'],
  237. atom_mask=result['final_atom_mask'],
  238. residue_index=features['residue_index'] + 1,
  239. chain_index=chain_index,
  240. b_factors=b_factors,
  241. )
  242. def from_feature(features: FeatureDict,
  243. b_factors: Optional[np.ndarray] = None) -> Protein:
  244. """Assembles a standard pdb from input atom positions & mask.
  245. Args:
  246. features: Dictionary holding model inputs.
  247. b_factors: (Optional) B-factors to use for the protein.
  248. Returns:
  249. A protein instance.
  250. """
  251. if 'asym_id' in features:
  252. chain_index = features['asym_id'] - 1
  253. else:
  254. chain_index = np.zeros_like((features['aatype']))
  255. if b_factors is None:
  256. b_factors = np.zeros_like(features['all_atom_mask'])
  257. return Protein(
  258. aatype=features['aatype'],
  259. atom_positions=features['all_atom_positions'],
  260. atom_mask=features['all_atom_mask'],
  261. residue_index=features['residue_index'] + 1,
  262. chain_index=chain_index,
  263. b_factors=b_factors,
  264. )