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

From Dynamo
Jump to navigation Jump to search
(31 intermediate revisions by 2 users not shown)
Line 1: Line 1:
This is the command line based version of the [[Creating Walkthrough on command line based tilt series alignment| GUI based alignment walkthrough]]
+
This is the command line based version of the [[Walkthrough on GUI based tilt series alignment | GUI based alignment walkthrough]]
 
 
  
 
= Create the workflow =
 
= Create the workflow =
Line 11: Line 10:
 
= Entering the data =
 
= Entering the data =
  
The <tt>u</tt> object contains several areas to interact with the workflow. They can be found by autocompletion using the <tt>tab</tt> key.
+
The <tt>u</tt> object contains several areas to interact with the workflow. They can be found by autocompletion using the <tt>tab</tt> key. Here, we will proceed step by step; remember that you can write all the command lines in a single <tt>.m</tt> script.
 +
 
 +
[[ File:alignCommandLine_autocompletionOfUserObject.png |thumb|center| 80px|autocompletion of user object]]
 +
[[ File:alignCommandLine_autocompletionOfFurtherFieldsOfUserObject.png |thumb|center| 200px|autocompletion of further fields of user object]]
  
 
== Basic data==
 
== Basic data==
  <nowiki>u.enter.tiltSeries(fileWorkshop);
+
The file of the tilt series can be defined somewhere else
u.enter.tiltAngles(-57:3:60);
+
<tt>file = '~/data/b001ts001.mrc'; </tt>
u.enter.discardedTiltIndices([1,2,40]); </nowiki>
+
 
 +
and then linked into the workflow:
 +
 
 +
  <nowiki>u.enter.tiltSeries(file);
 +
u.enter.tiltAngles(-57:3:60); </nowiki>
 +
 
 +
If you want to reject some of the hight tilts (or any other view that appears to have been damages):
 +
<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:
 
  <tt>u.enter.settingAcquisition.apix(2.7);</tt>
 
  <tt>u.enter.settingAcquisition.apix(2.7);</tt>
  
== Detection settings==
+
== Computation settings==
These here are the actual design decisions when running an alignment workflow:
+
We enable the use of parallel cores
 +
<tt>u.enter.settingComputing.parallelCPUUse(1); </tt>
 +
This will be using mainly for gold bead detection (cc computation) and reconstruction.
  
<code>u.enter.settingDetection.beadRadius(16);
+
To use all available physical cores use:
u.enter.settingDetection.maskRadius(28);</code>
+
<tt>u.enter.settingComputing.cpus('*');</tt>
  
 +
== Detection settings==
 +
These here are the actual design decisions when running an alignment workflow. The gold bead radius should have been  measured beforehand:
  
  <tt>u.enter.templateSidelength(64);</tt>
+
  <nowiki>u.enter.settingDetection.beadRadius(16);
 +
u.enter.settingDetection.maskRadius(28);
 +
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.
 
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>
 
  <tt>u.enter.settingDetection.detectionBinningFactor(1);</tt>
  
 +
== 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.
 +
<tt>u.area.indexing.step.tiltGapFiller.parameterSet.residualsThreshold(8);</tt>
 +
 +
= Running the workflow =
 +
We will skip CTF estimation and correction steps (both in the CTF area):
 +
<tt>u.run.all('noctf',1);</tt>
 +
 +
If should take around five minutes to run.
 +
 +
== Checking the results ==
 +
After completion, the <tt>view</tt> field allows you to check different output items for diagnose, as the tilt lines:
 +
 +
<tt>u.view.tiltLines.unaligned();</tt>
 +
 +
[[ File:alignCommandLine_viewOfTiltLines.png |thumb|center| 600px|view of tilt lines]]
 +
 +
or the reconstructions, as for example
 +
 +
<tt>u.view.reconstruction.binnedWBP</tt>
 +
 +
or
 +
 +
<tt>u.view.reconstruction.binnedSIRT</tt>
 +
 +
and aligned stacks themselves.
 +
 +
<tt>u.view.stack.alignedFull</tt>
 +
 +
 +
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;
 +
------------------------------------------------------------
 +
File with info:
 +
  workflows/workshop2.AWF/info/fitting.doc
 +
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 +
rms:  1.39
 +
psi: 85.10
 +
fit modus: normal
 +
