# *
# * Authors: Scipion Team
# *
# * 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-users@lists.sourceforge.net'
# *
# **************************************************************************
"""
This module contains the protocol for 3d refinement with Relion.
"""
import glob
from os import remove
from os.path import abspath, exists, getmtime
from pwem import ALIGN_PROJ
from pwem.objects import FSC, Integer
from pwem.protocols import ProtRefine3D
from pyworkflow import BETA
from relion.convert import Table, convert31
from reliontomo import Plugin
from tomo.objects import AverageSubTomogram
from tomo.protocols import ProtTomoBase
from reliontomo.protocols.protocol_base_tomo import ProtRelionBaseTomo
[docs]class ProtRelionSubtomoRefine3D(ProtRefine3D, ProtRelionBaseTomo, ProtTomoBase):
"""Protocol to refine a 3D map using Relion. Relion employs an empirical
Bayesian approach to refinement of (multiple) 3D reconstructions
or 2D class averages in electron cryo-microscopy (cryo-EM). Many
parameters of a statistical model are learned from the data,which
leads to objective and high-quality results.
"""
_label = '3D subtomogram auto-refine'
_devStatus = BETA
IS_CLASSIFY = False
PREFIXES = ['half1_', 'half2_']
def __init__(self, **args):
ProtRelionBaseTomo.__init__(self, **args)
def _initialize(self):
""" This function is mean to be called after the
working dir for the protocol have been set. (maybe after recovery from mapper)
"""
ProtRelionBaseTomo._initialize(self)
self.ClassFnTemplate = '%(ref)03d@%(rootDir)s/relion_it%(iter)03d_classes.mrcs'
self.numberOfClasses.set(1) # For refinement we only need just one "class"
# --------------------------- INSERT steps functions --------------------------------------------
def _setSamplingArgs(self, args):
""" Set sampling related params"""
# Sampling stuff
args['--auto_local_healpix_order'] = self.localSearchAutoSamplingDeg.get()
if not self.doContinue:
args['--healpix_order'] = self.angularSamplingDeg.get()
args['--offset_range'] = self.offsetSearchRangePix.get()
args['--offset_step'] = self.offsetSearchStepPix.get() * 2
args['--auto_refine'] = ''
args['--split_random_halves'] = ''
joinHalves = "--low_resol_join_halves"
if not joinHalves in self.extraParams.get():
args['--low_resol_join_halves'] = 40
if not Plugin.IS_30() and self.useFinerSamplingFaster:
args['--auto_ignore_angles'] = ''
args['--auto_resol_angles'] = ''
# --------------------------- STEPS functions --------------------------------------------
[docs] def createOutputStep(self):
subtomoSet = self._getInputParticles()
vol = AverageSubTomogram()
vol.setFileName(self._getExtraPath('relion_class001.mrc'))
vol.setSamplingRate(subtomoSet.getSamplingRate())
pattern = '*it*half%s_class*.mrc'
half1 = self._getLastFileName(self._getExtraPath(pattern % 1))
half2 = self._getLastFileName(self._getExtraPath(pattern % 2))
vol.setHalfMaps([half1, half2])
outSubtomoSet = self._createSetOfSubTomograms()
outSubtomoSet.copyInfo(subtomoSet)
self._fillDataFromIter(outSubtomoSet, self._lastIter())
self._defineOutputs(outputVolume=vol)
self._defineSourceRelation(subtomoSet, vol)
self._defineOutputs(outputParticles=outSubtomoSet)
self._defineTransformRelation(subtomoSet, outSubtomoSet)
fsc = FSC(objLabel=self.getRunName())
fn = self._getExtraPath("relion_model.star")
table = Table(fileName=fn, tableName='model_class_1')
resolution_inv = table.getColumnValues('rlnResolution')
frc = table.getColumnValues('rlnGoldStandardFsc')
fsc.setData(resolution_inv, frc)
self._defineOutputs(outputFSC=fsc)
self._defineSourceRelation(vol, fsc)
# --------------------------- INFO functions --------------------------------------------
def _validateNormal(self):
""" Should be overriden in subclasses to
return summary message for NORMAL EXECUTION.
"""
errors = []
if self.referenceVolume.get() is not None:
particlesDim = self._getInputParticles().getDim()
volumeDim = self.referenceVolume.get().getDim()
if particlesDim is None:
errors.append('Can not get dimensions from input particles!!!')
elif volumeDim is None:
errors.append('Can not get dimensions from reference volume!!!')
elif particlesDim[0] != volumeDim[0]:
errors.append('Volume and particles dimensions must be equal!!!')
return errors
def _validateContinue(self):
""" Should be overriden in subclasses to
return summary messages for CONTINUE EXECUTION.
"""
errors = []
continueRun = self.continueRun.get()
continueRun._initialize()
lastIter = continueRun._lastIter()
if self.continueIter.get() == 'last':
continueIter = lastIter
else:
continueIter = int(self.continueIter.get())
if continueIter > lastIter:
errors += ["The iteration from you want to continue must be %01d or less" % lastIter]
return errors
def _summaryNormal(self):
""" Should be overriden in subclasses to
return summary message for NORMAL EXECUTION.
"""
summary = []
it = self._lastIter()
if it:
if it >= 1:
modelStar = self._getFileName('half1_model', iter=it)
summary.append("Current resolution %s: *%0.2f* Å" %
self._getRlnCurrentResolution(modelStar))
return summary
def _summaryContinue(self):
""" Should be overriden in subclasses to
return summary messages for CONTINUE EXECUTION.
"""
summary = ["Continue from iteration %01d" % self._getContinueIter()]
return summary
# --------------------------- UTILS functions --------------------------------------------
def _createItemMatrix(self, item, row):
self.reader.setParticleTransform(item, row)
def _cleanUndesiredFiles(self):
"""Remove all files generated by relion_refine 3d excepting the ones which
correspond to the last iteration. Example for iteration 12:
relion_it012_half1_class001.mrc
relion_it012_half2_class001_angdist.bild
relion_it012_half2_class001.mrc
relion_it012_half2_model.star
relion_it012_optimiser.star
relion_it012_half1_class001_angdist.bild
relion_it012_data.star
relion_it012_sampling.star
relion_it012_half1_model.star
"""
itPref = 'relion_it'
clPref = 'class'
starExt = '.star'
mrcExt = '.mrc'
bildExt = '_angdist.bild'
halfPref = 'half'
# Refine calculations related files
calcFiles = ['data', '%s1_model' % halfPref, '%s2_model' % halfPref, 'optimiser', 'sampling']
for i in range(self._lastIter()):
for calcFile in calcFiles:
fn = abspath(self._getExtraPath('{}{:03d}_{}{}'.format(
itPref, i, calcFile, starExt)))
if exists(fn):
remove(fn)
# Classes related files
for itr in range(1, self.nClasses + 1):
for h in range(1, 3): # Two halves
fn = abspath(self._getExtraPath('{}{:03d}_{}{}_{}{:03d}'.format(
itPref, i, halfPref, h, clPref, itr)))
mrcFile = fn + mrcExt
bildFile = fn + bildExt
if exists(mrcFile):
remove(mrcFile)
if exists(bildFile):
remove(bildFile)
def _fillDataFromIter(self, subtomoClassesSet, iteration):
# tableName = 'particles@' if self.IS_GT30() else ''
dataStar = self._getFileName('data', iter=iteration)
# imgSet.setAlignmentProj()
for item in subtomoClassesSet:
item.setAlignmentProj()
# with open(dataStar) as fid:
# self.dataTable.readStar(fid)
self.reader = convert31.Reader(alignType=ALIGN_PROJ,
pixelSize=subtomoClassesSet.getSamplingRate())
mdIter = Table.iterRows(dataStar, key='rlnImageName')
subtomoClassesSet.copyItems(self._getInputParticles(),
doClone=False,
updateItemCallback=self._updateParticle,
itemDataIterator=mdIter)
def _updateParticle(self, particle, row):
self.reader.setParticleTransform(particle, row)
if not hasattr(particle, '_rlnRandomSubset'):
particle._rlnRandomSubset = Integer()
particle._rlnRandomSubset.set(row.rlnRandomSubset)
@staticmethod
def _getLastFileName(pattern):
files = glob.glob(pattern)
files.sort(key=getmtime)
return files[-1]