Source code for cistem.protocols.protocol_ctffind

# **************************************************************************
# *
# * Authors:     Josue Gomez BLanco (josue.gomez-blanco@mcgill.ca) [1]
# *              J.M. De la Rosa Trevin (delarosatrevin@scilifelab.se) [2]
# *              Grigory Sharov (gsharov@mrc-lmb.cam.ac.uk) [3]
# *
# * [1] Department of Anatomy and Cell Biology, McGill University
# * [2] SciLifeLab, Stockholm University
# * [3] MRC Laboratory of Molecular Biology (MRC-LMB)
# *
# * 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 3 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 pyworkflow.utils as pwutils
from pwem.protocols import ProtCTFMicrographs
from pwem.objects import CTFModel
from pwem import emlib

from .program_ctffind import ProgramCtffind


[docs]class CistemProtCTFFind(ProtCTFMicrographs): """ Estimate CTF for a set of micrographs with ctffind4. To find more information about ctffind4 visit: https://grigoriefflab.umassmed.edu/ctffind4 """ _label = 'ctffind4' def _defineParams(self, form): ProgramCtffind.defineInputParams(form) ProgramCtffind.defineProcessParams(form) self._defineStreamingParams(form) def _defineCtfParamsDict(self): ProtCTFMicrographs._defineCtfParamsDict(self) self._ctfProgram = ProgramCtffind(self) # -------------------------- STEPS functions ------------------------------ def _doCtfEstimation(self, mic, **kwargs): """ Run ctffind with required parameters. :param mic: input mic object :param kwargs: dict with arguments """ if self.usePowerSpectra: micFn = mic._powerSpectra.getFileName() powerSpectraPix = self.psSampling else: micFn = mic.getFileName() powerSpectraPix = None micDir = self._getTmpPath('mic_%06d' % mic.getObjId()) # Create micrograph dir pwutils.makePath(micDir) micFnMrc = os.path.join(micDir, pwutils.replaceBaseExt(micFn, 'mrc')) ih = emlib.image.ImageHandler() if not os.path.exists(micFn): raise FileNotFoundError("Missing input micrograph: %s" % micFn) if micFn.endswith('.mrc'): pwutils.createAbsLink(os.path.abspath(micFn), micFnMrc) else: ih.convert(micFn, micFnMrc, emlib.DT_FLOAT) try: program, args = self._ctfProgram.getCommand( micFn=micFnMrc, powerSpectraPix=powerSpectraPix, ctffindOut=self._getCtfOutPath(mic), ctffindPSD=self._getPsdPath(mic), **kwargs) self.runJob(program, args) pwutils.cleanPath(micDir) except Exception as e: self.error("ERROR: Ctffind has failed for %s. %s" % ( micFnMrc, self._getErrorFromCtffindTxt(mic, e))) def _getErrorFromCtffindTxt(self, mic, e): """ Parse output log for errors. :param mic: input mic object :return: the error string """ file = self._getCtfOutPath(mic) with open(file, "r") as fh: for line in fh.readlines(): if line.startswith("Error"): return line.replace("Error:", "") return e def _estimateCTF(self, mic, *args): """ Redefined func from the base class. """ self._doCtfEstimation(mic) def _reEstimateCTF(self, mic, ctf): """ Redefined func from the base class. """ self._doCtfEstimation(mic, **self._getRecalCtfParamsDict(ctf)) def _createCtfModel(self, mic, updateSampling=False): """ Redefined func from the base class. """ psd = self._getPsdPath(mic) ctfModel = self._ctfProgram.parseOutputAsCtf(self._getCtfOutPath(mic), psdFile=psd) ctfModel.setMicrograph(mic) return ctfModel def _createOutputStep(self): """ Do nothing in streaming case. """ pass # -------------------------- INFO functions ------------------------------- def _validate(self): errors = [] if self.inputType == 0: errors.append('Movie CTF estimation is not supported yet.') if self.usePowerSpectra: mic = self._getFirstMic() if not hasattr(mic, "_powerSpectra"): errors.append("Input micrographs do not have associated power spectra.") else: self.psSampling = mic._powerSpectra.getSamplingRate() if self.lowRes.get() > 50: errors.append("Minimum resolution cannot be > 50A.") valueStep = round(self.stepPhaseShift.get(), 2) valueMin = round(self.minPhaseShift.get(), 2) valueMax = round(self.maxPhaseShift.get(), 2) if not (self.minPhaseShift < self.maxPhaseShift and valueStep <= (valueMax - valueMin) and 0. <= valueMax <= 180.): errors.append('Wrong values for phase shift search.') return errors def _citations(self): return ["Mindell2003", "Rohou2015"] def _methods(self): methods = [] if self.inputMicrographs.get() is not None: methods.append("We calculated the CTF of %s using CTFFind. " % self.getObjectTag('inputMicrographs')) methods.append(self.methodsVar.get('')) if getattr(self, 'outputCTF'): methods.append('Output CTFs: %s' % self.getObjectTag('outputCTF')) return methods # -------------------------- UTILS functions ------------------------------ def _getRecalCtfParamsDict(self, ctfModel): """ Update values from user-adjusted params from the Java GUI. :param ctfModel: input CTF model object """ values = [float(x) for x in ctfModel.getObjComment().split()] sampling = self.inputMicrographs.get().getSamplingRate() return { 'step_focus': 500.0, 'lowRes': sampling / values[3], 'highRes': sampling / values[4], 'minDefocus': min([values[0], values[1]]), 'maxDefocus': max([values[0], values[1]]) } def _getMicExtra(self, mic, suffix): """ Return a file in extra direction with root of micFn. :param mic: input mic object :param suffix: file extension """ return self._getExtraPath(pwutils.removeBaseExt(os.path.basename( mic.getFileName())) + '_' + suffix) def _getPsdPath(self, mic): return self._getMicExtra(mic, 'ctf.mrc') def _getCtfOutPath(self, mic): return self._getMicExtra(mic, 'ctf.txt') def _parseOutput(self, filename): """ Try to find the output estimation parameters from filename. It searches for a first line without #. """ return self._ctfProgram.parseOutput(filename) def _getCTFModel(self, defocusU, defocusV, defocusAngle, psdFile): ctf = CTFModel() ctf.setStandardDefocus(defocusU, defocusV, defocusAngle) ctf.setPsdFile(psdFile) return ctf def _getFirstMic(self): """ Get first mic in the input set only once. """ return self.getInputMicrographs().getFirstItem()