total markers: 78
 +
total observations: 2590
 +
fit created at: 27-Aug-2019 18:01:40
 +
------------------------------------------------------------</nowiki>
 +
 +
or of the markers
 +
 +
<nowiki>>> u.info.markers;
 +
------------------------------------------------------------
 +
Info on working markers [used during detection, indexing refinement]
 +
  ... item saved and in memory
 +
Item in memory:
 +
Number of shapes          :  78
 +
Number of frames          :  40
 +
Number of observations    :  2590 (out of 3120)
 +
Number of empty frames    :  2
 +
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          :  78
 +
Number of frames          :  40
 +
Number of observations    :  2590 (out of 3120)
 +
Number of empty frames    :  2
 +
Number of complete frames :  0
 +
Number of complete shapes :  0
 +
 +
Item in disk:
 +
Markers in disk and currently in memory contain equal coordinates
 +
------------------------------------------------------------</nowiki>
 +
 +
= 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.
 +
 +
== Grep coordinates in the binned tomogram==
 +
 +
Open the binned reconstruction with <tt>dtmslice</tt>
 +
 +
<tt>dtmslice 'workflows/hivCommandLine.AWF/reconstruction/binnedReconstructionWBP.mrc';</tt>
 +
 +
[[ 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]]
 +
Now we create a new model of <tt>general</tt> type
 +
