algorithm.ts 24 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605
  1. /**
  2. * Copyright (c) 2020 mol* contributors, licensed under MIT, See LICENSE file for more info.
  3. *
  4. * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org>
  5. * @author Alexander Rose <alexander.rose@weirdbyte.de>
  6. */
  7. import { Structure, StructureElement, StructureProperties, Unit } from '../../mol-model/structure';
  8. import { Task, RuntimeContext } from '../../mol-task';
  9. import { CentroidHelper } from '../../mol-math/geometry/centroid-helper';
  10. import { AccessibleSurfaceAreaParams } from '../../mol-model-props/computed/accessible-surface-area';
  11. import { Vec3 } from '../../mol-math/linear-algebra';
  12. import { getElementMoleculeType } from '../../mol-model/structure/util';
  13. import { MoleculeType } from '../../mol-model/structure/model/types';
  14. import { AccessibleSurfaceArea } from '../../mol-model-props/computed/accessible-surface-area/shrake-rupley';
  15. import { ParamDefinition as PD } from '../../mol-util/param-definition';
  16. import { MembraneOrientation } from './prop';
  17. const LARGE_CA_THRESHOLD = 5000;
  18. const DEFAULT_UPDATE_INTERVAL = 10;
  19. const LARGE_CA_UPDATE_INTERVAL = 1;
  20. interface ANVILContext {
  21. structure: Structure,
  22. numberOfSpherePoints: number,
  23. stepSize: number,
  24. minThickness: number,
  25. maxThickness: number,
  26. asaCutoff: number,
  27. adjust: number,
  28. offsets: ArrayLike<number>,
  29. exposed: ArrayLike<number>,
  30. hydrophobic: ArrayLike<boolean>,
  31. centroid: Vec3,
  32. extent: number,
  33. large: boolean
  34. };
  35. export const ANVILParams = {
  36. numberOfSpherePoints: PD.Numeric(175, { min: 35, max: 700, step: 1 }, { description: 'Number of spheres/directions to test for membrane placement. Original value is 350.' }),
  37. stepSize: PD.Numeric(1, { min: 0.25, max: 4, step: 0.25 }, { description: 'Thickness of membrane slices that will be tested' }),
  38. minThickness: PD.Numeric(20, { min: 10, max: 30, step: 1}, { description: 'Minimum membrane thickness used during refinement' }),
  39. maxThickness: PD.Numeric(40, { min: 30, max: 50, step: 1}, { description: 'Maximum membrane thickness used during refinement' }),
  40. asaCutoff: PD.Numeric(40, { min: 10, max: 100, step: 1 }, { description: 'Relative ASA cutoff above which residues will be considered' }),
  41. adjust: PD.Numeric(14, { min: 0, max: 30, step: 1 }, { description: 'Minimum length of membrane-spanning regions (original values: 14 for alpha-helices and 5 for beta sheets). Set to 0 to not optimize membrane thickness.' })
  42. };
  43. export type ANVILParams = typeof ANVILParams
  44. export type ANVILProps = PD.Values<ANVILParams>
  45. /**
  46. * Implements:
  47. * Membrane positioning for high- and low-resolution protein structures through a binary classification approach
  48. * Guillaume Postic, Yassine Ghouzam, Vincent Guiraud, and Jean-Christophe Gelly
  49. * Protein Engineering, Design & Selection, 2015, 1–5
  50. * doi: 10.1093/protein/gzv063
  51. */
  52. export function computeANVIL(structure: Structure, props: ANVILProps) {
  53. return Task.create('Compute Membrane Orientation', async runtime => {
  54. return await calculate(runtime, structure, props);
  55. });
  56. }
  57. // avoiding namespace lookup improved performance in Chrome (Aug 2020)
  58. const v3add = Vec3.add;
  59. const v3clone = Vec3.clone;
  60. const v3create = Vec3.create;
  61. const v3distance = Vec3.distance;
  62. const v3dot = Vec3.dot;
  63. const v3magnitude = Vec3.magnitude;
  64. const v3normalize = Vec3.normalize;
  65. const v3scale = Vec3.scale;
  66. const v3scaleAndAdd = Vec3.scaleAndAdd;
  67. const v3set = Vec3.set;
  68. const v3squaredDistance = Vec3.squaredDistance;
  69. const v3sub = Vec3.sub;
  70. const v3zero = Vec3.zero;
  71. const centroidHelper = new CentroidHelper();
  72. async function initialize(structure: Structure, props: ANVILProps, accessibleSurfaceArea: AccessibleSurfaceArea): Promise<ANVILContext> {
  73. const l = StructureElement.Location.create(structure);
  74. const { label_atom_id, label_comp_id, x, y, z } = StructureProperties.atom;
  75. const asaCutoff = props.asaCutoff / 100;
  76. centroidHelper.reset();
  77. const offsets = new Array<number>();
  78. const exposed = new Array<number>();
  79. const hydrophobic = new Array<boolean>();
  80. const vec = v3zero();
  81. for (let i = 0, il = structure.units.length; i < il; ++i) {
  82. const unit = structure.units[i];
  83. const { elements } = unit;
  84. l.unit = unit;
  85. for (let j = 0, jl = elements.length; j < jl; ++j) {
  86. const eI = elements[j];
  87. l.element = eI;
  88. // consider only amino acids
  89. if (getElementMoleculeType(unit, eI) !== MoleculeType.Protein) {
  90. continue;
  91. }
  92. // only CA is considered for downstream operations
  93. if (label_atom_id(l) !== 'CA' && label_atom_id(l) !== 'BB') {
  94. continue;
  95. }
  96. // original ANVIL only considers canonical amino acids
  97. if (!MaxAsa[label_comp_id(l)]) {
  98. continue;
  99. }
  100. // while iterating use first pass to compute centroid
  101. v3set(vec, x(l), y(l), z(l));
  102. centroidHelper.includeStep(vec);
  103. // keep track of offsets and exposed state to reuse
  104. offsets.push(structure.serialMapping.getSerialIndex(l.unit, l.element));
  105. if (AccessibleSurfaceArea.getValue(l, accessibleSurfaceArea) / MaxAsa[label_comp_id(l)] > asaCutoff) {
  106. exposed.push(structure.serialMapping.getSerialIndex(l.unit, l.element));
  107. hydrophobic.push(isHydrophobic(label_comp_id(l)));
  108. }
  109. }
  110. }
  111. // calculate centroid and extent
  112. centroidHelper.finishedIncludeStep();
  113. const centroid = v3clone(centroidHelper.center);
  114. for (let k = 0, kl = offsets.length; k < kl; k++) {
  115. setLocation(l, structure, offsets[k]);
  116. v3set(vec, x(l), y(l), z(l));
  117. centroidHelper.radiusStep(vec);
  118. }
  119. const extent = 1.2 * Math.sqrt(centroidHelper.radiusSq);
  120. return {
  121. ...props,
  122. structure,
  123. offsets,
  124. exposed,
  125. hydrophobic,
  126. centroid,
  127. extent,
  128. large: offsets.length > LARGE_CA_THRESHOLD
  129. };
  130. }
  131. export async function calculate(runtime: RuntimeContext, structure: Structure, params: ANVILProps): Promise<MembraneOrientation> {
  132. // can't get away with the default 92 points here
  133. const asaProps = { ...PD.getDefaultValues(AccessibleSurfaceAreaParams), probeSize: 4.0, traceOnly: true, numberOfSpherePoints: 184 };
  134. const accessibleSurfaceArea = await AccessibleSurfaceArea.compute(structure, asaProps).runInContext(runtime);
  135. const ctx = await initialize(structure, params, accessibleSurfaceArea);
  136. const initialHphobHphil = HphobHphil.initial(ctx);
  137. const initialMembrane = (await findMembrane(runtime, 'Placing initial membrane...', ctx, generateSpherePoints(ctx, ctx.numberOfSpherePoints), initialHphobHphil))!;
  138. const refinedMembrane = (await findMembrane(runtime, 'Refining membrane placement...', ctx, findProximateAxes(ctx, initialMembrane), initialHphobHphil))!;
  139. let membrane = initialMembrane.qmax! > refinedMembrane.qmax! ? initialMembrane : refinedMembrane;
  140. if (ctx.adjust && !ctx.large) {
  141. membrane = await adjustThickness(runtime, 'Adjusting membrane thickness...', ctx, membrane, initialHphobHphil);
  142. }
  143. const normalVector = v3zero();
  144. const center = v3zero();
  145. v3sub(normalVector, membrane.planePoint1, membrane.planePoint2);
  146. v3normalize(normalVector, normalVector);
  147. v3add(center, membrane.planePoint1, membrane.planePoint2);
  148. v3scale(center, center, 0.5);
  149. const extent = adjustExtent(ctx, membrane, center);
  150. return {
  151. planePoint1: membrane.planePoint1,
  152. planePoint2: membrane.planePoint2,
  153. normalVector,
  154. centroid: center,
  155. radius: extent
  156. };
  157. }
  158. interface MembraneCandidate {
  159. planePoint1: Vec3,
  160. planePoint2: Vec3,
  161. stats: HphobHphil,
  162. normalVector?: Vec3,
  163. spherePoint?: Vec3,
  164. qmax?: number
  165. }
  166. namespace MembraneCandidate {
  167. export function initial(c1: Vec3, c2: Vec3, stats: HphobHphil): MembraneCandidate {
  168. return {
  169. planePoint1: c1,
  170. planePoint2: c2,
  171. stats
  172. };
  173. }
  174. export function scored(spherePoint: Vec3, planePoint1: Vec3, planePoint2: Vec3, stats: HphobHphil, qmax: number, centroid: Vec3): MembraneCandidate {
  175. const normalVector = v3zero();
  176. v3sub(normalVector, centroid, spherePoint);
  177. return {
  178. planePoint1,
  179. planePoint2,
  180. stats,
  181. normalVector,
  182. spherePoint,
  183. qmax
  184. };
  185. }
  186. }
  187. async function findMembrane(runtime: RuntimeContext, message: string | undefined, ctx: ANVILContext, spherePoints: Vec3[], initialStats: HphobHphil): Promise<MembraneCandidate | undefined> {
  188. const { centroid, stepSize, minThickness, maxThickness, large } = ctx;
  189. // best performing membrane
  190. let membrane: MembraneCandidate | undefined;
  191. // score of the best performing membrane
  192. let qmax = 0;
  193. // construct slices of thickness 1.0 along the axis connecting the centroid and the spherePoint
  194. const diam = v3zero();
  195. for (let n = 0, nl = spherePoints.length; n < nl; n++) {
  196. if (runtime.shouldUpdate && message && (n + 1) % (large ? LARGE_CA_UPDATE_INTERVAL : DEFAULT_UPDATE_INTERVAL) === 0) {
  197. await runtime.update({ message, current: (n + 1), max: nl });
  198. }
  199. const spherePoint = spherePoints[n];
  200. v3sub(diam, centroid, spherePoint);
  201. v3scale(diam, diam, 2);
  202. const diamNorm = v3magnitude(diam);
  203. const sliceStats = HphobHphil.sliced(ctx, stepSize, spherePoint, diam, diamNorm);
  204. const qvartemp = [];
  205. for (let i = 0, il = diamNorm - stepSize; i < il; i += stepSize) {
  206. const c1 = v3zero();
  207. const c2 = v3zero();
  208. v3scaleAndAdd(c1, spherePoint, diam, i / diamNorm);
  209. v3scaleAndAdd(c2, spherePoint, diam, (i + stepSize) / diamNorm);
  210. // evaluate how well this membrane slice embeddeds the peculiar residues
  211. const stats = sliceStats[Math.round(i / stepSize)];
  212. qvartemp.push(MembraneCandidate.initial(c1, c2, stats));
  213. }
  214. let jmax = Math.floor((minThickness / stepSize) - 1);
  215. for (let width = 0, widthl = maxThickness; width <= widthl;) {
  216. for (let i = 0, il = qvartemp.length - 1 - jmax; i < il; i++) {
  217. let hphob = 0;
  218. let hphil = 0;
  219. for (let j = 0; j < jmax; j++) {
  220. const ij = qvartemp[i + j];
  221. if (j === 0 || j === jmax - 1) {
  222. hphob += Math.floor(0.5 * ij.stats.hphob);
  223. hphil += 0.5 * ij.stats.hphil;
  224. } else {
  225. hphob += ij.stats.hphob;
  226. hphil += ij.stats.hphil;
  227. }
  228. }
  229. if (hphob !== 0) {
  230. const stats = { hphob, hphil };
  231. const qvaltest = qValue(stats, initialStats);
  232. if (qvaltest >= qmax) {
  233. qmax = qvaltest;
  234. membrane = MembraneCandidate.scored(spherePoint, qvartemp[i].planePoint1, qvartemp[i + jmax].planePoint2, stats, qmax, centroid);
  235. }
  236. }
  237. }
  238. jmax++;
  239. width = (jmax + 1) * stepSize;
  240. }
  241. }
  242. return membrane;
  243. }
  244. /** Adjust membrane thickness by maximizing the number of membrane segments. */
  245. async function adjustThickness(runtime: RuntimeContext, message: string | undefined, ctx: ANVILContext, membrane: MembraneCandidate, initialHphobHphil: HphobHphil): Promise<MembraneCandidate> {
  246. const { minThickness, large } = ctx;
  247. const step = 0.3;
  248. let maxThickness = v3distance(membrane.planePoint1, membrane.planePoint2);
  249. let maxNos = membraneSegments(ctx, membrane).length;
  250. let optimalThickness = membrane;
  251. let n = 0;
  252. const nl = Math.ceil((maxThickness - minThickness) / step);
  253. while (maxThickness > minThickness) {
  254. n++;
  255. if (runtime.shouldUpdate && message && n % (large ? LARGE_CA_UPDATE_INTERVAL : DEFAULT_UPDATE_INTERVAL) === 0) {
  256. await runtime.update({ message, current: n, max: nl });
  257. }
  258. const p = {
  259. ...ctx,
  260. maxThickness,
  261. stepSize: step
  262. };
  263. const temp = await findMembrane(runtime, void 0, p, [membrane.spherePoint!], initialHphobHphil);
  264. if (temp) {
  265. const nos = membraneSegments(ctx, temp).length;
  266. if (nos > maxNos) {
  267. maxNos = nos;
  268. optimalThickness = temp;
  269. }
  270. }
  271. maxThickness -= step;
  272. }
  273. return optimalThickness;
  274. }
  275. /** Report auth_seq_ids for all transmembrane segments. Will reject segments that are shorter than the adjust parameter specifies. Missing residues are considered in-membrane. */
  276. function membraneSegments(ctx: ANVILContext, membrane: MembraneCandidate): ArrayLike<{ start: number, end: number }> {
  277. const { offsets, structure, adjust } = ctx;
  278. const { normalVector, planePoint1, planePoint2 } = membrane;
  279. const { units } = structure;
  280. const { elementIndices, unitIndices } = structure.serialMapping;
  281. const testPoint = v3zero();
  282. const { auth_seq_id } = StructureProperties.residue;
  283. const d1 = -v3dot(normalVector!, planePoint1);
  284. const d2 = -v3dot(normalVector!, planePoint2);
  285. const dMin = Math.min(d1, d2);
  286. const dMax = Math.max(d1, d2);
  287. const inMembrane: { [k: string]: Set<number> } = Object.create(null);
  288. const outMembrane: { [k: string]: Set<number> } = Object.create(null);
  289. const segments: Array<{ start: number, end: number }> = [];
  290. let authAsymId;
  291. let lastAuthAsymId = null;
  292. let authSeqId;
  293. let lastAuthSeqId = units[0].model.atomicHierarchy.residues.auth_seq_id.value((units[0] as Unit.Atomic).chainIndex[0]) - 1;
  294. let startOffset = 0;
  295. let endOffset = 0;
  296. // collect all residues in membrane layer
  297. for (let k = 0, kl = offsets.length; k < kl; k++) {
  298. const unit = units[unitIndices[offsets[k]]];
  299. if (!Unit.isAtomic(unit)) throw 'Property only available for atomic models.';
  300. const elementIndex = elementIndices[offsets[k]];
  301. authAsymId = unit.model.atomicHierarchy.chains.auth_asym_id.value(unit.chainIndex[elementIndex]);
  302. if (authAsymId !== lastAuthAsymId) {
  303. if (!inMembrane[authAsymId]) inMembrane[authAsymId] = new Set<number>();
  304. if (!outMembrane[authAsymId]) outMembrane[authAsymId] = new Set<number>();
  305. lastAuthAsymId = authAsymId;
  306. }
  307. authSeqId = unit.model.atomicHierarchy.residues.auth_seq_id.value(unit.residueIndex[elementIndex]);
  308. v3set(testPoint, unit.conformation.x(elementIndex), unit.conformation.y(elementIndex), unit.conformation.z(elementIndex));
  309. if (_isInMembranePlane(testPoint, normalVector!, dMin, dMax)) {
  310. inMembrane[authAsymId].add(authSeqId);
  311. } else {
  312. outMembrane[authAsymId].add(authSeqId);
  313. }
  314. }
  315. for (let k = 0, kl = offsets.length; k < kl; k++) {
  316. const unit = units[unitIndices[offsets[k]]];
  317. if (!Unit.isAtomic(unit)) throw 'Property only available for atomic models.';
  318. const elementIndex = elementIndices[offsets[k]];
  319. authAsymId = unit.model.atomicHierarchy.chains.auth_asym_id.value(unit.chainIndex[elementIndex]);
  320. authSeqId = unit.model.atomicHierarchy.residues.auth_seq_id.value(unit.residueIndex[elementIndex]);
  321. if (inMembrane[authAsymId].has(authSeqId)) {
  322. // chain change
  323. if (authAsymId !== lastAuthAsymId) {
  324. segments.push({ start: startOffset, end: endOffset });
  325. lastAuthAsymId = authAsymId;
  326. startOffset = k;
  327. endOffset = k;
  328. }
  329. // sequence gaps
  330. if (authSeqId !== lastAuthSeqId + 1) {
  331. if (outMembrane[authAsymId].has(lastAuthSeqId + 1)) {
  332. segments.push({ start: startOffset, end: endOffset });
  333. startOffset = k;
  334. }
  335. lastAuthSeqId = authSeqId;
  336. endOffset = k;
  337. } else {
  338. lastAuthSeqId++;
  339. endOffset++;
  340. }
  341. }
  342. }
  343. segments.push({ start: startOffset, end: endOffset });
  344. const l = StructureElement.Location.create(structure);
  345. let startAuth;
  346. let endAuth;
  347. const refinedSegments: Array<{ start: number, end: number }> = [];
  348. for (let k = 0, kl = segments.length; k < kl; k++) {
  349. const { start, end } = segments[k];
  350. if (start === 0 || end === offsets.length - 1) continue;
  351. // evaluate residues 1 pos outside of membrane
  352. setLocation(l, structure, offsets[start - 1]);
  353. v3set(testPoint, l.unit.conformation.x(l.element), l.unit.conformation.y(l.element), l.unit.conformation.z(l.element));
  354. const d3 = -v3dot(normalVector!, testPoint);
  355. setLocation(l, structure, offsets[end + 1]);
  356. v3set(testPoint, l.unit.conformation.x(l.element), l.unit.conformation.y(l.element), l.unit.conformation.z(l.element));
  357. const d4 = -v3dot(normalVector!, testPoint);
  358. if (Math.min(d3, d4) < dMin && Math.max(d3, d4) > dMax) {
  359. // reject this refinement
  360. setLocation(l, structure, offsets[start]);
  361. startAuth = auth_seq_id(l);
  362. setLocation(l, structure, offsets[end]);
  363. endAuth = auth_seq_id(l);
  364. if (Math.abs(startAuth - endAuth) + 1 < adjust) {
  365. return [];
  366. }
  367. refinedSegments.push(segments[k]);
  368. }
  369. }
  370. return refinedSegments;
  371. }
  372. /** Filter for membrane residues and calculate the final extent of the membrane layer */
  373. function adjustExtent(ctx: ANVILContext, membrane: MembraneCandidate, centroid: Vec3): number {
  374. const { offsets, structure } = ctx;
  375. const { normalVector, planePoint1, planePoint2 } = membrane;
  376. const l = StructureElement.Location.create(structure);
  377. const testPoint = v3zero();
  378. const { x, y, z } = StructureProperties.atom;
  379. const d1 = -v3dot(normalVector!, planePoint1);
  380. const d2 = -v3dot(normalVector!, planePoint2);
  381. const dMin = Math.min(d1, d2);
  382. const dMax = Math.max(d1, d2);
  383. let extent = 0;
  384. for (let k = 0, kl = offsets.length; k < kl; k++) {
  385. setLocation(l, structure, offsets[k]);
  386. v3set(testPoint, x(l), y(l), z(l));
  387. if (_isInMembranePlane(testPoint, normalVector!, dMin, dMax)) {
  388. const dsq = v3squaredDistance(testPoint, centroid);
  389. if (dsq > extent) extent = dsq;
  390. }
  391. }
  392. return Math.sqrt(extent);
  393. }
  394. function qValue(currentStats: HphobHphil, initialStats: HphobHphil): number {
  395. if(initialStats.hphob < 1) {
  396. initialStats.hphob = 0.1;
  397. }
  398. if(initialStats.hphil < 1) {
  399. initialStats.hphil += 1;
  400. }
  401. const part_tot = currentStats.hphob + currentStats.hphil;
  402. return (currentStats.hphob * (initialStats.hphil - currentStats.hphil) - currentStats.hphil * (initialStats.hphob - currentStats.hphob)) /
  403. Math.sqrt(part_tot * initialStats.hphob * initialStats.hphil * (initialStats.hphob + initialStats.hphil - part_tot));
  404. }
  405. export function isInMembranePlane(testPoint: Vec3, normalVector: Vec3, planePoint1: Vec3, planePoint2: Vec3): boolean {
  406. const d1 = -v3dot(normalVector, planePoint1);
  407. const d2 = -v3dot(normalVector, planePoint2);
  408. return _isInMembranePlane(testPoint, normalVector, Math.min(d1, d2), Math.max(d1, d2));
  409. }
  410. function _isInMembranePlane(testPoint: Vec3, normalVector: Vec3, min: number, max: number): boolean {
  411. const d = -v3dot(normalVector, testPoint);
  412. return d > min && d < max;
  413. }
  414. /** Generates a defined number of points on a sphere with radius = extent around the specified centroid */
  415. function generateSpherePoints(ctx: ANVILContext, numberOfSpherePoints: number): Vec3[] {
  416. const { centroid, extent } = ctx;
  417. const points = [];
  418. let oldPhi = 0, h, theta, phi;
  419. for(let k = 1, kl = numberOfSpherePoints + 1; k < kl; k++) {
  420. h = -1 + 2 * (k - 1) / (2 * numberOfSpherePoints - 1);
  421. theta = Math.acos(h);
  422. phi = (k === 1 || k === numberOfSpherePoints) ? 0 : (oldPhi + 3.6 / Math.sqrt(2 * numberOfSpherePoints * (1 - h * h))) % (2 * Math.PI);
  423. const point = v3create(
  424. extent * Math.sin(phi) * Math.sin(theta) + centroid[0],
  425. extent * Math.cos(theta) + centroid[1],
  426. extent * Math.cos(phi) * Math.sin(theta) + centroid[2]
  427. );
  428. points[k - 1] = point;
  429. oldPhi = phi;
  430. }
  431. return points;
  432. }
  433. /** Generates sphere points close to that of the initial membrane */
  434. function findProximateAxes(ctx: ANVILContext, membrane: MembraneCandidate): Vec3[] {
  435. const { numberOfSpherePoints, extent } = ctx;
  436. const points = generateSpherePoints(ctx, 30000);
  437. let j = 4;
  438. let sphere_pts2: Vec3[] = [];
  439. const s = 2 * extent / numberOfSpherePoints;
  440. while (sphere_pts2.length < numberOfSpherePoints) {
  441. const dsq = (s + j) * (s + j);
  442. sphere_pts2 = [];
  443. for (let i = 0, il = points.length; i < il; i++) {
  444. if (v3squaredDistance(points[i], membrane.spherePoint!) < dsq) {
  445. sphere_pts2.push(points[i]);
  446. }
  447. }
  448. j += 0.2;
  449. }
  450. return sphere_pts2;
  451. }
  452. interface HphobHphil {
  453. hphob: number,
  454. hphil: number
  455. }
  456. namespace HphobHphil {
  457. export function initial(ctx: ANVILContext): HphobHphil {
  458. const { exposed, hydrophobic } = ctx;
  459. let hphob = 0;
  460. let hphil = 0;
  461. for (let k = 0, kl = exposed.length; k < kl; k++) {
  462. if (hydrophobic[k]) {
  463. hphob++;
  464. } else {
  465. hphil++;
  466. }
  467. }
  468. return { hphob, hphil };
  469. }
  470. const testPoint = v3zero();
  471. export function sliced(ctx: ANVILContext, stepSize: number, spherePoint: Vec3, diam: Vec3, diamNorm: number): HphobHphil[] {
  472. const { exposed, hydrophobic, structure } = ctx;
  473. const { units, serialMapping } = structure;
  474. const { unitIndices, elementIndices } = serialMapping;
  475. const sliceStats: HphobHphil[] = [];
  476. for (let i = 0, il = diamNorm - stepSize; i < il; i += stepSize) {
  477. sliceStats[sliceStats.length] = { hphob: 0, hphil: 0 };
  478. }
  479. for (let i = 0, il = exposed.length; i < il; i++) {
  480. const unit = units[unitIndices[exposed[i]]];
  481. const elementIndex = elementIndices[exposed[i]];
  482. v3set(testPoint, unit.conformation.x(elementIndex), unit.conformation.y(elementIndex), unit.conformation.z(elementIndex));
  483. v3sub(testPoint, testPoint, spherePoint);
  484. if (hydrophobic[i]) {
  485. sliceStats[Math.floor(v3dot(testPoint, diam) / diamNorm / stepSize)].hphob++;
  486. } else {
  487. sliceStats[Math.floor(v3dot(testPoint, diam) / diamNorm / stepSize)].hphil++;
  488. }
  489. }
  490. return sliceStats;
  491. }
  492. }
  493. /** ANVIL-specific (not general) definition of membrane-favoring amino acids */
  494. const HYDROPHOBIC_AMINO_ACIDS = new Set(['ALA', 'CYS', 'GLY', 'HIS', 'ILE', 'LEU', 'MET', 'PHE', 'SER', 'TRP', 'VAL']);
  495. /** Returns true if ANVIL considers this as amino acid that favors being embedded in a membrane */
  496. export function isHydrophobic(label_comp_id: string): boolean {
  497. return HYDROPHOBIC_AMINO_ACIDS.has(label_comp_id);
  498. }
  499. /** Accessible surface area used for normalization. ANVIL uses 'Total-Side REL' values from NACCESS, from: Hubbard, S. J., & Thornton, J. M. (1993). naccess. Computer Program, Department of Biochemistry and Molecular Biology, University College London, 2(1). */
  500. export const MaxAsa: { [k: string]: number } = {
  501. 'ALA': 69.41,
  502. 'ARG': 201.25,
  503. 'ASN': 106.24,
  504. 'ASP': 102.69,
  505. 'CYS': 96.75,
  506. 'GLU': 134.74,
  507. 'GLN': 140.99,
  508. 'GLY': 32.33,
  509. 'HIS': 147.08,
  510. 'ILE': 137.96,
  511. 'LEU': 141.12,
  512. 'LYS': 163.30,
  513. 'MET': 156.64,
  514. 'PHE': 164.11,
  515. 'PRO': 119.90,
  516. 'SER': 78.11,
  517. 'THR': 101.70,
  518. 'TRP': 211.26,
  519. 'TYR': 177.38,
  520. 'VAL': 114.28
  521. };
  522. function setLocation(l: StructureElement.Location, structure: Structure, serialIndex: number) {
  523. l.structure = structure;
  524. l.unit = structure.units[structure.serialMapping.unitIndices[serialIndex]];
  525. l.element = structure.serialMapping.elementIndices[serialIndex];
  526. return l;
  527. }