Difference between revisions of "Walkthrough on command line based tilt series alignment"

From Dynamo
Jump to navigation Jump to search
 
(35 intermediate revisions by the same user not shown)
Line 12: Line 12:
 
= Entering the data =
 
= Entering the data =
  
The variable <tt>u</tt> is an object that contains several areas to interact with the workflow. They can be found through autocompletion by using the <tt>tab</tt> key (or by pressing the enter key in the standalone version).
+
The variable <tt>u</tt> is an object that contains several ''fields'' to interact with the workflow. They can be found through autocompletion by using the <tt>tab</tt> key (or by pressing the enter key in the standalone version). In the following, these different fields will be used to pass our inputs to the workflow.
  
 
[[ File:alignCommandLine_autocompletionOfUserObject.png |thumb|center| 700px|autocompletion of user object to access different fields in matlab and standalone version]]
 
[[ File:alignCommandLine_autocompletionOfUserObject.png |thumb|center| 700px|autocompletion of user object to access different fields in matlab and standalone version]]
  
 
== Basic data==
 
== Basic data==
We use a binned version of the first tilt series from the EMPIAR entry 10164, depicting a set of virus like particles (VLP). Visit the course materials to access the tilt series or download it directly from EMPIAR. You can also follow this tutorial using your own data. The file of the tilt series can be defined as
+
Basic data can be entered using the ''enter'' field. Here we use a binned version of the first tilt series from the EMPIAR entry 10164, depicting a set of virus like particles (VLP). Visit the course materials to access the tilt series or download it directly from EMPIAR (you can also use your own data). The file of the tilt series and the angles are passed to the workflow:
  
<tt>file = 'b001ts001.mrc'; </tt>
+
  <nowiki>u.enter.tiltSeries('b001ts001.mrc');
 
 
and then linked into the workflow:
 
 
 
  <nowiki>u.enter.tiltSeries(file);
 
 
u.enter.tiltAngles(-57:3:60); </nowiki>
 
u.enter.tiltAngles(-57:3:60); </nowiki>
  
Optionally, you can reject tilts that are of too bad quality with the following command (in this tutorial we keep all tilts).
+
If necessary (not in this tutorial), you can use the ''discardedTiltIndices'' field to reject tilts that are of too bad quality using the following command (e.g., tilt 1 and 40):
 
  <tt>u.enter.discardedTiltIndices([1,40]); </tt>
 
  <tt>u.enter.discardedTiltIndices([1,40]); </tt>
  
 
== Acquisition settings==
 
== Acquisition settings==
Here you can enter parameters such as <tt>Cs</tt>, nominal defocus, etc. We will just provide the pixel size in Angstroms:
+
Here you can enter parameters such as <tt>Cs</tt>, nominal defocus, etc. For this tutorial we just provide the pixel size in Angstroms:
 
  <tt>u.enter.settingAcquisition.apix(2.7);</tt>
 
  <tt>u.enter.settingAcquisition.apix(2.7);</tt>
  
Line 42: Line 38:
  
 
== Detection settings==
 
== Detection settings==
These here are the actual design decisions when running an alignment workflow. The gold bead radius should have been  measured beforehand:
+
These here are the actual design decisions when running an alignment workflow. The gold bead radius has to be known beforehand ( see [[Walkthrough on GUI based tilt series alignment | GUI based alignment walkthrough]] on how to measure it). Binning is used internally only for quick and robust bead detection.
  
  <nowiki>u.enter.settingDetection.beadRadius(16);
+
  <nowiki>u.enter.settingDetection.detectionBinningFactor(1);
 +
u.enter.settingDetection.beadRadius(16);
 
u.enter.settingDetection.maskRadius(28);
 
u.enter.settingDetection.maskRadius(28);
u.enter.templateSidelength(64); </nowiki>
+
u.enter.templateSidelength(64);</nowiki>
 
 
The bin level is mainly used to accelerate the detection procedure. Needs to be chosen in a way that a binned gold bead still can be recognisable as such, with a radius of at least 4 pixels.
 
<tt>u.enter.settingDetection.detectionBinningFactor(1);</tt>
 
  
 
