# **************************************************************************
# *
# * 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'
INPUTARG = "-i %s"
OUTPUTARG = " -o %s"
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)
def _getFileName(self, key, **kwargs):
if key=="volume":
return self._getExtraPath("volume.mrc")
elif key=="reference_particles":
return self._getPath("reference_particles.xmd")
else:
return ""
# --------------------------- 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',
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)
def _getCommonParams(self):
params = self.INPUTARG % 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 = self.INPUTARG % self._getFileName("reference_particles")
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 = self.INPUTARG % self._getFileName("volume")
param += ' --params %s' % self._getExtraPath('params.txt')
param += self.OUTPUTARG % self._getFileName("reference_particles")
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 = self.INPUTARG % self._getPath('reference_particles.stk')
param += ' --mask circular %d' % R
self.runJob('xmipp_transform_mask', param, numberOfMpi=nproc,
numberOfThreads=nT)
[docs] def projectionLibraryStep(self, volDir, angularSampling):
# Generate projections from this reconstruction
nproc = self.numberOfMpi.get()
nT = self.numberOfThreads.get()
volName = self._getFileName("volume")
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 = self.INPUTARG % self._getExtraPath('scaled_particles.xmd')
elif anglesPath == 'ref_particles.xmd':
params = self.INPUTARG % self._getFileName("reference_particles")
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._getFileName("reference_particles")
aFile = self._getTmpPath('exp_particles.xmd')
aFileRef = self._getTmpPath('ref_particles.xmd')
aFileGallery = (volDir + '/gallery.doc')
params = self.INPUTARG % 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 += self.OUTPUTARG % 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 = self.INPUTARG % volName
params += ' --i2 %s' % neighbours
params += self.OUTPUTARG % (
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._getFileName("volume"))
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 _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