diff --git a/CHANGES.txt b/CHANGES.txt index 16f3b24..a92e81b 100644 --- a/CHANGES.txt +++ b/CHANGES.txt @@ -1,3 +1,10 @@ +Dec 22, 2023 - **v4.0.10** + Add new protocols: + * **3D Variability Analysis**: Protocol to compute the principle modes of variability with a dataset of aligned particles + * **3D Variability Display**: Protocol to create various versions of a 3D variability result that can be used for display + * **Blob Picker**: Automatically picks particles by searching for Gaussian signals. + * **Patch CTF Estimation**: Patch-based CTF estimation automatically estimates defocus variation for tilted, bent, deformed samples and is accurate for all particle sizes and types including flexible and membrane proteins. + Nov 17, 2023 - **v4.0.9** Compatibility with cryoSPARC v4.4.0 Handling aborted protocols/jobs diff --git a/README.rst b/README.rst index bfbb1fd..493ba2c 100644 --- a/README.rst +++ b/README.rst @@ -24,26 +24,30 @@ You will need to use `3.0.0 3, 0>4, 3>2, 2>1' is a valid connection" - " string. Each pair X>Y denotes that segment Y is " - "joined to segment X. The connections must form " - "a 'tree' structure and cannot have cycles. " - "The first pair X>Y must start with the " - "root of the tree as X. The connections " - "must be in breadth-first order of the tree. " - "When using Chimera Segger segmentation input, " - "the X>Y numbers should be region_ids from Segger (e.g., 948>947)") - - meshPreparationGroup.addParam('tetra_rigid_list', StringParam, - allowsNull=True, - default=None, - label="Rigid segments", - help="A comma separated list of segments to make " - "rigid. This is done by setting the rigidity " - "weight of tetra elements for this region to 20. ") - - - rigidityWeighting = form.addGroup('Rigidity weighting') - rigidityWeighting.addParam('rigidity_penalty_min', FloatParam, - default=0.5, - label="Min. rigidity weight", - help="Rigidity weights of tetra elements" - " are 1.0 in the most dense regions " - "of the input consensus map, and fall " - "off to this value (default 0.5) in the" - " least dense/empty regions. This helps " - "encourage the deformation model to expand/contract" - " empty space without distorting the protein density.") - - rigidityWeighting.addParam('rigidity_penalty_stiffen_low_density', BooleanParam, - default=False, - label="Stiffen low density regions", - help="Turning this on will cause the rigidity" - " weights of tetra elements at the periphery " - "of the input consensus density to be increased to 3.0 . " - "Empty regions will still have low rigidity, but non-empty " - "regions at the periphery of the structure will be rigidified." - " This helps to combat overfitting in smaller particles or" - " poor SNR data where otherwise low density peripheral " - "features start to 'fly around'. However, it can also cause " - "the deformations to be overly smooth and blur motion 'boundaries between domains. ") - - # --------------[Compute settings]--------------------------- form.addSection(label="Compute settings") addComputeSectionParams(form, allowMultipleGPUs=False, needGPU=False) @@ -234,11 +113,9 @@ def _defineParams(self, form): def _insertAllSteps(self): self._defineFileNames() self._defineParamsPrepareName() - self._defineParamsMeshName() self._initializeCryosparcProject() self._insertFunctionStep(self.convertInputStep) self._insertFunctionStep(self.dataPrepareStep) - self._insertFunctionStep(self.mergePrepareStep) self._insertFunctionStep(self.createOutputStep) # --------------------------- STEPS functions ------------------------------ @@ -246,10 +123,6 @@ def dataPrepareStep(self): self.info(pwutils.yellowStr("3D Flex Data Preparation started...")) self.doRun3DFlexDataPrepare() - def mergePrepareStep(self): - self.info(pwutils.yellowStr("3D Flex Mesh Preparation started...")) - self.doRun3DFlexMeshPrepare() - def createOutputStep(self): """ Create the protocol output. Convert cryosparc file to Relion file @@ -319,17 +192,7 @@ def _validate(self): def _defineParamsPrepareName(self): """ Define a list with 3D Flex Prepare Data parameters names""" self._paramsPrepareName = ['box_size_pix', 'bin_size_pix', 'alpha_min', - 'keep_num_particles'] - self.lane = str(self.getAttributeValue('compute_lane')) - - def _defineParamsMeshName(self): - """ Define a list with 3D Flex Mesh Prepare parameters names""" - self._maskMeshPrepareName = ['mask_in_lowpass_A', 'mask_in_threshold_level', - 'mask_dilate_A', 'mask_pad_A'] - self._paramsMeshName = ['tetra_num_cells', 'tetra_segments_path', - 'tetra_segments_fuse_list', 'tetra_rigid_list', - 'rigidity_penalty_min', - 'rigidity_penalty_stiffen_low_density'] + 'keep_num_particles'] self.lane = str(self.getAttributeValue('compute_lane')) def doRun3DFlexDataPrepare(self): @@ -358,45 +221,6 @@ def doRun3DFlexDataPrepare(self): "details.") clearIntermediateResults(self.projectName.get(), self.run3DFlexDataPrepJob.get()) - def doRun3DFlexMeshPrepare(self): - self._className = 'flex_meshprep' - params = {} - varDataPrepJob = str(self.run3DFlexDataPrepJob.get()) - input_group_connect = {"volume": str('%s.volume' % varDataPrepJob)} - - if self.refMask.get() is not None: - input_group_connect["mask"] = str(self.mask) - else: - for paramName in self._maskMeshPrepareName: - if self.getAttributeValue(paramName) is not None: - params[str(paramName)] = str(self.getAttributeValue(paramName)) - - for paramName in self._paramsMeshName: - if self.getAttributeValue(paramName) is not None: - params[str(paramName)] = str(self.getAttributeValue(paramName)) - - run3DFlexMeshPrepJob = enqueueJob(self._className, - self.projectName.get(), - self.workSpaceName.get(), - str(params).replace('\'', '"'), - str(input_group_connect).replace('\'', - '"'), - self.lane, False) - - self.run3DFlexMeshPrepJob = String(run3DFlexMeshPrepJob.get()) - self.currenJob.set(self.run3DFlexMeshPrepJob.get()) - self._store(self) - - waitForCryosparc(self.projectName.get(), - self.run3DFlexMeshPrepJob.get(), - "An error occurred in the 3D Flex Mesh Preparation process. " - "Please, go to cryoSPARC software for more " - "details.") - clearIntermediateResults(self.projectName.get(), - self.run3DFlexMeshPrepJob.get()) - - - diff --git a/cryosparc2/protocols/protocol_cryosparc_3D_flex_mesh_prepare.py b/cryosparc2/protocols/protocol_cryosparc_3D_flex_mesh_prepare.py new file mode 100644 index 0000000..12eab54 --- /dev/null +++ b/cryosparc2/protocols/protocol_cryosparc_3D_flex_mesh_prepare.py @@ -0,0 +1,318 @@ +# ************************************************************************** +# * +# * Authors: Yunior C. Fonseca Reyna (cfonseca@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 import BETA +from pyworkflow.protocol.params import (PointerParam, FloatParam, IntParam, + BooleanParam, FileParam, StringParam) +from . import ProtCryosparcBase +from ..convert import * +from ..utils import * +from ..constants import * + + +class ProtCryoSparc3DFlexMeshPrepare(ProtCryosparcBase): + """ + Prepares particles for use in 3DFlex training and reconstruction. At the same + way, Takes in a consensus (rigid) refinement density map, plus optionally + a segmentation and generates a tetrahedral mesh for 3DFlex. + """ + _label = '3D flex mesh prepare' + _devStatus = BETA + _protCompatibility = [V4_1_0, V4_1_1, V4_1_2, V4_2_0, V4_2_1, V4_3_1, + V4_4_0] + + # --------------------------- DEFINE param functions ---------------------- + def _defineFileNames(self): + """ Centralize how files are called within the protocol. """ + myDict = { + 'input_particles': self._getTmpPath('input_particles.star'), + 'out_particles': self._getPath() + '/output_particle.star', + 'stream_log': self._getPath() + '/stream.log' + } + self._updateFilenamesDict(myDict) + + def _defineParams(self, form): + form.addSection(label='Input') + form.addParam('dataPrepare', PointerParam, + pointerClass='ProtCryoSparc3DFlexDataPrepare', + label="Data prepare protocol", + important=True, + help='Prepared particles and consensus volume from the 3D flex data prepare protocol.') + form.addParam('refMask', PointerParam, pointerClass='VolumeMask', + default=None, + label='Input Mask', + allowsNull=True, + help='Mask raw data') + + form.addSection(label='Mesh Prepare') + + solventMaskGroup = form.addGroup('Solvent mask preparation', + condition="refMask is None") + solventMaskGroup.addParam('mask_in_lowpass_A', IntParam, default=10, + condition="refMask is None", + label="Filter input volume res. (A)", + help="Filter the input consensus reconstruction " + "volume to this resolution (A) before thresholding " + "to create the outer solvent mask. Solvent mask is" + " only generated if it is not already input.") + + solventMaskGroup.addParam('mask_in_threshold_level', FloatParam, default=0.5, + condition="refMask is None", + label="Mask threshold", + help="Threshold the input consensus reconstruction " + "volume (after filtering) at this absolute " + "density level to make the solvent mask. " + "Solvent mask is only generated if it is " + "not already input.") + + solventMaskGroup.addParam('mask_dilate_A', IntParam, + default=2, + condition="refMask is None", + label="Mask dilation (A)", + help="After thresholding, dilate by this much" + " (A) to create solvent mask") + + solventMaskGroup.addParam('mask_pad_A', IntParam, + default=5, + condition="refMask is None", + label="Mask soft padding (A)", + help="After thresholding and dilation, soft " + "pad by this much (A) to create solvent mask") + + meshPreparationGroup = form.addGroup('Mesh preparation') + meshPreparationGroup.addParam('tetra_num_cells', IntParam, + default=20, + label="Base num. tetra cells", + help="Number of tetrahedral cells that fit " + "across the extent of the box. Use this " + "to set the size of mesh elements. If " + "set to e.g. 20, the base tetramesh will" + " have elements that are the right size " + "to create a spacing of 20 tetra elements " + "across the box extent in each x,y,z " + "direction. A higher number makes a finer mesh.") + + meshPreparationGroup.addParam('tetra_segments_path', FileParam, + allowsNull=True, + default=None, + label="Segmentation file path", + help="Absolute path to a segmentation file " + "in either .seg format from UCSF " + "Chimera Segger tool or else .mrc format " + "(see CryoSPARC guide for details), " + "defining subdomain regions that " + "should each have a submesh. " + "Submeshes are fused to make final " + "mesh using the segment connections list.") + + meshPreparationGroup.addParam('tetra_segments_fuse_list', StringParam, + allowsNull=True, + default=None, + label="Segment connections", + help="A comma and '>' separated list of " + "connections between segments to use when " + "fusing sub-meshes to make the final mesh. " + "See CryoSPARC guide for full explanation. " + "For example, '0>3, 0>4, 3>2, 2>1' is a valid connection" + " string. Each pair X>Y denotes that segment Y is " + "joined to segment X. The connections must form " + "a 'tree' structure and cannot have cycles. " + "The first pair X>Y must start with the " + "root of the tree as X. The connections " + "must be in breadth-first order of the tree. " + "When using Chimera Segger segmentation input, " + "the X>Y numbers should be region_ids from Segger (e.g., 948>947)") + + meshPreparationGroup.addParam('tetra_rigid_list', StringParam, + allowsNull=True, + default=None, + label="Rigid segments", + help="A comma separated list of segments to make " + "rigid. This is done by setting the rigidity " + "weight of tetra elements for this region to 20. ") + + + rigidityWeighting = form.addGroup('Rigidity weighting') + rigidityWeighting.addParam('rigidity_penalty_min', FloatParam, + default=0.5, + label="Min. rigidity weight", + help="Rigidity weights of tetra elements" + " are 1.0 in the most dense regions " + "of the input consensus map, and fall " + "off to this value (default 0.5) in the" + " least dense/empty regions. This helps " + "encourage the deformation model to expand/contract" + " empty space without distorting the protein density.") + + rigidityWeighting.addParam('rigidity_penalty_stiffen_low_density', BooleanParam, + default=False, + label="Stiffen low density regions", + help="Turning this on will cause the rigidity" + " weights of tetra elements at the periphery " + "of the input consensus density to be increased to 3.0 . " + "Empty regions will still have low rigidity, but non-empty " + "regions at the periphery of the structure will be rigidified." + " This helps to combat overfitting in smaller particles or" + " poor SNR data where otherwise low density peripheral " + "features start to 'fly around'. However, it can also cause " + "the deformations to be overly smooth and blur motion 'boundaries between domains. ") + + + # --------------[Compute settings]--------------------------- + form.addSection(label="Compute settings") + addComputeSectionParams(form, allowMultipleGPUs=False, needGPU=False) + + def _insertAllSteps(self): + self._defineFileNames() + self._defineParamsPrepareName() + self._defineParamsMeshName() + self._initializeCryosparcProject() + self._insertFunctionStep(self.convertInputStep) + self._insertFunctionStep(self.mergePrepareStep) + self._insertFunctionStep(self.createOutputStep) + + # --------------------------- STEPS functions ------------------------------ + def mergePrepareStep(self): + self.info(pwutils.yellowStr("3D Flex Mesh Preparation started...")) + self.doRun3DFlexMeshPrepare() + + def createOutputStep(self): + """ + Create the protocol output. Convert cryosparc file to Relion file + """ + pass + # print(pwutils.yellowStr("Creating the output..."), flush=True) + # csOutputFolder = os.path.join(self.projectDir.get(), + # self.run3DFlexDataPrepJob.get()) + # csParticlesName = "%s_passthrough_particles.cs" % self.run3DFlexDataPrepJob.get() + # fnVolName = "%s_map.mrc" % self.run3DFlexDataPrepJob.get() + # + # # Copy the CS output volume and half to extra folder + # copyFiles(csOutputFolder, self._getExtraPath(), files=[csParticlesName, + # fnVolName]) + # + # csFile = os.path.join(self._getExtraPath(), csParticlesName) + # + # outputStarFn = self._getFileName('out_particles') + # argsList = [csFile, outputStarFn] + # convertCs2Star(argsList) + # + # fnVol = os.path.join(self._getExtraPath(), fnVolName) + # imgSet = self._getInputParticles() + # vol = Volume() + # fixVolume([fnVol]) + # vol.setFileName(fnVol) + # vol.setSamplingRate(calculateNewSamplingRate(vol.getDim(), + # imgSet.getSamplingRate(), + # imgSet.getDim())) + # outImgSet = self._createSetOfParticles() + # outImgSet.copyInfo(imgSet) + # self._fillDataFromIter(outImgSet) + # + # self._defineOutputs(outputVolume=vol) + # self._defineSourceRelation(self.inputParticles.get(), vol) + # self._defineOutputs(outputParticles=outImgSet) + + # ------------------------- Utils methods ---------------------------------- + + def _fillDataFromIter(self, imgSet): + outImgsFn = 'particles@' + self._getFileName('out_particles') + imgSet.setAlignmentProj() + imgSet.copyItems(self._getInputParticles(), + updateItemCallback=self._createItemMatrix, + itemDataIterator=emtable.Table.iterRows(fileName=outImgsFn)) + + def _createItemMatrix(self, particle, row): + createItemMatrix(particle, row, align=ALIGN_PROJ) + setCryosparcAttributes(particle, row, RELIONCOLUMNS.rlnRandomSubset.value) + + def _validate(self): + validateMsgs = cryosparcValidate() + return validateMsgs + + def _defineParamsPrepareName(self): + """ Define a list with 3D Flex Prepare Data parameters names""" + self._paramsPrepareName = ['box_size_pix', 'bin_size_pix', 'alpha_min', + 'keep_num_particles'] + self.lane = str(self.getAttributeValue('compute_lane')) + + def _defineParamsMeshName(self): + """ Define a list with 3D Flex Mesh Prepare parameters names""" + self._maskMeshPrepareName = ['mask_in_lowpass_A', 'mask_in_threshold_level', + 'mask_dilate_A', 'mask_pad_A'] + self._paramsMeshName = ['tetra_num_cells', 'tetra_segments_path', + 'tetra_segments_fuse_list', 'tetra_rigid_list', + 'rigidity_penalty_min', + 'rigidity_penalty_stiffen_low_density'] + self.lane = str(self.getAttributeValue('compute_lane')) + + def doRun3DFlexMeshPrepare(self): + self._className = 'flex_meshprep' + params = {} + varDataPrepJob = str(self.dataPrepare.get().run3DFlexDataPrepJob) + input_group_connect = {"volume": str('%s.volume' % varDataPrepJob)} + + if self.refMask.get() is not None: + input_group_connect["mask"] = str(self.mask) + else: + for paramName in self._maskMeshPrepareName: + if self.getAttributeValue(paramName) is not None: + params[str(paramName)] = str(self.getAttributeValue(paramName)) + + for paramName in self._paramsMeshName: + if self.getAttributeValue(paramName) is not None: + params[str(paramName)] = str(self.getAttributeValue(paramName)) + + run3DFlexMeshPrepJob = enqueueJob(self._className, + self.projectName.get(), + self.workSpaceName.get(), + str(params).replace('\'', '"'), + str(input_group_connect).replace('\'', + '"'), + self.lane, False) + + self.run3DFlexMeshPrepJob = String(run3DFlexMeshPrepJob.get()) + self.currenJob.set(self.run3DFlexMeshPrepJob.get()) + self._store(self) + + waitForCryosparc(self.projectName.get(), + self.run3DFlexMeshPrepJob.get(), + "An error occurred in the 3D Flex Mesh Preparation process. " + "Please, go to cryoSPARC software for more " + "details.") + clearIntermediateResults(self.projectName.get(), + self.run3DFlexMeshPrepJob.get()) + + + + + + + + + + diff --git a/cryosparc2/protocols/protocol_cryosparc_blob_picker.py b/cryosparc2/protocols/protocol_cryosparc_blob_picker.py new file mode 100644 index 0000000..db192d1 --- /dev/null +++ b/cryosparc2/protocols/protocol_cryosparc_blob_picker.py @@ -0,0 +1,324 @@ +# ************************************************************************** +# * +# * Authors: Yunior C. Fonseca Reyna (cfonseca@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' +# * +# ************************************************************************** +import os +import emtable + +from pwem.objects import Coordinate, CTFModel +import pyworkflow.utils as pwutils +from pyworkflow import NEW +from pyworkflow.protocol.params import (PointerParam, FloatParam, + BooleanParam, IntParam, + String) + +from .protocol_base import ProtCryosparcBase +from .. import RELIONCOLUMNS +from ..convert import convertCs2Star +from ..utils import (addComputeSectionParams, cryosparcValidate, enqueueJob, waitForCryosparc, clearIntermediateResults, + copyFiles) + + +class ProtCryoSparcBlobPicker(ProtCryosparcBase): + """ + Automatically picks particles by searching for Gaussian signals. + """ + _label = 'blob_picker' + _className = "blob_picker_gpu" + _devStatus = NEW + + def _defineParams(self, form): + form.addSection(label='Input') + form.addParam('inputMicrographs', PointerParam, important=True, + label=pwutils.Message.LABEL_INPUT_MIC, + pointerClass='SetOfMicrographs') + + form.addParam('diameter', IntParam, default=None, + label='Minimum particle diameter (px)', + help='Minimum particle diameter (px)') + + form.addParam('diameter_max', IntParam, default=None, + label='Maximum particle diameter (px)', + help='Maximum particle diameter (px)') + + form.addParam('use_circle', BooleanParam, default=True, + label='Use circular blob', + help='Use three circular blobs, at minimum, average, and maximum particle diameters based on parameters above.') + + form.addParam('use_ellipse', BooleanParam, default=False, + label='Use elliptical blob', + help='Use an elliptical blob with minor diameter equal to minimum particle diamater param above, and major diameter equal to maximum particle diameter above. You may want to turn off circular blob with this on.') + + form.addParam('use_ring', BooleanParam, default=False, + label='Use ring blob', + help='Use a ring-shaped blob with inner diameter equal to minimum particle diameter param above, and outer diameter equal to maximum particle diameter above. You may want to turn off circular and elliptical blob with this on.') + + form.addParam('estimate_ctf', BooleanParam, default=False, + label='Estimate CTF before pick?', + help='Estimate CTF using cryoSPARC Patch CTF algorithm') + + # form.addParam('lowpass_res_template', IntParam, default=20, + # label='Lowpass filter to apply to templates (A)', + # help='Lowpass filter to apply to templates, (A)s') + # + # form.addParam('lowpass_res', IntParam, default=20, + # label='Lowpass filter to apply to micrographs (A)', + # help='Lowpass filter to apply to micrographs, (A)s') + # + # form.addParam('angular_spacing_deg', IntParam, default=5, + # label='Angular sampling (degrees)', + # help='Angular sampling of templates in degrees. Lower value will mean finer rotations.') + + form.addParam('min_distance', FloatParam, default=1.0, + label='Min. separation dist (diameters)', + help='Minimum distance between particles in units of particle diameter (min diameter for blob picker). The lower this value, the more and closer particles it picks.') + + form.addParam('num_process', IntParam, default=None, + allowsNull=True, + label='Number of mics to process', + help='Number of micrographs to process. None means all.') + + form.addParam('max_num_hits', IntParam, default=4000, + label='Maximum number of local maxima to consider', + help='Maximum number of local maxima (peaks) considered.') + + """ + job.param_add('template', "num_plot", base_value=10, title="Number of mics to plot", param_type="number", hidden=False, advanced=False, desc='Number of micrographs to plot.') + job.param_add('template', "recenter_templates", base_value=True, title="Recenter templates", param_type="boolean", hidden=False, advanced=True, desc='Whether or not to recenter the input templates.') + """ + + # --------------[Compute settings]--------------------------- + form.addSection(label="Compute settings") + addComputeSectionParams(form, allowMultipleGPUs=False) + + # --------------------------- INSERT steps functions ----------------------- + + def _insertAllSteps(self): + self._defineParamsName() + self._initializeCryosparcProject() + self._insertFunctionStep(self.convertInputStep) + self._insertFunctionStep(self.processStep) + self._insertFunctionStep(self.createOutputStep) + + # --------------------------- STEPS functions ------------------------------ + + def processStep(self): + if self.estimate_ctf.get(): + self.info(pwutils.yellowStr("Patch CTF estimate started...")) + self.doPatchCTFEstimate() + self.micrographs = String(str(self.runPatchCTF.get()) + '.exposures') + + self.info(pwutils.yellowStr("Blob picker started...")) + self.doBlobPicker() + + def createOutputStep(self): + """ + Create the protocol output. Convert cryosparc file to star file + """ + self.info(pwutils.yellowStr("Create output started...")) + self._initializeUtilsVariables() + micSetPtr = self._getInputMicrographs() + + micList = {os.path.basename(mic.getFileName()): mic.clone() for mic in micSetPtr} + + csOutputFolder = os.path.join(self.projectDir.get(), + self.runBlobPicker.get()) + # Copy the CS output coordinates to extra folder + outputPath = os.path.join(self._getExtraPath(), self.runBlobPicker.get()) + copyFiles(csOutputFolder, outputPath) + csPickedParticlesName = 'picked_particles.cs' + + csFile = os.path.join(outputPath, csPickedParticlesName) + outputStarFn = self._getExtraPath('output_coordinates.star') + argsList = [csFile, outputStarFn] + convertCs2Star(argsList) + + outputCoords = self._fillSetOfCoordinates(micSetPtr, outputStarFn, micList) + + # Copy the CTF output to extra folder + if self.estimate_ctf.get(): + csOutputFolder = os.path.join(self.projectDir.get(), + self.runPatchCTF.get()) + outputPath = os.path.join(self._getExtraPath(), self.runPatchCTF.get()) + copyFiles(csOutputFolder, outputPath) + + ctfEstimatedFileName = 'exposures_ctf_estimated.cs' + csFile = os.path.join(outputPath, ctfEstimatedFileName) + outputStarFn = self._getExtraPath('ctf.star') + argsList = [csFile, outputStarFn] + convertCs2Star(argsList) + + outputCtfSet = self._fillSetOfCTF(outputStarFn, micList) + + self._defineOutputs(outputCTF=outputCtfSet) + self._defineSourceRelation(micSetPtr, outputCtfSet) + + self._defineOutputs(outputCoordinates=outputCoords) + self._defineSourceRelation(micSetPtr, outputCoords) + + def _fillSetOfCoordinates(self, micSetPtr, outputStarFn, micList): + + outputCoords = self._createSetOfCoordinates(micSetPtr) + boxSixe = (self.diameter.get() + self.diameter_max.get()) / 2 + outputCoords.setBoxSize(int(boxSixe)) + + coord = Coordinate() + mdFileName = '%s@%s' % ('particles', outputStarFn) + table = emtable.Table(fileName=outputStarFn) + + for row in table.iterRows(mdFileName): + coord.setObjId(None) + micName = os.path.basename(row.get(RELIONCOLUMNS.rlnMicrographName.value)) + splitMicName = micName.split('_') + if len(splitMicName) > 1: + micName = '_'.join(splitMicName[1:]) + else: + micName = splitMicName[-1] + coord.setMicrograph(micList[micName]) + x = row.get(RELIONCOLUMNS.rlnCoordinateX.value) + y = row.get(RELIONCOLUMNS.rlnCoordinateY.value) + dim = micList[micName].getDimensions() + flipY = dim[1] - y + coord.setPosition(x, flipY) + # Add it to the set + outputCoords.append(coord) + + return outputCoords + + def _fillSetOfCTF(self, outputCTFFn, micList): + + inputMics = self._getInputMicrographs() + outputCtfSet = self._createSetOfCTF() + outputCtfSet.setMicrographs(inputMics) + mics = list(micList.values()) + + ctf = CTFModel() + mdFileName = '%s@%s' % ('micrograph', outputCTFFn) + table = emtable.Table(fileName=outputCTFFn) + + for mic, row in enumerate(table.iterRows(mdFileName)): + ctf.setDefocusU(row.get(RELIONCOLUMNS.rlnDefocusU.value)) + ctf.setDefocusV(row.get(RELIONCOLUMNS.rlnDefocusV.value)) + ctf.setPhaseShift(row.get(RELIONCOLUMNS.rlnPhaseShift.value)) + ctf.setResolution(row.get(RELIONCOLUMNS.rlnCtfMaxResolution.value)) + ctf.setDefocusAngle(row.get(RELIONCOLUMNS.rlnDefocusAngle.value)) + ctf.setMicrograph(mics[mic]) + outputCtfSet.append(ctf) + + return outputCtfSet + + def _defineParamsName(self): + """ Define a list with all protocol parameters names""" + self._paramsName = ['diameter', 'diameter_max', 'use_circle', + 'use_ellipse', 'use_ring', 'min_distance', + 'num_process', 'max_num_hits'] + + self.lane = str(self.getAttributeValue('compute_lane')) + + # --------------------------- INFO functions ------------------------------- + def _validate(self): + """ Should be overwritten in subclasses to + return summary message for NORMAL EXECUTION. + """ + validateMsgs = cryosparcValidate() + + # if not validateMsgs: + # micrographs = self._getInputMicrographs() + # if micrographs is not None and not micrographs.hasCTF(): + # validateMsgs.append("The micrographs has not associated a CTF model") + + return validateMsgs + + def _summary(self): + summary = [] + return summary + + def doPatchCTFEstimate(self): + input_group_connect = {"exposures": self.micrographs.get()} + params = {'classic_mode': 'False'} + className = 'patch_ctf_estimation_multi' + try: + gpusToUse = self.getGpuList() + except Exception: + gpusToUse = False + + runPatchCTFJob = enqueueJob(className, + self.projectName.get(), + self.workSpaceName.get(), + str(params).replace('\'', '"'), + str(input_group_connect).replace('\'', '"'), + self.lane, gpusToUse) + + self.runPatchCTF = String(runPatchCTFJob.get()) + self.currenJob.set(runPatchCTFJob.get()) + self._store(self) + + waitForCryosparc(self.projectName.get(), + self.runPatchCTF.get(), + "An error occurred in the ctf estimation process. " + "Please, go to cryoSPARC software for more " + "details.") + + def doBlobPicker(self): + + input_group_connect = {"micrographs": self.micrographs.get()} + params = {} + micSetPtr = self._getInputMicrographs() + samplingRate = micSetPtr.getSamplingRate() + for paramName in self._paramsName: + if (paramName != 'diameter' and paramName != 'diameter_max' and + paramName != 'num_process'): + params[str(paramName)] = str(self.getAttributeValue(paramName)) + elif paramName == 'diameter' and self.diameter.get() is not None: + params[str(paramName)] = str(int(self.diameter.get()*samplingRate)) + elif paramName == 'diameter_max' and self.diameter_max.get() is not None: + params[str(paramName)] = str(int(self.diameter_max.get()*samplingRate)) + elif paramName == 'num_process' and self.num_process.get() is not None: + params[str(paramName)] = str(self.num_process.get()) + + # Determinate the GPUs to use (in dependence of + # the cryosparc version) + try: + gpusToUse = self.getGpuList() + except Exception: + gpusToUse = False + + runBlobPickerJob = enqueueJob(self._className, + self.projectName.get(), + self.workSpaceName.get(), + str(params).replace('\'', '"'), + str(input_group_connect).replace('\'', '"'), + self.lane, gpusToUse) + + self.runBlobPicker = String(runBlobPickerJob.get()) + self.currenJob.set(runBlobPickerJob.get()) + self._store(self) + + waitForCryosparc(self.projectName.get(), + self.runBlobPicker.get(), + "An error occurred in the particles picking process. " + "Please, go to cryoSPARC software for more " + "details.") + clearIntermediateResults(self.projectName.get(), self.runBlobPicker.get()) \ No newline at end of file diff --git a/cryosparc2/protocols/protocol_cryosparc_helical_refinement.py b/cryosparc2/protocols/protocol_cryosparc_helical_refinement.py index 8c01ccb..b0c70bb 100644 --- a/cryosparc2/protocols/protocol_cryosparc_helical_refinement.py +++ b/cryosparc2/protocols/protocol_cryosparc_helical_refinement.py @@ -52,7 +52,6 @@ class ProtCryoSparcHelicalRefine3D(ProtCryoSparc3DHomogeneousRefine): other cryoSPARC refinement jobs. """ _label = '3D helical refinement' - _devStatus = BETA _fscColumns = 4 _protCompatibility = [V3_3_1, V3_3_2, V4_0_0, V4_0_1, V4_0_2, V4_0_3, V4_1_0, V4_1_1, V4_1_2, V4_2_0, V4_2_1, V4_3_1, V4_4_0] diff --git a/cryosparc2/protocols/protocol_cryosparc_new_3D_classification.py b/cryosparc2/protocols/protocol_cryosparc_new_3D_classification.py index 4e3b7b0..7474127 100644 --- a/cryosparc2/protocols/protocol_cryosparc_new_3D_classification.py +++ b/cryosparc2/protocols/protocol_cryosparc_new_3D_classification.py @@ -60,7 +60,6 @@ class ProtCryoSparcNew3DClassification(ProtCryosparcBase): """ _label = '3D Classification' _className = "class_3D" - _devStatus = BETA _protCompatibility = [V3_3_1, V3_3_2, V4_0_0, V4_0_1, V4_0_2, V4_0_3, V4_1_0, V4_1_1, V4_1_2, V4_2_0, V4_2_1, V4_3_1, V4_4_0] diff --git a/cryosparc2/protocols/protocol_cryosparc_new_local_refine.py b/cryosparc2/protocols/protocol_cryosparc_new_local_refine.py index df39f91..08615d1 100644 --- a/cryosparc2/protocols/protocol_cryosparc_new_local_refine.py +++ b/cryosparc2/protocols/protocol_cryosparc_new_local_refine.py @@ -56,7 +56,6 @@ class ProtCryoSparcLocalRefine(ProtCryosparcBase, ProtOperateParticles): Subtract projections of a masked volume from particles. """ _label = 'local refinement' - _devStatus = NEW _protCompatibility = [V3_3_1, V3_3_2, V4_0_0, V4_0_1, V4_0_2, V4_0_3, V4_1_0, V4_1_1, V4_1_2, V4_2_0, V4_2_1, V4_3_1, V4_4_0] _className = "new_local_refine" diff --git a/cryosparc2/protocols/protocol_cryosparc_patch_ctf_estimation.py b/cryosparc2/protocols/protocol_cryosparc_patch_ctf_estimation.py new file mode 100644 index 0000000..9973857 --- /dev/null +++ b/cryosparc2/protocols/protocol_cryosparc_patch_ctf_estimation.py @@ -0,0 +1,222 @@ +# ************************************************************************** +# * +# * Authors: Yunior C. Fonseca Reyna (cfonseca@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' +# * +# ************************************************************************** + + +import os +import emtable +import numpy + +from pwem.objects import CTFModel +import pyworkflow.utils as pwutils +from pyworkflow import NEW +from pyworkflow.protocol.params import (PointerParam, FloatParam, + BooleanParam, IntParam, + String) + +from .protocol_base import ProtCryosparcBase +from .. import RELIONCOLUMNS +from ..convert import convertCs2Star +from ..utils import (addComputeSectionParams, cryosparcValidate, enqueueJob, waitForCryosparc, + copyFiles) + + +class ProtCryoSparcPatchCTFEstimate(ProtCryosparcBase): + """ + Patch-based CTF estimation automatically estimates defocus variation for tilted, bent, + deformed samples and is accurate for all particle sizes and types including flexible and membrane proteins. + """ + _label = 'ctf_estimation' + _className = "patch_ctf_estimation_multi" + _devStatus = NEW + + def _defineParams(self, form): + form.addSection(label='Input') + form.addParam('inputMicrographs', PointerParam, important=True, + label=pwutils.Message.LABEL_INPUT_MIC, + pointerClass='SetOfMicrographs') + + form.addParam('amp_contrast', FloatParam, default=0.1, + label='Amplitude Contrast', + help='Amplitude constrast to use. Typically 0.07 or 0.1 for cryo-EM data.') + + form.addParam('res_min_align', IntParam, default=25, + label='Minimum resolution (A)', + help='Minimum resolution (in A) to consider when estimating CTF.') + + form.addParam('res_max_align', IntParam, default=4, + label='Maximum resolution (A)', + help='Maximum resolution (in A) to consider when estimating CTF.') + + form.addParam('df_search_min', IntParam, default=1000, + label='Maximum resolution (A)', + help='Defocus range for gridsearch.') + + form.addParam('df_search_max', IntParam, default=40000, + label='Maximum resolution (A)', + help='Defocus range for gridsearch.') + + form.addParam('phase_shift_min', IntParam, default=0, + label='Min. search phase-shift (rad)', + help='Phase-shift range for gridsearch.') + + form.addParam('phase_shift_max', FloatParam, default=numpy.pi, + label='Min. search phase-shift (rad)', + help='Phase-shift range for gridsearch.') + + form.addParam('do_phase_shift_refine_only', BooleanParam, default=False, + label='Do phase refine only', + help='Whether to carry out refinement over phase shift only') + + """job.param_add('ctf_settings', "override_K_Y", base_value=None, title="Override knots Y", param_type="number", + hidden=False, advanced=True, + desc='Override automatically selected spline order for Y dimension (vertical)') + job.param_add('ctf_settings', "override_K_X", base_value=None, title="Override knots X", param_type="number", + hidden=False, advanced=True, + desc='Override automatically selected spline order for X dimension (horizontal)') + + job.param_add_section('compute_settings', title='Compute settings', desc='') + job.param_add('compute_settings', "compute_num_gpus", base_value=1, title="Number of GPUs to parallelize", + param_type="number", hidden=False, advanced=False, + desc='Number of GPUs over which to parallelize computation.')""" + + # --------------[Compute settings]--------------------------- + form.addSection(label="Compute settings") + addComputeSectionParams(form, allowMultipleGPUs=False) + + # --------------------------- INSERT steps functions ----------------------- + + def _insertAllSteps(self): + self._defineParamsName() + self._initializeCryosparcProject() + self._insertFunctionStep(self.convertInputStep) + self._insertFunctionStep(self.processStep) + self._insertFunctionStep(self.createOutputStep) + + # --------------------------- STEPS functions ------------------------------ + + def processStep(self): + self.info(pwutils.yellowStr("Patch CTF estimate started...")) + self.doPatchCTFEstimate() + self.micrographs = String(str(self.runPatchCTF.get()) + '.exposures') + + def createOutputStep(self): + """ + Create the protocol output. Convert cryosparc file to star file + """ + self.info(pwutils.yellowStr("Create output started...")) + self._initializeUtilsVariables() + micSetPtr = self._getInputMicrographs() + + micList = {os.path.basename(mic.getFileName()): mic.clone() for mic in micSetPtr} + + # Copy the CTF output to extra folder + csOutputFolder = os.path.join(self.projectDir.get(), + self.runPatchCTF.get()) + outputPath = os.path.join(self._getExtraPath(), self.runPatchCTF.get()) + copyFiles(csOutputFolder, outputPath) + + ctfEstimatedFileName = 'exposures_ctf_estimated.cs' + csFile = os.path.join(outputPath, ctfEstimatedFileName) + outputStarFn = self._getExtraPath('ctf.star') + argsList = [csFile, outputStarFn] + convertCs2Star(argsList) + + outputCtfSet = self._fillSetOfCTF(outputStarFn, micList) + + self._defineOutputs(outputCTF=outputCtfSet) + self._defineSourceRelation(micSetPtr, outputCtfSet) + + def _fillSetOfCTF(self, outputCTFFn, micList): + + inputMics = self._getInputMicrographs() + outputCtfSet = self._createSetOfCTF() + outputCtfSet.setMicrographs(inputMics) + mics = list(micList.values()) + + ctf = CTFModel() + mdFileName = '%s@%s' % ('micrograph', outputCTFFn) + table = emtable.Table(fileName=outputCTFFn) + + for mic, row in enumerate(table.iterRows(mdFileName)): + ctf.setDefocusU(row.get(RELIONCOLUMNS.rlnDefocusU.value)) + ctf.setDefocusV(row.get(RELIONCOLUMNS.rlnDefocusV.value)) + ctf.setPhaseShift(row.get(RELIONCOLUMNS.rlnPhaseShift.value)) + ctf.setResolution(row.get(RELIONCOLUMNS.rlnCtfMaxResolution.value)) + ctf.setDefocusAngle(row.get(RELIONCOLUMNS.rlnDefocusAngle.value)) + ctf.setMicrograph(mics[mic]) + outputCtfSet.append(ctf) + + return outputCtfSet + + def _defineParamsName(self): + """ Define a list with all protocol parameters names""" + + self._paramsName = ['amp_contrast', 'res_min_align', 'res_max_align', 'df_search_min', + 'df_search_max', 'phase_shift_min', 'phase_shift_max', 'do_phase_shift_refine_only'] + + self.lane = str(self.getAttributeValue('compute_lane')) + + # --------------------------- INFO functions ------------------------------- + def _validate(self): + """ Should be overwritten in subclasses to + return summary message for NORMAL EXECUTION. + """ + validateMsgs = cryosparcValidate() + return validateMsgs + + def _summary(self): + summary = [] + return summary + + def doPatchCTFEstimate(self): + input_group_connect = {"exposures": self.micrographs.get()} + params = {'classic_mode': 'False'} + try: + gpusToUse = self.getGpuList() + except Exception: + gpusToUse = False + + for paramName in self._paramsName: + params[str(paramName)] = str(self.getAttributeValue(paramName)) + + runPatchCTFJob = enqueueJob(self._className, + self.projectName.get(), + self.workSpaceName.get(), + str(params).replace('\'', '"'), + str(input_group_connect).replace('\'', '"'), + self.lane, gpusToUse) + + self.runPatchCTF = String(runPatchCTFJob.get()) + self.currenJob.set(runPatchCTFJob.get()) + self._store(self) + + waitForCryosparc(self.projectName.get(), + self.runPatchCTF.get(), + "An error occurred in the ctf estimation process. " + "Please, go to cryoSPARC software for more " + "details.") + diff --git a/cryosparc2/utils.py b/cryosparc2/utils.py index 38d94d2..a1ca483 100644 --- a/cryosparc2/utils.py +++ b/cryosparc2/utils.py @@ -335,7 +335,7 @@ def getProjectInformation(project_uid, info='project_dir'): """ Get information about a single project :param project_uid: the id of the project - :return: the information related to the project thats stored in the database + :return: the information related to the project that's stored in the database """ import ast getProject_cmd = (getCryosparcProgram() + @@ -456,6 +456,41 @@ def doImportVolumes(protocol, refVolumePath, refVolume, volType, msg): return importedVolume +def doImportMicrographs(protocol): + print(pwutils.yellowStr("Importing micrographs..."), flush=True) + className = "import_micrographs" + micrographs = protocol._getInputMicrographs() + acquisition = micrographs.getAcquisition() + micList = list(micrographs.getFiles()) + + micFolder = os.path.join(protocol._getExtraPath('micrographs')) + os.makedirs(micFolder, exist_ok=True) + + for micPath in micList: + micName = os.path.basename(micPath) + micLink = os.path.join(micFolder, micName) + os.symlink(os.path.abspath(micPath), micLink) + micExt = '*%s' % os.path.splitext(micList[0])[1] + + params = {"blob_paths": str(os.path.join(os.getcwd(), micFolder, micExt)), + "psize_A": str(micrographs.getSamplingRate()), + "accel_kv": str(acquisition.getVoltage()), + "cs_mm": str(acquisition.getSphericalAberration()), + "total_dose_e_per_A2": str(0.1), + "output_constant_ctf": "True" + } + + import_particles = enqueueJob(className, protocol.projectName, protocol.workSpaceName, + str(params).replace('\'', '"'), '{}', protocol.lane) + + waitForCryosparc(protocol.projectName.get(), import_particles.get(), + "An error occurred importing particles. " + "Please, go to cryoSPARC software for more " + "details.") + + return import_particles + + def doJob(jobType, projectName, workSpaceName, params, input_group_connect): """ do_job(job_type, puid='P1', wuid='W1', uuid='devuser', params={},