== Changing generic parameters ==
 
== Changing generic parameters ==
Yo can find handles to the parameters of the individuals steps through autocompletion on the <tt>area</tt>, then <tt>step</tt> items.
+
You can find all parameters of the individuals steps using autocompletion on the <tt>area</tt> field and then by adding <tt>step</tt> in the end. For example, the residual threshold of the tilt gap filler can be changed as follows:
 
  <tt>u.area.indexing.step.tiltGapFiller.parameterSet.residualsThreshold(8);</tt>
 
  <tt>u.area.indexing.step.tiltGapFiller.parameterSet.residualsThreshold(8);</tt>
  
 
= Running the workflow =
 
= Running the workflow =
We will skip CTF estimation and correction steps (both in the CTF area):
+
After setting the input parameters, the complete workflow can be run using the following command (we skip the CTF correction here). It should take about 5 minutes to run.
 
  <tt>u.run.all('noctf',1);</tt>
 
  <tt>u.run.all('noctf',1);</tt>
 
If should take around five minutes to run.
 
  
 
== Checking the results ==
 
== Checking the results ==
After completion, the <tt>view</tt> field allows you to check different output items for diagnose, as the tilt lines:
+
After completion, the <tt>view</tt> field allows you to check different output items for diagnose, such as the tilt lines:
  
 
  <tt>u.view.tiltLines.unaligned();</tt>
 
  <tt>u.view.tiltLines.unaligned();</tt>
Line 68: Line 60:
 
[[ File:alignCommandLine_viewOfTiltLines.png |thumb|center| 600px|view of tilt lines]]
 
[[ File:alignCommandLine_viewOfTiltLines.png |thumb|center| 600px|view of tilt lines]]
  
or the reconstructions, as for example
+
or the different reconstructions  
 
 
<tt>u.view.reconstruction.binnedWBP</tt>
 
 
 
or
 
  
  <tt>u.view.reconstruction.binnedSIRT</tt>
+
  <nowiki>u.view.reconstruction.binnedWBP();
 +
u.view.reconstruction.binnedSIRT();</nowiki>
  
and aligned stacks themselves.
+
or the aligned stacks themselves.
  
  <tt>u.view.stack.alignedFull</tt>
+
  <tt>u.view.stack.alignedFull();</tt>
  
 
+
You can also use the autocompletion of the <tt>info</tt> field (as shown in the beginning) to find and display different logs:
You can also use the autocompletion (in the Matlab version) of the <tt>info</tt> field in order to read the logs on the fit.
 
  
 
  <nowiki>>> u.info.fit;
 
  <nowiki>>> u.info.fit;
 
------------------------------------------------------------
 
------------------------------------------------------------
 
File with info:
 
File with info:
   workflows/workshop2.AWF/info/fitting.doc
+
   workflows/hivCommandLine.AWF/info/fitting.doc
 
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  
 
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -  
rms:  1.39
+
rms:  1.41
psi: 85.10
+
psi: 85.20
 
fit modus: normal
 
fit modus: normal
total markers: 78
+
total markers: 79
total observations: 2590
+
total observations: 2702
fit created at: 27-Aug-2019 18:01:40
+
fit created at: 11-Aug-2021 10:56:24
 
------------------------------------------------------------</nowiki>
 
------------------------------------------------------------</nowiki>
 
or of the markers
 
  
 
  <nowiki>>> u.info.markers;
 
  <nowiki>>> u.info.markers;
Line 103: Line 89:
 
   ... item saved and in memory
 
   ... item saved and in memory
 
Item in memory:
 
Item in memory:
Number of shapes          :  78
+
Number of shapes          :  79
 
Number of frames          :  40
 
Number of frames          :  40
Number of observations    :  2590 (out of 3120)
+
Number of observations    :  2702 (out of 3160)
Number of empty frames    :  2
+
Number of empty frames    :  1
 
Number of complete frames :  0
 
Number of complete frames :  0
 
Number of complete shapes :  0
 
Number of complete shapes :  0
Line 118: Line 104:
 
   ... item saved and in memory
 
   ... item saved and in memory
 
