diff --git a/src/mol-model-props/computed/membrane-orientation.ts b/src/mol-model-props/computed/membrane-orientation.ts index bcbb939d03b3b4a44f684acc7d57b328da10bbe6..bd765d3be76be56bb28130eb8e1cf89f74bf4618 100644 --- a/src/mol-model-props/computed/membrane-orientation.ts +++ b/src/mol-model-props/computed/membrane-orientation.ts @@ -10,8 +10,9 @@ import { Structure } from '../../mol-model/structure'; import { CustomStructureProperty } from '../common/custom-structure-property'; import { CustomProperty } from '../common/custom-property'; import { CustomPropertyDescriptor } from '../../mol-model/custom-property'; -import { ANVILParams, MembraneOrientation, ANVILProps } from './membrane-orientation/ANVIL'; +import { ANVILParams, ANVILProps, computeANVIL } from './membrane-orientation/ANVIL'; import { AccessibleSurfaceAreaProvider } from './accessible-surface-area'; +import { MembraneOrientation } from '../../mol-model/structure/model/properties/membrane-orientation'; function getMembraneOrientationParams(data?: Structure) { let defaultType = 'anvil' as 'anvil' | 'opm'; // TODO flip - OPM should be default if PDB identifier is known @@ -49,7 +50,7 @@ export const MembraneOrientationProvider: CustomStructureProperty.Provider<Membr async function computeAnvil(ctx: CustomProperty.Context, data: Structure, props: ANVILProps): Promise<MembraneOrientation> { await AccessibleSurfaceAreaProvider.attach(ctx, data); const p = { ...PD.getDefaultValues(MembraneOrientationParams), ...props }; - return await MembraneOrientation.compute(data, p).runInContext(ctx.runtime); + return await computeANVIL(data, p).runInContext(ctx.runtime); } async function computeOpm(structure: Structure): Promise<MembraneOrientation> { diff --git a/src/mol-model-props/computed/membrane-orientation/ANVIL.ts b/src/mol-model-props/computed/membrane-orientation/ANVIL.ts index 5e0ccc1336940c8627dc91892e2de14ef1f72964..500a345ef969e5d7cdc2734d97b9d40a7be98109 100644 --- a/src/mol-model-props/computed/membrane-orientation/ANVIL.ts +++ b/src/mol-model-props/computed/membrane-orientation/ANVIL.ts @@ -15,6 +15,7 @@ import { Vec3 } from '../../../mol-math/linear-algebra'; import { AccessibleSurfaceArea } from '../accessible-surface-area/shrake-rupley'; import { AccessibleSurfaceAreaProvider } from '../accessible-surface-area'; import { ANVILContext } from './anvil/common'; +import { MembraneOrientation } from '../../../mol-model/structure/model/properties/membrane-orientation'; export const ANVILParams = { numberOfSpherePoints: PD.Numeric(120, { min: 35, max: 700, step: 1 }, { description: 'Number of spheres/directions to test for membrane placement. Original value is 350.' }), @@ -27,379 +28,331 @@ export const ANVILParams = { export type ANVILParams = typeof ANVILParams export type ANVILProps = PD.Values<ANVILParams> -export { MembraneOrientation }; -interface MembraneOrientation extends Array<Vec3> {} - -namespace MembraneOrientation { - /** - * Implements: - * Membrane positioning for high- and low-resolution protein structures through a binary classification approach - * Guillaume Postic, Yassine Ghouzam, Vincent Guiraud, and Jean-Christophe Gelly - * Protein Engineering, Design & Selection, 2015, 1–5 - * doi: 10.1093/protein/gzv063 - */ - export function compute(structure: Structure, props: Partial<ANVILProps> = {}) { - const p = { ...PD.getDefaultValues(ANVILParams), ...props }; - return Task.create('Compute Membrane Topology', async runtime => { - return await calculate(runtime, structure, p); - }); - } +/** + * Implements: + * Membrane positioning for high- and low-resolution protein structures through a binary classification approach + * Guillaume Postic, Yassine Ghouzam, Vincent Guiraud, and Jean-Christophe Gelly + * Protein Engineering, Design & Selection, 2015, 1–5 + * doi: 10.1093/protein/gzv063 + */ +export function computeANVIL(structure: Structure, props: Partial<ANVILProps> = {}) { + const p = { ...PD.getDefaultValues(ANVILParams), ...props }; + return Task.create('Compute Membrane Topology', async runtime => { + return await calculate(runtime, structure, p); + }); +} + +const l = StructureElement.Location.create(void 0); +const centroidHelper = new CentroidHelper(); +function initialize(structure: Structure, params: ANVILProps): ANVILContext { + const { label_atom_id, x, y, z } = StructureProperties.atom; + const elementCount = structure.polymerResidueCount; + centroidHelper.reset(); + l.structure = structure; + + let offsets = new Int32Array(elementCount); + let exposed: boolean[] = new Array<boolean>(elementCount); + + const accessibleSurfaceArea = structure && AccessibleSurfaceAreaProvider.get(structure); + const asa = accessibleSurfaceArea.value!; - const l = StructureElement.Location.create(void 0); - const centroidHelper = new CentroidHelper(); const vec = Vec3(); - function initialize(structure: Structure, params: ANVILProps): ANVILContext { - const { label_atom_id, x, y, z } = StructureProperties.atom; - const elementCount = structure.polymerResidueCount; - centroidHelper.reset(); - l.structure = structure; - - let offsets = new Int32Array(elementCount); - let exposed: boolean[] = new Array<boolean>(elementCount); - - // ensure ASA - const accessibleSurfaceArea = structure && AccessibleSurfaceAreaProvider.get(structure); - const asa = accessibleSurfaceArea.value!; - - let m = 0; - for (let i = 0, il = structure.units.length; i < il; ++i) { - const unit = structure.units[i]; - const { elements } = unit; - l.unit = unit; - - for (let j = 0, jl = elements.length; j < jl; ++j) { - const eI = elements[j]; - l.element = eI; - - // consider only amino acids - if (getElementMoleculeType(unit, eI) !== MoleculeType.Protein) { - continue; - } - - // only CA is considered for downstream operations - if (label_atom_id(l) !== 'CA') { - continue; - } - - // while iterating use first pass to compute centroid - Vec3.set(vec, x(l), y(l), z(l)); - centroidHelper.includeStep(vec); - - // keep track of offsets and exposed state to reuse - offsets[m] = l.element; - exposed[m] = AccessibleSurfaceArea.getValue(l, asa) > params.afilter; - - m++; + let m = 0; + for (let i = 0, il = structure.units.length; i < il; ++i) { + const unit = structure.units[i]; + const { elements } = unit; + l.unit = unit; + + for (let j = 0, jl = elements.length; j < jl; ++j) { + const eI = elements[j]; + l.element = eI; + + // consider only amino acids + if (getElementMoleculeType(unit, eI) !== MoleculeType.Protein) { + continue; } - } - // omit potentially empty tail - offsets = offsets.slice(0, m); - exposed = exposed.slice(0, m); + // only CA is considered for downstream operations + if (label_atom_id(l) !== 'CA') { + continue; + } - // calculate centroid and extent - centroidHelper.finishedIncludeStep(); - const centroid = centroidHelper.center; - for (let k = 0, kl = offsets.length; k < kl; k++) { - setLocation(l, structure, offsets[k]); + // while iterating use first pass to compute centroid Vec3.set(vec, x(l), y(l), z(l)); - centroidHelper.radiusStep(vec); - } - const extent = 1.2 * Math.sqrt(centroidHelper.radiusSq); + centroidHelper.includeStep(vec); - return { - structure: structure, - numberOfSpherePoints: params.numberOfSpherePoints, - stepSize: params.stepSize, - minThickness: params.maxThickness, - maxThickness: params.minThickness, - afilter: params.afilter, - membranePointDensity: params.membranePointDensity, - - offsets: offsets, - exposed: exposed, - centroid: centroid, - extent: extent - }; - } + // keep track of offsets and exposed state to reuse + offsets[m] = l.element; + exposed[m] = AccessibleSurfaceArea.getValue(l, asa) > params.afilter; - export async function calculate(runtime: RuntimeContext, structure: Structure, params: ANVILProps): Promise<MembraneOrientation> { - const { label_comp_id } = StructureProperties.residue; - - const ctx = initialize(structure, params); - const initialHphobHphil = HphobHphil.filtered(ctx, label_comp_id); - - const initialMembrane = findMembrane(ctx, generateSpherePoints(ctx, ctx.numberOfSpherePoints), initialHphobHphil, label_comp_id); - const alternativeMembrane = findMembrane(ctx, findProximateAxes(ctx, initialMembrane), initialHphobHphil, label_comp_id); - - const membrane = initialMembrane.qmax! > alternativeMembrane.qmax! ? initialMembrane : alternativeMembrane; - - return createMembraneLayers(ctx, membrane); + m++; + } } - function createMembraneLayers(ctx: ANVILContext, membrane: MembraneCandidate): Vec3[] { - const { extent, membranePointDensity } = ctx; - const out: Vec3[] = []; - const radius = extent * extent; - const normalVector = membrane.normalVector!; - - createMembraneLayer(ctx, out, membrane.planePoint1, normalVector, membranePointDensity, radius); - createMembraneLayer(ctx, out, membrane.planePoint2, normalVector, membranePointDensity, radius); - - return out; + // omit potentially empty tail + offsets = offsets.slice(0, m); + exposed = exposed.slice(0, m); + + // calculate centroid and extent + centroidHelper.finishedIncludeStep(); + const centroid = centroidHelper.center; + for (let k = 0, kl = offsets.length; k < kl; k++) { + setLocation(l, structure, offsets[k]); + Vec3.set(vec, x(l), y(l), z(l)); + centroidHelper.radiusStep(vec); } - - function createMembraneLayer(ctx: ANVILContext, out: Vec3[], point: Vec3, normalVector: Vec3, density: number, radius: number) { - const { x, y, z } = StructureProperties.atom; - const { offsets, structure } = ctx; - const test = Vec3(); - - const d = -Vec3.dot(normalVector, point); - for (let i = -1000, il = 1000; i < il; i += density) { - for (let j = -1000, jl = 1000; j < jl; j += density) { - const rep = Vec3.create(i, j, -(d + i * normalVector[0] + j * normalVector[1]) / normalVector[2]); - if (Vec3.squaredDistance(rep, point) < radius) { - // it has a nice effect to ensure that membrane points don't overlap with structure representation - let valid = true; - for (let k = 0, kl = ctx.offsets.length; k < kl; k++) { - setLocation(l, structure, offsets[k]); - Vec3.set(test, x(l), y(l), z(l)); - const dsq = Vec3.squaredDistance(test, rep); - if (dsq < 16) { - valid = false; - } - } - if (valid) { - out.push(rep); - } - } - } - } + const extent = 1.2 * Math.sqrt(centroidHelper.radiusSq); + + return { + structure: structure, + numberOfSpherePoints: params.numberOfSpherePoints, + stepSize: params.stepSize, + minThickness: params.maxThickness, + maxThickness: params.minThickness, + afilter: params.afilter, + membranePointDensity: params.membranePointDensity, + + offsets: offsets, + exposed: exposed, + centroid: centroid, + extent: extent + }; +} + +export async function calculate(runtime: RuntimeContext, structure: Structure, params: ANVILProps): Promise<MembraneOrientation> { + const { label_comp_id } = StructureProperties.residue; + + const ctx = initialize(structure, params); + const initialHphobHphil = HphobHphil.filtered(ctx, label_comp_id); + + const initialMembrane = findMembrane(ctx, generateSpherePoints(ctx, ctx.numberOfSpherePoints), initialHphobHphil, label_comp_id); + const alternativeMembrane = findMembrane(ctx, findProximateAxes(ctx, initialMembrane), initialHphobHphil, label_comp_id); + + const membrane = initialMembrane.qmax! > alternativeMembrane.qmax! ? initialMembrane : alternativeMembrane; + + return MembraneOrientation(membrane.planePoint1, membrane.planePoint2, membrane.normalVector!, ctx.extent * ctx.extent); +} + +interface MembraneCandidate { + planePoint1: Vec3, + planePoint2: Vec3, + stats: HphobHphil, + normalVector?: Vec3, + spherePoint?: Vec3, + qmax?: number, + membraneAtoms?: Vec3[] +} + +namespace MembraneCandidate { + export function initial(c1: Vec3, c2: Vec3, stats: HphobHphil): MembraneCandidate { + return { + planePoint1: c1, + planePoint2: c2, + stats: stats + }; } - interface MembraneCandidate { - planePoint1: Vec3, - planePoint2: Vec3, - stats: HphobHphil, - normalVector?: Vec3, - spherePoint?: Vec3, - centroid?: Vec3, - qmax?: number, - membraneAtoms?: Vec3[] + export function scored(spherePoint: Vec3, c1: Vec3, c2: Vec3, stats: HphobHphil, qmax: number, centroid: Vec3): MembraneCandidate { + const diam_vect = Vec3(); + Vec3.sub(diam_vect, centroid, spherePoint); + return { + planePoint1: c1, + planePoint2: c2, + stats: stats, + normalVector: diam_vect, + spherePoint: spherePoint, + qmax: qmax, + membraneAtoms: [] + }; } - - namespace MembraneCandidate { - export function initial(c1: Vec3, c2: Vec3, stats: HphobHphil): MembraneCandidate { - return { - planePoint1: c1, - planePoint2: c2, - stats: stats - }; - } - - export function scored(spherePoint: Vec3, c1: Vec3, c2: Vec3, stats: HphobHphil, qmax: number, centroid: Vec3): MembraneCandidate { - const diam_vect = Vec3(); - Vec3.sub(diam_vect, centroid, spherePoint); - return { - planePoint1: c1, - planePoint2: c2, - stats: stats, - normalVector: diam_vect, - spherePoint: spherePoint, - centroid: centroid, - qmax: qmax, - membraneAtoms: [] - }; +} + +function findMembrane(ctx: ANVILContext, spherePoints: Vec3[], initialStats: HphobHphil, label_comp_id: StructureElement.Property<string>): MembraneCandidate { + const { centroid, stepSize, minThickness, maxThickness } = ctx; + // best performing membrane + let membrane: MembraneCandidate; + // score of the best performing membrane + let qmax = 0; + + // construct slices of thickness 1.0 along the axis connecting the centroid and the spherePoint + const diam = Vec3(); + for (let i = 0, il = spherePoints.length; i < il; i++) { + const spherePoint = spherePoints[i]; + Vec3.sub(diam, centroid, spherePoint); + Vec3.scale(diam, diam, 2); + const diamNorm = Vec3.magnitude(diam); + const qvartemp = []; + + for (let i = 0, il = diamNorm - stepSize; i < il; i += stepSize) { + const c1 = Vec3(); + const c2 = Vec3(); + Vec3.scaleAndAdd(c1, spherePoint, diam, i / diamNorm); + Vec3.scaleAndAdd(c2, spherePoint, diam, (i + stepSize) / diamNorm); + + // evaluate how well this membrane slice embeddeds the peculiar residues + const stats = HphobHphil.filtered(ctx, label_comp_id, (testPoint: Vec3) => isInMembranePlane(testPoint, diam, c1, c2)); + qvartemp.push(MembraneCandidate.initial(c1, c2, stats)); } - } - - function findMembrane(ctx: ANVILContext, spherePoints: Vec3[], initialStats: HphobHphil, label_comp_id: StructureElement.Property<string>): MembraneCandidate { - const { centroid, stepSize, minThickness, maxThickness } = ctx; - // best performing membrane - let membrane: MembraneCandidate; - // score of the best performing membrane - let qmax = 0; - - // construct slices of thickness 1.0 along the axis connecting the centroid and the spherePoint - const diam = Vec3(); - for (let i = 0, il = spherePoints.length; i < il; i++) { - const spherePoint = spherePoints[i]; - Vec3.sub(diam, centroid, spherePoint); - Vec3.scale(diam, diam, 2); - const diamNorm = Vec3.magnitude(diam); - const qvartemp = []; - - for (let i = 0, il = diamNorm - stepSize; i < il; i += stepSize) { - const c1 = Vec3(); - const c2 = Vec3(); - Vec3.scaleAndAdd(c1, spherePoint, diam, i / diamNorm); - Vec3.scaleAndAdd(c2, spherePoint, diam, (i + stepSize) / diamNorm); - - // evaluate how well this membrane slice embeddeds the peculiar residues - const stats = HphobHphil.filtered(ctx, label_comp_id, (testPoint: Vec3) => isInMembranePlane(testPoint, diam, c1, c2)); - qvartemp.push(MembraneCandidate.initial(c1, c2, stats)); - } - let jmax = (minThickness / stepSize) - 1; - - for (let width = 0, widthl = maxThickness; width < widthl;) { - const imax = qvartemp.length - 1 - jmax; - - for (let i = 0, il = imax; i < il; i++) { - const c1 = qvartemp[i].planePoint1; - const c2 = qvartemp[i + jmax].planePoint2; - - let hphob = 0; - let hphil = 0; - let total = 0; - for (let j = 0; j < jmax; j++) { - const ij = qvartemp[i + j]; - if (j === 0 || j === jmax - 1) { - hphob += 0.5 * ij.stats.hphob; - hphil += 0.5 * ij.stats.hphil; - } else { - hphob += ij.stats.hphob; - hphil += ij.stats.hphil; - } - total += ij.stats.total; + let jmax = (minThickness / stepSize) - 1; + + for (let width = 0, widthl = maxThickness; width < widthl;) { + const imax = qvartemp.length - 1 - jmax; + + for (let i = 0, il = imax; i < il; i++) { + const c1 = qvartemp[i].planePoint1; + const c2 = qvartemp[i + jmax].planePoint2; + + let hphob = 0; + let hphil = 0; + let total = 0; + for (let j = 0; j < jmax; j++) { + const ij = qvartemp[i + j]; + if (j === 0 || j === jmax - 1) { + hphob += 0.5 * ij.stats.hphob; + hphil += 0.5 * ij.stats.hphil; + } else { + hphob += ij.stats.hphob; + hphil += ij.stats.hphil; } + total += ij.stats.total; + } - const stats = HphobHphil.of(hphob, hphil, total); + const stats = HphobHphil.of(hphob, hphil, total); - if (hphob !== 0) { - const qvaltest = qValue(stats, initialStats); - if (qvaltest > qmax) { - qmax = qvaltest; - membrane = MembraneCandidate.scored(spherePoint, c1, c2, HphobHphil.of(hphob, hphil, total), qmax, centroid); - } + if (hphob !== 0) { + const qvaltest = qValue(stats, initialStats); + if (qvaltest > qmax) { + qmax = qvaltest; + membrane = MembraneCandidate.scored(spherePoint, c1, c2, HphobHphil.of(hphob, hphil, total), qmax, centroid); } } - jmax++; - width = (jmax + 1) * stepSize; } + jmax++; + width = (jmax + 1) * stepSize; } - - return membrane!; } - function qValue(currentStats: HphobHphil, initialStats: HphobHphil): number { - if(initialStats.hphob < 1) { - initialStats.hphob = 0.1; - } - - if(initialStats.hphil < 1) { - initialStats.hphil += 1; - } + return membrane!; +} - const part_tot = currentStats.hphob + currentStats.hphil; - return (currentStats.hphob * (initialStats.hphil - currentStats.hphil) - currentStats.hphil * (initialStats.hphob - currentStats.hphob)) / - Math.sqrt(part_tot * initialStats.hphob * initialStats.hphil * (initialStats.hphob + initialStats.hphil - part_tot)); +function qValue(currentStats: HphobHphil, initialStats: HphobHphil): number { + if(initialStats.hphob < 1) { + initialStats.hphob = 0.1; } - function isInMembranePlane(testPoint: Vec3, normalVector: Vec3, planePoint1: Vec3, planePoint2: Vec3): boolean { - const d1 = -Vec3.dot(normalVector, planePoint1); - const d2 = -Vec3.dot(normalVector, planePoint2); - const d = -Vec3.dot(normalVector, testPoint); - return d > Math.min(d1, d2) && d < Math.max(d1, d2); + if(initialStats.hphil < 1) { + initialStats.hphil += 1; } - // generates a defined number of points on a sphere with radius = extent around the specified centroid - function generateSpherePoints(ctx: ANVILContext, numberOfSpherePoints: number): Vec3[] { - const { centroid, extent } = ctx; - const points = []; - let oldPhi = 0, h, theta, phi; - for(let k = 1, kl = numberOfSpherePoints + 1; k < kl; k++) { - h = -1 + 2 * (k - 1) / (numberOfSpherePoints - 1); - theta = Math.acos(h); - phi = (k === 1 || k === numberOfSpherePoints) ? 0 : (oldPhi + 3.6 / Math.sqrt(numberOfSpherePoints * (1 - h * h))) % (2 * Math.PI); - - const point = Vec3.create( - extent * Math.sin(phi) * Math.sin(theta) + centroid[0], - extent * Math.cos(theta) + centroid[1], - extent * Math.cos(phi) * Math.sin(theta) + centroid[2] - ); - points[k - 1] = point; - oldPhi = phi; - } - - return points; + const part_tot = currentStats.hphob + currentStats.hphil; + return (currentStats.hphob * (initialStats.hphil - currentStats.hphil) - currentStats.hphil * (initialStats.hphob - currentStats.hphob)) / + Math.sqrt(part_tot * initialStats.hphob * initialStats.hphil * (initialStats.hphob + initialStats.hphil - part_tot)); +} + +function isInMembranePlane(testPoint: Vec3, normalVector: Vec3, planePoint1: Vec3, planePoint2: Vec3): boolean { + const d1 = -Vec3.dot(normalVector, planePoint1); + const d2 = -Vec3.dot(normalVector, planePoint2); + const d = -Vec3.dot(normalVector, testPoint); + return d > Math.min(d1, d2) && d < Math.max(d1, d2); +} + +// generates a defined number of points on a sphere with radius = extent around the specified centroid +function generateSpherePoints(ctx: ANVILContext, numberOfSpherePoints: number): Vec3[] { + const { centroid, extent } = ctx; + const points = []; + let oldPhi = 0, h, theta, phi; + for(let k = 1, kl = numberOfSpherePoints + 1; k < kl; k++) { + h = -1 + 2 * (k - 1) / (numberOfSpherePoints - 1); + theta = Math.acos(h); + phi = (k === 1 || k === numberOfSpherePoints) ? 0 : (oldPhi + 3.6 / Math.sqrt(numberOfSpherePoints * (1 - h * h))) % (2 * Math.PI); + + const point = Vec3.create( + extent * Math.sin(phi) * Math.sin(theta) + centroid[0], + extent * Math.cos(theta) + centroid[1], + extent * Math.cos(phi) * Math.sin(theta) + centroid[2] + ); + points[k - 1] = point; + oldPhi = phi; } - // generates sphere points close to that of the initial membrane - function findProximateAxes(ctx: ANVILContext, membrane: MembraneCandidate): Vec3[] { - const { numberOfSpherePoints, extent } = ctx; - const points = generateSpherePoints(ctx, 30000); - let j = 4; - let sphere_pts2: Vec3[] = []; - while (sphere_pts2.length < numberOfSpherePoints) { - const d = 2 * extent / numberOfSpherePoints + j; - const dsq = d * d; - sphere_pts2 = []; - for (let i = 0, il = points.length; i < il; i++) { - if (Vec3.squaredDistance(points[i], membrane.spherePoint!) < dsq) { - sphere_pts2.push(points[i]); - } + return points; +} + +// generates sphere points close to that of the initial membrane +function findProximateAxes(ctx: ANVILContext, membrane: MembraneCandidate): Vec3[] { + const { numberOfSpherePoints, extent } = ctx; + const points = generateSpherePoints(ctx, 30000); + let j = 4; + let sphere_pts2: Vec3[] = []; + while (sphere_pts2.length < numberOfSpherePoints) { + const d = 2 * extent / numberOfSpherePoints + j; + const dsq = d * d; + sphere_pts2 = []; + for (let i = 0, il = points.length; i < il; i++) { + if (Vec3.squaredDistance(points[i], membrane.spherePoint!) < dsq) { + sphere_pts2.push(points[i]); } - j += 0.2; } - return sphere_pts2; + j += 0.2; } + return sphere_pts2; +} - interface HphobHphil { - hphob: number, - hphil: number, - total: number - } - - namespace HphobHphil { - export function of(hphob: number, hphil: number, total?: number) { - return { - hphob: hphob, - hphil: hphil, - total: !!total ? total : hphob + hphil - }; - } +interface HphobHphil { + hphob: number, + hphil: number, + total: number +} - const testPoint = Vec3(); - export function filtered(ctx: ANVILContext, label_comp_id: StructureElement.Property<string>, filter?: (test: Vec3) => boolean): HphobHphil { - const { offsets, exposed, structure } = ctx; - const { x, y, z } = StructureProperties.atom; - let hphob = 0; - let hphil = 0; - for (let k = 0, kl = offsets.length; k < kl; k++) { - // ignore buried residues - if (!exposed[k]) { - continue; - } +namespace HphobHphil { + export function of(hphob: number, hphil: number, total?: number) { + return { + hphob: hphob, + hphil: hphil, + total: !!total ? total : hphob + hphil + }; + } - setLocation(l, structure, offsets[k]); - Vec3.set(testPoint, x(l), y(l), z(l)); + const testPoint = Vec3(); + export function filtered(ctx: ANVILContext, label_comp_id: StructureElement.Property<string>, filter?: (test: Vec3) => boolean): HphobHphil { + const { offsets, exposed, structure } = ctx; + const { x, y, z } = StructureProperties.atom; + let hphob = 0; + let hphil = 0; + for (let k = 0, kl = offsets.length; k < kl; k++) { + // ignore buried residues + if (!exposed[k]) { + continue; + } - // testPoints have to be in putative membrane layer - if (filter && !filter(testPoint)) { - continue; - } + setLocation(l, structure, offsets[k]); + Vec3.set(testPoint, x(l), y(l), z(l)); - if (isHydrophobic(label_comp_id(l))) { - hphob++; - } else { - hphil++; - } + // testPoints have to be in putative membrane layer + if (filter && !filter(testPoint)) { + continue; } - return of(hphob, hphil); - } - // ANVIL-specific (not general) definition of membrane-favoring amino acids - const HYDROPHOBIC_AMINO_ACIDS = ['ALA', 'CYS', 'GLY', 'HIS', 'ILE', 'LEU', 'MET', 'PHE', 'SER', 'THR', 'VAL']; - function isHydrophobic(label_comp_id: string): boolean { - return HYDROPHOBIC_AMINO_ACIDS.indexOf(label_comp_id) !== -1; + if (isHydrophobic(label_comp_id(l))) { + hphob++; + } else { + hphil++; + } } + return of(hphob, hphil); } - function setLocation(l: StructureElement.Location, structure: Structure, serialIndex: number) { - l.structure = structure; - l.unit = structure.units[structure.serialMapping.unitIndices[serialIndex]]; - l.element = structure.serialMapping.elementIndices[serialIndex]; - return l; + // ANVIL-specific (not general) definition of membrane-favoring amino acids + const HYDROPHOBIC_AMINO_ACIDS = ['ALA', 'CYS', 'GLY', 'HIS', 'ILE', 'LEU', 'MET', 'PHE', 'SER', 'THR', 'VAL']; + function isHydrophobic(label_comp_id: string): boolean { + return HYDROPHOBIC_AMINO_ACIDS.indexOf(label_comp_id) !== -1; } +} + +function setLocation(l: StructureElement.Location, structure: Structure, serialIndex: number) { + l.structure = structure; + l.unit = structure.units[structure.serialMapping.unitIndices[serialIndex]]; + l.element = structure.serialMapping.elementIndices[serialIndex]; + return l; } \ No newline at end of file diff --git a/src/mol-model/structure/model/properties/membrane-orientation.ts b/src/mol-model/structure/model/properties/membrane-orientation.ts new file mode 100644 index 0000000000000000000000000000000000000000..eb15083465fd683453b1746892cedb6729c7f652 --- /dev/null +++ b/src/mol-model/structure/model/properties/membrane-orientation.ts @@ -0,0 +1,24 @@ +/** + * Copyright (c) 2020 mol* contributors, licensed under MIT, See LICENSE file for more info. + * + * @author Sebastian Bittrich <sebastian.bittrich@rcsb.org> + */ + +import { Vec3 } from '../../../../mol-math/linear-algebra'; + +interface MembraneOrientation { + // point in membrane boundary + readonly p1: Vec3, + // point in opposite side of membrane boundary + readonly p2: Vec3, + // normal vector of membrane layer + readonly normal: Vec3, + // the radius of the membrane layer + readonly radius: number +} + +function MembraneOrientation(p1: Vec3, p2: Vec3, normal: Vec3, radius: number): MembraneOrientation { + return { p1, p2, normal, radius }; +} + +export { MembraneOrientation }; \ No newline at end of file diff --git a/src/tests/browser/render-structure.ts b/src/tests/browser/render-structure.ts index 45857da4208ae2ee79e9d14d90bfee9930de94d3..ed599fcc4a8d4d3b6184b32c15f98333fd91f9c0 100644 --- a/src/tests/browser/render-structure.ts +++ b/src/tests/browser/render-structure.ts @@ -31,7 +31,8 @@ import { SpheresBuilder } from '../../mol-geo/geometry/spheres/spheres-builder'; import { Spheres } from '../../mol-geo/geometry/spheres/spheres'; import { Color } from '../../mol-util/color'; import { createRenderObject } from '../../mol-gl/render-object'; -import { MembraneOrientation } from '../../mol-model-props/computed/membrane-orientation/ANVIL'; +import { MembraneOrientation } from '../../mol-model/structure/model/properties/membrane-orientation'; +import { Vec3 } from '../../mol-math/linear-algebra'; const parent = document.getElementById('app')!; parent.style.width = '100%'; @@ -122,12 +123,12 @@ function getGaussianSurfaceRepr() { return GaussianSurfaceRepresentationProvider.factory(reprCtx, GaussianSurfaceRepresentationProvider.getParams); } -function getMembraneRepr(membrane: MembraneOrientation) { - // TODO is a representation provider the right place for this? - const spheresBuilder = SpheresBuilder.create(membrane.length, 1); - for (let i = 0, il = membrane.length; i < il; i++) { - spheresBuilder.add(membrane[i][0], membrane[i][1], membrane[i][2], 0); - } +function getMembraneRepr(structure: Structure, membrane: MembraneOrientation) { + const density = 1; + const radius = membrane.radius; + const spheresBuilder = SpheresBuilder.create(); + createMembraneLayer(spheresBuilder, membrane.p1, membrane.normal, density, radius); + createMembraneLayer(spheresBuilder, membrane.p2, membrane.normal, density, radius); const spheres = spheresBuilder.getSpheres(); const values = Spheres.Utils.createValuesSimple(spheres, {}, Color(0xCCCCCC), 1); @@ -138,6 +139,19 @@ function getMembraneRepr(membrane: MembraneOrientation) { return repr; } +function createMembraneLayer(spheresBuilder: SpheresBuilder, point: Vec3, normalVector: Vec3, density: number, radius: number) { + const d = -Vec3.dot(normalVector, point); + const rep = Vec3(); + for (let i = -1000, il = 1000; i < il; i += density) { + for (let j = -1000, jl = 1000; j < jl; j += density) { + Vec3.set(rep, i, j, -(d + i * normalVector[0] + j * normalVector[1]) / normalVector[2]); + if (Vec3.squaredDistance(rep, point) < radius) { + spheresBuilder.add(rep[0], rep[1], rep[2], 0); + } + } + } +} + async function init() { const ctx = { runtime: SyncRuntimeContext, assetManager: new AssetManager() }; @@ -172,7 +186,7 @@ async function init() { const ballAndStickRepr = getBallAndStickRepr(); const molecularSurfaceRepr = getMolecularSurfaceRepr(); const gaussianSurfaceRepr = getGaussianSurfaceRepr(); - const membraneRepr = getMembraneRepr(MembraneOrientationProvider.get(structure).value!); + const membraneRepr = getMembraneRepr(structure, MembraneOrientationProvider.get(structure).value!); if (show.cartoon) { cartoonRepr.setTheme({