# -*- coding: utf-8 -*-
# **************************************************************************
# *
# * Authors: Estrella Fernandez Gimenez (me.fernandez@cnb.csic.es)
# *
# *
# * Unidad de Bioinformatica of Centro Nacional de Biotecnologia , CSIC
# *
# * This program is free software; you can redistribute it and/or modify
# * it under the terms of the GNU General Public License as published by
# * the Free Software Foundation; either version 2 of the License, or
# * (at your option) any later version.
# *
# * This program is distributed in the hope that it will be useful,
# * but WITHOUT ANY WARRANTY; without even the implied warranty of
# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# * GNU General Public License for more details.
# *
# * You should have received a copy of the GNU General Public License
# * along with this program; if not, write to the Free Software
# * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
# * 02111-1307 USA
# *
# * All comments concerning this program package may be sent to the
# * e-mail address 'scipion@cnb.csic.es'
# *
# **************************************************************************
from pyworkflow.protocol.params import PointerParam, BooleanParam, IntParam
from pwem.convert import headers
from pwem.objects import Volume, Transform
from pwem.protocols import EMProtocol
from pyworkflow import BETA, UPDATED, NEW, PROD
[docs]class XmippProtLocalVolAdj(EMProtocol):
"""Protocol to adjust locally volume intensity to a reference volume.
Occupancy volume is saved in protocol folder. Based on
https://www.sciencedirect.com/science/article/pii/S1047847723000874
AI Generated
## Overview
The Volume Local Adjustment protocol locally adjusts the intensity of one
volume so that it matches a reference volume.
This protocol is useful when two aligned maps have local differences in
intensity or occupancy. Instead of applying a single global scale factor to the
whole map, the protocol estimates local adjustments in neighborhoods around
each voxel. This allows the input volume to be adapted locally to the reference
volume.
The protocol can also perform a real-space subtraction after local adjustment.
In that case, the output represents the reference volume minus the locally
adjusted input volume.
The main output is an adjusted volume, or a subtraction volume if subtraction
is enabled.
## Inputs and General Workflow
The protocol requires:
- a reference volume;
- an input volume to be adjusted;
- a mask defining the region of interest;
- a neighborhood size.
The reference volume and input volume are passed to the Xmipp local-volume
adjustment program. The mask defines where the adjustment is performed, and the
neighborhood parameter defines the local region used to estimate the
adjustment.
The protocol writes the main result as `output_volume.mrc`. It also saves
additional files in the protocol folder, including the occupancy or local
adjustment information generated by the underlying program.
## Reference Volume
The **Reference volume** parameter defines the map used as the target for local
adjustment.
The input volume is adjusted to match this reference volume. The reference
volume also provides the sampling rate and, when available, the origin metadata
for the output volume.
The reference should be the map whose local intensity scale the user wants to
preserve.
## Input Volume
The **Input volume** parameter defines the map that will be locally adjusted to
the reference.
This volume should already be aligned with the reference volume and should
represent the same structure or comparable density. The protocol does not
perform map alignment.
If the input volume and reference volume are shifted, rotated, or sampled
differently, the local adjustment will be unreliable.
## Mask for Reference Volume
The **Mask for reference volume** parameter defines the region of interest for
the adjustment.
The mask should contain signal regions as white voxels, with value 1, and
background regions as black voxels, with value 0.
The mask focuses the local adjustment on the molecular region and prevents
background voxels from dominating the local estimates.
The mask must have the same pixel size as the input volumes. The protocol
validates this requirement.
## Neighborhood
The **Neighborhood** parameter defines the side length, in angstroms, of the
local region used for adjustment.
For each local area, the protocol estimates how the input volume should be
scaled or adjusted to match the reference volume. A smaller neighborhood makes
the adjustment more local and potentially more sensitive to fine variations. A
larger neighborhood makes the adjustment smoother and more stable but less
spatially specific.
The protocol summary reports the neighborhood both in angstroms and in pixels,
using the sampling rate of the reference volume.
The default value is 5 Å.
## Perform Subtraction
The **Perform subtraction?** option controls whether the protocol only adjusts
the input volume or also subtracts it from the reference.
If this option is disabled, the output is the locally adjusted version of the
input volume.
If this option is enabled, the protocol performs:
\[
\text{output} = \text{reference volume} - \text{locally adjusted input volume}
\]
This is useful when the user wants to highlight residual density after locally
matching one map to another.
## Output Volume
The main output is **outputVolume**.
This volume is written as:
`output_volume.mrc`
It is assigned the sampling rate of the reference volume. If the reference
volume is an MRC file with origin information, the protocol copies that origin
to the output.
The output should be interpreted according to the subtraction setting:
- without subtraction, it is the locally adjusted input volume;
- with subtraction, it is the residual difference between the reference volume
and the locally adjusted input volume.
## Additional Saved Files
The underlying Xmipp program is run with a save option pointing to the protocol
folder.
This means that additional diagnostic or support files may be written besides
the main output volume. These files can include local adjustment or occupancy
information generated during the procedure.
These files are useful for advanced inspection of how the local adjustment was
performed.
## Validation Rules
The protocol checks that the reference volume and input volume have the same
pixel size.
It also checks that the mask and the volumes have the same pixel size.
If these conditions are not met, the protocol reports validation errors.
Although the validation focuses on pixel size, the user should also ensure that
the volumes and mask have compatible box sizes, origins, and coordinate frames.
## Interpreting the Result
The output volume should be interpreted as a locally intensity-adjusted map, or
as a locally adjusted subtraction result.
Local adjustment can compensate for regional differences in intensity or
occupancy between two maps. This can make comparisons more meaningful when a
global scaling factor is not sufficient.
However, the result depends on the mask, neighborhood size, alignment, and map
quality. If the maps are not well aligned, or if the neighborhood is too small,
the adjustment may fit noise or create local artifacts. If the neighborhood is
too large, real local differences may be smoothed out.
## Practical Recommendations
Use this protocol only with volumes that are already aligned.
Make sure the reference volume, input volume, and mask have the same sampling
rate and compatible coordinate frames.
Use a mask that includes the molecular region of interest and avoids excessive
background.
Start with the default neighborhood size. Increase it if the adjustment appears
too noisy; decrease it if broad smoothing hides local differences.
Use subtraction when the goal is to identify residual density after matching
the input volume locally to the reference.
Inspect both the output volume and any additional saved diagnostic files.
Do not interpret residual density automatically as biological signal. Check
whether it could arise from misalignment, masking, local overfitting, or map
quality differences.
## Final Perspective
Volume Local Adjustment is a local map-matching protocol.
For biological users, its main value is that it adjusts one volume to a
reference in a spatially adaptive way. This is useful when local intensity or
occupancy differences make global adjustment insufficient.
The optional subtraction mode turns the protocol into a residual-map generator,
helping highlight differences between two locally matched maps. The result is
most reliable when the input maps are well aligned, consistently sampled, and
analyzed with an appropriate mask and neighborhood size.
"""
_label = 'volume local adjustment'
_possibleOutputs = Volume
_devStatus = NEW
# --------------------------- DEFINE param functions --------------------------------------------
def _defineParams(self, form):
form.addSection(label='Input')
form.addParam('vol1', PointerParam, pointerClass='Volume', label="Reference volume",
help='Specify a volume to be used as reference volume.')
form.addParam('vol2', PointerParam, pointerClass='Volume', label="Input volume",
help='Specify a volume which will be adjusted to the reference volume.')
form.addParam('mask', PointerParam, pointerClass='VolumeMask', label="Mask for reference volume",
help='Specify a mask to define region of interest (which is signal in white (1s) and background in '
'black (0s))')
form.addParam('neighborhood', IntParam, label="Neighborhood (A)", default=5,
help='Side length (in Angstroms) of a square which will define the region of adjustment')
form.addParam('subtract', BooleanParam, label="Perform subtraction?", default=False,
help='Perform subtraction of reference volume minus input volume in real space')
# --------------------------- INSERT steps functions --------------------------------------------
def _insertAllSteps(self):
self._insertFunctionStep('adjustStep')
self._insertFunctionStep('createOutputStep')
# --------------------------- STEPS functions --------------------------------------------
[docs] def adjustStep(self):
vol1 = self.vol1.get().clone()
fnVol1 = vol1.getFileName()
vol2 = self.vol2.get().getFileName()
if fnVol1.endswith('.mrc'):
fnVol1 += ':mrc'
if vol2.endswith('.mrc'):
vol2 += ':mrc'
program = "xmipp_local_volume_adjust"
args = '--i1 %s --i2 %s -o %s --mask %s --neighborhood %d --sampling %s --save %s' % \
(fnVol1, vol2, self._getExtraPath("output_volume.mrc"), self.mask.get().getFileName(),
self.neighborhood.get(), vol1.getSamplingRate(), self._getExtraPath())
if self.subtract.get():
args += ' --sub'
self.runJob(program, args)
[docs] def createOutputStep(self):
vol1 = self.vol1.get()
volume = Volume()
volume.setSamplingRate(vol1.getSamplingRate())
if vol1.getFileName().endswith('mrc'):
origin = Transform()
ccp4header = headers.Ccp4Header(vol1.getFileName(), readHeader=True)
shifts = ccp4header.getOrigin()
origin.setShiftsTuple(shifts)
volume.setOrigin(origin)
volume.setFileName(self._getExtraPath("output_volume.mrc"))
filename = volume.getFileName()
if filename.endswith('.mrc') or filename.endswith('.map'):
volume.setFileName(filename + ':mrc')
self._defineOutputs(outputVolume=volume)
# --------------------------- INFO functions --------------------------------------------
def _summary(self):
neighborhood = self.neighborhood.get()
vol1 = self.vol1.get()
summary = ["Volume 1: %s\nVolume 2: %s\nInput mask 1: %s\n\nNeighborhood: %d Å (%d px)" %
(vol1.getFileName(), self.vol2.get().getFileName(), self.mask.get().getFileName(),
neighborhood, round(neighborhood/vol1.getSamplingRate()))]
if self.subtract.get():
summary.append("\nSubtraction performed")
return summary
def _methods(self):
methods = []
if not hasattr(self, 'outputVolume'):
methods.append("Output volume not ready yet.")
else:
methods.append("Volume %s adjusted to volume %s" % (self.vol2.get().getFileName(),
self.vol1.get().getFileName()))
return methods
def _validate(self):
errors = []
if self.vol1.get().getSamplingRate() != self.vol2.get().getSamplingRate():
errors.append('Input volumes should have same pixel size')
if self.vol1.get().getSamplingRate() != self.mask.get().getSamplingRate():
errors.append('\nInput mask and volumes should have same pixel size')
return errors