Item in memory:
 
Item in memory:
Number of shapes          :  78
+
Number of shapes          :  79
 
Number of frames          :  40
 
Number of frames          :  40
Number of observations    :  2590 (out of 3120)
+
Number of observations    :  2702 (out of 3160)
Number of empty frames    :  2
+
Number of empty frames    :  1
 
Number of complete frames :  0
 
Number of complete frames :  0
 
Number of complete shapes :  0
 
Number of complete shapes :  0
Line 128: Line 114:
 
Markers in disk and currently in memory contain equal coordinates
 
Markers in disk and currently in memory contain equal coordinates
 
------------------------------------------------------------</nowiki>
 
------------------------------------------------------------</nowiki>
 +
 +
If you are following this tutorial during a ''Dynamo'' workshop, the exercise ends here.
  
 
= On the fly reconstruction of particles=
 
= On the fly reconstruction of particles=
  
In this section of the a, we show how to avoid reconstructing ''full'' tomograms of ''full'' resolution if only a set of subtomograms need to be reconstructed at high resolution. We will use the gold beads as proxies for particles of interest. We will thus locate a bunch of them in the '''binned''' tomogram, and then use their coordinates in that tomogram to directly reconstruct them in full resolution from the full resolution stack.
+
In this section we show how to avoid having to reconstruct a '''full''' tomogram at '''full''' resolution for the case where only a set of subtomograms are need. For this exercise we will use the gold beads as proxies that represent any particles of interest. We will first locate a few of them in the '''binned''' tomogram, and then use those coordinates to directly reconstruct the subtomograms in full resolution and skipping the need to reconstruct the whole full sized tomogram..
  
== Grep coordinates in the binned tomogram==
+
== Define coordinates in the binned tomogram==
  
Open the binned reconstruction with <tt>dtmslice</tt>
+
First, we open the binned reconstruction tht we created before using the visualization tool <tt>dtmslice</tt> and adjust the contrast if necessary (arrow):
  
<tt>dtmslice 'workflows/hivCommandLine.AWF/reconstruction/binnedReconstructionWBP.mrc';</tt>
+
<tt>dtmslice 'workflows/hivCommandLine.AWF/reconstruction/binnedReconstructionWBP.mrc';</tt>
  
 
[[ File:alignCommandLine_dtmsliceOfBinnedReconstruction.png |thumb|center| 600px|dtmslice of binned reconstruction]]
 
[[ File:alignCommandLine_dtmsliceOfBinnedReconstruction.png |thumb|center| 600px|dtmslice of binned reconstruction]]
You will need to adjust the contrast:
 
[[ File:alignCommandLine_dtmsliceContrastControl.png |thumb|center| 200px|dtmslice contrast control]]
 
 
[[ File:alignCommandLine_dtmsliceAfterContrastCorrection.png |thumb|center| 600px|dtmslice after contrast correction]]
 
[[ File:alignCommandLine_dtmsliceAfterContrastCorrection.png |thumb|center| 600px|dtmslice after contrast correction]]
Now we create a new model of <tt>general</tt> type
+
To manually define the particle coordinates, we first open a new model of type <tt>general</tt>:
 
[[ File:alignCommandLine_createANewModel.png |thumb|center| 600px|create a new model]]
 
[[ File:alignCommandLine_createANewModel.png |thumb|center| 600px|create a new model]]
We use the projections in order to quickly find the 3d positions of a set of tomograms
+
To quickly find the 3d positions of a set of particles, we use a projection of the tomogram:
 
[[ File:alignCommandLine_showAFullProjection.png |thumb|center| 600px|show a full projection]]
 
[[ File:alignCommandLine_showAFullProjection.png |thumb|center| 600px|show a full projection]]
 
[[ File:alignCommandLine_fullProjectionOfBinnedReconstruction.png |thumb|center| 600px|full projection of binned reconstruction]]
 
[[ File:alignCommandLine_fullProjectionOfBinnedReconstruction.png |thumb|center| 600px|full projection of binned reconstruction]]
You can right click on the (x,y) position of a marker and get it depicted on its cutting xz and yz slices
+
Right click on the (x,y) position of a particle opens two small windows showing the x-z and y-z planes:
 
