# **************************************************************************
# *
# * Authors: J.M. De la Rosa Trevin (delarosatrevin@scilifelab.se) [1]
# * Federico P. de Isidro Gomez (fp.deisidro@cnb.csic.es) [2]
# *
# * [1] SciLifeLab, Stockholm University
# * [2] Centro Nacional de Biotecnologia, CSIC, Spain
# *
# * 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 pyworkflow as pw
from pyworkflow import BETA
import pyworkflow.protocol.params as params
import pyworkflow.utils.path as path
from pwem.protocols import EMProtocol
import tomo.objects as tomoObj
from tomo.protocols import ProtTomoBase
from tomo.convert import writeTiStack
import tomo.constants as constants
from imod import Plugin
from imod import utils
from pwem.emlib.image import ImageHandler
[docs]class ProtImodEtomo(EMProtocol, ProtTomoBase):
"""
Simple wrapper around etomo to manually reconstruct a Tomogram.
More info:
https://bio3d.colorado.edu/imod/doc/etomoTutorial.html
"""
_label = 'Etomo interactive'
_devStatus = BETA
# -------------------------- DEFINE param functions -----------------------
def _defineParams(self, form):
form.addSection('Input')
form.addParam('inputSetOfTiltSeries',
params.PointerParam,
pointerClass='SetOfTiltSeries',
important=True,
label='Input set of Tilt-Series',
help='Input set of tilt-series to be processed with etomo.')
form.addParam('excludeList',
params.StringParam,
default='',
label='Exclusion list',
help='Provide tilt images IDs (usually starting at 1) that you want to exclude from the '
'processing separated by blank spaces.')
form.addParam('markersDiameter',
params.FloatParam,
default=10,
label='Fiducial markers diameter (nm)',
help='Diameter of gold beads in nanometers.')
# -------------------------- INSERT steps functions ---------------------
# Overwrite the following function to prevent streaming from base class
def _stepsCheck(self):
pass
def _insertAllSteps(self):
self.inputTiltSeries = None
self._insertFunctionStep(self.runEtomoStep, interactive=True)
# --------------------------- STEPS functions ----------------------------
[docs] def runEtomoStep(self):
from imod.viewers import ImodGenericViewer
setOftiltSeries = self.inputSetOfTiltSeries.get()
view = ImodGenericViewer(None, self, setOftiltSeries,
displayAllButton=False, isInteractive=True,
itemDoubleClick=True)
view.show()
self.createOutput()
[docs] def runAllSteps(self, obj):
for item in self.inputSetOfTiltSeries.get():
if item.getTsId() == obj.getTsId():
self.runEtomo(item)
break
[docs] def getFilePath(self, ts, suffix="", extension=""):
tsId = ts.getTsId()
extraPrefix = self._getExtraPath(tsId)
return os.path.join(extraPrefix, ts.getFirstItem().parseFileName(suffix=suffix,
extension=extension))
[docs] def runEtomo(self, ts):
tsId = ts.getTsId()
edfFilePath = self._getExtraPath(os.path.join(tsId, ts.getFirstItem().parseFileName(extension=".edf")))
if not os.path.exists(edfFilePath):
self.convertInputStep(ts)
if ts is not None:
extraPrefix = self._getExtraPath(tsId)
args = '--fg '
args += ts.getFirstItem().parseFileName(extension=".edf")
Plugin.runImod(self, 'etomo', args, cwd=extraPrefix)
[docs] def createOutput(self):
outputPrealiSetOfTiltSeries = None
outputAliSetOfTiltSeries = None
outputSetOfLandmarkModelsNoGaps = None
outputSetOfCoordinates3D = None
outputSetOfFullTomograms = None
outputSetOfPostProcessTomograms = None
setOfTiltSeries = self.inputSetOfTiltSeries.get()
for ts in setOfTiltSeries:
self.inputTiltSeries = ts
tsId = ts.getTsId()
"""Prealigned tilt-series"""
prealiFilePath = self.getFilePath(ts, extension=".preali")
if os.path.exists(prealiFilePath):
if outputPrealiSetOfTiltSeries is None:
outputPrealiSetOfTiltSeries = self._createSetOfTiltSeries(suffix='Preali')
outputPrealiSetOfTiltSeries.copyInfo(setOfTiltSeries)
outputPrealiSetOfTiltSeries.setDim(setOfTiltSeries.getDim())
self._defineOutputs(outputPrealignedSetOfTiltSeries=outputPrealiSetOfTiltSeries)
self._defineSourceRelation(self.inputTiltSeries,
outputPrealiSetOfTiltSeries)
else:
outputPrealiSetOfTiltSeries.enableAppend()
newTs = ts.clone()
newTs.copyInfo(ts)
outputPrealiSetOfTiltSeries.append(newTs)
ih = ImageHandler()
for index, tiltImage in enumerate(ts.iterItems(iterate=False)):
newTi = tiltImage.clone()
newTi.copyInfo(tiltImage, copyId=True, copyTM=False)
newTi.setAcquisition(tiltImage.getAcquisition())
newTi.setLocation(index + 1, prealiFilePath)
index += 1
xPreali, _, _, _ = ih.getDimensions(newTi.getFileName()+":mrc")
newTi.setSamplingRate(self.getPixSizeFromDimensions(xPreali))
newTs.append(newTi)
xPreali, yPreali, zPreali, _ = ih.getDimensions(newTs.getFirstItem().getFileName()+":mrc")
newTs.setDim((xPreali, yPreali, zPreali))
newTs.write(properties=False)
outputPrealiSetOfTiltSeries.setSamplingRate(self.getPixSizeFromDimensions(xPreali))
outputPrealiSetOfTiltSeries.write()
self._store(outputPrealiSetOfTiltSeries)
"""Aligned tilt-series"""
aligFilePath = self.getFilePath(ts, extension=".ali")
if os.path.exists(aligFilePath):
if outputAliSetOfTiltSeries is None:
outputAliSetOfTiltSeries = self._createSetOfTiltSeries(suffix='Ali')
outputAliSetOfTiltSeries.copyInfo(setOfTiltSeries)
outputAliSetOfTiltSeries.setDim(setOfTiltSeries.getDim())
self._defineOutputs(outputAlignedSetOfTiltSeries=outputAliSetOfTiltSeries)
self._defineSourceRelation(self.inputSetOfTiltSeries,
outputAliSetOfTiltSeries)
else:
outputAliSetOfTiltSeries.enableAppend()
newTs = ts.clone()
newTs.copyInfo(ts)
outputAliSetOfTiltSeries.append(newTs)
ih = ImageHandler()
tltFilePath = self.getFilePath(ts, suffix='_fid',
extension=".tlt")
if os.path.exists(tltFilePath):
tltList = utils.formatAngleList(tltFilePath)
else:
tltList = None
for index, tiltImage in enumerate(ts.iterItems(iterate=False)):
newTi = tiltImage.clone()
newTi.copyInfo(tiltImage, copyId=True, copyTM=False)
newTi.setAcquisition(tiltImage.getAcquisition())
newTi.setLocation(index + 1, aligFilePath)
if tltList is not None:
newTi.setTiltAngle(float(tltList[index]))
xAli, _, _, _ = ih.getDimensions(newTi.getFileName()+":mrc")
newTi.setSamplingRate(self.getPixSizeFromDimensions(xAli))
newTs.append(newTi)
xAli, yAli, zAli, _ = ih.getDimensions(newTs.getFirstItem().getFileName() + ":mrc")
newTs.setDim((xAli, yAli, zAli))
newTs.write(properties=False)
outputAliSetOfTiltSeries.setSamplingRate(self.getPixSizeFromDimensions(xAli))
outputAliSetOfTiltSeries.write()
self._store(outputAliSetOfTiltSeries)
"""Output set of coordinates 3D (associated to the aligned tilt-series)"""
coordFilePath = self.getFilePath(ts, suffix='fid', extension=".xyz")
if os.path.exists(coordFilePath):
if outputSetOfCoordinates3D is None:
outputSetOfCoordinates3D = self._createSetOfCoordinates3D(volSet=outputAliSetOfTiltSeries,
suffix='LandmarkModel')
outputSetOfCoordinates3D.setSamplingRate(self.inputTiltSeries.getSamplingRate())
outputSetOfCoordinates3D.setPrecedents(outputAliSetOfTiltSeries)
self._defineOutputs(outputSetOfCoordinates3D=outputSetOfCoordinates3D)
self._defineSourceRelation(self.inputSetOfTiltSeries,
outputSetOfCoordinates3D)
else:
outputSetOfCoordinates3D.enableAppend()
coordList = utils.format3DCoordinatesList(coordFilePath)
for element in coordList:
newCoord3D = tomoObj.Coordinate3D()
newCoord3D.setVolume(ts)
newCoord3D.setX(element[0], constants.BOTTOM_LEFT_CORNER)
newCoord3D.setY(element[1], constants.BOTTOM_LEFT_CORNER)
newCoord3D.setZ(element[2], constants.BOTTOM_LEFT_CORNER)
newCoord3D.setVolId(ts.getObjId())
outputSetOfCoordinates3D.append(newCoord3D)
outputSetOfCoordinates3D.update(newCoord3D)
outputSetOfCoordinates3D.write()
self._store(outputSetOfCoordinates3D)
"""Landmark models with no gaps"""
modelFilePath = self.getFilePath(ts, suffix="_nogaps", extension=".fid")
residFilePath = self.getFilePath(ts, extension=".resid")
if os.path.exists(modelFilePath) and os.path.exists(residFilePath):
modelFilePathTxt = self.getFilePath(ts, suffix="_nogaps_fid", extension=".txt")
paramsNoGapPoint2Model = {
'inputFile': modelFilePath,
'outputFile': modelFilePathTxt
}
argsNoGapPoint2Model = "-InputFile %(inputFile)s " \
"-OutputFile %(outputFile)s"
Plugin.runImod(self, 'model2point', argsNoGapPoint2Model % paramsNoGapPoint2Model)
if outputSetOfLandmarkModelsNoGaps is None:
outputSetOfLandmarkModelsNoGaps = self._createSetOfLandmarkModels(suffix='NoGaps')
outputSetOfLandmarkModelsNoGaps.copyInfo(setOfTiltSeries)
self._defineOutputs(outputSetOfLandmarkModelsNoGaps=outputSetOfLandmarkModelsNoGaps)
self._defineSourceRelation(self.inputSetOfTiltSeries,
outputSetOfLandmarkModelsNoGaps)
else:
outputSetOfLandmarkModelsNoGaps.enableAppend()
fiducialNoGapList = utils.formatFiducialList(modelFilePathTxt)
landmarkModelNoGapsFilePath = self.getFilePath(ts, suffix="_nogaps", extension=".sfid")
fiducialNoGapsResidList = utils.formatFiducialResidList(residFilePath)
landmarkModelNoGaps = tomoObj.LandmarkModel(tsId,
landmarkModelNoGapsFilePath,
modelFilePath)
prevTiltIm = 0
chainId = 0
indexFake = 0
for fiducial in fiducialNoGapList:
if int(float(fiducial[2])) <= prevTiltIm:
chainId += 1
prevTiltIm = int(float(fiducial[2]))
if (indexFake < len(fiducialNoGapsResidList) and
fiducial[2] == fiducialNoGapsResidList[indexFake][2]):
landmarkModelNoGaps.addLandmark(xCoor=fiducial[0],
yCoor=fiducial[1],
tiltIm=fiducial[2],
chainId=chainId,
xResid=fiducialNoGapsResidList[indexFake][3],
yResid=fiducialNoGapsResidList[indexFake][4])
indexFake += 1
else:
landmarkModelNoGaps.addLandmark(xCoor=fiducial[0],
yCoor=fiducial[1],
tiltIm=fiducial[2],
chainId=chainId,
xResid='0',
yResid='0')
outputSetOfLandmarkModelsNoGaps.append(landmarkModelNoGaps)
outputSetOfLandmarkModelsNoGaps.write()
self._store(outputSetOfLandmarkModelsNoGaps)
"""Full reconstructed tomogram"""
reconstructTomoFilePath = self.getFilePath(ts, suffix="_full",
extension=".rec")
if os.path.exists(reconstructTomoFilePath):
if outputSetOfFullTomograms is None:
outputSetOfFullTomograms = self._createSetOfTomograms(suffix='Full')
outputSetOfFullTomograms.copyInfo(setOfTiltSeries)
self._defineOutputs(outputSetOfFullTomograms=outputSetOfFullTomograms)
self._defineSourceRelation(self.inputSetOfTiltSeries,
outputSetOfFullTomograms)
else:
outputSetOfFullTomograms.enableAppend()
newTomogram = tomoObj.Tomogram()
newTomogram.setLocation(reconstructTomoFilePath)
newTomogram.setTsId(tsId)
newTomogram.setSamplingRate(ts.getSamplingRate())
# Set default tomogram origin
newTomogram.setOrigin(newOrigin=False)
outputSetOfFullTomograms.append(newTomogram)
outputSetOfFullTomograms.write()
self._store(outputSetOfFullTomograms)
"""Post-processed reconstructed tomogram"""
posprocessedRecTomoFilePath = self.getFilePath(ts, extension=".rec")
if os.path.exists(posprocessedRecTomoFilePath):
if outputSetOfPostProcessTomograms is None:
outputSetOfPostProcessTomograms = self._createSetOfTomograms()
outputSetOfPostProcessTomograms.copyInfo(setOfTiltSeries)
self._defineOutputs(outputSetOfPostProcessTomograms=outputSetOfPostProcessTomograms)
self._defineSourceRelation(self.inputSetOfTiltSeries,
outputSetOfPostProcessTomograms)
else:
outputSetOfPostProcessTomograms.enableAppend()
newTomogram = tomoObj.Tomogram()
ih = ImageHandler()
outputDim, _, _, _ = ih.getDimensions(posprocessedRecTomoFilePath)
newTomogram.setLocation(posprocessedRecTomoFilePath)
newTomogram.setTsId(tsId)
newTomogram.setSamplingRate(self.getPixSizeFromDimensions(outputDim))
# Set default tomogram origin
newTomogram.setOrigin(newOrigin=False)
outputSetOfPostProcessTomograms.append(newTomogram)
outputSetOfPostProcessTomograms.write()
self._store(outputSetOfPostProcessTomograms)
self.closeMappers()
# --------------------------- UTILS functions ----------------------------
@staticmethod
def _writeEtomoEdf(fn, paramsDict):
template = """
#%(date)s // Generated by Scipion, version: %(version)s
Setup.DataSource=CCD
ReconstructionState.InvalidEdgeFunctionsA=no result
Setup.BackupDirectory=
ReconstructionState.InvalidEdgeFunctionsB=no result
ProcessTrack.PostProcessing=Not started
Setup.Combine.ManualCleanup=false
Setup.AxisA.TiltAngle.Type=File
Setup.FiducialessAlignmentB=false
Setup.FiducialessAlignmentA=false
ProcessTrack.FinalAlignedStack-A=Not started
ProcessTrack.FinalAlignedStack-B=Not started
Setup.tiltalign.TargetPatchSizeXandY=700,700
Setup.AxisB.TiltAngle.RangeMin=%(minTilt)f
Setup.Combine.TempDirectory=
Setup.Combine.UseList=
Setup.Combine.PatchBoundaryYMax=0
Setup.WholeTomogramSampleA=false
Setup.DatasetName=%(name)s
Setup.FiducialDiameter=%(markerDiameter)f
Setup.WholeTomogramSampleB=false
Setup.SetFEIPixelSize=true
ReconstructionState.A.AdjustOrigin=true
Setup.Setup.OrigImageStackExt=.st
ProcessTrack.RevisionNumber=2.0
Setup.ProjectLog.FrameSize.Width=683
Setup.Combine.PatchBoundaryYMin=0
Setup.Combine.MaxPatchBoundaryZMax=0
ProcessTrack.CoarseAlignment-A=Not started
Setup.ViewType=Single View
ProcessTrack.CoarseAlignment-B=Not started
ReconstructionState.TrimvolFlipped=no result
Setup.Track.B.TrackMethod=Seed
Setup.Track.A.TrackMethod=Seed
Setup.A.SizeToOutputInXandY=/
ProcessTrack.TomogramGeneration-A=Not started
ProcessTrack.TomogramGeneration-B=Not started
ProcessTrack.CleanUp=Not started
Setup.AxisA.TiltAngle.RangeMin=%(minTilt)f
Setup.AxisB.TiltAngle.TiltAngleFilename=
Setup.Pos.A.NewDialog=true
Setup.Squeezevol.LinearInterpolation=false
Setup.Stack.B.Is.Twodir=false
Setup.Combine.RevisionNumber=1.2
Setup.Stack.B.CTF.AutoFit.RangeAndStep=-Infinity,-Infinity
Setup.UseLocalAlignmentsB=true
Setup.UseLocalAlignmentsA=true
Setup.AxisB.TiltAngle.RangeStep=1.0
Setup.Combine.FiducialMatchListA=
Setup.Combine.FiducialMatchListB=
Setup.B.SizeToOutputInXandY=/
Setup.Binning=1
Setup.Combine.ModelBased=false
ReconstructionState.MadeZFactorsB=no result
ReconstructionState.MadeZFactorsA=no result
ReconstructionState.SqueezevolFlipped=no result
ProcessTrack.FiducialModel-B=Not started
ProcessTrack.FiducialModel-A=Not started
Setup.Combine.PatchBoundaryXMax=0
ProcessTrack.Setup=Complete
ProcessTrack.PreProcessing-A=Not started
ProcessTrack.PreProcessing-B=Not started
ProcessTrack.FineAlignment-B=Not started
ProcessTrack.FineAlignment-A=Not started
Setup.AxisA.TiltAngle.TiltAngleFilename=
Setup.AxisA.TiltAngle.RangeStep=1.0
Setup=-Infinity,-Infinity
Setup.Combine.PatchBoundaryXMin=0
Setup.MagGradientFile=
Setup.RevisionNumber=1.12
Setup.Track.B.SeedModel.Transfer=true
Setup.Track.A.Raptor.UseRawStack=false
Setup.PixelSize=%(pixelSize)f
ReconstructionState.B.AdjustOrigin=true
ReconstructionState.Combine.ScriptsCreated=no result
Setup.Combine.Transfer=true
Setup.DistortionFile=
ReconstructionState.UsedLocalAlignmentsA=no result
ReconstructionState.UsedLocalAlignmentsB=no result
Setup.ProjectLog.FrameSize.Height=230
Setup.Stack.B.Twodir=0.0
Setup.Combine.FiducialMatch=BothSides
Setup.AxisType=Single Axis
Setup.ImageRotationA=%(rotationAngle)f
ProcessTrack.TomogramPositioning-A=Not started
Setup.ImageRotationB=
Setup.Stack.A.Is.Twodir=false
Setup.Pos.B.NewDialog=true
ProcessTrack.TomogramPositioning-B=Not started
Setup.Combine.PatchBoundaryZMax=0
Setup.DefaultGpuProcessing=false
Setup.Track.A.SeedModel.Auto=true
Setup.Combine.PatchSize=M
Setup.AxisB.TiltAngle.Type=Extract
Setup.Combine.PatchBoundaryZMin=0
Setup.ProjectLog.FrameLocation.Y=55
Setup.ProjectLog.FrameLocation.X=95
Setup.DefaultParallel=false
Setup.ProjectLog.Visible=true
Setup.tiltalign.NumberOfLocalPatchesXandY=5,5
Setup.Combine.PatchRegionModel=
ReconstructionState.NewstFiducialessAlignmentA=no result
ReconstructionState.NewstFiducialessAlignmentB=no result
Setup.Stack.A.Twodir=0.0
ProcessTrack.TomogramCombination=Not started
"""
with open(fn, 'w') as f:
f.write(template % paramsDict)
[docs] def getPixSizeFromDimensions(self, outputDim):
ih = ImageHandler()
originalDim, _, _, _ = ih.getDimensions(self.inputTiltSeries.getFirstItem().getFileName())
return self.inputTiltSeries.getSamplingRate() * round(originalDim/outputDim)
[docs] def getResizeFactorFromDimensions(self, outputDim):
ih = ImageHandler()
originalDim, _, _, _ = ih.getDimensions(self.inputTiltSeries.get().getFirstItem().getFileName())
return round(outputDim / originalDim)
# --------------------------- INFO functions ----------------------------
def _summary(self):
summary = ["The following outputs have been generated from the "
"operations performed over the input tilt-series:"]
if hasattr(self, 'outputPrealignedSetOfTiltSeries'):
summary.append("- Tilt-series prealignment.")
if hasattr(self, 'outputAlignedSetOfTiltSeries'):
summary.append("- Tilt-series alignment.")
if hasattr(self, 'outputSetOfCoordinates3D'):
summary.append("- Landmark 3D coordinates have been extracted.")
if hasattr(self, 'outputSetOfLandmarkModelsGaps'):
summary.append("- Landmark model with gaps has been generated.")
if hasattr(self, 'outputSetOfLandmarkModelsNoGaps'):
summary.append("- Landmark model without gaps has been generated.")
if hasattr(self, 'outputSetOfFullTomograms'):
summary.append("- Full raw reconstructed tomogram.")
if hasattr(self, 'outputSetOfPostProcessTomograms'):
summary.append("- Post processed reconstructed tomogram.")
if summary == ["The following operations has been performed over the input tilt-series:"]:
summary = ["Output classes not ready yet."]
return summary
def _methods(self):
methods = ["The following outputs have been generated from the "
"operations performed over the input tilt-series:"]
if hasattr(self, 'outputPrealignedSetOfTiltSeries'):
methods.append("- Tilt-series prealignment.")
if hasattr(self, 'outputAlignedSetOfTiltSeries'):
methods.append("- Tilt-series alignment.")
if hasattr(self, 'outputSetOfCoordinates3D'):
methods.append("- Landmark 3D coordinates have been extracted.")
if hasattr(self, 'outputSetOfLandmarkModelsGaps'):
methods.append("- Landmark model with gaps has been generated.")
if hasattr(self, 'outputSetOfLandmarkModelsNoGaps'):
methods.append("- Landmark model without gaps has been generated.")
if hasattr(self, 'outputSetOfFullTomograms'):
methods.append("- Full raw reconstructed tomogram.")
if hasattr(self, 'outputSetOfPostProcessTomograms'):
methods.append("- Post processed reconstructed tomogram.")
if methods == ["The following operations has been performed over the input tilt-series:"]:
methods = ["Output classes not ready yet."]
return methods