wip adjust thickness

src/extensions/anvil/algorithm.ts

@@ -24,6 +24,7 @@ interface ANVILContext {
     minThickness: number,
     maxThickness: number,
     asaCutoff: number,
+    adjust: number,
 
     offsets: ArrayLike<number>,
     exposed: ArrayLike<boolean>,
@@ -32,11 +33,12 @@ interface ANVILContext {
 };
 
 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.' }),
+    numberOfSpherePoints: PD.Numeric(350, { min: 35, max: 700, step: 1 }, { description: 'Number of spheres/directions to test for membrane placement. Original value is 350.' }),
     stepSize: PD.Numeric(1, { min: 0.25, max: 4, step: 0.25 }, { description: 'Thickness of membrane slices that will be tested' }),
     minThickness: PD.Numeric(20, { min: 10, max: 30, step: 1}, { description: 'Minimum membrane thickness used during refinement' }),
     maxThickness: PD.Numeric(40, { min: 30, max: 50, step: 1}, { description: 'Maximum membrane thickness used during refinement' }),
-    asaCutoff: PD.Numeric(40, { min: 10, max: 100, step: 1 }, { description: 'Absolute ASA cutoff above which residues will be considered' })
+    asaCutoff: PD.Numeric(40, { min: 10, max: 100, step: 1 }, { description: 'Relative ASA cutoff above which residues will be considered' }),
+    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.' })
 };
 export type ANVILParams = typeof ANVILParams
 export type ANVILProps = PD.Values<ANVILParams>