[[ File:alignCommandLine_rightClickOnAPoint.png |thumb|center| 600px|right click on a point]]
 
[[ File:alignCommandLine_rightClickOnAPoint.png |thumb|center| 600px|right click on a point]]
Left clicking on it will provide a z coordinate to it
+
Use the y-z depiction to left clicking on the marker. After clicking, the marker coordinates are saved in the model:
 
[[ File:alignCommandLine_xViewOfClickedPoint.png |thumb|center| 600px|x view of clicked point]]
 
[[ File:alignCommandLine_xViewOfClickedPoint.png |thumb|center| 600px|x view of clicked point]]
 
[[ File:alignCommandLine_mainClickToFixThe3dCoordinates.png |thumb|center| 600px|main click to fix the 3d coordinates]]
 
[[ File:alignCommandLine_mainClickToFixThe3dCoordinates.png |thumb|center| 600px|main click to fix the 3d coordinates]]
Keep clicking till you have a small set (no more than ten)
+
Go back to the projection view and repeat the procedure until you have a few more markers (but no more than ten in this exercise). The x-z and y-z planes will be updated automatically when clicking on new markers on the projection and do not need to be closed every time. When you are done, you can close all the plane views and also the projection view.
 
[[ File:alignCommandLine_clickAroundTenMarkers.png |thumb|center| 600px|click around ten markers]]
 
[[ File:alignCommandLine_clickAroundTenMarkers.png |thumb|center| 600px|click around ten markers]]
We export the coordinates of the model into the Matlab workspace already in table format, so that we will be able to use it in <tt>dtcrop</tt> on the binned tomogram.
+
To save the coordinates we first export them into the Matlab workspace in table format:
 
[[ File:alignCommandLine_exportATableIntoTheMatlabWorkspace.png |thumb|center| 600px|export a table into the matlab workspace]]
 
[[ File:alignCommandLine_exportATableIntoTheMatlabWorkspace.png |thumb|center| 600px|export a table into the matlab workspace]]
The Matlab shell should report the variable name used to store the table.
+
The Matlab shell should report the variable name used to store the table:
 
[[ File:alignCommandLine_theCommandLineShowsTheVariableNameAssignedToTheCreatedTable.png |thumb|center| 400px|the command line shows the variable name assigned to the created table]]
 
[[ File:alignCommandLine_theCommandLineShowsTheVariableNameAssignedToTheCreatedTable.png |thumb|center| 400px|the command line shows the variable name assigned to the created table]]
  
Now we can save the table as a file:
+
Now we can save the table with the coordinates as a file in the current directory (this table will be used to crop the subtomograms in the next steps):
<tt>write(temp_table,'binnedGoldbeads.tbl');</tt>
+
<tt>dwrite(temp_table,'binnedGoldbeads.tbl');</tt>
  
 
=== Cropping the gold bead particles ===
 
=== Cropping the gold bead particles ===
  
We can now crop the particles using the table file
+
We can now crop the particles using the table file that we generated before:
 
  <nowiki>tomo = 'workflows/hivCommandLine.AWF/reconstruction/binnedReconstructionWBP.mrc';
 
  <nowiki>tomo = 'workflows/hivCommandLine.AWF/reconstruction/binnedReconstructionWBP.mrc';
 
tbl  = 'binnedGoldbeads.tbl';
 
tbl  = 'binnedGoldbeads.tbl';
Line 172: Line 158:
 
o = dtcrop(tomo,tbl, targetBinnedParticles, sidelengthBinned);</nowiki>
 
o = dtcrop(tomo,tbl, targetBinnedParticles, sidelengthBinned);</nowiki>
  
And then we visualise a projection along <tt>z</tt> of each cropped particle
+
To monitor if everything worked correctly, we visualize a projection along <tt>z</tt> of each cropped particle using the command:
 
  <tt>dslices(targetBinnedParticles,'proj','c');shg;</tt>
 
  <tt>dslices(targetBinnedParticles,'proj','c');shg;</tt>
  
