Source code for xmipp3.protocols.protocol_multireference_alignability

# **************************************************************************
# *
# * Authors:         Javier Vargas (jvargas@cnb.csic.es) (2016)
# *
# * 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 os.path import join, isfile
from shutil import copyfile
import os

from pyworkflow.object import Float, String
from pyworkflow.protocol.params import (PointerParam, FloatParam,
                                        STEPS_PARALLEL,
                                        StringParam, BooleanParam, IntParam,
                                        LEVEL_ADVANCED, USE_GPU, GPU_LIST)

from pyworkflow.utils.path import moveFile, makePath, cleanPattern
from pyworkflow.gui.plotter import Plotter

from pwem.objects import Volume
import pwem.emlib.metadata as md
from pwem.protocols import ProtAnalysis3D

from xmipp3.constants import CUDA_ALIGN_SIGNIFICANT
from xmipp3.convert import writeSetOfParticles, writeSetOfVolumes, \
    getImageLocation
from xmipp3.base import isXmippCudaPresent

[docs]class XmippProtMultiRefAlignability(ProtAnalysis3D): """ Performs soft alignment validation of a set of particles confronting them against a given 3DEM map. This protocol produces particle alignment precision and accuracy parameters. """ _label = 'multireference alignability' def __init__(self, *args, **kwargs): ProtAnalysis3D.__init__(self, *args, **kwargs) # --------------------------- DEFINE param functions -------------------------------------------- def _defineParams(self, form): form.addHidden(USE_GPU, BooleanParam, default=True, label="Use GPU for execution", help="This protocol has both CPU and GPU implementation.\ Select the one you want to use.") form.addHidden(GPU_LIST, StringParam, default='0', expertLevel=LEVEL_ADVANCED, label="Choose GPU IDs", help="Add a list of GPU devices that can be used") form.addSection(label='Input') form.addParam('inputVolumes', PointerParam, pointerClass='Volume', label="Input volume", help='Select the input volume(s).') form.addParam('inputParticles', PointerParam, pointerClass='SetOfParticles', pointerCondition='hasAlignment', label="Input particles", important=True, help='Select the input projection images.') form.addParam('symmetryGroup', StringParam, default='c1', label="Symmetry group", help='See [[Xmipp Symmetry][http://www2.mrc-lmb.cam.ac.uk/Xmipp/index.php/Conventions_%26_File_formats#Symmetry]] page ' 'for a description of the symmetry format accepted by Xmipp') form.addParam('angularSampling', FloatParam, default=5, expertLevel=LEVEL_ADVANCED, label="Angular Sampling (degrees)", help='Angular distance (in degrees) between neighboring projection points ') form.addParam('numOrientations', FloatParam, default=7, expertLevel=LEVEL_ADVANCED, label="Number of Orientations for particle", help='Parameter to define the number of most similar volume \n' ' projected images for each projection image') form.addParam('doNotUseWeights', BooleanParam, default=False, expertLevel=LEVEL_ADVANCED, label="Do not use the weights", help='Do not use the weights in the clustering calculation') form.addParam('pseudoSymmetryGroup', StringParam, default='', expertLevel=LEVEL_ADVANCED, label="Pseudo symmetry group", help='Add only in case the map is close to a symmetry different and more restrict than the one reported in the parameter Symmetry group.' 'See [[Xmipp Symmetry][http://www2.mrc-lmb.cam.ac.uk/Xmipp/index.php/Conventions_%26_File_formats#Symmetry]] page ' 'for a description of the symmetry format accepted by Xmipp') form.addParam('minTilt', FloatParam, default=0, expertLevel=LEVEL_ADVANCED, label="Minimum allowed tilt angle", help='Tilts below this value will not be considered for the alignment') form.addParam('maxTilt', FloatParam, default=180, expertLevel=LEVEL_ADVANCED, label="Maximum allowed tilt angle without mirror check", help='Tilts above this value will not be considered for the alignment without mirror check') form.addSection(label='Preprocess') form.addParam('doWiener', BooleanParam, default='True', label="CTF correction", help='Perform CTF correction by Wiener filtering.') form.addParam('isIsotropic', BooleanParam, default='True', label="Isotropic Correction", condition='doWiener', help='If true, Consider that there is not astigmatism and then it is performed an isotropic correction.') form.addParam('padding_factor', IntParam, default=2, expertLevel=LEVEL_ADVANCED, label="Padding factor", condition='doWiener', help='Padding factor for Wiener correction ') form.addParam('wiener_constant', FloatParam, default=-1, expertLevel=LEVEL_ADVANCED, label="Wiener constant", condition='doWiener', help=' Wiener-filter constant (if < 0: use FREALIGN default)') form.addParam('correctEnvelope', BooleanParam, default='False', expertLevel=LEVEL_ADVANCED, label="Correct for CTF envelope", condition='doWiener', help=' Only in cases where the envelope is well estimated correct for it') form.addParam('targetResolution', FloatParam, default=8, label='Target resolution (A)', help='Low pass filter the particles to this resolution. This usually helps a lot obtaining good alignment. You should have a good' ' reason to modify this value outside the range [8-10] A') form.addParallelSection(threads=1, mpi=1) # --------------------------- INSERT steps functions -------------------------------------------- def _insertAllSteps(self): convertId = self._insertFunctionStep('convertInputStep', self.inputParticles.get().getObjId()) deps = [] # store volumes steps id to use as dependencies for last step commonParams = self._getCommonParams() commonParamsRef = self._getCommonParamsRef() sym = self.symmetryGroup.get() for i, vol in enumerate(self._iterInputVols()): volName = getImageLocation(vol) volDir = self._getVolDir(i + 1) pmStepId = self._insertFunctionStep('projectionLibraryStep', volName, volDir, self.angularSampling.get(), prerequisites=[convertId]) sigStepId1 = self._insertFunctionStep('significantStep', volName, volDir, 'exp_particles.xmd', commonParams, prerequisites=[pmStepId]) phanProjStepId = self._insertFunctionStep('phantomProject', prerequisites=[ sigStepId1]) sigStepId2 = self._insertFunctionStep('significantStep', volName, volDir, 'ref_particles.xmd', commonParamsRef, prerequisites=[ phanProjStepId]) if (not (self.pseudoSymmetryGroup.get() == '')): sym = self.pseudoSymmetryGroup.get() volStepId = self._insertFunctionStep('alignabilityStep', volName, volDir, sym, prerequisites=[sigStepId2]) deps.append(volStepId) self._insertFunctionStep('createOutputStep', prerequisites=deps)
[docs] def convertInputStep(self, particlesId): """ Write the input images as a Xmipp metadata file. particlesId: is only need to detect changes in input particles and cause restart from here. """ writeSetOfParticles(self.inputParticles.get(), self._getPath('input_particles.xmd')) if self.doWiener.get(): params = ' -i %s' % self._getPath('input_particles.xmd') params += ' -o %s' % self._getExtraPath( 'corrected_ctf_particles.stk') params += ' --save_metadata_stack %s' % self._getExtraPath( 'corrected_ctf_particles.xmd') params += ' --pad %s' % self.padding_factor.get() params += ' --wc %s' % self.wiener_constant.get() params += ' --sampling_rate %s' % self.inputParticles.get().getSamplingRate() if self.inputParticles.get().isPhaseFlipped(): params += ' --phase_flipped ' if self.correctEnvelope: params += ' --correct_envelope ' nproc = self.numberOfMpi.get() nT = self.numberOfThreads.get() self.runJob('xmipp_ctf_correct_wiener2d', params) newTs, newXdim = self._getModifiedSizeAndSampling() if self.doWiener.get(): params = ' -i %s' % self._getExtraPath('corrected_ctf_particles.xmd') else : params = ' -i %s' % self._getPath('input_particles.xmd') params += ' -o %s' % self._getExtraPath('scaled_particles.stk') params += ' --save_metadata_stack %s' % self._getExtraPath('scaled_particles.xmd') params += ' --fourier %d' % newXdim self.runJob('xmipp_image_resize',params) from pwem.emlib.image import ImageHandler img = ImageHandler() img.convert(self.inputVolumes.get(), self._getExtraPath("volume.vol")) Xdim = self.inputVolumes.get().getDim()[0] if Xdim != newXdim: self.runJob("xmipp_image_resize", "-i %s --dim %d" % \ (self._getExtraPath("volume.vol"), newXdim), numberOfMpi=1)
def _getCommonParams(self): params = ' -i %s' % self._getExtraPath('scaled_particles.xmd') params += ' --sym %s' % self.symmetryGroup.get() params += ' --dontReconstruct' params += ' --useForValidation %0.3f' % (self.numOrientations.get() - 1) params += ' --dontCheckMirrors' return params def _getCommonParamsRef(self): params = ' -i %s' % self._getPath('reference_particles.xmd') params += ' --sym %s' % self.symmetryGroup.get() params += ' --dontReconstruct' params += ' --useForValidation %0.3f' % (self.numOrientations.get() - 1) params += ' --dontCheckMirrors' return params def _getModifiedSizeAndSampling(self): Xdim = self.inputParticles.get().getDimensions()[0] Ts = self.inputParticles.get().getSamplingRate() newTs = self.targetResolution.get() * 0.4 newTs = max(Ts, newTs) newXdim = Xdim * Ts / newTs return newTs, newXdim
[docs] def phantomProject(self): nproc = self.numberOfMpi.get() nT = self.numberOfThreads.get() newTs, newXdim = self._getModifiedSizeAndSampling() pathParticles = self._getExtraPath('scaled_particles.xmd') R = -int(newXdim / 2) f = open(self._getExtraPath('params.txt'), 'w') f.write("""# XMIPP_STAR_1 * # data_block1 _dimensions2D '%d %d' _projAngleFile %s _ctfPhaseFlipped %d _ctfCorrected %d _applyShift 0 _noisePixelLevel '0 0'""" % (newXdim, newXdim, pathParticles, self.inputParticles.get().isPhaseFlipped(), self.doWiener.get())) f.close() param = ' -i %s' % self._getExtraPath("volume.vol") param += ' --params %s' % self._getExtraPath('params.txt') param += ' -o %s' % self._getPath('reference_particles.xmd') param += ' --sampling_rate % 0.3f' % newTs param += ' --method fourier' #while (~isfile(self._getExtraPath('params'))): # print('No created') self.runJob('xmipp_phantom_project', param, numberOfMpi=1,numberOfThreads=1) param = ' -i %s' % self._getPath('reference_particles.stk') param += ' --mask circular %d' % R self.runJob('xmipp_transform_mask', param, numberOfMpi=nproc, numberOfThreads=nT)
[docs] def projectionLibraryStep(self, volName, volDir, angularSampling): # Generate projections from this reconstruction nproc = self.numberOfMpi.get() nT = self.numberOfThreads.get() volName = self._getExtraPath("volume.vol") makePath(volDir) fnGallery = (volDir + '/gallery.stk') params = '-i %s -o %s --sampling_rate %f --sym %s --method fourier 1 0.25 bspline --compute_neighbors --angular_distance %f --experimental_images %s --max_tilt_angle %f --min_tilt_angle %f' \ % ( volName, fnGallery, angularSampling, self.symmetryGroup.get(), -1, self._getExtraPath('scaled_particles.xmd'), self.maxTilt.get(), self.minTilt.get()) self.runJob("xmipp_angular_project_library", params, numberOfMpi=nproc, numberOfThreads=nT)
[docs] def significantStep(self, volName, volDir, anglesPath, params): nproc = self.numberOfMpi.get() nT = self.numberOfThreads.get() fnGallery = (volDir + '/gallery.doc') if not self.useGpu.get(): params += ' --initgallery %s' % fnGallery params += ' --odir %s' % volDir params += ' --iter %d' % 1 self.runJob('xmipp_reconstruct_significant', params, numberOfMpi=nproc, numberOfThreads=nT) copyfile(volDir + '/angles_iter001_00.xmd', self._getTmpPath(anglesPath)) else: count=0 GpuListCuda='' if self.useQueueForSteps() or self.useQueue(): GpuList = os.environ["CUDA_VISIBLE_DEVICES"] GpuList = GpuList.split(",") for elem in GpuList: GpuListCuda = GpuListCuda+str(count)+' ' count+=1 else: GpuList = ' '.join([str(elem) for elem in self.getGpuList()]) GpuListAux = '' for elem in self.getGpuList(): GpuListCuda = GpuListCuda+str(count)+' ' GpuListAux = GpuListAux+str(elem)+',' count+=1 os.environ["CUDA_VISIBLE_DEVICES"] = GpuListAux if anglesPath == 'exp_particles.xmd': params = ' -i %s' % self._getExtraPath('scaled_particles.xmd') elif anglesPath == 'ref_particles.xmd': params = ' -i %s' % self._getPath('reference_particles.xmd') params += ' --keepBestN %f' % (self.numOrientations.get() - 1) params += ' -r %s' % fnGallery params += ' -o %s' % self._getTmpPath(anglesPath) params += ' --dev %s ' % GpuListCuda self.runJob(CUDA_ALIGN_SIGNIFICANT, params, numberOfMpi=1)
[docs] def alignabilityStep(self, volName, volDir, sym): nproc = self.numberOfMpi.get() nT = self.numberOfThreads.get() makePath(volDir) inputFile = self._getPath('input_particles.xmd') inputFileRef = self._getPath('reference_particles.xmd') aFile = self._getTmpPath('exp_particles.xmd') aFileRef = self._getTmpPath('ref_particles.xmd') aFileGallery = (volDir + '/gallery.doc') params = ' -i %s' % inputFile params += ' -i2 %s' % inputFileRef params += ' --volume %s' % volName params += ' --angles_file %s' % aFile params += ' --angles_file_ref %s' % aFileRef params += ' --gallery %s' % aFileGallery params += ' --odir %s' % volDir params += ' --sym %s' % sym params += ' --check_mirrors' if self.doNotUseWeights: params += ' --dontUseWeights' self.runJob('xmipp_multireference_aligneability', params, numberOfMpi=nproc, numberOfThreads=nT)
[docs] def neighbourhoodDirectionStep(self, volName, volDir, sym): aFileGallery = (volDir + '/gallery.doc') neighbours = (volDir + '/neighbours.xmd') params = ' --i1 %s' % self._getPath('input_particles.xmd') params += ' --i2 %s' % aFileGallery params += ' -o %s' % neighbours params += ' --dist %s' % (self.angDist.get() + 1) params += ' --sym %s' % sym self.runJob('xmipp_angular_neighbourhood', params, numberOfMpi=1, numberOfThreads=1)
[docs] def angularAccuracyStep(self, volName, volDir, indx): nproc = self.numberOfMpi.get() nT = self.numberOfThreads.get() neighbours = (volDir + '/neighbours.xmd') params = ' -i %s' % volName params += ' --i2 %s' % neighbours params += ' -o %s' % ( volDir + '/pruned_particles_alignability_accuracy.xmd') self.runJob('xmipp_angular_accuracy_pca', params, numberOfMpi=nproc, numberOfThreads=nT)
[docs] def createOutputStep(self): outputVols = self._createSetOfVolumes() for i, vol in enumerate(self._iterInputVols()): volDir = self._getVolDir(i + 1) volume = vol.clone() volPrefix = 'vol%03d_' % (i + 1) m_pruned = md.MetaData() m_pruned.read(volDir + '/pruned_particles_alignability.xmd') prunedMd = self._getExtraPath( volPrefix + 'pruned_particles_alignability.xmd') moveFile(join(volDir, 'pruned_particles_alignability.xmd'), prunedMd) m_volScore = md.MetaData() m_volScore.read(volDir + '/validationAlignability.xmd') validationMd = self._getExtraPath( volPrefix + 'validation_alignability.xmd') moveFile(join(volDir, 'validationAlignability.xmd'), validationMd) imgSet = self.inputParticles.get() outImgSet = self._createSetOfParticles(volPrefix) outImgSet.copyInfo(imgSet) outImgSet.copyItems(imgSet, updateItemCallback=self._setWeight, itemDataIterator=md.iterRows(prunedMd, sortByLabel=md.MDL_ITEM_ID)) mdValidatoin = md.getFirstRow(validationMd) weight = mdValidatoin.getValue(md.MDL_WEIGHT_PRECISION_ALIGNABILITY) volume.weightAlignabilityPrecision = Float(weight) weight = mdValidatoin.getValue(md.MDL_WEIGHT_ACCURACY_ALIGNABILITY) volume.weightAlignabilityAccuracy = Float(weight) weight = mdValidatoin.getValue(md.MDL_WEIGHT_PRECISION_MIRROR) volume.weightMirror = Float(weight) volume.cleanObjId() # clean objects id to assign new ones inside the set outputVols.append(volume) self._defineOutputs(outputParticles=outImgSet) self.createPlot2D(volPrefix, m_pruned) outputVols.setSamplingRate(volume.getSamplingRate()) self._defineOutputs(outputVolumes=outputVols) cleanPattern(self._getPath("reference_particles.*")) cleanPattern(self._getExtraPath("scaled_particles.*")) cleanPattern(self._getExtraPath("reference_particles.*")) cleanPattern(self._getExtraPath("corrected_ctf_particles.*")) cleanPattern(self._getExtraPath("volume.vol")) cleanPattern(self._getExtraPath("params.txt"))
# --------------------------- INFO functions -------------------------------------------- def _validate(self): validateMsgs = [] # if there are Volume references, it cannot be empty. if self.inputVolumes.get() and not self.inputVolumes.hasValue(): validateMsgs.append('Please provide an input reference volume.') if self.inputParticles.get() and not self.inputParticles.hasValue(): validateMsgs.append('Please provide input particles.') if self.useGpu and not isXmippCudaPresent(CUDA_ALIGN_SIGNIFICANT): validateMsgs.append("You have asked to use GPU, but I cannot find the Xmipp GPU programs in the path") return validateMsgs def _summary(self): summary = [ "Input particles: %s" % self.inputParticles.get().getNameId()] summary.append("-----------------") if self.inputVolumes.get(): for i, vol in enumerate(self._iterInputVols()): summary.append("Input volume(s)_%d: [%s]" % (i + 1, vol)) summary.append("-----------------") if (not hasattr(self, 'outputVolumes')): summary.append("Output volumes not ready yet.") else: for i, vol in enumerate(self._iterInputVols()): VolPrefix = 'vol%03d_' % (i + 1) mdVal = md.MetaData(self._getExtraPath( VolPrefix + 'validation_alignability.xmd')) weightAccuracy = mdVal.getValue( md.MDL_WEIGHT_ACCURACY_ALIGNABILITY, mdVal.firstObject()) weightPrecision = mdVal.getValue( md.MDL_WEIGHT_PRECISION_ALIGNABILITY, mdVal.firstObject()) weightAlignability = mdVal.getValue(md.MDL_WEIGHT_ALIGNABILITY, mdVal.firstObject()) summary.append("ALIGNABILITY ACCURACY parameter_%d : %f" % ( i + 1, weightAccuracy)) summary.append("ALIGNABILITY PRECISION parameter_%d : %f" % ( i + 1, weightPrecision)) summary.append( "ALIGNABILITY ACCURACY & PRECISION parameter_%d : %f" % ( i + 1, weightAlignability)) summary.append("-----------------") return summary def _methods(self): messages = [] if (hasattr(self, 'outputVolumes')): messages.append('The quality parameter(s) has been obtained using ' 'the approach [Vargas2014a] with angular sampling ' 'of %f and number of orientations of %f' % ( self.angularSampling.get(), self.numOrientations.get())) return messages def _citations(self): return ['Vargas2014a'] # --------------------------- UTILS functions -------------------------------------------- def _getVolDir(self, volIndex): return self._getTmpPath('vol%03d' % volIndex) def _iterInputVols(self): """ In this function we will encapsulate the logic to iterate through the input volumes. This give the flexibility of having Volumes, SetOfVolumes or a combination of them as input and the protocol code remain the same. """ inputVols = self.inputVolumes.get() if isinstance(inputVols, Volume): yield inputVols else: for vol in inputVols: yield vol def _defineMetadataRootName(self, mdrootname, volId): if mdrootname == 'P': VolPrefix = 'vol%03d_' % (volId) return self._getExtraPath(VolPrefix + 'clusteringTendency.xmd') if mdrootname == 'Volume': VolPrefix = 'vol%03d_' % (volId) return self._getExtraPath(VolPrefix + 'validation.xmd') def _definePName(self): fscFn = self._defineMetadataRootName('P') return fscFn def _defineVolumeName(self, volId): fscFn = self._defineMetadataRootName('Volume', volId) return fscFn def _setWeight(self, item, row): item._xmipp_scoreAlignabilityPrecision = Float( row.getValue(md.MDL_SCORE_BY_ALIGNABILITY_PRECISION)) item._xmipp_scoreAlignabilityAccuracy = Float( row.getValue(md.MDL_SCORE_BY_ALIGNABILITY_ACCURACY)) item._xmipp_scoreMirror = Float(row.getValue(md.MDL_SCORE_BY_MIRROR)) item._xmipp_weight = Float( float(item._xmipp_scoreAlignabilityAccuracy)*float(item._xmipp_scoreAlignabilityPrecision))
[docs] def createPlot2D(self,volPrefix,md): from pwem import emlib figurePath = self._getExtraPath(volPrefix + 'softAlignmentValidation2D.png') figureSize = (8, 6) # alignedMovie = mic.alignMetaData plotter = Plotter(*figureSize) figure = plotter.getFigure() ax = figure.add_subplot(111) ax.grid() ax.set_title('Soft alignment validation map') ax.set_xlabel('Angular Precision') ax.set_ylabel('Angular Accuracy') for objId in md: x = md.getValue(emlib.MDL_SCORE_BY_ALIGNABILITY_PRECISION, objId) y = md.getValue(emlib.MDL_SCORE_BY_ALIGNABILITY_ACCURACY, objId) ax.plot(x, y, 'r.',markersize=1) ax.grid(True, which='both') ax.autoscale_view(True, True, True) plotter.savefig(figurePath) plotter.show() return plotter