@@ -58,18 +60,20 @@ export function computeANVIL(structure: Structure, props: ANVILProps) {
 const centroidHelper = new CentroidHelper();
 function initialize(structure: Structure, props: ANVILProps): ANVILContext {
     const l = StructureElement.Location.create(structure);
-    const { label_atom_id, x, y, z } = StructureProperties.atom;
+    const { label_atom_id, label_comp_id, x, y, z } = StructureProperties.atom;
     const elementCount = structure.polymerResidueCount;
     centroidHelper.reset();
 
-    let offsets = new Int32Array(elementCount);
+    let offsets = new Array<number>(elementCount);
     let exposed = new Array<boolean>(elementCount);
 
     const accessibleSurfaceArea = structure && AccessibleSurfaceAreaProvider.get(structure);
     const asa = accessibleSurfaceArea.value!;
+    const asaCutoff = props.asaCutoff / 100;
 
     const vec = Vec3();
     let m = 0;
+    let e = 0, b = 0;
     for (let i = 0, il = structure.units.length; i < il; ++i) {
         const unit = structure.units[i];
         const { elements } = unit;
@@ -95,13 +99,19 @@ function initialize(structure: Structure, props: ANVILProps): ANVILContext {
 
             // keep track of offsets and exposed state to reuse
             offsets[m] = structure.serialMapping.getSerialIndex(l.unit, l.element);
-            exposed[m] = AccessibleSurfaceArea.getValue(l, asa) > props.asaCutoff;
-
+            exposed[m] = AccessibleSurfaceArea.getValue(l, asa) / MaxAsa[label_comp_id(l)] > asaCutoff;
+            if (AccessibleSurfaceArea.getValue(l, asa) / MaxAsa[label_comp_id(l)] > asaCutoff) {
+                e++;
+            } else {
+                b++;
+            }
             m++;
         }
     }
+    console.log('CAs = ' + m);
+    console.log('exposed ' + e + ' - buried ' + b);
 
-    // omit potentially empty tail1
+    // omit potentially empty tail
     offsets = offsets.slice(0, m);
     exposed = exposed.slice(0, m);
 
@@ -114,28 +124,52 @@ function initialize(structure: Structure, props: ANVILProps): ANVILContext {
         centroidHelper.radiusStep(vec);
     }
     const extent = 1.2 * Math.sqrt(centroidHelper.radiusSq);
+    console.log(`center: ${centroid}, radius: ${Math.sqrt(centroidHelper.radiusSq)}`);
 
     return {
         ...props,
-        structure: structure,
+        structure,
 
-        offsets: offsets,
-        exposed: exposed,
-        centroid: centroid,
-        extent: extent
+        offsets,
+        exposed,
+        centroid,
+        extent
     };
 }
 
 export async function calculate(runtime: RuntimeContext, structure: Structure, params: ANVILProps): Promise<MembraneOrientation> {
-    const { label_comp_id } = StructureProperties.atom;
-
     const ctx = initialize(structure, params);
-    const initialHphobHphil = HphobHphil.filtered(ctx, label_comp_id);
+    const initialHphobHphil = HphobHphil.filtered(ctx);
+    console.log(`init: ${initialHphobHphil.hphob} - ${initialHphobHphil.hphil}`);
+
+    if (runtime.shouldUpdate) {
+        await runtime.update({ message: 'Placing initial membrane...' });
+    }
+    console.time('ini-membrane');
+    const initialMembrane = findMembrane(ctx, generateSpherePoints(ctx, ctx.numberOfSpherePoints), initialHphobHphil)!;
+    console.timeEnd('ini-membrane');
+    console.log(`initial: ${initialMembrane.qmax}`);
+
+    if (runtime.shouldUpdate) {
+        await runtime.update({ message: 'Refining membrane placement...' });
+    }
+    console.time('ref-membrane');
+    const refinedMembrane = findMembrane(ctx, findProximateAxes(ctx, initialMembrane), initialHphobHphil)!;
+    console.timeEnd('ref-membrane');
+    console.log(`refined: ${refinedMembrane.qmax}`);
+    let membrane = initialMembrane.qmax! > refinedMembrane.qmax! ? initialMembrane : refinedMembrane;
+
+    if (ctx.adjust) {
+        if (runtime.shouldUpdate) {
+            await runtime.update({ message: 'Adjusting membrane thickness...' });
+        }
+        console.time('adj-thickness');
+        membrane = adjustThickness(ctx, membrane, initialHphobHphil);
+        console.timeEnd('adj-thickness');
+        console.log('Membrane width: ' + Vec3.distance(membrane.planePoint1, membrane.planePoint2));
+    }
 
-    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 {
         planePoint1: membrane.planePoint1,
@@ -166,7 +200,7 @@ namespace MembraneCandidate {
 
     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);
+        Vec3.sub(diam_vect, spherePoint, centroid);
         return {
             planePoint1: c1,
             planePoint2: c2,
@@ -178,36 +212,40 @@ namespace MembraneCandidate {
     }
 }
 
-function findMembrane(ctx: ANVILContext, spherePoints: Vec3[], initialStats: HphobHphil, label_comp_id: StructureElement.Property<string>): MembraneCandidate {
+function findMembrane(ctx: ANVILContext, spherePoints: Vec3[], initialStats: HphobHphil): MembraneCandidate | undefined {
     const { centroid, stepSize, minThickness, maxThickness } = ctx;
     // best performing membrane
-    let membrane: MembraneCandidate;
+    let membrane: MembraneCandidate | undefined;
     // 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];
+    for (let n = 0, nl = spherePoints.length; n < nl; n++) {
+        const spherePoint = spherePoints[n];
         Vec3.sub(diam, centroid, spherePoint);
         Vec3.scale(diam, diam, 2);
         const diamNorm = Vec3.magnitude(diam);
-        const qvartemp = [];
+        const qvartemp = []; // TODO use fixed length array
 
         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);
+            const d1 = -Vec3.dot(diam, c1);
+            const d2 = -Vec3.dot(diam, c2);
+            const dMin = Math.min(d1, d2);
+            const dMax = Math.max(d1, d2);
 
             // 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));
+            const stats = HphobHphil.filtered(ctx, (testPoint: Vec3) => _isInMembranePlane(testPoint, diam, dMin, dMax));
             qvartemp.push(MembraneCandidate.initial(c1, c2, stats));
         }
 
-        let jmax = (minThickness / stepSize) - 1;
+        let jmax = Math.floor((minThickness / stepSize) - 1);
 
-        for (let width = 0, widthl = maxThickness; width < widthl;) {
+        for (let width = 0, widthl = maxThickness; width <= widthl;) {
             const imax = qvartemp.length - 1 - jmax;
 
             for (let i = 0, il = imax; i < il; i++) {
@@ -220,7 +258,7 @@ function findMembrane(ctx: ANVILContext, spherePoints: Vec3[], initialStats: Hph
                 for (let j = 0; j < jmax; j++) {
                     const ij = qvartemp[i + j];
                     if (j === 0 || j === jmax - 1) {
-                        hphob += 0.5 * ij.stats.hphob;
+                        hphob += Math.floor(0.5 * ij.stats.hphob);
                         hphil += 0.5 * ij.stats.hphil;
                     } else {
                         hphob += ij.stats.hphob;
@@ -229,11 +267,10 @@ function findMembrane(ctx: ANVILContext, spherePoints: Vec3[], initialStats: Hph
                     total += ij.stats.total;
                 }
 
-                const stats = HphobHphil.of(hphob, hphil, total);
-
                 if (hphob !== 0) {
+                    const stats = HphobHphil.of(hphob, hphil, total);
                     const qvaltest = qValue(stats, initialStats);
-                    if (qvaltest > qmax) {
+                    if (qvaltest >= qmax) {
                         qmax = qvaltest;
                         membrane = MembraneCandidate.scored(spherePoint, c1, c2, HphobHphil.of(hphob, hphil, total), qmax, centroid);
                     }
@@ -244,7 +281,141 @@ function findMembrane(ctx: ANVILContext, spherePoints: Vec3[], initialStats: Hph
         }
     }
 
-    return membrane!;
+    return membrane;
+}
+
+function adjustThickness(ctx: ANVILContext, membrane: MembraneCandidate, initialHphobHphil: HphobHphil): MembraneCandidate {
+    const { minThickness } = ctx;
+    const step = 0.3;
+    let maxThickness = Vec3.distance(membrane.planePoint1, membrane.planePoint2);
+
+    let maxNos = membraneSegments(ctx, membrane).length;
+    let optimalThickness = membrane;
+
+    while (maxThickness > minThickness) {
+        const p = {
+            ...ctx,
+            maxThickness,
+            stepSize: step
+        };
+        const temp = findMembrane(p, [membrane.spherePoint!], initialHphobHphil);
+        if (temp) {
+            const nos = membraneSegments(ctx, temp).length;
+            if (nos > maxNos) {
+                maxNos = nos;
+                optimalThickness = temp;
+            }
+        }
+        maxThickness -= step;
+    }
+
+    console.log('Number of TM segments: ' + maxNos);
+    return optimalThickness;
+}
+
+const testPoint = Vec3();
+function membraneSegments(ctx: ANVILContext, membrane: MembraneCandidate): ArrayLike<{ start: number, end: number }> {
+    const { offsets, structure, adjust } = ctx;
+    const { normalVector, planePoint1, planePoint2 } = membrane;
+    const l = StructureElement.Location.create(structure);
+    const { auth_asym_id } = StructureProperties.chain;
+    const { auth_seq_id } = StructureProperties.residue;
+    const { x, y, z } = StructureProperties.atom;
+
+    const d1 = -Vec3.dot(normalVector!, planePoint1);
+    const d2 = -Vec3.dot(normalVector!, planePoint2);
+    const dMin = Math.min(d1, d2);
+    const dMax = Math.max(d1, d2);
+
+    const inMembrane: { [k: string]: Set<number> } = Object.create(null);
+    const outMembrane: { [k: string]: Set<number> } = Object.create(null);
+    const segments: Array<{ start: number, end: number }> = [];
+    setLocation(l, structure, offsets[0]);
+    let authAsymId;
+    let lastAuthAsymId = null;
+    let authSeqId;
+    let lastAuthSeqId = auth_seq_id(l) - 1;
+    let startOffset = 0;
+    let endOffset = 0;
+
+    // collect all residues in membrane layer
+    for (let k = 0, kl = offsets.length; k < kl; k++) {
+        setLocation(l, structure, offsets[k]);
+        authAsymId = auth_asym_id(l);
+        if (authAsymId !== lastAuthAsymId) {
+            if (!inMembrane[authAsymId]) inMembrane[authAsymId] = new Set<number>();
+            if (!outMembrane[authAsymId]) outMembrane[authAsymId] = new Set<number>();
+            lastAuthAsymId = authAsymId;
+        }
+
+        authSeqId = auth_seq_id(l);
+        Vec3.set(testPoint, x(l), y(l), z(l));
+        if (_isInMembranePlane(testPoint, normalVector!, dMin, dMax)) {
+            inMembrane[authAsymId].add(authSeqId);
+        } else {
+            outMembrane[authAsymId].add(authSeqId);
+        }
+    }
+
+    for (let k = 0, kl = offsets.length; k < kl; k++) {
+        setLocation(l, structure, offsets[k]);
+        authAsymId = auth_asym_id(l);
+        authSeqId = auth_seq_id(l);
+        if (inMembrane[authAsymId].has(authSeqId)) {
+            // chain change
+            if (authAsymId !== lastAuthAsymId) {
+                segments.push({ start: startOffset, end: endOffset });
+                lastAuthAsymId = authAsymId;
+                startOffset = k;
+                endOffset = k;
+            }
+
+            // sequence gaps
+            if (authSeqId !== lastAuthSeqId + 1) {
+                if (outMembrane[authAsymId].has(lastAuthSeqId + 1)) {
+                    segments.push({ start: startOffset, end: endOffset });
+                    startOffset = k;
+                }
+                lastAuthSeqId = authSeqId;
+                endOffset = k;
+            } else  {
+                lastAuthSeqId++;
+                endOffset++;
+            }
+        }
+    }
+    segments.push({ start: startOffset, end: endOffset });
+
+    let startAuth;
+    let endAuth;
+    const refinedSegments: Array<{ start: number, end: number }> = [];
+    for (let k = 0, kl = segments.length; k < kl; k++) {
+        const { start, end } = segments[k];
+        if (start === 0 || end === offsets.length - 1) continue;
+
+        // evaluate residues 1 pos outside of membrane
+        setLocation(l, structure, offsets[start - 1]);
+        Vec3.set(testPoint, x(l), y(l), z(l));
+        const d3 = -Vec3.dot(normalVector!, testPoint);
+
+        setLocation(l, structure, offsets[end + 1]);
+        Vec3.set(testPoint, x(l), y(l), z(l));
+        const d4 = -Vec3.dot(normalVector!, testPoint);
+
+        if (Math.min(d3, d4) < dMin && Math.max(d3, d4) > dMax) {
+            // reject this refinement
+            setLocation(l, structure, offsets[start]);
+            startAuth = auth_seq_id(l);
+            setLocation(l, structure, offsets[end]);
+            endAuth = auth_seq_id(l);
+            if (Math.abs(startAuth - endAuth) + 1 < adjust) {
+                return [];
+            }
+            refinedSegments.push(segments[k]);
+        }
+    }
+
+    return refinedSegments;
 }
 
 function qValue(currentStats: HphobHphil, initialStats: HphobHphil): number {
@@ -264,8 +435,12 @@ function qValue(currentStats: HphobHphil, initialStats: HphobHphil): number {
 export function isInMembranePlane(testPoint: Vec3, normalVector: Vec3, planePoint1: Vec3, planePoint2: Vec3): boolean {
     const d1 = -Vec3.dot(normalVector, planePoint1);
     const d2 = -Vec3.dot(normalVector, planePoint2);
+    return _isInMembranePlane(testPoint, normalVector, Math.min(d1, d2), Math.max(d1, d2));
+}
+
+function _isInMembranePlane(testPoint: Vec3, normalVector: Vec3, min: number, max: number): boolean {
     const d = -Vec3.dot(normalVector, testPoint);
-    return d > Math.min(d1, d2) && d < Math.max(d1, d2);
+    return d > min && d < max;
 }
 
 // generates a defined number of points on a sphere with radius = extent around the specified centroid
@@ -296,9 +471,9 @@ function findProximateAxes(ctx: ANVILContext, membrane: MembraneCandidate): Vec3
     const points = generateSpherePoints(ctx, 30000);
     let j = 4;
     let sphere_pts2: Vec3[] = [];
+    const s = 2 * extent / numberOfSpherePoints;
     while (sphere_pts2.length < numberOfSpherePoints) {
-        const d = 2 * extent / numberOfSpherePoints + j;
-        const dsq = d * d;
+        const dsq = (s + j) * (s + j);
         sphere_pts2 = [];
         for (let i = 0, il = points.length; i < il; i++) {
             if (Vec3.squaredDistance(points[i], membrane.spherePoint!) < dsq) {
@@ -326,8 +501,9 @@ namespace HphobHphil {
     }
 
     const testPoint = Vec3();
-    export function filtered(ctx: ANVILContext, label_comp_id: StructureElement.Property<string>, filter?: (test: Vec3) => boolean): HphobHphil {
+    export function filtered(ctx: ANVILContext, filter?: (test: Vec3) => boolean): HphobHphil {
         const { offsets, exposed, structure } = ctx;
+        const { label_comp_id } = StructureProperties.atom;
         const l = StructureElement.Location.create(structure);
         const { x, y, z } = StructureProperties.atom;
         let hphob = 0;
@@ -339,11 +515,13 @@ namespace HphobHphil {
             }
 
             setLocation(l, structure, offsets[k]);
-            Vec3.set(testPoint, x(l), y(l), z(l));
 
             // testPoints have to be in putative membrane layer
-            if (filter && !filter(testPoint)) {
-                continue;
+            if (filter) {
+                Vec3.set(testPoint, x(l), y(l), z(l));
+                if (!filter(testPoint)) {
+                    continue;
+                }
             }
 
             if (isHydrophobic(label_comp_id(l))) {
@@ -357,11 +535,35 @@ namespace HphobHphil {
 }
 
 // ANVIL-specific (not general) definition of membrane-favoring amino acids
-const HYDROPHOBIC_AMINO_ACIDS = new Set(['ALA', 'CYS', 'GLY', 'HIS', 'ILE', 'LEU', 'MET', 'PHE', 'SER', 'THR', 'VAL']);
+const HYDROPHOBIC_AMINO_ACIDS = new Set(['ALA', 'CYS', 'GLY', 'HIS', 'ILE', 'LEU', 'MET', 'PHE', 'SER', 'TRP', 'VAL']);
 export function isHydrophobic(label_comp_id: string): boolean {
     return HYDROPHOBIC_AMINO_ACIDS.has(label_comp_id);
 }
 
+/** 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). */
+export const MaxAsa: { [k: string]: number } = {
+    'ALA': 69.41,
+    'ARG': 201.25,
+    'ASN': 106.24,
+    'ASP': 102.69,
+    'CYS': 96.75,
+    'GLU': 134.74,
+    'GLN': 140.99,
+    'GLY': 32.33,
+    'HIS': 147.08,
+    'ILE': 137.96,
+    'LEU': 141.12,
+    'LYS': 163.30,
+    'MET': 156.64,
+    'PHE': 164.11,
+    'PRO': 119.90,
+    'SER': 78.11,
+    'THR': 101.70,
+    'TRP': 211.26,
+    'TYR': 177.38,
+    'VAL': 114.28
+};
+
 function setLocation(l: StructureElement.Location, structure: Structure, serialIndex: number) {
     l.structure = structure;
     l.unit = structure.units[structure.serialMapping.unitIndices[serialIndex]];

src/extensions/anvil/prop.ts

@@ -76,7 +76,8 @@ export const MembraneOrientationProvider: CustomStructureProperty.Provider<Membr
 });
 
 async function computeAnvil(ctx: CustomProperty.Context, data: Structure, props: Partial<ANVILProps>): Promise<MembraneOrientation> {
-    await AccessibleSurfaceAreaProvider.attach(ctx, data);
+    // can't get away with the default 92 points here
+    await AccessibleSurfaceAreaProvider.attach(ctx, data, { probeSize: 4.0, alphaOnly: true, numberOfSpherePoints: 184 });
     const p = { ...PD.getDefaultValues(ANVILParams), ...props };
     return await computeANVIL(data, p).runInContext(ctx.runtime);
 }