[[ 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
 +
[[ File:alignCommandLine_showAFullProjection.png |thumb|center| 600px|show a full projection]]
 +
[[ 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
 +
[[ File:alignCommandLine_rightClickOnAPoint.png |thumb|center| 600px|right click on a point]]
 +
Left clicking on it will provide a z coordinate to it
 +
[[ 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]]
 +
Keep clicking till you have a small set (no more than ten)
 +
[[ 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.
 +
[[ 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.
 +
[[ 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:
 +
<tt>write(temp_table,'binnedGoldbeads.tbl');</tt>
 +
 +
=== Cropping the gold bead particles ===
 +
 +
We can now crop the particles using the table file
 +
<nowiki>tomo = 'workflows/hivCommandLine.AWF/reconstruction/binnedReconstructionWBP.mrc';
 +
tbl  = 'binnedGoldbeads.tbl';
 +
targetBinnedParticles = 'binnedGoldBeads.Data';
 +
sidelengthBinned = 32;
 +
o = dtcrop(tomo,tbl, targetBinnedParticles, sidelengthBinned);</nowiki>
 +
 +
And then we visualise a projection along <tt>z</tt> of each cropped particle
 +
<tt>dslices(targetBinnedParticles,'proj','c');shg;</tt>
 +
 +
[[ File:alignCommandLine_dslicesOnBinnedGoldBeads.png |thumb|center| 600px|dslices on binned gold beads]]</nowiki>
 +
 +
== 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.
 +
 +
<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>:
 +
<nowiki>alignmentBinLevel = 2;
 +
tblFull = dynamo_table_rescale(tbl,'factor',2^alignmentBinLevel); </nowiki>
 +
 +
Now we fetch the full sized aligned stack stored inside the alignment folder:
 +
<tt>source  = 'workflows/hivCommandLine.AWF/align/alignedFullStack.mrc';</tt>
  
== Detection settings==
+
We create a data folder to contain the particles to be reconstructed on the fly:
<tt>u.enter.settingComputing.parallelCPUUse(1); </tt>
+
<tt>targetDirectory = 'tempCrop.Data'</tt>;
  
== Changing generic parameters ==
+
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.
Yo can find handles to the parameters of the individuals steps through auto ompletion on the <tt>area</tt>, then <tt>step</tt> items.
+
  <tt>angles =  'workflows/hivCommandLine.AWF/align/reconstructionTiltAngles.tlt';</tt>
  <tt>u.area.indexing.step.tiltGapFiller.parameterSet.residualsThreshold(8);</tt>
 
  
% computing parameters
+
The (binned) coordinates were measured in a binned tomogram, whose dimensionality is requested by <tt>dynamo_table_rec</tt>. We need to provide <tt>dynamo_table_rec</tt> with
 +
<nowiki>visualizationTomogram =  'workflows/hivCommandLine.AWF/reconstruction/binnedReconstructionWBP.mrc';
 +
sizeTomogram = dynamo_read_size(visualizationTomogram);
 +
sizeTomogramFull = sizeTomogram*2^alignmentBinLevel;</nowiki>
  
 +
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.
 +
<tt>shiftTomogramCenter = [0,0,0];</tt>
 +
We proceed now to the reconstruction. Our sid length should be 4 times the size we used before:
 +
<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.
 +
<nowiki>dynamo_table_rec(tblFull,source,angles,targetDirectory, sidelengthFullSize,....
 +
    'applyRampFilter',1,....
 +
    'sizeTomogram',sizeTomogramFull,....
 +
    'shiftTomogramCenter',shiftTomogramCenter);</nowiki>
  
= Running the workflow =
+
We can now check that the particles have been reconstructed in full size:
  <tt>u.run.all('noctf',1);</tt>
+
  <tt>figure;dslices(targetDirectory,'proj','c');shg;</tt>
 +
[[ File:alignCommandLine_dslicesOnFullSizedGoldBeads.png |thumb|center| 600px|dslices on full sized gold beads]]

Revision as of 17:16, 14 July 2021

This is the command line based version of the GUI based alignment walkthrough

Create the workflow

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

Here, we instruct Dynamo to skip opening the GUI. We also create the workflow in a different folder.

Entering the data

The u object contains several areas to interact with the workflow. They can be found by autocompletion using the tab key. Here, we will proceed step by step; remember that you can write all the command lines in a single .m script.

autocompletion of user object
autocompletion of further fields of user object

Basic data

The file of the tilt series can be defined somewhere else

file = '~/data/b001ts001.mrc'; 

and then linked into the workflow:

u.enter.tiltSeries(file);
u.enter.tiltAngles(-57:3:60); 

If you want to reject some of the hight tilts (or any other view that appears to have been damages):

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

Acquisition settings

Here you can enter parameters such as Cs, nominal defocus, etc. We will 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 should have been measured beforehand:

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

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.

u.enter.settingDetection.detectionBinningFactor(1);

Changing generic parameters

Yo can find handles to the parameters of the individuals steps through autocompletion on the area, then step items.

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

Running the workflow

We will skip CTF estimation and correction steps (both in the CTF area):

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

If should take around five minutes to run.

Checking the results

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

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

or the reconstructions, as for example

u.view.reconstruction.binnedWBP

or

u.view.reconstruction.binnedSIRT

and aligned stacks themselves.

u.view.stack.alignedFull


You can also use the autocompletion (in the Matlab version) of the info field in order to read the logs on the fit.

>> u.info.fit;
------------------------------------------------------------
File with info:
  workflows/workshop2.AWF/info/fitting.doc
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
rms:  1.39
psi: 85.10
fit modus: normal
total markers: 78
total observations: 2590
fit created at: 27-Aug-2019 18:01:40
------------------------------------------------------------

or of the markers

>> u.info.markers;
------------------------------------------------------------
Info on working markers [used during detection, indexing refinement]
  ... item saved and in memory
Item in memory:
Number of shapes          :  78
Number of frames          :  40
Number of observations    :  2590 (out of 3120)
Number of empty frames    :  2
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          :  78
Number of frames          :  40
Number of observations    :  2590 (out of 3120)
Number of empty frames    :  2
Number of complete frames :  0
Number of complete shapes :  0
 
Item in disk:
Markers in disk and currently in memory contain equal coordinates
------------------------------------------------------------

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.

Grep coordinates in the binned tomogram

Open the binned reconstruction with dtmslice

dtmslice 'workflows/hivCommandLine.AWF/reconstruction/binnedReconstructionWBP.mrc';

dtmslice of binned reconstruction

You will need to adjust the contrast:

dtmslice contrast control
dtmslice after contrast correction

Now we create a new model of general type

create a new model

We use the projections in order to quickly find the 3d positions of a set of tomograms

show a full projection
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 a point

Left clicking on it will provide a z coordinate to it

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

Keep clicking till you have a small set (no more than ten)

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 dtcrop on the binned tomogram.

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 as a file: write(temp_table,'binnedGoldbeads.tbl');

Cropping the gold bead particles

We can now crop the particles using the table file

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

And then we visualise a projection along z of each cropped particle

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

</nowiki>

Reconstruct particles on the fly

The command that reconstructs full sized particles from an aligned stack is dynama_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