Source code for gctf.protocols.protocol_gctf

# **************************************************************************
# *
# * Authors:     Grigory Sharov (gsharov@mrc-lmb.cam.ac.uk)
# *
# * 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 pyworkflow.constants import PROD
from pwem import emlib
from pwem.objects import CTFModel
from pwem.protocols import ProtCTFMicrographs


from .program_gctf import ProgramGctf


[docs]class ProtGctf(ProtCTFMicrographs): """ Estimates CTF on a set of micrographs using Gctf. To find more information about Gctf go to: https://www2.mrc-lmb.cam.ac.uk/research/locally-developed-software/zhang-software/#gctf """ _label = 'ctf estimation' _devStatus = PROD def _defineCtfParamsDict(self): ProtCTFMicrographs._defineCtfParamsDict(self) self._gctfProgram = ProgramGctf(self) def _defineParams(self, form): ProgramGctf.defineInputParams(form) ProgramGctf.defineProcessParams(form) self._defineStreamingParams(form) # -------------------------- STEPS functions ------------------------------ def _estimateCTF(self, mic, *args): self._estimateCtfList([mic], *args) def _estimateCtfList(self, micList, *args, **kwargs): """ Estimate several micrographs at once, probably a bit more efficient. """ try: micPath = self._getMicrographDir(micList[0]) if len(micList) > 1: micPath += ('-%04d' % micList[-1].getObjId()) pwutils.makePath(micPath) ih = emlib.image.ImageHandler() for mic in micList: micFn = mic.getFileName() # We convert the input micrograph on demand if not in .mrc downFactor = self.ctfDownFactor.get() micFnMrc = os.path.join(micPath, pwutils.replaceBaseExt(micFn, 'mrc')) if downFactor != 1: # Replace extension by 'mrc' cause there are some formats # that cannot be written (such as dm3) ih.scaleFourier(micFn, micFnMrc, downFactor) sps = self._params['scannedPixelSize'] * downFactor kwargs['scannedPixelSize'] = sps else: ih.convert(micFn, micFnMrc, emlib.DT_FLOAT) program, args = self._gctfProgram.getCommand(**kwargs) args += ' %s/*.mrc' % micPath self.runJob(program, args) # , env=gctf.Plugin.getEnviron()) def _getFile(micBase, suffix): return os.path.join(micPath, micBase + suffix) for mic in micList: micFn = mic.getFileName() micBase = pwutils.removeBaseExt(micFn) micFnMrc = _getFile(micBase, '.mrc') # Let's clean the temporary mrc micrograph pwutils.cleanPath(micFnMrc) # move output from tmp to extra micFnCtf = _getFile(micBase, self._gctfProgram.getExt()) micFnCtfLog = _getFile(micBase, '_gctf.log') micFnCtfFit = _getFile(micBase, '_EPA.log') micFnCtfOut = self._getPsdPath(micFn) micFnCtfLogOut = self._getCtfOutPath(micFn) micFnCtfFitOut = self._getCtfFitOutPath(micFn) pwutils.moveFile(micFnCtf, micFnCtfOut) pwutils.moveFile(micFnCtfLog, micFnCtfLogOut) pwutils.moveFile(micFnCtfFit, micFnCtfFitOut) pwutils.cleanPath(micPath) except: print("ERROR: Gctf has failed on %s/*.mrc" % micPath) import traceback traceback.print_exc() def _reEstimateCTF(self, mic, ctf): """ Re-run gctf with required parameters """ self._estimateCtfList([mic], **self._getRecalCtfParamsDict(ctf)) def _createCtfModel(self, mic, updateSampling=True): # When downsample option is used, we need to update the # sampling rate of the micrograph associated with the CTF # since it could be downsampled if updateSampling: newSampling = mic.getSamplingRate() * self.ctfDownFactor.get() mic.setSamplingRate(newSampling) micFn = mic.getFileName() ctf = self._gctfProgram.parseOutputAsCtf( self._getCtfOutPath(micFn), psdFile=self._getPsdPath(micFn)) ctf.setMicrograph(mic) return ctf def _createOutputStep(self): pass # -------------------------- INFO functions ------------------------------- def _validate(self): errors = [] nprocs = max(self.numberOfMpi.get(), self.numberOfThreads.get()) if nprocs < len(self.getGpuList()): errors.append("Multiple GPUs can not be used by a single process. " "Make sure you specify more processors than GPUs. ") return errors def _methods(self): if self.inputMicrographs.get() is None: return ['Input micrographs not available yet.'] methods = "We calculated the CTF of " methods += self.getObjectTag('inputMicrographs') methods += " using Gctf [Zhang2016]. " methods += self.methodsVar.get('') if self.hasAttribute('outputCTF'): methods += 'Output CTFs: %s' % self.getObjectTag('outputCTF') return [methods] # -------------------------- UTILS functions ------------------------------ def _getRecalCtfParamsDict(self, ctfModel): values = [float(x) for x in ctfModel.getObjComment().split()] sampling = ctfModel.getMicrograph().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 _getPsdPath(self, micFn): micFnBase = pwutils.removeBaseExt(micFn) return self._getExtraPath(micFnBase + '_ctf.mrc') def _getCtfOutPath(self, micFn): micFnBase = pwutils.removeBaseExt(micFn) return self._getExtraPath(micFnBase + '_ctf.log') def _getCtfFitOutPath(self, micFn): micFnBase = pwutils.removeBaseExt(micFn) return self._getExtraPath(micFnBase + '_ctf_EPA.log') def _parseOutput(self, filename): """ Try to find the output estimation parameters from filename. It search for lines containing: Final Values and Resolution limit. """ return self._gctfProgram.parseOutput(filename) def _getCTFModel(self, defocusU, defocusV, defocusAngle, psdFile): ctf = CTFModel() ctf.setStandardDefocus(defocusU, defocusV, defocusAngle) ctf.setPsdFile(psdFile) return ctf