Source code for xmipp3.protocols.protocol_ctf_consensus

# **************************************************************************
# *
# * Authors:     J.M. De la Rosa Trevin (jmdelarosa@cnb.csic.es)
# *              Roberto Marabini (roberto@cnb.csic.es)
# *              Tomas Majtner (tmajtner@cnb.csic.es)  -- streaming version
# *              Amaya Jimenez (ajimenez@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
from datetime import datetime
from cmath import rect, phase
from math import radians, degrees

from pyworkflow import VERSION_3_0
from pwem.objects import SetOfCTF, SetOfMicrographs
from pyworkflow.object import Set, Integer, Pointer
import pyworkflow.protocol.params as params
import pyworkflow.utils as pwutils

from pwem.protocols import ProtCTFMicrographs
from pwem.emlib.metadata import Row
from pyworkflow.protocol.constants import (STATUS_NEW)
from pyworkflow import BETA, UPDATED, NEW, PROD

from pwem import emlib
import xmipp3
from xmipp3.convert import setXmippAttribute, getScipionObj

ACCEPTED = 'Accepted'
DISCARDED = 'Discarded'
INPUT1 = 1
INPUT2 = 2

[docs]class XmippProtCTFConsensus(ProtCTFMicrographs): """ Protocol to make a selection of meaningful CTFs in basis of the defocus values, the astigmatism, the resolution, other Xmipp parameters, and the agreement with a secondary CTF for the same set of micrographs. """ _label = 'ctf consensus' _devStatus = PROD _lastUpdateVersion = VERSION_3_0 def __init__(self, **args): ProtCTFMicrographs.__init__(self, **args) self._freqResol = {} self.stepsExecutionMode = params.STEPS_PARALLEL def _defineParams(self, form): form.addSection(label='Input') form.addParam('inputCTF', params.PointerParam, pointerClass='SetOfCTF', label="Input CTF", important=True, help='Select the estimated CTF to evaluate') form.addParam('useDefocus', params.BooleanParam, default=True, label='Use Defocus for selection', help='Use this button to decide if carry out the ' 'selection taking into account or not the defocus ' 'values.') line = form.addLine('Defocus (A)', condition="useDefocus", help='Minimum and maximum values for defocus in ' 'Angstroms.\nMicrographs out of this range ' 'will left out of the output.') line.addParam('minDefocus', params.FloatParam, default=4000, label='Min') line.addParam('maxDefocus', params.FloatParam, default=40000, label='Max') form.addParam('useAstigmatism', params.BooleanParam, default=False, label='Use Astigmatism for selection', help='Use this button to decide if carry out the ' 'selection taking into account or not the ' 'astigmatism value.') form.addParam('astigmatism', params.FloatParam, default=1000, label='Astigmatism (A)', condition="useAstigmatism", help='Maximum value allowed for astigmatism in ' 'Angstroms. If the evaluated CTF has a ' 'larger Astigmatism, it will be discarded.') form.addParam('useAstigmatismPercentage', params.BooleanParam, default=True, label='Use Astigmatism percentage for selection', help='Use this button to decide if carry out the ' 'selection taking into account or not the ' 'astigmatism value.') form.addParam('astigmatismPer', params.FloatParam, default=0.10, label='Astigmatism percentage', condition="useAstigmatismPercentage", help='Maximum value allowed for astigmatism ' 'percentage (|defocus_U-defocus_V|/mean_defocus). If the evaluated CTF has a ' 'larger Astigmatism Percentage, it will be discarded.') form.addParam('useResolution', params.BooleanParam, default=True, label='Use Resolution for selection', help='Use this button to decide if carry out the ' 'selection taking into account or not the ' 'resolution value.') form.addParam('resolution', params.FloatParam, default=17, label='Resolution (A)', condition="useResolution", help='Minimum value for resolution in Angstroms. ' 'If the evaluated CTF has not reached that minimum, ' 'it will be discarded.') form.addSection(label='Xmipp criteria') form.addParam('useCritXmipp', params.BooleanParam, default=False, label='Use Xmipp criteria for selection', help='Use this button to decide if carrying out the ' 'selection taking into account the Xmipp parameters.\n' 'Only available when Xmipp CTF estimation was used ' 'for the _Input CTF_ or for the _Secondary CTF_.') form.addParam('critFirstZero', params.FloatParam, default=5, condition="useCritXmipp", label='Minimum 1st zero', help='Minimun value of CritFirstZero') line = form.addLine('First zero astigmatism', condition="useCritXmipp", help='Minimum and maximum values for ' 'CritFirstZeroRatio') line.addParam('minCritFirstZeroRatio', params.FloatParam, default=0.9, label='Min') line.addParam('maxCritFirstZeroRatio', params.FloatParam, default=1.1, label='Max') form.addParam('critCorr', params.FloatParam, default=0.05, condition="useCritXmipp", label='Correlation 1st-3rd zero', help='Minimum value of correlation between 1st and 3rd zeros') form.addParam('critCtfMargin', params.FloatParam, default=1.5, condition="useCritXmipp", label='CTF Margin', help='Minimum value of CritCtfMargin') form.addParam('critIceness', params.FloatParam, default=1, condition="useCritXmipp", label='CritIceness', help='Minimum value of the iceness.') line = form.addLine('Non astigmatic validity', condition="useCritXmipp", help='Minimum and maximum values for ' 'CritNonAstigmaticValidity') line.addParam('minCritNonAstigmaticValidity', params.FloatParam, default=0.3, label='Min') line.addParam('maxCritNonAstigmaticValidity', params.FloatParam, default=9, label='Max') form.addSection(label='Consensus') form.addParam('calculateConsensus', params.BooleanParam, default=False, label='Calculate Consensus Resolution', help='Option for calculating consensus resolution. ' 'The algorithm assumes that two CTF are ' 'consistent if the phase (wave aberration function) ' 'of the two CTFs are closer than 90 degrees.\n' 'The reported consensusResolution is the resolution ' 'at which the two CTF phases differ in 90 degrees.') form.addParam('inputCTF2', params.PointerParam, pointerClass='SetOfCTF', condition="calculateConsensus", label="Secondary CTF", help='CTF to be compared with reference CTF') form.addParam('minConsResol', params.FloatParam, condition="calculateConsensus", default=15.0, label='Minimum consensus resolution (A).', help="Minimum value for the consensus resolution in " "Angstroms.\nIf there are noticeable discrepancies " "between the two estimations below this resolution, " "it will be discarded.") form.addParam('averageDefocus', params.BooleanParam, condition="calculateConsensus", default=False, label='Average equivalent metadata?', help='If *Yes*, making an average of those metadata present ' 'in both CTF estimations (defocus, astigmatism angle...)\n ' 'If *No*, the primary estimation metadata will persist.') form.addParam('includeSecondary', params.BooleanParam, condition="calculateConsensus", default=False, label='Include all secondary metadata?', help='If *Yes*, all metadata in the *Secondary CTF* will ' 'be included in the resulting CTF.\n ' 'If *No*, only the primary metadata (plus consensus ' 'scores) will be in the resulting CTF.') form.addParallelSection(threads=4, mpi=1) # --------------------------- INSERT steps functions ------------------------- def _insertAllSteps(self): self.initializeParams() if self.calculateConsensus: self.ctfFn2 = self.inputCTF2.get().getFileName() self.allCtf2 = {} self._insertFunctionStep('createOutputStep', prerequisites=[], wait=True)
[docs] def createOutputStep(self): self._closeOutputSet()
[docs] def initializeParams(self): self.finished = False self.insertedDict = {} self.initializeRejDict() self.setSecondaryAttributes() self.ctfFn1 = self.inputCTF.get().getFileName() self.allCtf1 = {} pwutils.makePath(self._getExtraPath('DONE'))
def _getFirstJoinStepName(self): # This function will be used for streaming, to check which is # the first function that need to wait for all ctfs # to have completed, this can be overriden in subclasses # (e.g., in Xmipp 'sortPSDStep') return 'createOutputStep' def _getFirstJoinStep(self): for s in self._steps: if s.funcName == self._getFirstJoinStepName(): return s return None def _insertNewCtfsSteps(self, newIDs1, newIDs2, insertedDict): deps = [] newIDs = list(set(newIDs1).intersection(set(newIDs2))) md1 = emlib.MetaData() md2 = emlib.MetaData() for ctfID in newIDs: if ctfID not in insertedDict: ctf1 = self.allCtf1.get(ctfID) ctf2 = self.allCtf2.get(ctfID) try: self._ctfToMd(ctf1, md1) self._ctfToMd(ctf2, md2) self._freqResol[ctfID] = emlib.errorMaxFreqCTFs2D(md1, md2) except TypeError as exc: print("Error reading ctf for id:%s. %s" % (ctfID, exc)) self._freqResol[ctfID] = 9999 stepId = self._insertFunctionStep('selectCtfStep', ctfID, prerequisites=[]) deps.append(stepId) insertedDict[ctfID] = stepId return deps def _insertNewSelectionSteps(self, insertedDict, newIDs): deps = [] # For each ctf insert the step to process it for ctfID in newIDs: if ctfID not in insertedDict: stepId = self._insertFunctionStep('selectCtfStep', ctfID, prerequisites=[]) deps.append(stepId) insertedDict[ctfID] = stepId return deps def _stepsCheck(self): self._checkNewInput() self._checkNewOutput() def _checkNewInput(self): if self.calculateConsensus: # Check if there are new ctf to process from the input set self.lastCheck = getattr(self, 'lastCheck', datetime.now()) mTime = max(datetime.fromtimestamp(os.path.getmtime(self.ctfFn1)), datetime.fromtimestamp(os.path.getmtime(self.ctfFn2))) self.debug('Last check: %s, modification: %s' % (pwutils.prettyTime(self.lastCheck), pwutils.prettyTime(mTime))) # If the input movies.sqlite have not changed since our last check, # it does not make sense to check for new input data if self.lastCheck > mTime and (hasattr(self, 'outputCTF') or hasattr(self, "outputCTFDiscarded")): return None ctfsSet1 = self._loadInputCtfSet(self.ctfFn1) ctfsSet2 = self._loadInputCtfSet(self.ctfFn2) ctfDict1 = {ctf.getObjId(): ctf.clone() for ctf in ctfsSet1.iterItems()} ctfDict2 = {ctf.getObjId(): ctf.clone() for ctf in ctfsSet2.iterItems()} newIds1 = [idCTF1 for idCTF1 in ctfDict1.keys() if idCTF1 not in self.insertedDict] self.allCtf1.update(ctfDict1) newIds2 = [idCTF2 for idCTF2 in ctfDict2.keys() if idCTF2 not in self.insertedDict] self.allCtf2.update(ctfDict2) self.lastCheck = datetime.now() self.isStreamClosed = ctfsSet1.isStreamClosed() and \ ctfsSet2.isStreamClosed() ctfsSet1.close() ctfsSet2.close() outputStep = self._getFirstJoinStep() if len(set(self.allCtf1)) > len(set(self.insertedDict)) and \ len(set(self.allCtf2)) > len(set(self.insertedDict)): fDeps = self._insertNewCtfsSteps(newIds1, newIds2, self.insertedDict) if outputStep is not None: outputStep.addPrerequisites(*fDeps) self.updateSteps() else: now = datetime.now() self.lastCheck = getattr(self, 'lastCheck', now) mTime = datetime.fromtimestamp(os.path.getmtime(self.ctfFn1)) self.debug('Last check: %s, modification: %s' % (pwutils.prettyTime(self.lastCheck), pwutils.prettyTime(mTime))) # If the input ctfs.sqlite have not changed since our last check, # it does not make sense to check for new input data if self.lastCheck > mTime and (hasattr(self, 'outputCTF') or hasattr(self, "outputCTFDiscarded")): return None # Open input ctfs.sqlite and close it as soon as possible ctfSet = self._loadInputCtfSet(self.ctfFn1) ctfDict = {ctf.getObjId(): ctf.clone() for ctf in ctfSet.iterItems()} newIds = [idCTF for idCTF in ctfDict.keys() if idCTF not in self.insertedDict] self.allCtf1.update(ctfDict) self.isStreamClosed = ctfSet.isStreamClosed() ctfSet.close() self.lastCheck = now newCtf = any(ctf.getObjId() not in self.insertedDict for ctf in self.allCtf1.values()) outputStep = self._getFirstJoinStep() if newCtf: fDeps = self._insertNewSelectionSteps(self.insertedDict, newIds) if outputStep is not None: outputStep.addPrerequisites(*fDeps) self.updateSteps() def _checkNewOutput(self): """ Check for already selected CTF and update the output set. """ # Load previously done items (from text file) doneListDiscarded = self._readCertainDoneList(DISCARDED) doneListAccepted = self._readCertainDoneList(ACCEPTED) # Check for newly done items ctfListIdAccepted = self._readtCtfId(True) ctfListIdDiscarded = self._readtCtfId(False) newDoneAccepted = [ctfId for ctfId in ctfListIdAccepted if ctfId not in doneListAccepted] newDoneDiscarded = [ctfId for ctfId in ctfListIdDiscarded if ctfId not in doneListDiscarded] firstTimeAccepted = len(doneListAccepted) == 0 firstTimeDiscarded = len(doneListDiscarded) == 0 allDone = len(doneListAccepted) + len(doneListDiscarded) +\ len(newDoneAccepted) + len(newDoneDiscarded) # We have finished when there is not more input ctf (stream closed) # and the number of processed ctf is equal to the number of inputs if self.calculateConsensus: maxCtfSize = min(len(self.allCtf1), len(self.allCtf2)) else: maxCtfSize = len(self.allCtf1) self.finished = (self.isStreamClosed and allDone == maxCtfSize) streamMode = Set.STREAM_CLOSED if self.finished else Set.STREAM_OPEN def readOrCreateOutputs(doneList, newDone, label=''): if len(doneList) > 0 or len(newDone) > 0: cSet = self._loadOutputSet(SetOfCTF, 'ctfs'+label+'.sqlite') mSet = self._loadOutputSet(SetOfMicrographs, 'micrographs'+label+'.sqlite') label = ACCEPTED if label == '' else DISCARDED self.fillOutput(cSet, mSet, newDone, label) return cSet, mSet return None, None ctfSet, micSet = readOrCreateOutputs(doneListAccepted, newDoneAccepted) ctfSetDiscarded, micSetDiscarded = readOrCreateOutputs(doneListDiscarded, newDoneDiscarded, DISCARDED) if not self.finished and not newDoneDiscarded and not newDoneAccepted: # If we are not finished and no new output have been produced # it does not make sense to proceed and updated the outputs # so we exit from the function here return def updateRelationsAndClose(cSet, mSet, first, label=''): if os.path.exists(self._getPath('ctfs'+label+'.sqlite')): micsAttrName = 'outputMicrographs'+label self._updateOutputSet(micsAttrName, mSet, streamMode) # Set micrograph as pointer to protocol to prevent pointee end up as another attribute (String, Booelan,...) # that happens somewhere while scheduling. cSet.setMicrographs(Pointer(self, extended=micsAttrName)) self._updateOutputSet('outputCTF'+label, cSet, streamMode) if first: self._defineTransformRelation(self.inputCTF.get().getMicrographs(), mSet) # self._defineTransformRelation(cSet, mSet) self._defineTransformRelation(self.inputCTF, cSet) self._defineCtfRelation(mSet, cSet) mSet.close() cSet.close() updateRelationsAndClose(ctfSet, micSet, firstTimeAccepted) updateRelationsAndClose(ctfSetDiscarded, micSetDiscarded, firstTimeDiscarded, DISCARDED) if self.finished: # Unlock createOutputStep if finished all jobs outputStep = self._getFirstJoinStep() if outputStep and outputStep.isWaiting(): outputStep.setStatus(STATUS_NEW)
[docs] def fillOutput(self, ctfSet, micSet, newDone, label): if newDone: inputCtfSet = self._loadInputCtfSet(self.ctfFn1) if self.calculateConsensus: inputCtfSet2 = self._loadInputCtfSet(self.ctfFn2) for ctfId in newDone: ctf = inputCtfSet[ctfId].clone() mic = ctf.getMicrograph().clone() ctf.setEnabled(self._getEnable(ctfId)) mic.setEnabled(self._getEnable(ctfId)) if self.calculateConsensus: ctf2 = inputCtfSet2[ctfId] conRes = self._freqResol[ctfId] setAttribute(ctf, '_consensus_resolution', conRes) setAttribute(ctf, '_ctf2_defocus_diff', max(abs(ctf.getDefocusU()-ctf2.getDefocusU()), abs(ctf.getDefocusV()-ctf2.getDefocusV()))) setAttribute(ctf, '_ctf2_defocusAngle_diff', anglesDifference(ctf.getDefocusAngle(), ctf2.getDefocusAngle())) if ctf.hasPhaseShift() and ctf2.hasPhaseShift(): setAttribute(ctf, '_ctf2_phaseShift_diff', anglesDifference(ctf.getPhaseShift(), ctf2.getPhaseShift())) setAttribute(ctf, '_ctf2_resolution', ctf2.getResolution()) setAttribute(ctf, '_ctf2_fitQuality', ctf2.getFitQuality()) if ctf2.hasAttribute('_xmipp_ctfmodel_quadrant'): # To check CTF in Xmipp _quadrant is the best copyAttribute(ctf2, ctf, '_xmipp_ctfmodel_quadrant') else: setAttribute(ctf, '_ctf2_psdFile', ctf2.getPsdFile()) if self.averageDefocus: newDefocusU = 0.5*(ctf.getDefocusU() + ctf2.getDefocusU()) newDefocusV = 0.5*(ctf.getDefocusV() + ctf2.getDefocusV()) newDefocusAngle = averageAngles(ctf.getDefocusAngle(), ctf2.getDefocusAngle()) ctf.setStandardDefocus(newDefocusU, newDefocusV, newDefocusAngle) if ctf.hasPhaseShift() and ctf2.hasPhaseShift(): newPhaseShift = averageAngles(ctf.getPhaseShift(), ctf2.getPhaseShift()) ctf.setPhaseShift(newPhaseShift) else: setAttribute(ctf, '_ctf2_defocusRatio', ctf2.getDefocusRatio()) setAttribute(ctf, '_ctf2_astigmatism', abs(ctf2.getDefocusU() - ctf2.getDefocusV())) if self.includeSecondary: for attr in self.secondaryAttributes: copyAttribute(ctf2, ctf, attr) # main _astigmatism always but after consensus if so setAttribute(ctf, '_astigmatism', abs(ctf.getDefocusU() - ctf.getDefocusV())) # percentage _astigmatism always but after consensus if so astigmatismPer = abs(ctf.getDefocusU() - ctf.getDefocusV())/(0.5 * (ctf.getDefocusU() + ctf.getDefocusV())) setAttribute(ctf, '_astigmatismPercentage', astigmatismPer) ctfSet.append(ctf) micSet.append(mic) self._writeCertainDoneList(ctfId, label) inputCtfSet.close() if self.calculateConsensus: inputCtfSet2.close()
[docs] def setSecondaryAttributes(self): if self.calculateConsensus and self.includeSecondary: item = self.inputCTF.get().getFirstItem() ctf1Attr = set(item.getObjDict().keys()) item = self.inputCTF2.get().getFirstItem() ctf2Attr = set(item.getObjDict().keys()) self.secondaryAttributes = ctf2Attr - ctf1Attr else: self.secondaryAttributes = set()
def _loadOutputSet(self, SetClass, baseName): """ Load the output set if it exists or create a new one. """ setFile = self._getPath(baseName) if os.path.exists(setFile): outputSet = SetClass(filename=setFile) if (outputSet.__len__() == 0): pwutils.path.cleanPath(setFile) if os.path.exists(setFile): outputSet = SetClass(filename=setFile) outputSet.loadAllProperties() outputSet.enableAppend() else: outputSet = SetClass(filename=setFile) outputSet.setStreamState(outputSet.STREAM_OPEN) micSet = self.inputCTF.get().getMicrographs() if isinstance(outputSet, SetOfMicrographs): outputSet.copyInfo(micSet) elif isinstance(outputSet, SetOfCTF): outputSet.setMicrographs(micSet) return outputSet def _ctfToMd(self, ctf, ctfMd): """ Write the proper metadata for Xmipp from a given CTF """ ctfMd.clear() ctfRow = Row() xmipp3.convert.ctfModelToRow(ctf, ctfRow) xmipp3.convert.micrographToRow(ctf.getMicrograph(), ctfRow, alignType=xmipp3.convert.ALIGN_NONE) ctfRow.addToMd(ctfMd)
[docs] def initializeRejDict(self): self.discDict = {'defocus': 0, 'astigmatism': 0, 'astigmatismPer': 0, 'singleResolution': 0, '_xmipp_ctfCritFirstZero': 0, '_xmipp_ctfCritfirstZeroRatio': 0, '_xmipp_ctfCritCorr13': 0, '_xmipp_ctfCritIceness': 0, '_xmipp_ctfCritCtfMargin': 0, '_xmipp_ctfCritNonAstigmaticValidty': 0, 'consensusResolution': 0 } for k in self.discDict: setattr(self, "rejBy"+k, Integer(0)) self._store()
[docs] def selectCtfStep(self, ctfId): # Depending on the flags selected by the user, we set the values of # the params to compare with doneFn = self._getCTFDone(ctfId) if self.isContinued() and self._isCTFDone(ctfId): self.info("Skipping CTF with ID: %s, seems to be done" % ctfId) return # Clean old finished files pwutils.cleanPath(doneFn) def compareValue(ctf, label, comp, crit): """ Returns True if the ctf.label NOT complain the crit by comp """ if hasattr(ctf, label): if comp == 'lt': discard = getattr(ctf, label).get() < crit elif comp == 'bt': discard = getattr(ctf, label).get() > crit else: raise Exception("'comp' must be either 'lt' or 'bt'.") else: print("%s not found. Skipping evaluation on that." % label) return False if discard: self.discDict[label] += 1 return discard minDef, maxDef = self._getDefociValues() maxAstig = self._getMaxAstisgmatism() maxAstigPer = self._getMaxAstigmatismPer() minResol = self._getMinResol() ctf = self.allCtf1.get(ctfId) defocusU = ctf.getDefocusU() defocusV = ctf.getDefocusV() astigm = abs(defocusU - defocusV) astigmPer = abs(defocusU - defocusV)/((defocusU+defocusV)/2) resol = self._getCtfResol(ctf) defRangeCrit = (defocusU < minDef or defocusU > maxDef or defocusV < minDef or defocusV > maxDef) if defRangeCrit: self.discDict['defocus'] += 1 astigCrit = astigm > maxAstig if astigCrit: self.discDict['astigmatism'] += 1 astigPer = (astigmPer > maxAstigPer) if astigPer: self.discDict['astigmatismPer'] += 1 singleResolCrit = resol > minResol if singleResolCrit: self.discDict['singleResolution'] += 1 firstCondition = defRangeCrit or astigCrit or singleResolCrit or astigPer consResolCrit = False if self.calculateConsensus: # FIXME: this conditional is to avoid an error, but it shouldn't happen!!! if ctfId in self._freqResol: consResolCrit = self.minConsResol < self._freqResol[ctfId] if consResolCrit: self.discDict['consensusResolution'] += 1 else: consResolCrit = True self.discDict['consensusResolution'] += 1 self._freqResol[ctfId] = 9999 secondCondition = False if self.useCritXmipp: firstZero = self._getCritFirstZero() minFirstZero, maxFirstZero = self._getCritFirstZeroRatio() corr = self._getCritCorr() iceness = self._getIceness() ctfMargin = self._getCritCtfMargin() minNonAstigmatic, maxNonAstigmatic = \ self._getCritNonAstigmaticValidity() if self.xmippCTF == INPUT1: ctfX = ctf else: ctfX = self.allCtf2.get(ctfId) secondCondition = ( compareValue(ctfX, '_xmipp_ctfCritFirstZero', 'lt', firstZero) or compareValue(ctfX, '_xmipp_ctfCritfirstZeroRatio', 'lt', minFirstZero) or compareValue(ctfX, '_xmipp_ctfCritfirstZeroRatio', 'bt', maxFirstZero) or compareValue(ctfX, '_xmipp_ctfCritCorr13', 'lt', corr) or compareValue(ctfX, '_xmipp_ctfCritIceness', 'bt', iceness) or compareValue(ctfX, '_xmipp_ctfCritCtfMargin', 'lt', ctfMargin) or compareValue(ctfX, '_xmipp_ctfCritNonAstigmaticValidty', 'lt', minNonAstigmatic) or compareValue(ctfX, '_xmipp_ctfCritNonAstigmaticValidty', 'bt', maxNonAstigmatic)) """ Write to a text file the items that have been done. """ if firstCondition or consResolCrit or secondCondition: fn = self._getCtfSelecFileDiscarded() with open(fn, 'a') as f: f.write('%d F\n' % ctf.getObjId()) else: if (ctf.isEnabled()): fn = self._getCtfSelecFileAccepted() with open(fn, 'a') as f: f.write('%d T\n' % ctf.getObjId()) else: fn = self._getCtfSelecFileAccepted() with open(fn, 'a') as f: f.write('%d F\n' % ctf.getObjId()) for k, v in self.discDict.items(): setattr(self, "rejBy"+k, Integer(v)) self._store() # Mark this ctf as finished open(doneFn, 'w').close()
def _readDoneList(self): """ Read from a file the id's of the items that have been done. """ doneFile = self._getAllDone() doneList = [] # Check what items have been previously done if os.path.exists(doneFile): with open(doneFile) as f: doneList += [int(line.strip()) for line in f] return doneList def _getAllDone(self): return self._getExtraPath('DONE_all.TXT') def _writeDoneList(self, idList): """ Write to a text file the items that have been done. """ with open(self._getAllDone(), 'a') as f: for id in idList: f.write('%d\n' % id) def _isCTFDone(self, id): """ A mic is done if the marker file exists. """ return os.path.exists(self._getCTFDone(id)) def _getCTFDone(self, id): """ Return the file that is used as a flag of termination. """ return self._getExtraPath('DONE', 'ctf_%06d.TXT' % id) def _citations(self): return ['Marabini2014a'] def _summary(self): if not (hasattr(self, "outputCTF") or hasattr(self, "outputCTFDiscarded")): return ['No CTF processed, yet.'] acceptedSize = (self.outputCTF.getSize() if hasattr(self, "outputCTF") else 0) discardedSize = (self.outputCTFDiscarded.getSize() if hasattr(self, "outputCTFDiscarded") else 0) message = ["%d/%d CTF processed (%d accepted and %d discarded)." % (acceptedSize+discardedSize, self.inputCTF.get().getSize(), acceptedSize, discardedSize)] def addDiscardedStr(label): obj = getattr(self, "rejBy%s" % label, Integer(0)) number = obj.get() return "" if number == 0 else " (%d discarded)" % number if any([self.useDefocus, self.useAstigmatism, self.useResolution, self.useAstigmatismPercentage]): message.append("*General Criteria*:") if self.useDefocus: message.append(" - _Defocus_. Range: %.0f - %.0f %s" % (self.minDefocus, self.maxDefocus, addDiscardedStr('defocus'))) if self.useAstigmatism: message.append(" - _Astigmatism_. Threshold: %.0f %s" % (self.astigmatism, addDiscardedStr('astigmatism'))) if self.useAstigmatismPercentage: message.append(" - _AstigmatismPer_. Threshold: %.3f %s" % (self.astigmatismPer, addDiscardedStr('astigmatismPer'))) if self.useResolution: message.append(" - _Resolution_. Threshold: %.0f %s" % (self.resolution, addDiscardedStr('singleResolution'))) if self.useCritXmipp: message.append("*Xmipp criteria*:") message.append(" - _First zero_. Threshold: %.0f %s" % (self.critFirstZero, addDiscardedStr('_xmipp_ctfCritFirstZero'))) message.append(" - _First zero astigmatism_. Range: %.2f - %.2f %s" % (self.minCritFirstZeroRatio, self.maxCritFirstZeroRatio, addDiscardedStr('_xmipp_ctfCritfirstZeroRatio'))) message.append(" - _Correlation Experimental-Estimated_. Threshold: " "%.2f %s" % (self.critCorr, addDiscardedStr('_xmipp_ctfCritCorr13'))) message.append(" - _CTF margin_. Threshold: %.2f %s" % (self.critCtfMargin, addDiscardedStr('_xmipp_ctfCritCtfMargin'))) message.append(" - _Iceness_. Threshold: %.2f %s" % (self.critIceness, addDiscardedStr('_xmipp_ctfCritIceness'))) message.append(" - _Non Astigmatic validation range_: %.2f - %.2f %s" % (self.minCritNonAstigmaticValidity, self.maxCritNonAstigmaticValidity, addDiscardedStr('_xmipp_ctfCritNonAstigmaticValidty'))) if self.calculateConsensus: def getProtocolInfo(inCtf): protocol = self.getMapper().getParent(inCtf) runName = protocol.getRunName() classLabel = protocol.getClassLabel() if runName == classLabel: infoStr = runName else: infoStr = "%s (%s)" % (runName, classLabel) return infoStr message.append("*CTF consensus*:") message.append(" - _Consensus resolution. Threshold_: %.0f %s" % (self.minConsResol, addDiscardedStr('consensusResolution'))) message.append(" > _Primary CTF_: %s" % getProtocolInfo(self.inputCTF.get())) message.append(" > _Reference CTF_: %s" % getProtocolInfo(self.inputCTF2.get())) return message def _validate(self): """ The function of this hook is to add some validation before the protocol is launched to be executed. It should return a list of errors. If the list is empty the protocol can be executed. """ # same micrographs in both CTF?? errors = [] if self.useCritXmipp.get() and not self.calculateConsensus.get(): if self.usingXmipp(self.inputCTF.get().getFirstItem()): self.xmippCTF = INPUT1 else: errors.append("The primary CTF input ( _Input CTF_ ) must be " "estimated using the _Xmipp - CTF estimation_ " "protocol.") if self.useCritXmipp.get() and self.calculateConsensus.get(): if self.usingXmipp(self.inputCTF.get().getFirstItem()): self.xmippCTF = INPUT1 elif self.usingXmipp(self.inputCTF2.get().getFirstItem()): self.xmippCTF = INPUT2 else: errors.append("One of the CTF inputs ( _Input CTF_ or " "_Secundary CTF_) must be estimated using the " "_Xmipp - CTF estimation_ protocol.") return errors
[docs] def usingXmipp(self, ctf): return ctf.hasAttribute('_xmipp_ctfCritFirstZero')
def _getCtfResol(self, ctf): resolution = ctf.getResolution() if resolution is not None: return resolution else: return 0 def _readCertainDoneList(self, label): """ Read from a text file the id's of the items that have been done. """ doneFile = self._getCertainDone(label) doneList = [] # Check what items have been previously done if os.path.exists(doneFile): with open(doneFile) as f: doneList += [int(line.strip()) for line in f] return doneList def _writeCertainDoneList(self, ctfId, label): """ Write to a text file the items that have been done. """ doneFile = self._getCertainDone(label) with open(doneFile, 'a') as f: f.write('%d\n' % ctfId) def _getCertainDone(self, label): return self._getExtraPath('DONE_'+label+'.TXT') def _getCtfSelecFileAccepted(self): return self._getExtraPath('selection-ctf-accepted.txt') def _getCtfSelecFileDiscarded(self): return self._getExtraPath('selection-ctf-discarded.txt') def _readtCtfId(self, accepted): if accepted: fn = self._getCtfSelecFileAccepted() else: fn = self._getCtfSelecFileDiscarded() ctfList = [] # Check what items have been previously done if os.path.exists(fn): with open(fn) as f: ctfList += [int(line.strip().split()[0]) for line in f] return ctfList def _getEnable(self, ctfId): fn = self._getCtfSelecFileAccepted() # Check what items have been previously done if os.path.exists(fn): with open(fn) as f: for line in f: if ctfId == int(line.strip().split()[0]): if line.strip().split()[1] == 'T': return True else: return False def _loadInputCtfSet(self, ctfFn): self.debug("Loading input db: %s" % ctfFn) ctfSet = SetOfCTF(filename=ctfFn) ctfSet.loadAllProperties() return ctfSet def _getDefociValues(self): if not self.useDefocus: return 0, 1000000 else: return self.minDefocus.get(), self.maxDefocus.get() def _getMaxAstigmatismPer(self): if not self.useAstigmatismPercentage: return 0.1 else: return self.astigmatismPer.get() def _getMaxAstisgmatism(self): if not self.useAstigmatism: return 10000000 else: return self.astigmatism.get() def _getMinResol(self): if not self.useResolution: return 1000000 else: return self.resolution.get() def _getIceness(self): return self.critIceness.get() def _getCritFirstZero(self): return self.critFirstZero.get() def _getCritFirstZeroRatio(self): return self.minCritFirstZeroRatio.get(),\ self.maxCritFirstZeroRatio.get() def _getCritCorr(self): return self.critCorr.get() def _getCritCtfMargin(self): return self.critCtfMargin.get() def _getCritNonAstigmaticValidity(self): return (self.minCritNonAstigmaticValidity.get(), self.maxCritNonAstigmaticValidity.get())
[docs]def averageAngles(angle1, angle2): c1 = rect(1, radians(angle1*2)) c2 = rect(1, radians(angle2*2)) return degrees(phase((c1 + c2)*0.5))/2
[docs]def anglesDifference(angle1, angle2): if (angle1 > angle2) == (abs(angle2 - angle1) > 90): aux = angle1 angle1 = angle2 angle2 = aux return (angle1 - angle2) % 180
[docs]def setAttribute(obj, label, value): if value is None: return setattr(obj, label, getScipionObj(value))
[docs]def copyAttribute(src, dst, label, default=None): setAttribute(dst, label, getattr(src, label, default))