[[ File:alignCommandLine_dslicesOnBinnedGoldBeads.png |thumb|center| 600px|dslices on binned gold beads]]</nowiki>
+
[[ File:alignCommandLine_dslicesOnBinnedGoldBeads.png |thumb|center| 600px|dslices on binned gold beads]]
 
 
 
== Reconstruct particles on the fly ==
 
== Reconstruct particles on the fly ==
  
The command that reconstructs full sized particles from an aligned stack is <tt>dynama_table_trec</tt>. This is an independent function, not part of the workflow, implying that we will need to perform some format conversion.
+
The command that reconstructs full sized particles from an aligned stack is <tt>dynamo_table_trec</tt>. This is an independent function, not part of the workflow, implying that we will need to perform some format conversion.
  
<tt>tbl    = dread('binnedGoldbeads.tbl');</tt>
+
<tt>tbl    = dread('binnedGoldbeads.tbl');</tt>
  
 
We scale the table. "Bin level 2" means that each pixel in the binned tomogram is 2 times twice its size in the full resolution tomogram. This is the factor that we input into <tt>dynamo_table_rescale</tt>:
 
We scale the table. "Bin level 2" means that each pixel in the binned tomogram is 2 times twice its size in the full resolution tomogram. This is the factor that we input into <tt>dynamo_table_rescale</tt>:
Line 206: Line 191:
 
  <tt>sidelengthFullSize = 128;</tt>
 
  <tt>sidelengthFullSize = 128;</tt>
 
We use the <tt>applyRampFilter</tt> option in order to explicitely estate the kind of filter that we want during the back projection.
 
We use the <tt>applyRampFilter</tt> option in order to explicitely estate the kind of filter that we want during the back projection.
  <nowiki>dynamo_table_rec(tblFull,source,angles,targetDirectory, sidelengthFullSize,....
+
  <nowiki>dynamo_table_rec(tblFull,source,angles,targetDirectory, sidelengthFullSize, 'applyRampFilter',1, 'sizeTomogram',sizeTomogramFull, 'shiftTomogramCenter',shiftTomogramCenter);</nowiki>
    'applyRampFilter',1,....
 
    'sizeTomogram',sizeTomogramFull,....
 
    'shiftTomogramCenter',shiftTomogramCenter);</nowiki>
 
  
 
We can now check that the particles have been reconstructed in full size:
 
We can now check that the particles have been reconstructed in full size:
 
  <tt>figure;dslices(targetDirectory,'proj','c');shg;</tt>
 
  <tt>figure;dslices(targetDirectory,'proj','c');shg;</tt>
 
[[ File:alignCommandLine_dslicesOnFullSizedGoldBeads.png |thumb|center| 600px|dslices on full sized gold beads]]
 
[[ File:alignCommandLine_dslicesOnFullSizedGoldBeads.png |thumb|center| 600px|dslices on full sized gold beads]]

Latest revision as of 11:15, 11 August 2021

Dynamo includes a package for automated alignment and reconstruction of tilt series. This walkthrough guides users through a series of steps on how to use it in command line based mode. The GUI based alignment walkthrough shows how to operate this procedure through the GUI. Here, we go through each command one by one (for advanced applications, all commands can also be summarized in a Matlab script and run at once).

Create the workflow

The following commands create a new alignment workflow in a separate folder. Make sure Dynamo is running and type the following commands into the Matlab command line (or the Dynamo console, if you work on the standalone version):

name = 'hivCommandLine';
folder = 'workflows';
u = dtsa(name,'--nogui','-path',folder,'fp',1); 

In the Matlab workspace you should now see the variables that you just created with these commands. If you work with the standalone version you can simply type the word whos to see a list of all variables in the workspace.

Entering the data

The variable u is an object that contains several fields to interact with the workflow. They can be found through autocompletion by using the tab key (or by pressing the enter key in the standalone version). In the following, these different fields will be used to pass our inputs to the workflow.

autocompletion of user object to access different fields in matlab and standalone version

Basic data

Basic data can be entered using the enter field. Here we use a binned version of the first tilt series from the EMPIAR entry 10164, depicting a set of virus like particles (VLP). Visit the course materials to access the tilt series or download it directly from EMPIAR (you can also use your own data). The file of the tilt series and the angles are passed to the workflow:

u.enter.tiltSeries('b001ts001.mrc');
u.enter.tiltAngles(-57:3:60); 

If necessary (not in this tutorial), you can use the discardedTiltIndices field to reject tilts that are of too bad quality using the following command (e.g., tilt 1 and 40):

u.enter.discardedTiltIndices([1,40]); 

Acquisition settings

Here you can enter parameters such as Cs, nominal defocus, etc. For this tutorial we just provide the pixel size in Angstroms:

u.enter.settingAcquisition.apix(2.7);

Computation settings

We enable the use of parallel cores

u.enter.settingComputing.parallelCPUUse(1); 

This will be using mainly for gold bead detection (cc computation) and reconstruction.

To use all available physical cores use:

u.enter.settingComputing.cpus('*');

Detection settings

These here are the actual design decisions when running an alignment workflow. The gold bead radius has to be known beforehand ( see GUI based alignment walkthrough on how to measure it). Binning is used internally only for quick and robust bead detection.

u.enter.settingDetection.detectionBinningFactor(1);
u.enter.settingDetection.beadRadius(16);
u.enter.settingDetection.maskRadius(28);
u.enter.templateSidelength(64);

Changing generic parameters

You can find all parameters of the individuals steps using autocompletion on the area field and then by adding step in the end. For example, the residual threshold of the tilt gap filler can be changed as follows:

u.area.indexing.step.tiltGapFiller.parameterSet.residualsThreshold(8);

Running the workflow

After setting the input parameters, the complete workflow can be run using the following command (we skip the CTF correction here). It should take about 5 minutes to run.

u.run.all('noctf',1);

Checking the results

After completion, the view field allows you to check different output items for diagnose, such as the tilt lines:

u.view.tiltLines.unaligned();
view of tilt lines

or the different reconstructions

u.view.reconstruction.binnedWBP();
u.view.reconstruction.binnedSIRT();

or the aligned stacks themselves.

u.view.stack.alignedFull();

You can also use the autocompletion of the info field (as shown in the beginning) to find and display different logs:

>> u.info.fit;
------------------------------------------------------------
File with info:
  workflows/hivCommandLine.AWF/info/fitting.doc
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
rms:  1.41
psi: 85.20
fit modus: normal
total markers: 79
total observations: 2702
fit created at: 11-Aug-2021 10:56:24
------------------------------------------------------------
>> u.info.markers;
------------------------------------------------------------
Info on working markers [used during detection, indexing refinement]
  ... item saved and in memory
Item in memory:
Number of shapes          :  79
Number of frames          :  40
Number of observations    :  2702 (out of 3160)
Number of empty frames    :  1
Number of complete frames :  0
Number of complete shapes :  0
 
Item in disk:
Markers in disk and currently in memory contain equal coordinates
 
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
 
Info on final markers [used from alignment step]
  ... item saved and in memory
Item in memory:
Number of shapes          :  79
Number of frames          :  40
Number of observations    :  2702 (out of 3160)
Number of empty frames    :  1
Number of complete frames :  0
Number of complete shapes :  0
 
Item in disk:
Markers in disk and currently in memory contain equal coordinates
------------------------------------------------------------

If you are following this tutorial during a Dynamo workshop, the exercise ends here.

On the fly reconstruction of particles

In this section we show how to avoid having to reconstruct a full tomogram at full resolution for the case where only a set of subtomograms are need. For this exercise we will use the gold beads as proxies that represent any particles of interest. We will first locate a few of them in the binned tomogram, and then use those coordinates to directly reconstruct the subtomograms in full resolution and skipping the need to reconstruct the whole full sized tomogram..

Define coordinates in the binned tomogram

First, we open the binned reconstruction tht we created before using the visualization tool dtmslice and adjust the contrast if necessary (arrow):

dtmslice 'workflows/hivCommandLine.AWF/reconstruction/binnedReconstructionWBP.mrc';
dtmslice of binned reconstruction
dtmslice after contrast correction

To manually define the particle coordinates, we first open a new model of type general:

create a new model

To quickly find the 3d positions of a set of particles, we use a projection of the tomogram:

show a full projection
full projection of binned reconstruction

Right click on the (x,y) position of a particle opens two small windows showing the x-z and y-z planes:

right click on a point

Use the y-z depiction to left clicking on the marker. After clicking, the marker coordinates are saved in the model:

x view of clicked point
main click to fix the 3d coordinates

Go back to the projection view and repeat the procedure until you have a few more markers (but no more than ten in this exercise). The x-z and y-z planes will be updated automatically when clicking on new markers on the projection and do not need to be closed every time. When you are done, you can close all the plane views and also the projection view.

click around ten markers

To save the coordinates we first export them into the Matlab workspace in table format:

export a table into the matlab workspace

The Matlab shell should report the variable name used to store the table:

the command line shows the variable name assigned to the created table

Now we can save the table with the coordinates as a file in the current directory (this table will be used to crop the subtomograms in the next steps):

dwrite(temp_table,'binnedGoldbeads.tbl');

Cropping the gold bead particles

We can now crop the particles using the table file that we generated before:

tomo = 'workflows/hivCommandLine.AWF/reconstruction/binnedReconstructionWBP.mrc';
tbl  = 'binnedGoldbeads.tbl';
targetBinnedParticles = 'binnedGoldBeads.Data';
sidelengthBinned = 32;
o = dtcrop(tomo,tbl, targetBinnedParticles, sidelengthBinned);

To monitor if everything worked correctly, we visualize a projection along z of each cropped particle using the command:

dslices(targetBinnedParticles,'proj','c');shg;
dslices on binned gold beads

Reconstruct particles on the fly

The command that reconstructs full sized particles from an aligned stack is dynamo_table_trec. This is an independent function, not part of the workflow, implying that we will need to perform some format conversion.

tbl     = dread('binnedGoldbeads.tbl');

We scale the table. "Bin level 2" means that each pixel in the binned tomogram is 2 times twice its size in the full resolution tomogram. This is the factor that we input into dynamo_table_rescale:

alignmentBinLevel = 2;
tblFull = dynamo_table_rescale(tbl,'factor',2^alignmentBinLevel); 

Now we fetch the full sized aligned stack stored inside the alignment folder:

source  = 'workflows/hivCommandLine.AWF/align/alignedFullStack.mrc';

We create a data folder to contain the particles to be reconstructed on the fly:

targetDirectory = 'tempCrop.Data';

Angles actually used during the reconstruction are contained in the reconstructionTiltAngles' item of the output folder. Note that they might differ from the nominal angles, as some views might have been eliminated, both explicitely by the user, of by the program not being able to align all tilts.

angles =  'workflows/hivCommandLine.AWF/align/reconstructionTiltAngles.tlt';

The (binned) coordinates were measured in a binned tomogram, whose dimensionality is requested by dynamo_table_rec. We need to provide dynamo_table_rec with

visualizationTomogram =  'workflows/hivCommandLine.AWF/reconstruction/binnedReconstructionWBP.mrc';
sizeTomogram = dynamo_read_size(visualizationTomogram);
sizeTomogramFull = sizeTomogram*2^alignmentBinLevel;

The tomogram center can be shifted in order to account for a possible shift in the centre in the binned reconstruction that was used to crop the particles.

shiftTomogramCenter = [0,0,0];

We proceed now to the reconstruction. Our sid length should be 4 times the size we used before:

sidelengthFullSize = 128;

We use the applyRampFilter option in order to explicitely estate the kind of filter that we want during the back projection.

dynamo_table_rec(tblFull,source,angles,targetDirectory, sidelengthFullSize, 'applyRampFilter',1, 'sizeTomogram',sizeTomogramFull, 'shiftTomogramCenter',shiftTomogramCenter);

We can now check that the particles have been reconstructed in full size:

figure;dslices(targetDirectory,'proj','c');shg;
dslices on full sized gold beads