View source for High resolution walkthrough with HIV1 (Tsinghua Spring School 2025)

[[Category:WalkthroughsSpecific]]
= The compact version of the HIV1 walkthrough  =

The compact version of the HIV1 walkthrough is meant to be a faster and less resource requiring version of the [[High resolution walkthrough with HIV1 | HIV1 walkthrough]].

The main scripts of the workshop will not be run, but only inspected. The script have already been run and their output has already been generated. As such, only the inspection of the results will be done. This will render the precess faster since no computing time will be required. It also means that the softwares IMOD, CTFFIND4 and RELION4 are optional.

<pre style="color: blue">As such, do NOT do the step highlighted in red. Only do the steps highlighted in blue</pre>


Chimera will also be used but it is best to use a local Chimera installation and download the the density maps you wish to inspect on the local machine.

= Goal of the walktrough =

This walktrough was designed to describe and discuss the steps necessary to go from the raw microscope data of a cryo-electron tomography experiment to a high-resolution average with the relative particles position and orientations.

In practice the EMPIAR dataset entry 10164 (https://www.ebi.ac.uk/empiar/EMPIAR-10164/) is used to analyse the structure of the GAG coat of immature HIV1 virus-like particles.

A set of common cryoem software tools will be used such as: MotionCor2, CTFFIND 4, IMOD, Dynamo, Chimera, RELION and Bsoft.

One of the goals is to show how these different softwares can be used in a synergistic approach to recover high resolution data. Additionally, the step-by-step description of the workflow is meant to allow users to be able to achieve the final result even by using other similar softwares that were not introduced in this walktrough.

[[File:General_pipeline.png|centre|thumb|800px| Cryoem general pipeline shown from grid preparation to high-resolution final average. The steps presented in the figure are, from left to right: grid preparation, data acquisition in the microscope, creation of the tilt series, contrast transfer function extraction, tilt series alignment, CTF correction, tomogram generation and modelling, subtomogram averaging and finally recovery of the high resolution average.]]

= Structure of the walkthrough =

The full HIV1 walkthrough pipeline has been split into multiple sequential steps. Each of these steps corresponds to a subsection of this wiki page. It is possible to reach an average with about ~4 A of resolution starting from 5 tilt series by following all of the steps described in this walkthrough. The only exception is the process to create a tilt series from the microscope movies, which is only discussed in the "Motion correction and dose weighting" subsection. The walkthrough is paired with a presentation giving context on each step. The presentation can be found in: <!--https://drive.google.com/file/d/1em6wWQUaB3mwIyL9LTnDyPG4Ks43quBK/view?usp=sharing-->{{Template:DownloadLinkHighResHIV1presentation}}<br>
<br>
<br>
However, depending on the hardware available, a full execution of this pipeline can easily take a week of computing time. Additionally, the generated data consumes considerable amounts of memory resources. As such, this walkthrough is has been constructed to be scalable. The full version of the walkthrough code package is divided into two sections: scripts and data. The following figure represents the package tree. The data part includes all of the pre-computed intermediate files.



[[File:Folder_structure_HIV1_v2.png|centre|thumb|800px| Structure of package containing the data and scripts for the HIV1 walktrough. The full version of the 'data' part of the package is 2.2 TB. The 'scripts' part is divided into steps. Note the correspondence between the folder names and the subsection names in this wiki page. Each step folder is divided into three parts: inputs, outputs and scripts. The input contains a set of paths to the where the input data is meant to be found in the 'data' part of the package. The output contains a set of path of where the generated files of this step should be moved to. The scripts includes the MATLAB codes that will be executed following the practical notes in the wiki page.]]

There are two main approaches on the implementation of walktrough with different scales of memory and computing resource requirements: with pre-computations or full execution. The two approaches are described in the following tables.

{| class="wikitable" style="width: 85%;"
! with pre-computations
|-
| Requires an already filled 'data' section.
|-
| Does not require large computing resources.
|-
| Intended for in-person workshop.
|-
| The 'data' section has been made available from a physical memory.
|-
| The 'script' section can be downloaded from <!--https://drive.google.com/file/d/1sM8OBQ7LF7OfhMPU0WwfJ31hVeXwTbCh/view?usp=drive_link-->{{Template:DownloadLinkHighResHIV1scripts}}.
|-
| The script generally do not need modification.
|-
| By default most scripts process only one tilt series. This is done to save computing time.
|-
| The scripts access their inputs data from the 'data' section.
|}

{| class="wikitable" style="width: 85%;"
! full execution
|-
| Does not require an already filled 'data' section.
|-
| Requires large computing resources.
|-
| Intended for individual practical use.
|-
| The 'data' section starts by being mostly empty. The 'data' section can be downloaded from <!--https://drive.google.com/file/d/1k4bLLswsacYdIGcNRc76YcbAiycxyPuX/view?usp=drive_link-->{{Template:DownloadLinkHighResHIV1lightData}}.
|-
| The 'script' section can be downloaded from <!--https://drive.google.com/file/d/1sM8OBQ7LF7OfhMPU0WwfJ31hVeXwTbCh/view?usp=drive_link-->{{Template:DownloadLinkHighResHIV1scripts}}.
|-
| The scripts will need modifications.
|-
| The scripts need to be modified to process all of the tilt series. The scripts are already structured to simplify the needed modifications.
|-
| The scripts access their inputs data from the 'data' section. As such, the 'data' section needs to be filled with the outputs of the scripts after each step.
|}

As mentioned previously, the "full-execution" approach requires some limited script modifications. The MATLAB scripts have been generally structured so as to present three parts. The first part consists of the paths to the input and output files. The second part is the section of parameters. The third part is a for loop on the input tilt series. As such, the typical modifications necessary for the "full-execution" approach are limited to the first part, ensuring that the inputs for all five tilt series are given.<br>
<br>
Additionally, in the "full-execution" approach, the results of the scripts are used to populate the empty 'data' section so that the following scripts may be able to access their input files. This means that the files need to be positioned in the correct folder tree location and with a consistent file name. To help with this process, each step has an output text file in its folder with the path and file names that need to be applied to ensure consistency.
<br>

The approach used in this workshop will be the one with pre-computations. In this configuration, two new folders are provided. A reduced version of the data is stored in a folder called '''lightweight_data'''. This folder follows the same structure as the data folder described in the previous figure, but it contains only a subset of the files so as to reduce the memory requirements. It is designed so that the task described in this walkthrough remain possible. A second complementary folder is '''scripts_pre_executed'''. This contains the same files as the 'script' section, but all of these files have been completed and executed. Additionally, the output files of the tasks are available for inspection.

<pre style="color: blue">Pairing '''lightweight_data''' and '''scripts_pre_executed''' is sufficient to complete the tasks of this walkthrough that have been written in blue.</pre>

<br>
<br>
{| class="wikitable" style="margin-left: auto; margin-right: auto; border: none;"
|+ Time requirements for each step
|-
! Step !! Computing time !! Interactive time !! Workshop time !! Notes
|-
! Recovery of acquisition parameters 
| seconds 
|colspan=2 style="text-align: center;"| 5 min 
| 
|-
! Motion correction 
| 1 min per tilt series with GPU 
|colspan=2 style="text-align: center;"| minutes 
| Not done in this walkthrough
|-
! CTF estimation 
| 5 min for all tilt series 
|colspan=2 style="text-align: center;"| minutes 
| Using CTFFIND4
|-
! Automatic tilt series alignment: alignment 
| 8 min per tilt series 
|colspan=2 style="text-align: center;"| minutes 
| Executed by default only on the first tilt series
|-
! Automatic tilt series alignment: dose weighting 
| 3 min 
|colspan=2 style="text-align: center;"| minutes 
| Executed by default on all tilt series
|-
! Automatic tilt series alignment: CTF correction 
| 1 hour 30 min per tilt series without GPU, 5 min per tilt series with 1 GPU 
|colspan=2 style="text-align: center;"| minutes 
| Using IMOD, executed by default only on the first tilt series
|-
! Automatic tilt series alignment: tomogram reconstruction 
| 12 min per tilt series 
|colspan=2 style="text-align: center;"| minutes 
| Executed by default only on the first tilt series
|-
! DMMPP: preparation and management of data 
| 5 min 
| 3 min 
| 15 min 
|
|-
! DMMPP: particle picking 
| - 
| 1 hour 30 min 
| 10 min 
|
|-
! DMMPP: particle extraction 
| 1 hour 30 min 
| - 
| 15 min 
|
|-
! StA: creation of initial reference 
| 1 hour 30 min 
| - 
| 15 min 
|
|-
! StA: extraction of first per tomogram average 
| 15 min (CPU) + 1 hour 30 min (8 GPUs) per tilt series 
| 15 min 
| 15 min 
| Use of UCSF Chimera
|-
! StA: computing an average of averages 
| 15 min (CPU) 
| - 
| 15 min 
|
|-
! StA: extraction of second per tomogram average 
| 30 min (CPU) + 3 hours 30 min (8 GPUs) per tilt series 
| 25 min 
| 15 min 
| Use of UCSF Chimera
|-
! StA: extraction of coarse half maps 
| 22 min (CPU) + 48 hours (8 GPUs) 
| - 
| 15 min 
|
|-
! StA: extraction of refined half maps 
| 1 hour (CPU) + 24 hours (8 GPUs) + 1 hour (CPU) per tilt series 
| - 
| 15 min 
|
|-
! FSC computation 
| seconds 
|colspan=2 style="text-align: center;"| 15 min 
| Using RELION
|-
! RELION conversion 
| seconds 
|colspan=2 style="text-align: center;"| 5 min 
| Using RELION
|-
! RELION refinement 
| >40 hours using 8 GPU 
|colspan=2 style="text-align: center;"| minutes 
| Using RELION, split into generation of a RELION template (~1h), extraction of RELION particles (~1h) and then the execution of the RELION refinement
|-
|}
<br>
<br>
<br>
----

= Recovery of acquisition parameters =

A number of important characteristics of a cryo-electron microscopy experiment are necessary to properly process the data.<br>
As such, this first step of the walkthrough is a guided process to recover the relevant information for this datasets.<br>
In general, other that the sample movies themselves, the required information can be split in two parts: the characteristics of the microscope and the characteristics of the experiment.<br>

The needed microscope characteristics are:<br>
- voltage.<br>
- spherical aberration.<br>
- amplitude contrast.<br>
- gain reference.<br>
The needed experiment characteristics are:<br>
- pixel size.<br>
- image size of the tilt series.<br>
- dose.<br>
- tilt scheme.<br>
<br><br>
The walkthrough is based on the EMPIAR dataset entry 10164. The corresponding paper is:      Florian K. M. Schur et al., An atomic model of HIV-1 capsid-SP1 reveals structures regulating assembly and maturation.Science353,506-508(2016). DOI:10.1126/science.aaf9620. and the EMD deposition is EMD-4015, https://www.ebi.ac.uk/emdb/EMD-4015?tab=experiment.<br><br>
The microscope parameters are describe in the paper. The microscope that was used is a FEI TITAN KRIOS with an acceleration voltage of 300 kV. The nominal spherical aberration is of 2.7 mm. The amplitude contrast in cryo-electron tomography is ~10%. In the following the value of 7% is picked.<br>
The gain reference is not available for this experiment. However, the gain reference is only useful during motion correction of the movie. In this walkthrough we will start with already motion corrected and stacked tilt series. As such, the gain reference is not required for this walkthrough. In any case, more information on motion correction and on the gain reference is provided in the next section.<br><br><br>
The experiment characteristics still need to be found. Firstly:


<pre style="color: blue">access the EMPIAR dataset entry 10164 by accessing the link https://www.ebi.ac.uk/empiar/EMPIAR-10164/</pre>
<pre style="color: blue">Go to the "Image set" section of the EMPIAR page.</pre>

[[File:image_set_section.png|centre|thumb|600px| Screenshot of the EMPIAR "Image set" section.]]

The subset of tilt series that will be used in the walktrough are TS_01, TS_03, TS_43, TS_45 and TS_54.<br><br>
The dataset in the EMPIAR consists in movies. As such, the '''Frames per image''' of 8 means that each movie is a set of 8 frames. In addition, we can note that '''Image size''' is (7420,7676). This refers to the size of a movie. In summary, if one of the movies in this EMPIAR dataset would be stored as a matrix, it would have the shape 7420x7676x8.<br><br>

The pixel size can be found by inspecting '''Pixel spacing''', which get a pixel size of 1.35 A, and in the '''Details''' section the information of the dose can be found.
<pre style="color: blue">Note the dose and pixel size of the tilt series TS_01, TS_03, TS_43, TS_45 and TS_54.</pre>
The size of the image size of the tilt series can now be found.
<pre style="color: blue">Update the variable root to the path the the data folder.</pre>
This can be done by running <tt>root = '/GLOBAL/PATH/TO/DATA/lightweight_data/';</tt>

'''Change your MATLAB working folder to the Experimental_parameters subfolder of scripts.'''

Then run:
<pre style="color: blue">dynamo_read_mrc_size([root,'data/TS_not_dose_weighted/01.mrc'])</pre>

This will return the dimensions of the tilt series TS_01. The order is size in x, size in y and size in z. In a tilt series the z dimension is the tilt. Now run:
<pre style="color: blue">dynamo_read_mrc_size([root,'data/TS_not_dose_weighted/43.mrc'])</pre>
You will note that the tilt series TS_01 is one tilt shorter than the tilt series TS_43. By inspecting the tilt angles files:
<pre style="color: blue">open([root,'data/angles_files/43.rawtlt'])</pre>
<pre style="color: blue">open([root,'data/angles_files/01.rawtlt'])</pre>
you can see what tilt angle is has been removed in the tilt series TS_01. Additionally you can:
<pre style="color: blue">Annotate the image size of the tilt series.</pre>
You can note that the x and y sizes of the tilt series are half of the size of the movies that was indicated in the EMPIAR
'''Image size''' section. This is because during motion correction the micrographs are cropped to avoid aliasing. More information on this can be found in the paper: Electron counting and beam-induced motion correction enable near-atomic-resolution single-particle cryo-EM, Li, Xueming, Mooney, Paul, et al. Nature Methods, 2013, https://doi.org/10.1038/nmeth.2472<br><br>

The final experiment characteristic needed is the tilt scheme. This information is generally already known, but it can also be recovered from the microscope SerialEM logs. These logs are provided for this dataset in the EMPIAR, but only partially: the logs refer to a binned output and do not present information on the dose. This is why the dose information was presented in the '''Details''' section instead.<br><br>
One of the these serialEM logs has been made available in the data folder. You can open it:
<pre style="color: blue">open([root,'data/SERIALEM_doc/TS_01.mrc.mdoc'])</pre>

[[File:serial_log.png|centre|thumb|600px| Screenshot of the SerialEM binned log.]]

The log file starts with an header of six lines with some general characteristics. The binned SerialEM pixel size is consistent with eh pixel size noted in the EMPIAR.

<pre style="color: blue">Note the Tilt axis angle</pre>
The tilt axis angle in SerialEM is meant to be the angle between the axis of rotation of the sample during tilting and the axis of the detector.

The log file is then split in section by the line [ZValue=????]. This line separates the different microscope acquisitions (i.e. each movie generation). The number of ZValue indicates the order of the acquisitions. Taking the screen shot as an example: the first block is about the first movie taken because the Zvalue is 0. The tilt angle was 0, that is the first acquisition was take with a flat sample. The second block is about the second movie because the Zvalue is 1. The second tilt angle was +3.<br><br>
This means that by inspection the log the tilt scheme can be recovered. The custom function parse_metadata will do this task. Run:<br><br>

<pre style="color: blue">[tiltAngles, indices] = parse_metadata([root,'data/SERIALEM_doc/TS_01.mrc.mdoc']);</pre>

And you can inspect the acquisition scheme by running:
<pre style="color: blue">figure;
scatter(indices,tiltAngles,'.');
xlabel('order of exposure');
ylabel('tilt angle');
title('Exposition scheme');</pre>
[[File:Exposition_scheme_HIV1.png|centre|thumb|600px| Dose-symmetric tilt scheme used for the HIV1 dataset.]]

We can note that the tilt scheme is a dose-symmetric with steps of 3 degrees like: 0,+3,-3,-6,+6,+9,-9,...<br><br>
A copy of the tilt scheme in an appropriate format for the following steps is already computed and can be inspected by:

<pre style="color: blue">n_of_preexposures = csvread([root,'data/Microscope_parameters/order_list.csv']);</pre>

= Motion correction and dose weighting =
== Initial comments ==

The output of a microscope after a tilt series acquisition will be a gain image and a movie for each tilt. In practice, often a single gain file is used for all the tilt series experiments done in a microscope session. The following figure present an example of a gain reference file from a different dataset.

[[File:an_example_of_a_gain_reference.png|centre|thumb|600px| An example of a gain reference file of a K3 detector inspected with the 3dmod command of IMOD.]]

Some more information on reference gains can be found at the links:<br>
https://zmb.dozuki.com/Guide/TEM+-+TFS+Titan+Krios+G3i+-+2a:+Collect+Gain+References+for+the+K3/183

https://relion.readthedocs.io/en/release-3.1/Reference/MovieCompression.html

The first steps that can be done are the correction of the sample motions caused by the electron beam exposure and the gain correction.

In this case, MotionCor2 will be used (https://emcore.ucsf.edu/ucsf-software).
MotionCor2 will use the movies and the gain to generate a motion-corrected micrograph per each movie. MotionCor2 should not be used to do the dose weighting at the same time as the motion correction. This is because the implementation of the dose weighting procedure in Motioncorr2 has been designed for single particle experiments, not cryo-electron tomography.<br> 

Dose weighting is meant to generate a more strict filtering of high-frequency components of the power spectra in frames that have been exposed to large doses of electrons per square angstroms. This means that the order of exposition of the tilts is crucial. This process will be done in this walkthrough after the alignment step.

Some steps in the following walkthrough are more adapted to either the dose weighted micrographs or the not dose weighted micrographs, so it is often good to keep both available.

Finally, the micrographs should be stacked in order of tilt (generally from negative angles to positive angles).

== Input ==

- a set of microscope movies.<br>
- a gain.<br>
- microscope parameters: pixel size, voltage, spherical aberration and amplitude contrast. (SPECIFY)<br>
- the tilt angle of each movie.<br>
- the order of tilt exposition.<br>

== Practical notes ==

Here is an example of a MotionCor2 call that will generate both a dose weighted and a not dose weighted micrograph:
<code>
'motioncorr2 -InMrc movie.mrc -OutMrc tilt.mrc -Patch [5 5] -Iter 30 -bft 200 -Tol 0.5 -Kv 300 -PixSize 2.73 -FmDose FMDOSE -InitDose INITDOSE -FtBin 1'
</code>

where INITDOSE will depend on the number of previous tilts in the exposition scheme. As mentioned before, the dose weighting results of Motioncorr2 are meant to be used in single particle experiments. As such, in this case, the dose weighting will instead be done during the tilt series alignment step using a MATLAB script.<br><br>
Note that the previous motioncorr2 call did not use a gain. If you want to use the gain you can add it with the <code>-GAIN PATHTOGAIN</code> flag. Sometimes the gain needs to be flipped or rotated depending on what where the setting of SerialEM during the acquisition and this can be done with the <code>-RotGain</code> and <code>-FlipGain</code> flags. Note also that if Motioncorr2 cannot find or use the gain file provided, it will not terminate the motion correction and will instead continue computing without gain correction.<br>

More practically, dynamo has a set of common preprocessing scripts that can be downloaded at: https://github.com/C-CINA/TomographyTools/tree/master/matlab_files

With a particular note on the script preprocessing_script.m, motioncor2Wrapper.m and applyExposureFilter.m.

The script applyExposureFilter.m can be used to apply dose weighting on a micrograph that was already motion-corrected.

== Outputs ==

- the dose weighted tilt series.<br>
- the not dose weighted tilt series.<br>

= CTF estimation and correction =
== Initial comments ==

The CTF of each micrograph can be estimated and corrected in different ways. In this case, CTFFIND4 (https://grigoriefflab.umassmed.edu/ctffind4) was used to estimate the CTF parameters and IMOD was used to correct it using the ctfphaseflip function of IMOD (https://bio3d.colorado.edu/imod/doc/man/ctfphaseflip.html).

CTF can be tricky in cryoem because of the low dose per micrograph and the fact that the sample is tilted, causing the CTF to be different in different regions of the micrograph. In addition, the thickness of the sample means that even the CTF estimated in a region of the micrograph is in actuality caused by the entire thickness of the sample.

These problems can apply both to the CTF estimation and the CTF correction. Some strategies that can be used to limit these problems are going to be presented in the following.

CTF estimation:
It can be noted that the CTF can be done either on the dose weighted or on the not dose weighted tilt series. It should be more proper to use the not dose weighted micrographs, but in practice using the dose weighted micrographs can be more straightforward. In our case, CTFFIND4 was used on the not dose weighted tilt series.

CTF correction:
ctfphaseflip locally corrects the CTF (that was estimated over the entire micrograph) by correcting the local defocus in each area using the geometric Z-offset generated by the tilt. However, this means that ctfphaseflip requires an aligned tilt series. Indeed, the automatic fiducial-based tilt series aligner of dynamo can be used to correct the CTF after generating the aligned tilt series (see following steps).

[[File:CTF_problems_in_cryoem.png|centre|thumb|800px| Scheme showing some difficulties in handling the CTF in Cryo-electron tomography. Both the tilt of the sample and the sample thickness cause each particle to be affected by a locally different CTF. In the scheme each local CTF is noted as an 'overall' CTF(z) in z<sub>R</sub> in addition dz corrections. Note that z<sub>R</sub> is a priori unknown.]]

== Input ==

for the CTF estimation:<br>
- the not dose weighted tilt series.<br><br>


for the CTF correction:<br>
- the aligned, dose-weighted tilt series.<br>
- the CTF parameters of the tilt series.<br>
- the tilt angles.<br>
- the microscope parameters.<br>

== Practical notes ==

At this time CTF estimation part will be executed. The CTF correction will be done during the alignment and reconstruction step.

The relevant information will be found in the 'CTF_estimation' folder.

The script <code>use_CTFFIND4_to_estimate_CTF.m</code> will use <code>ctffind4Wrapper.m</code> to estimate the CTF parameters using CTFFIND4.

Note that this is done by making a system call for the CTFFIND4 software. As such, it is necessary to know where the software has been installed.
Line 178 of <code>ctffind4Wrapper.m</code> should be changed to be able to reach the CTFFIND4 executable. An easy way to check if the executable can be reached is to try  <code>which ctffind4</code> in the terminal. Note that the CTFFIND4 wrapper has been created for CTFfind version 4.1.5. Using a different CTFfind version may require some limited modifications of the script.

<pre style="color: blue">Update the variable root to the path the the data folder.</pre>

'''Change your MATLAB working folder to the CTF_estimation/scripts subfolder.'''

Find the CTFFIND4 installation:

<pre style="color: red">which ctffind4</pre>

Once you have the path to the CTFFIND4 installation you can update the script.

<pre style="color: red">change line 178 of ctffind4Wrapper from:

    fprintf(fileID,'ctffind4 <<EOF\n');    % path to ctffind4 installation (adapt if necessary)

to:

    fprintf(fileID,'PATHTOMYCTFFIND4EXECUTABLE <<EOF\n');

where PATHTOMYCTFFIND4EXECUTABLE was just identified.
</pre>

Now the CTF microscope parameters need to be given to the script.

<pre style="color: red">
Update the script use_CTFFIND4_to_estimate_CTF.m in the parameters section.
</pre>

Now the script can be run:
<pre style="color: red"> run use_CTFFIND4_to_estimate_CTF.m</pre>

CTFFIND4 will generate three files for each tilt series: a 'XXX_diag.mrc' file, a 'XXX_diag.txt' file and a 'XXX_diag_avrot.txt' file. Some notes on the role these files:<br>

- 'XXX_diag.mrc': this file will present visually the estimated CTF on the power spectra of each tilt. The dynamo command <code>dshow</code> can be convenient to quickly inspect the results<br>
- 'XXX_diag.txt': this file will annotate the CTF parameters estimated for each tilt.<br>

Finally, the CTFFIND4 result are converted in a '.defocus' file compatible with IMOD. This is done because IMOD will be used later to correct the CTF. This conversion is done in <code>use_CTFFIND4_to_estimate_CTF.m</code> by calling the function <code>ctffind2imod.m</code><br><br>

Inspect a CTFFIND4 result with <code>dshow</code>:

<pre style="color: blue">dshow('01_diag.mrc')</pre>

You will get something like this figure:


[[File:Ctffind4_on_unaligned_TS.png|centre|thumb|600px| Inspection using <code>dshow</code> of a 'XXX_diag.mrc' file generated by CTFFIND4 on an unaligned tilt series. Note the positioning of the calculated CTF zeros on the power spectra.]]

The Z dimension scroller correspond to moving between different tilts, and correspondingly to different relative sample thickness as well as accumulated doses.<br><br>
The parameters found by CTFFIND4 can be inspected by opening the XXX_diag.txt file:
<pre style="color: blue">open('01_diag.txt')</pre>

== Outputs ==

for the CTF estimation:<br>
- CTF parameters of the tilt series.<br><br>


for the CTF correction:<br>
- the aligned, dose-weighted, CTF-corrected tilt series.<br>

= Alignment and reconstruction =
== Initial comments ==

Each tilt series needs to be aligned to then be able to reconstruct them in tomograms. In this approach a fiducial-based automatic alignment is used. The general concept of fiducial-based alignment is shown in the following figure. In particular, the alignment procedure will describe the alignment as a general rotation angle for the entire tilt series and then additional shifts in x and y directions for each tilt. Generally, the rotation angle can be seen as the angular difference between the axes of the detector and the axis of rotation of the sample. The x and y shifts try to model the imprecision in the sample movement during each tilt acquisition.



[[File:Concept_of_fiducial_alignment.png|centre|thumb|800px| Basic concept of fiducial-based tilt series alignment. Tilt series become naturally misaligned during the process of tilt series acquisition in the microscope (left). However, the detection of key fiducials in the micrographs can be used to determine their position in the 3D sample (center). These fiducial positions can the be used to correct the tilt series alignment (right). An automatic and robust approach to achieve this goal has been implemented in dynamo.]]


Once the tilt series has been aligned, the ctfphaseflip correction can be applied as described in the previous section.

Finally, the aligned, dose-weighted, CTF-corrected tilt series can be used to reconstruct the tomograms. Weighted backprojection is the preferred method to generate the tomograms that will be used to generate subtomograms for the STA. SIRT can be used to generate tomograms with better contrast to more easily annotate the data with dynamo models.


The automatic fiducial-based workflow (AWF) in dynamo can be used to do: the alignment, the estimation of CTF parameter, the CTF correction and finally the tomogram reconstruction. The AWF can be used either through the GUI or programmatically.

From the dynamo wiki:

https://www.dynamo-em.org/w/index.php?title=Walkthrough_on_GUI_based_tilt_series_alignment

https://www.dynamo-em.org/w/index.php?title=Walkthrough_on_command_line_based_tilt_series_alignment

https://www.dynamo-em.org/w/index.php?title=Programmatic_control_of_alignment_and_reconstruction_workflows
<br>

The AWF will generate an AWF folder in which the results will be saved.

https://www.dynamo-em.org/w/index.php?title=Tilt_series_alignment_workflow_folder
<br>

The quality of an alignment is difficult to measure, but some general considerations can be made:

https://www.dynamo-em.org/w/index.php?title=Considerations_for_tilt_series_alignment_in_IMOD<br>
https://www.dynamo-em.org/w/index.php?title=Tomogram_quality_assessment
<br>

In this practical setup, a script will be used to automatically and programmatically align all tilt series, but the CTF handling and the reconstruction will be done separately. This choice was made to present in mode detail how to achieve a CTF correction and to mention how other softwares (like IMOD) can be used to reconstruct TS that were aligned using dynamo. 

Note the tilt angle conventions with CTF:<br>
https://www.dynamo-em.org//w/index.php?title=Tilt_angle_convention

== Input ==

- (not dose-weighted) tilt series.<br>
- gold beads parameters.<br>
- CTF parameters.<br>

== Practical notes ==

The relevant information will be found in the 'Alignment_and_reconstruction' folder. 

The work will be split in four parts: the alignment, the dose weighting, the CTF correction and finally the reconstruction.

'''Change your MATLAB working folder to the Alignment_and_reconstruction/scripts subfolder.'''

=== Practical notes: alignment ===
The alignment will be generated using the script <code>make_AWF.m</code>. This script will programmatically initialize and run the dynamo automatic fiducial-based tilt series alignment but skipping the CTF handling and the reconstruction.

The algorithm only requires as parameters the gold bead radius, the gold bead mask, the patch sidelength and the binning to apply. These parameters are determined by inspecting a tilt series. Open the first tilt series by running:

<pre style="color: blue">dshow([root,'data/TS_not_dose_weighted/01.mrc']);</pre>

This will open the <code>dshow</code> GUI representation of the tilt series. 

[[File:dshow_HIV1_raw.png|centre|thumb|600px| <code>dshow</code> GUI representation of the tilt series]]

Using the GUI tools:

<pre style="color: blue">Use the GUI binning tool to render the gold beads easier to see</pre>

The binning applied is purely visual. The dimensions of the tilt series does not change and this means that the measurement with the distance tool are in full-scale pixels.

<pre style="color: blue">Use the measuring tool to estimate the radius of the gold beads</pre>

This procedure allows to estimate the radius of the gold bead and the binning to apply. This is because the AWF will initially detect the gold beads using template matching on the visually binned tilt series just as it was done manually now.<br>

The next figure present visually the relationship between parameters gold bead radius, the gold bead mask radius and the patch sidelength.

[[File:AWF_GB_definitions.png|centre|thumb|600px| Schematic representation of the AWF gold bead parameters.]]

Generally we suggest to set the mask radius parameter as about 50% bigger than the beads radius and then set the patch sidelength as the smallest that can contain the mask diameter.

<pre style="color: blue">Use this rule of thumb to estimate the gold bead radius and sidelength parameters</pre>

The required data to set up an AWF are: the unaligned TS, the tilt angles, the pixel size and the gold bead parameters. In addition, the dimensions of the of the reconstructed tomogram are used in the AWF. However, in this case, these values will act as fill-ins because the actual reconstruction will be done separately.

<pre style="color: red">Update the parameters of the script make_AWF.m (except zFullSize_a, as said before)</pre>

The script will do the alignment only for the first tilt series because of time limitations. Precomputed alignment are proveided in the data/AWF_projects folder.
Once the alignment as been completed, one (or all five if the for loop at line 33 has been changed to <code>for 1:n_tomos</code>) .AWF project will have been generated. In particular, the aligned tilt series can be found in '.AWF/align/alignedFullStack.mrc'. The actual alignment parameters can be found in the files 'stackAligner.dyn' and 'stackAligner.star'. The information in the two files is essentially the same, but in a human-readable format in the .star file and a dynamo-readable format in the .dyn file. Inspect the alignment parameter by running:

<pre style="color: blue">open('HIV_AWF_fortestsdata/TS_not_dose_weighted/01.AWF/align/stackAligner.star')</pre>

You will see the alignment parameters shift in x, shift in y and psis angular correction. Psis is the angle between the tilt axis of the sample and the axis of the camera.

<pre style="color: blue">compare the psis value found in the alignment with the SerialEM tilt axis angle</pre>

The difference of angle should be less than one degree. By default AWF scans the psis angle at intervals of 0.1 degrees. It is sometimes possible that there might be an offset of 180 degrees in the estimation of the tilt axis angle in different softwares depending on how the axis are defined.<br><br>
The overall alignment can be inspected in multiple ways. If you are using the pre-computed results, you can load the alignment pipeline object with <code>u = dread('HIV_AWF_fortestsdata/TS_not_dose_weighted/01.AWF');</code>, otherwise it will already be in your MATLAB workspace. First, to see the information about the markers found, run:
<pre style="color: blue">u.info.markers;</pre>
You should see in the log that there are no empty frames. The number of empty frames is the number of frames removed by the AWF and those frames will not appear in the aligned tilt series.<br><br>
Then you can inspect the fitting log. Run:
<pre style="color: blue">u.info.fit;</pre>
You will see the same psi as in the one in <code>stackAligner.star</code>. The value of the rms is also provided.<br><br>
A visual inspection of the aligned tilt series can alos be done. Run:
<pre style="color: blue">u.view.stack.alignedFull();</pre>
The gold beads should appear like they are moving horizontally between tilt.

[[File:ali_ts1.png|centre|thumb|600px| Aligned tilt series]]

<pre style="color: blue">Increase the thickness value of the GUI to more than the number of tilt to have a superposition of all aligned micrographs</pre>

[[File:av_ali_ts1.png|centre|thumb|600px| Aligned tilt series using the superposition trick]]

You can see the original tilt lines with:
<pre style="color: blue">u.view.tiltLines.unaligned;</pre>

[[File:tilt_lines.png|centre|thumb|600px| Tilt lines before alignment]]
It can be noted that in this tilt series the tilt lines are mostly vertical before alignment and horizontal after alignment.

=== Practical notes: dose weighting ===

The dose weighting will be done using the script <code>apply_DW.m</code>. This script will dose weight a TS by calling the function <code>applyExposureFilter.m</code>. This function implement a dose filtering based on the paper "Measuring the optimal exposure for single particle cryo-EM using a 2.6 Å reconstruction of rotavirus VP6", Grant & Grigorieff, elife 2015.

The required data to do the dose filtering are the aligned tilt series, the dose per tilt and the tilt scheme. The tilt scheme information has been provided as a .csv file and merely notes the order of exposition of all tilts.

<pre style="color: red">Update the parameters of the script apply_DW.m</pre>

The <code>apply_DW.m</code> script will generate .mrc files of dose weighted tilt series. Run the script:
<pre style="color: red">run apply_DW.m</pre>

You can compare the tilt series before and after dose weighting by running:
<pre style="color: blue">dtmshow({'01_DW.mrc',[root,'data/AWF_projects/HIV_ali_AWF/TS_01/01.AWF/align/alignedFullStack.mrc']});</pre>
This will open two synchronized representation of the tilt series before and after dose weighting. The difference will be more obvious in tilt that have been subjected to larger accumulated doses.

[[File:notDW_vs_DW.png|centre|thumb|600px| <code>dshow</code> Before and after dose weighting]]

=== Practical notes: CTF correction ===

The CTF correction will be done using the script <code>correct_CTF.m</code>. This script do a system call to the <code>ctfphaseflip</code> function of IMOD. The tested version of IMOD is 4.11.12.

<pre style="color: red">Update the line 26 of the script correct_CTF.m from:

    system(['ctfphaseflip -InputStack ',path_to_dwTS{tiltserie},...

to:

    system(['PATHTOYOURIMOD/ctfphaseflip -InputStack ',path_to_dwTS{tiltserie},...
</pre>

a way to check that this worked is to type:

<pre style="color: red">system(['PATHTOYOURIMOD/ctfphaseflip'])</pre>

if the path is correct you should get a description of the function inputs.<br><br>
ctfphaseflip is much faster with a GPU. If you have one available you can simply add the -gpu 0 flag.<br>

The required data to do the dose filtering are the aligned and dose-weighted tilt series, the tilt angles, the CTF parameters (in .defocus format) and the microscope parameters.

<pre style="color: red">Update the parameters of the script correct_CTF.m</pre>

In general, it is possible that alignment procedure may have discarded some frames (for example, black frames). In this case, it is important to ensure that the tilt series .mrc file, the CTF .defocus and the angle file match.

The <code>correct_CTF.m</code> script will generate .mrc files of CTF-corrected tilt series. For time constrains, this is done only on the first tilt series. As usual, precomputed results can be found in the folder data/CTF_corrections.<br><br>

You can inspect the tilt series before and after CTF correction by typing:
<pre style="color: blue">dtmshow({'stackCTF_DW_ali_01.mrc',[root,'data/TS_dose_weighted_after_alignment/alignedFullStack_TS01_DW.mrc']});</pre>
An easy way to check if CTF correction was applied to inspect the black areas cause by the alignment. A CTF-corrected tilt series will present a wave pattern in these areas.

[[File:DW_DW_and_CTF.png|centre|thumb|600px| <code>dshow</code> Before and after dose CTF correction]]

=== Practical notes: reconstruction ===

The reconstruction will be done using the script <code>reconstruct_tomograms.m</code>. This script will use the dynamo function that reconstructs a tomogram using the weighted back-projection approach. It is the same function that would have been used had we wanted to do the reconstruction directly in the AWF project.

The required data to do the reconstructions are the dose-weighted, CTF corrected aligned TS, the tilt angles and the requested Z thickness of the tomogram.

<pre style="color: red">Update the parameters of the script reconstruct_tomograms.m with a thickness of 2000</pre>

The reconstruction is equivalent to an IMOD reconstruction using the <code>tilt .... --ROTATE 90</code> command of a tilt series aligned with the IMOD-compatible alignment parameters. These alignment parameters were generated in the AWF folder as the files 'stackAlignerImod.dyn' and 'stackAlignerImod.star'.<br><br>
<pre style="color: red">Note that the reconstruction of a full-size tomogram is resource intensive, especially in memory. In this case, the expected size of a single tomogram is about 100 GB.</pre> All the tomograms have been pre-computed and stored in the data folder. If you wish to backproject a tomogram it may be done by running:
<pre style="color: red">run reconstruct_tomograms.m</pre>
For time constrains, this is done only on the first tilt series. As usual, precomputed results can be found in the folder data/tomograms.<br><br>
You can check the dimensions of the tomogram with: 
<pre style="color: red">dynamo_read_mrc_size(['recCTF_DW_ali_01.mrc'])</pre>
and you can inspect the tomogram with: <code>dshow(['recCTF_DW_ali_01.mrc'])</code>.
However, this is a large full-scale tomogram, so it will take a long time to load. As such, it might be best to inspect only binned tomograms. This will be done in the next section.

== Outputs ==

- tomograms

= Dynamo models and particles =
== Initial comments ==
After the tomographic reconstruction, the user needs to identify the target of the biological investigation. In general, one would like to have a vectorial representation that captures the geometry of the particle distribution. This process is known as particle picking (sometimes referred to as volume annotation).

Here, the focus is on the main structural component of HIV-1: the polyprotein Gag. Since, for instance, its ordered cleavage causes a dramatic structural rearrangement of the virus to the mature, infectious virion. In the immature virus particle, Gag is radially arranged (left, figure below). The particle segmentation model capturing such a surface is reasonably spherical/ellipsoidal. Therefore, the choice should be a dipole. This Dynamo model class is very convenient since it allows the user on-screen interaction to be reduced in directly clicking the centers that roughly correspond to the presence of individual particles.

[[File:Schematic-representations-of-typical-immature-viral-morphologies-dmmpp.png|centre|thumb|600px| Schematic representations of typical immature viral morphologies. Main components of the polyprotein Gag, consisting of the MA (matrix), CA (capsid), NC (nucleocapsid), and p6 domains as well as the spacer peptides SP1 and SP2. Source:  https://www.pnas.org/doi/abs/10.1073/pnas.1811237115]]

Once the dipole sets for the particle's global region are ready, we proceed with "seed oversampling". This technique involves creating cut points along the models to extract smaller sub-volumes. Understanding the underlying distribution, especially the degree of homogeneity and symmetry, is key to selecting a sensible oversampling configuration; it is key to success on the road to high resolution. Although, at this stage of the project the positions of these cropping points are not supposed to be realistically centered on a copy of the polyprotein of interest.

At the end of the extraction process, you should have a crop table, a set of particles, and a template (i.e. the average of the subvolumes).

== Input ==

- tomograms

== Practical notes ==
Prior to a more practical depiction, it is essential to pause and become aware of the data organisation and naming convention. In this regard, the user is advised to have a highly consistent data structure. The project is about to grow in size. Therefore, errors probably mean lots of hours of work lost.


[[File:motivation-file-naming-convention-dmmpp.png|centre|thumb|300px| Motivation for establishing a file naming convention. A file naming convention is a framework for naming your files in a wav that describes what they contain and how the relate to other files.]]



For more information:

https://www.dynamo-em.org//w/index.php?title=Practical_Suggestions_for_Tomographic_Reconstruction

=== Preparation and management of data ===
The organization of folder can be created by running:


<pre style="color: blue">
hiv1Dirs = {'catalogue', 'particles', 'projects'} ;
cellfun(@(x) mkdir(x), hiv1Dirs);</pre>

the data structure should look as
     .
     │
     └── workplace
        ├── catalogue
        ├── particles 
        └── projects   

Next, move the current directory to catalogue 
<pre style="color: blue">
cd catalogue </pre>

Then,  

<pre style="color: red">Create a '.vll' file inside the catalogue folder containing absolute paths to the tomograms (one path per line). </pre>

A '.doc' file that also contains the path (absolute or relative, it is recommended to use absolute paths) to the tomograms and their volume ID needs should also be stored in the catalogue folder. That is, exactly the same content as your '.vll' with an ID in the beginning of each tomogram path. Remember, in MATLAB, indexing start from 1.

<pre style="color: red">Create a '.doc' file inside the catalogue folder containing and ID followed by the absolute paths to the tomograms (one ID and path per line).</pre>

The  line  
<pre style="color: red">
dcm -create c000 -fromvll tomograms.vll
</pre>

will help to create a catalogue, as a centralised manager of tomograms. 


The GUI will pop up as soon as the catalog 'c000.ctlg' is called via   
<pre style="color: red">dcm c000</pre>

Now,

<pre style="color: red">
Create the binned version of the tomogram selecting a factor of 2. </pre>

[[File:dcm-gui-create-pre-binned-dmmpp.png|centre|thumb|600px| Creation of tomograms binned versions using <code>dcm GUI</code>.]]


Note, these versions will be saved inside the folder of the tomograms, not inside the catalogue folder.

=== Particle picking === 

As commented before, this part involves on-screen interaction. For this, the main tool to work is the ''dcm''. Just make sure to be inside the right directory ('./catalogue’). '''As the lightweight version of this walkthrough does not contain tomograms, most of the tasks in this section ill be skipped. However, similar walkthroughs on modelling, particle picking and catalogue usage can be found in 'Walkthrough on hexagonal lattice coherence through neighborhood analysis' and 'Spherical geometries' '''

Then,

<pre style="color: red">Open one tomogram through the catalogue in ''dtmslice'' (load the previously binned versions, else, there is risk of potential delay and slow performance).</pre>


[[File:dcm-gui-open-pre-binned-tomoslice-dmmpp.png|centre|thumb|600px| Openning pre-binned tomograms via <code>dcm GUI</code> onto tomoslice.]]

The first thing one will probably want to do is to change the contrast, bandpass it and, if not binned, play with a certain thickness. 

<pre style="color: red">Adjust the tomogram appearance.</pre>

 

Next, 

<pre style="color: red">Interactively create a dipole model set.</pre>

  
[[File:tomoslice-create-dipole-set-dmmpp.png|centre|thumb|600px| Creation of dipoleSet model in <code>tomoslice</code>.]]


<pre style="color: red">
To generate the surface parametrisation of the VLPs, every visible target is marked with only 2 clicks: one on its center and one on its surface (press keyboard letter '''c''' and '''n''', respectively). Pressing '''enter''' saves the current dipole in the set and activates the next dipole.  
</pre>

The annotation does not need to be very accurate at this point, since the following alignment projects are designed to cope with inaccuracies introduced at labelling (up to 40 pixels are tolerable). The VLP quality is not relevant, thus try to maximise the number of annotations (approx. 8 to 9 VLPs per tomogram). The reason is that bad/junk particles will be excluded in a dedicated step later on. Although, depending on your specific pipeline and requirements, data curation might be important and, for other projects, you might want to have only good quality targets. Especially if you have a lot of tomograms. 

[[File:3d-silce-along-with-active-pool-model-dmmpp.png|centre|thumb|400px| Annotated volume 3D slice along with active pool model.]]


'''Just make sure to save the dipoles to disk after the labelling of each tomogram is done.'''

[[File:reminder-save-models-to-disk-dmmpp.png|centre|thumb|400px| Saving the dipole set from memory to disk  inside the active mode menu.]]

The catalogue will manage for you folder creation and naming (e.g., './c000/tomograms/volume_1/models/mdipoleSet.omd').

<pre style="color: red">Repeat the annotation process for the remaining tomograms.</pre>

=== Particle extraction  ===

As commented in the opening remarks, the challenge in this section is to set the oversampling configuration. You will need to define a distribution of boxes based on the expected distribution of physical particles. To this end, open one of the binned tomograms:

<pre style="color: red">dtmshow -otf pathToTomogram </pre>

In order to get the support needed w.r.t. the decision making process:

<pre style="color: red">
Use the ruler tool to measure distance.
Decide the optimal relation of:
- physical distance between the particles / sampling distance in pixel of dipole model surfaces
- physical size of the particle / side length that you ask for the particle cropping.
</pre>


[[File:support-measurement-tool-for-oversampling-configuration-dmmpp.png|centre|thumb|400px| Zoom onto a tomogram VLP via <code>dtmshow</code> along with a line measurement of distance to assist in the parameter selection for VLP oversampling. ]]


Once you have decided the oversampling configuration, the time for scripting has come. Get out the catalogue directory to the workspace and:
<pre style="color: blue">mkdir scripts</pre>

Then, simply move or copy the downloaded wild card script and open it via 

<pre style="color: blue">
movefile([pathToDownloads '/Dynamo_models_and_particles/scripts_wc/model_workflow.m'], './scripts')
open model_workflow.m 
</pre>

Remember, you will find a '''[!!]''' and/or an '''NaN''' wherever scripts need to be completed. 

Is time for to fill in the quantities missing to:

<pre style="color: red">
Configure geometry related inputs ( i.e., sampling distance and the box size for particle extraction ). 
Add the separation of crop points onto the vesicle model. 
</pre>

Additionally, the following task will need to be completed

<pre style="color: red">
Create crop points on VLP surface (oversampling).  
Add the radius to vesicle. 
Create crop table from vesicle model. 
</pre>


Finally, each dipole on the annotated tomogram is processed by running the <code>model_workflow.m</code> script. 

Thus,

<pre style="color: red">
run model_workflow.m
</pre>

This will create a regular lattice of coordinates (about 800 to 1,000 per tomogram) on the surface. These coordinates will be the center of the extracted subvolumes. For each tomogram and the associated set of dipoles a regular lattice of cropped subvolumens, the output will be a crop table, and a template of the initial average of those oriented particles (subvolumes).  Technically, the so-called particle extraction process involves the creation of a merged table with <code>dynamo_table_merge</code>. Then, the cropping of sub-volumes with <code>dtcrop</code>. It will not only crop the table but it will also manage the sub-volume tagging. Finally, templates generate averaging all the particles with <code>daverage</code>. While the user can implement its own version, the recommendation is to use the in-house version given that it accounts for missing wedge compensation.



For more details about particle extraction, Dynamo average, and Fourier compensation, the reader is referred to:

https://www.dynamo-em.org/w/index.php?title=Particle_extraction

https://www.dynamo-em.org//w/index.php?title=Daverage

https://www.dynamo-em.org/w/index.php?title=Fourier_compensation_during_averaging


''How to find out if everything went well?'' It is time to do a '''sanity check.'''


Display final table positions of the seed oversampling along with their orientations.
'''NOTE: for example use <tt>pathToTargetParticleFolder = 'PATH_TO_PRECOMPUTED_DATA/particles/pa_TS01_s256';</tt> to inspect the data from the first tomogram.'''
<pre style="color: blue">
finalTable = dread([pathToTargetParticleFolder '/crop.tbl']);
figure; title([num2str(size(finalTable,1)) ' particles.']);
h = dpktbl.plots.sketch(finalTable,'haxis',gca()); 
h.zlength.value=100;
</pre>

[[File:table-from-seed-oversampling-dmmpp.png|centre|thumb|500px|Particles pose from crop table.]]


A simple average of the extracted particles, without alignment, should already reveal the VLP curvature of the surface. To check:
<pre style="color: blue">
finalAverage = dread([pathToTargetParticleFolder '/template.em']);
dshow(finalAverage)
</pre>




[[File:average-bin-2-from-seed-oversampling-dmmpp.png|centre|thumb|500px| Average of cropped particles (projected in x-view). The curvature of the VLP surface is visible. ]]

== Outputs ==

- dipole models.<br>
- initial cropping table.<br>
- particles.<br>

= Subtomogram averaging =
== Initial comments ==
Once multiple copies of individual protein complexes are extracted from the tomograms, the next step is to mutually align and average these sub-volumes. The ultimate goal: to yield a signal-enhanced 3D structure up to sub-nanometer or even near-atomic resolution. This technique is known as subtomogram averaging (StA).

The general StA approach seeks to recover the common signal in the collected particles by an iterative process that alternates alignment and averaging steps. In simple words,  aligning subtomograms is the process of computing for each available particle the set of shifts and angles that will bring them to a common reference. When all the particles have been aligned, they can be averaged together. This process will run in an iterative way in order to gradually improve the attained resolution. In practice, the full dataset is split into two independent datasets that are aligned and averaged separately. The half-sets are commonly referred to as even and odd and are divided based on the particle indexing number.   The 2 half-maps are needed to validate the final map via  gold standard Fourier shell correlation (FSC) computation.


[[File:alignment-project-main-block-iterative-workflow-sta.png|centre|thumb|600px| Mechanics of StA. Each iteration of StA consists of three main phases: preparation, alignment, and averaging. For the first phase, the preparation phase, the particles, reference(s), and mask(s) are loaded to memory. For the alignment phase, the reference and the mask are rotated through a predefined list of search angles. At each of the tested angles, the missing wedge in Fourier space and the rotated mask in real space are applied to the newly rotated reference. The cross-correlation between the reference and every particle is calculated, after which the shifts corresponding to the highest cross-correlation are extracted. For the averaging step, the particles selected for averaging are aligned and averaged using the angles and shifts that gave the best cross-correlation values for those particles. The average is weighted in Fourier space according to the number of times each voxel was measured. At the end of each iteration, the updated references and alignment parameters serve as input for the next iteration. ]]

However, reaching a point where the quality is sufficient to opt for gold-standard refinement is neither simple nor straightforward. With Dynamo as the core tool for high-demand computation, and Chimera for interactive particle selection within the density maps, the protocol to be followed is as outlined here:

'''Creation of initial reference.''' <br> 
- Data preparation and management.<br> 
- Alignment project for references. <br> 

'''Extraction of first per tomogram average.''' <br> 
- Data preparation and management.<br> 
- Particles centering and re-averaging. <br> 
- Alignment project per tomogram. <br> 

'''Computing an average of averages.'''<br> 
- Data preparation and management.<br> 
- Alignment project for average of averages. <br> 

'''Extraction of second per tomogram average.''' <br> 
- Data preparation and management. <br> 
- Subboxing. <br> 
- Alignment project per tomogram.  <br> 
 
'''Extraction of coarse half maps.''' <br> 
- Data preparation and management. <br> 
- Particle filtering, registration, and re-cropping. <br> 
- Gold standard alignment. <br> 

'''Extraction of refined half maps.''' <br> 
- Data preparation and management.<br> 
- Particle filtering and new initial reference computation. <br> 
- Gold standard alignment. <br>

== Input ==

- tomograms.<br>
- initial cropping table.<br>

== Practical notes ==

  '''If the reader is already familiar with <code>dcp</code>, they may choose to skip ahead to the next section.'''

The alignment project is one of the most important objects in Dynamo. It is used to design the operations carried on the particles, to select and configure the environment where the actual computations will occur, to start the computations and to access and visualise the results. For more information about the elements involved:

https://www.dynamo-em.org//w/index.php?title=Alignment_project


The command line–based workflow is convenient since allows for the projects to be streamlined and automatised for increased throughput. To create the project, it is enough to give its name (a string) to <code>dcp.new</code> (for instance, <code>dcp.new(theNameOfMyProject)</code>). 

Regarding file management, it is ideal to work with folders for particle storage. In some cases, it will be sufficient to initialise the new project with this folder, the table, and the mask(s). This is not possible when more than one tomogram is used. Each set of particles will be inside a different folder. In this case, access is centralised by creating a '.star' file.  In MATLAB, the user must store everything in generic containers (cell array). Create a dynamo object  <code>dpkdata.containers.ParticleListFile.mergeDataFolders</code>, and use the <code>writeFile</code> to actually write the file to disk. Dynamo objects are versatile, and, for example, <code>dcp</code>, offers the flexibility to use the same functions with different file formats. 


Once all the files are ready, these can be linked to the existing project through <code>dvput</code> followed by the name of the project, the flag, and the element (for example, <code>dvput(theNameOfMyProject, 'data', theFolderOfMydata)</code> or <code>dvput(theNameOfMyProject, 'data', starFileName)</code>). Just make sure to create default masks in case you have not created ones ''ad hoc'' for your target, <code>dvput(theNameOfMyProject, 'masks', default)</code>. 

The user can choose to work fully interactively. If so, the '''GUI dcp''' is the main tool. The Wizard will open unless the flag '''show''' is passed to 0 at creation time (e.g., <code>dcp.new(theNameOfMyProject, 'show', 0)</code>). It will indicate with a different color what should the user input onto the project. 


[[File:empty-sta-project-dcp-gui-sta.png|centre|thumb|400px| Empty StA project <code>dcp</code> GUI.]]


Analogous to the CLI version, if applicable, do not forget to click the option for default masks use. 


[[File:default-mask-project-dcp-gui-sta.png|centre|thumb|500px| Selection of default mask creation through  <code>dcp</code> GUI.]]



Note, that in contrast to other softwares, the data (particles) and metadata (table(s)) are explicitly different files. Therefore, at this stage of the project, it is highly recommended for the user to get acquainted with table convention that can be found here:

https://www.dynamo-em.org/w/index.php?title=Table_column_convention



Next, numerical parameters will be configured. They define the mathematical part of the project. Among others, in terms of sets of angles to scan, shifts to allow, used bandpass, angular refinement policy or binning options.
It might feel overwhelming at first. Expertise about each parameter of such configuration is required to select the optimal configuration in one-shot. 


A good start for analysis is the direction of the proteins based on the particle model ( defined by the vertical axis first two Euler angles). For instance, if you know the proteins are more or less orthogonal to the membrane, there is no point in scanning 360 degrees. On the other hand, the sampling of the cone must be coupled to the status of the project. That is to say, in an initial state, when there is uncertainty about the orientation, it makes no sense to sample every degree. However, a good estimate could be every 10 degrees. This reasoning also applies if you have reached your subvolumes by template matching, since at most, you have a directional error of 30 degrees. 


As rule of thumb, the third Euler angle (defining the inplane rotation) can be assumed to be biased. Either for the missing wedge, or for some arbitrary sampling of the particles. Therefore, azimuth angle randomization shall be performed. The function needed is  <code>dynamo_table_randomize_azimuth</code>. Do not forget to create a template with the new table, <code>daverage</code>, and save to disk, <code>dwrite </code>. Again, depending on your target, a global scanning might be necessary. In contrast, if you have some premise that 90 degrees is sufficient, as is the case for HVI-1, go ahead and set it up. 

For those users who are still not confident with the orientation parameters (depicted in blue in the figure of numerical parameters), the recommendation is to complete the advance starters guide. In this walkthrough, the orientation is unravelled into separate projects. The links to the tutorial and data are as follows:

https://www.dynamo-em.org/w/index.php?title=Advanced_starters_guide

https://wiki.dynamo.biozentrum.unibas.ch/w/doc/data/fhv/crop.rec




To configure the position parameters, think about how likely are the particles to be shifted in x,y, and z. In most cases, after a couple of trials cropping, more precision can be expected in x,y. Whereas the z, is known to be harder to fit. Thus, it seems reasonable to permit more translational freedom in that direction.  The modus of this area of search, unless having a good reason, should be set to 1. Meaning that the limits are understood as the center of the particle. Note that mode 0 is dangerous, and accumulated shifts can lead to very large displacements. These two values are squared in gray in the figure.


Regarding filtering, try to be realistic. The default value for high pass is okay. It will prevent noise from ruining the alignment. However, the low pass filter needs to be a meaningful value. In any case, even if you are not sure, do not push Nyquist. Instead, support this decision by visualising the binned template, for instance, you might use <code>dshow</code> and calculate the pixel size with  <code> dynamo_resolution </code>. Although, if you have half maps available, it is worth it to compute the FSC curve. 

Symmetrization is a powerful technique. However, it should initially be avoided since it can induce bias. As the project evolves, if the target happens to be clearly symmetric, exploit this fact. Just be aware that symmetries that are not a subset of the sampling scheme are not recommended. Namely, any other family but circular families ('''cn''').  

The quality, signal-to-noise ratio, and homogeneity of the target will determine the number of iterations needed. It is a trade-off between computational time and the improvement being obtained. To do this, you can open another instance of MATLAB and buy the results. 

The logical trend is that with higher resolution, the parameters can be more restrictive, accompanied by a larger particle size, and the low pass filter cutoff frequency reduced. Generally, the higher the level of detail, the more iterations will be needed for convergence. 


[[File:numerical-parameters-gui-sta.png|centre|thumb|600px| Numerical parameters configuration overview via <code>dcp</code> GUI. ]]

The computational environment will be determined by the resources you have available. Try to maximise the number of cores and GPUs as this is directly proportional to the computational time. If it is a long project, you probably do not want to use ''gpu 0''. So that the pc used to visualise the results is not affected by the computations. In case you are sharing resources, make sure that there is no mismatch between the system ids and MATLAB (remember that the indexing for your system starts at 0 and for MATLAB at 1). 

[[File:environment-parameters-gui-sta.png|centre|thumb|400px| Environment  parameters GUI configuration overview via <code>dcp</code> GUI. ]]


The project is now ready to be checked, unfolded, and launched. Do not worry if some path is broken or not found, <code>dcp</code>  will warn and prevent the project's launching unless it can be done safely. For this, you might want to run: 

	   
        <code>outputCheck = dynamo_vpr_check(theNameOfMyProject);</code>

        <code>if isempty(outputCheck)</code>
          <code>dvrun(theNameOfMyProject,'unfold',true);</code>
        <code>else</code>
          <code>disp('Project not safe.')</code>
        <code>end</code>


=== Creation of initial reference ===

Place yourself in the './scripts' directory. Then, simply move or copy the downloaded already prepared card script and open it via

<pre style="color: blue"> movefile([pathToDownloads '/Subtomogram_averaging/scripts_wc/*'], './scripts')
 open aligment_for_reference.m </pre>

Note that all the scripts needed for the StA section have been moved already to the scripts folder.

Within the script, ensure you have updated the environment parameters and paths as necessary. Remember, a '''[!!]''' and/or '''NaN''' will precede the lines that require your input. Consider establishing a global variable for consistent use throughout the subsequent scripts.

<pre style="color: blue"> Change paths.</pre>

The last table azimut needs to be randomised to avoid missing wedge bias: 
<pre style="color: red"> Create first template with randomized azimuth. </pre>


Everything is ready to create a new project.  Below, the numerical parameters.

[[File:initial-reference-project-configuration-sta.png|centre|thumb|600px| Numerical parameters configuration for initial reference alignment project.]]


<pre style="color: blue"> Fill in the input values, computational, and environment parameters for pr_0. </pre>

Launch it,

<pre style="color: red"> run aligment_for_reference.m </pre>


The last iteration average should look something qualitatively close to the one shown in the figure. To extract and plot the last average,

<pre style="color: red"> aPath  =  ddb([ pr_0 ':a:ite=last']);
 a      = dread(aPath);
 dshow(a)
</pre>

'''To do the same but with the pre-computed results, do instead:'''
<pre style="color: blue"> pr_0 = [root,'data/STA/projects/pr_TS01_0'];
 aPath  =  [pr_0,'/results/ite_0006/averages/average_ref_001_ite_0006.em'];
 a      = dread(aPath);
 dshow(a)
</pre>


At this stage, with only one tomogram, the wild average (shifting mode 2 in the second round allows for this connotation) already reveals a strong ''c6'' symmetry. 

[[File:last-average-from-initial-reference-project-sta.png|centre|thumb|600px| Last average from initial reference alignment project (visualisation with <code>dshow</code>).]]

=== Extraction of first per tomogram average ===

Let us start by preparing the last average by running the <code>prepare_average_for_chimera.m</code> script.

<pre style="color: red"> movefile([pathToDownloads '/Subtomogram_averaging/scripts/*'], './scripts')
 prepare_average_for_chimera(pr_0)</pre> 

'''To do the same but with the pre-computed results, do instead:'''

<pre style="color: blue"> movefile([pathToDownloads '/Subtomogram_averaging/scripts/*'], './scripts')
 theNameOfMyMap = 'result_pr_TS01_0_INVERTED.em';
 dwrite(dynamo_bandpass(a,[1 23])*(-1),theNameOfMyMap);</pre> 


This function will simply filter via <code>dynamo_bandpass</code>, and invert the map.

To open the file into Chimera,

<pre style="color: blue"> dynamo_chimera(theNameOfMyMap)</pre> 

If for some reason is not able to find the application, just give Dynamo the path,
<pre style="color: blue"> dynamo_chimera(theNameOfMyMap,'path',thePathToChimera) </pre>

The volumen viewer window will pop out. 

<pre style="color: blue"> Adjust the density threshold and, if helps, change the color. </pre>


[[File:step-one-preprocesing-chimera-sta.png|centre|thumb|350px| Adjusting the density threshold of Chimera volume viewer. ]]


The second pre-processing step is to remove noise. 

<pre style="color: blue">Go to Tools > Volume Data > Hide Dust.
Adjust the percentage until one unique blob is in the scene. </pre>


[[File:noise-configuration-chimera-sta.png|centre|thumb|375px| Adjusting the noise removal percentage of Chimera hide dust.  ]]

The density map is ready to be annotated: 

<pre style="color: blue">Open the annotation window inside Tools > Volume Data > Volume Tracer. 
Create a new marker set via File > New marker set. 
Select a proper configuration for marker placement.  </pre>

[[File:chimera-volume-tracer-sta.png|centre|thumb|275px| Configuration for marker annotation of Chimera volume tracer. ]]

Next, rotate the average upside down: 

<pre style="color: blue">Select a proper perspective with the help of Tool > Viewing Control > Side View (check the figure below and pay attention to the slicing tool position). </pre> 

Marking the tip of the CA-N terminal domain,

<pre style="color: blue">Define the center of unit cell with a marker.</pre> 

Save them to disk:

<pre style="color: blue">Volume tracer: File > Save current marker set >  'reference_center.cmm' </pre>


[[File:chimera-reference-center-gui-sta.png|centre|thumb|600px|Center of unit cell defined in Chimera UCSF. The center of a unit cell (yellow dot) depicts the tip of the CA-N terminal domain. Marking via volume tracer tool and viewing controls adjusted via side view tool. ]]


From here, copy and open the <code>first_aligment_per_tomogram.m</code> script.

<pre style="color: blue"> movefile([pathToDownloads '/Subtomogram_averaging/scripts_wc/*'], './scripts')
 open first_aligment_per_tomogram.m</pre> 


As usual, change the paths and set up the available hardware resources,

<pre style="color: blue">
Change paths.
Specify hardware resources. </pre> 

Using the 'reference_center.cmm', the task you will need to complete is the centering and re-average of the particles.  
<pre style="color: red">Compute the vector pointing from box center to the new coordinate.
Create transformation object, and transform the table per se. 
Re-average the particles (Do not forget Fourier compensation).
</pre>

Of course, before launching an alignment project, you should perform a sanity check. To do this, visualise the volume before and after the centre and click on the corresponding coordinate.

<pre style="color: red">dshow -otf aBeforeCentering 
dshow -otf aAfterCentering 

Pin the coordinate and compare. </pre> 

'''To do the same but with the pre-computed results, do instead:'''

<pre style="color: blue">dshow(aPath) ;
recenterPath = [root,'data/STA/projects/result_pr_ts001_0_recentered.em'];
dshow(recenterPath); 

Pin the coordinate and compare. </pre> 


[[File:sanity_check_centering-sta.png|centre|thumb|600px| Results of pr_0 before, and after entering and re-avering. Left and right, respectively.  Coordinate of the reference center unit is displayed for comparison purposes.  Coordinate picking tool on the bottom. ]]

From here, running the initial per tomogram alignment  is straightforward. The configuration of the project is identical to the alignment for reference but from here on, imposing a ''c6'' symmetry. Thus,

<pre style="color: red">Configure the c6 symmetry on both rounds of the project.</pre>


[[File:first-per-tomogram-and-averages-of-averages-project-configuration-sta.png|centre|thumb|600px| Numerical parameters configuration for first per tomogram alignment project.]]

Each project will take about 1 hour 30 min (estimation with 8 gpus). To directly run the script:

<pre style="color: red">run first_aligment_per_tomogram.m</pre>

The results of the last iteration symmetrised average should look something qualitatively close to the ones shown in the figure. Visualise the results and compare against the one you have obtained. 

<pre style="color: red">aTsAll = [];
for idx = 1:nTomo
    aPath = ddb([pr_1{idx} ':as:ite=last']);
    aTsAll = [aTsAll; dread(aPath)];
end
dshow(aTsAll)
</pre>

'''To do the same but with the pre-computed results, do instead:'''
<pre style="color: blue">aTsAll = [];
nTomo = 5;
path_to_av = [root,'data/STA/projects/pr_ts_all_1/results/ite_0006/averages/'];
for idx = 1:nTomo
    aPath = [path_to_av,'average_ref_00',num2str(idx),'_ite_0006.em'];
    aTsAll = [aTsAll; dread(aPath)];
end
dshow(aTsAll)
</pre>

The first thing to notice is that each average is centred. Thanks to this operation, the periphery reveals more proteic units. The biologically most relevant feature is the hexameric intra-structure of Gag, which can be observed in each of the units for all tomograms. The differences found in each tomogram can be attributed to the inherent quality of the tomogram, the structural variety of the biological target (size, imperfections, degree of immaturity of the VLP), and the randomness of the inherent process that takes place during the alignment project.

[[File:last-average-first-per-tomogram-average-sta.png|centre|thumb|600px| Last average from initial reference alignment project. ]]

=== Computing an average of averages ===
To complete this small project, open the corresponding <code>aligment_for_average.m</code> and change the paths inside,

<pre style="color: blue">open aligment_for_average.m 
Change paths.</pre>


The configuration of the project is maintained from  first per tomogram alignment. The difference lies in the data used. Here, as the names implies, all the tomograms averages will be merged onto one single average:

<pre style="color: red">Create corresponding merged table, and set:
- No missing wedge compensation needed.
- Particle tag number = tomogram number
- Centers.</pre>

The motivation behind this computation is to minimise manual labour in the next step. It will run quickly (near 15 minutes 8 workers), since it only contains 5 particles. Launch the project by typing:

<pre style="color: red">run aligment_for_average.m</pre> 

To check the final average,

<pre style="color: red">aPath = ddb([pr_a ':as:ite=last']);
dtmshow(aPath)</pre>

'''To do the same but with the pre-computed results, do instead:'''

<pre style="color: blue">pr_a = [root,'data/STA/projects/pr_a'];
aPath  =  [pr_a,'/results/ite_0006/averages/average_symmetrized_ref_001_ite_0006.em'];
dtmshow(aPath)</pre>

The symmetrised average should look something qualitatively close to the one shown in the figure. 
[[File:Final-average-from-average-of-averages-sta.png|centre|thumb|600px|  Last average from ''average of averages'' project. ]]

=== Extraction of second per tomogram average ===

The '''average of averages''' is prepared as previously:
<pre style="color: red">prepare_average_for_chimera(pr_a)</pre>

'''To do the same but with the pre-computed results, do instead:'''

<pre style="color: blue">
 theNameOfMyMap = 'result_pr_a_INVERTED.em';
 a = dread(aPath);
 dwrite(dynamo_bandpass(a,[1 23])*(-1),theNameOfMyMap);</pre> 


To open the prepared map inside Chimera,

<pre style="color: blue">dynamo_chimera(theNameOfMyMap)</pre>

This time, every 18-meric assembly of the capsid protein p24 are marked, so:

<pre style="color: blue">Volume tracer: File > New marker set
Mark the centers of every unit cell of the hexameric lattice 
Save the  current marker set under 'particle_centers.cmm'.</pre> 

[[File:chimera-particles-center-gui.png|centre|thumb|600px|All unit cells defined in Chimera UCSF.  The ''average of averages'' procedure, allows for candidate coordinates of all unit cells can be determined for the full dataset by manually labelling only 1 single density map.  ]]

The rest of this is implemented inside the script second_aligment_per_tomogram.m There are a lot of things going on for data preparation before running the alignment. Thus, it is worth diving into what is going on in detail.

We can use our previous results to set up an approach that centers the alignment on the individual. Each one of the particles in our new data set will be a subbox (a lattice unit cell) that we have just defined in the averages. The so-called '''subboxing technique'' will be used to redefine the area of interest inside a previously defined average. 

In this example, the candidate coordinates will be determined using a two-step subboxing procedure. Note that, in principle, one subboxing step would suffice. However, this would have to be done for each tomogram individually. Therefore, the motivation behind the average of average computations is to minimise manual labour. Since the average is composed of multiple subvolumes with individual orientations, these coordinates can be then translated (or mapped) onto each of those subvolumes and finally also onto the tomograms themselves. 


Let us now begin with the hands-on portion of the script:

<pre style="color: red">open second_aligment_per_tomogram</pre>


The first step consists in mapping coordinates back onto the averages. Similarly to what was done for re-entering, 
the vector from box center to new coordinate for each coordinate is computed. To complete this process,

<pre style="color: red">Create a mini table for each average that contains the corresponding coordinates.
Merge the mini tables.</pre> 

The second step involves mapping coordinates back onto the tomograms. Now, the vector from each average box center to the new subbox center is computed. Finally, as above, 
<pre style="color: red">Obtain the last subboxing table per tomogram.
Combine the tables into one.</pre> 

The coordinates describing the same volume should be removed. Therefore, 
<pre style="color: red">Exclude the repeated unit cell. </pre> 


Before cropping the particles, let us visualize the merged table generated containing the candidate positions before and after subboxing:

<pre style="color: red">Add a break point in the line where ckpt variable is declared.
Run second_aligment_per_tomogram.m script section</pre>


About 5000-9000 coordinates per tomogram should have been defined. One table should suffice as a means of comparison:

<pre style="color: red">docFilePath = [root 'data/dynamo_models/tomograms.doc'];
fileID  = fopen(docFilePath); D = textscan(fileID,'%d %s'); fclose(fileID); tomoID  = D{1,1}';
idx = 1;

tSubboxPath = ['subbox_' D{1,2}{idx,1}(end-7:end-4) '.tbl'];
tOriPath = ddb([pr_1{idx} ':rt']);

refs = {'Before', 'After'};
tbls = {tOriPath{1},tSubboxPath};

figure;
for i = 1:length(refs)
subplot(1,2,i);title([refs{i} ' subboxing: ' num2str(size(dread(tbls{i}),1)) ' particles.']); h = dpktbl.plots.sketch(tbls{i},'haxis',gca()); h.zlength.value=100;
end</pre>


'''To do the same but with the pre-computed results, do instead:'''

<pre style="color: blue">docFilePath = [root 'data/dynamo_models/tomograms.doc'];
fileID  = fopen(docFilePath); D = textscan(fileID,'%d %s'); fclose(fileID); tomoID  = D{1,1}';
idx = 1;

tSubboxPath = [root,'data/STA/projects/subbox_' D{1,2}{idx,1}(end-7:end-4) '.tbl'];

path_to_tables = [root,'data/STA/projects/pr_ts_all_1/results/ite_0006/averages/'];
tOriPath = [path_to_av,'refined_table_ref_00',num2str(idx),'_ite_0006.tbl'];

refs = {'Before', 'After'};
tbls = {tOriPath,tSubboxPath};

figure;
for i = 1:length(refs)
subplot(1,2,i);title([refs{i} ' subboxing: ' num2str(size(dread(tbls{i}),1)) ' particles.']); h = dpktbl.plots.sketch(tbls{i},'haxis',gca()); h.zlength.value=100;
end</pre>

[[File:subboxing-tbl-vs-last-ite-per-tomogram-project-tbl-sta.png|centre|thumb|600px| Resulting coordinates from particle picking via <code>dtplot</code>. Last table from first per tomogram alignment project, left; Table after subboxing, right. The original shapes of each VLP surface are preserved.]]

We can complete the rest of the script. So, next, we will use those coordinates to extract a new set of subvolumes with a smaller sidelength. 

<pre style="color: red">Re-crop subboxed particles onto corresponding folder with a box size of 192. </pre>

Prior to running the second per tomogram project, the subboxed particles are re-averaged, and save them to disk.  Since we will need input templates for the incoming project. Its configuration is found below. If we compare this second configuration with the previous, the angular search space and shift limit is reduced. This is possible because particle positions are now more accurate. Recall that the hexagonal lattice structure inherent to the immature HIV-1 CA-SP1 is the target of the present localised alignment. Nevertheless, high resolution is not of interest yet. This is the reason why subvolumes are still binned on-the-flight.

[[File:second-per-tomogram-project-configuration-sta.png|centre|thumb|600px| Numerical parameters configuration for second per tomogram alignment project.]]


This section, per tomogram, is estimated to take 30 minutes for the subboxing block (8 workers), and 3.5 hours for the alignment project (8 gpus). To actually compute the cropping of the particles and the corresponding alignment projects:
<pre style="color: red">run second_aligment_per_tomogram.m</pre>


The results after running the project demonstrate how localised alignments are capable of driving the alignment robustly. To replicate the figure:

<pre style="color: red">aTsAll = [];
for idx = 1:nTomo
    aPath = ddb([pr_2{idx} ':a']);
    aTsAll = [aTsAll; dread(aPath)];
end
dshow(aTsAll)
</pre>

'''To do the same but with the pre-computed results, do instead:'''

<pre style="color: blue">aTsAll = [];
nTomo = 5;
path_to_av = [root,'data/STA/projects/pr_ts_all_2/results/ite_0003/averages/'];
for idx = 1:nTomo
    aPath = [path_to_av,'average_symmetrized_ref_00',num2str(idx),'_ite_0003.em'];
    aTsAll = [aTsAll; dread(aPath)];
end
dshow(aTsAll)
</pre>



The higher contrast is associated with a stronger signal. This is important since a potential downside of this technique is the reduction of the amount of signal available to drive the alignment due to smaller box sizes. Therefore, is an indicative of a better estimation of their pose; a better concordance between particles (e.g. in the cavities between MA - CA , there is no density. In previous stages of the project, in these areas, there was uniformly distributed noise). Furthermore, in the  orthoslice through the CA-SP1 region (last row), each unit of lattice locally morphology was previously fuzzy.  Whereas now, the density is ordered and some helical patterns can be grasped. 


[[File:last-average-second-per-tomogram-average-sta.png|centre|thumb|600px|  Last average from second per tomogram alignment project. ]]

=== Extraction of coarse half maps ===
This block is 100% automated. As we have been doing so far, some functionalities have to be implemented in order to be able to launch it. 

Open the corresponding script:

<pre style="color: blue">open coarse_gold_standard_alignment.m</pre> 

From the last project, the tables will be gathered. Column 10 stores the CC between each aligned subvolume and the final average. The exclusion process is done by fitting a Gaussian mixture model (GMM) to the such  distribution. The MATLAB function fitgmdist will help to classify them into a combination of two populations. The one with the lower CC score is assumed to be the “bad” class (i.e., particularly bad quality or subvolumes containing pure noise). Therefore, it  needs to be excluded from the processing. The threshold is then defined by taking the minimum between the 2 Gaussian peaks and by adding an empirically defined constant of 0.01. 


While the table get filtered, to visualise the results add the following lines in the loop:

<pre style="color: red">figure;
hold on;
histogram(t(:,10), 'Normalization', 'pdf', 'EdgeColor', 'none')
plot(x, pdf(gmdist, x'))
xline(thresh)
xlim([0 0.5])
legend('CC distribution','GMM fit','Threshold')
stackName{idx} = D{1,2}{idx,1}(end-7:end-4);
title(['Stack: ' stackName{idx}])</pre>


'''Tips:'''

<pre style="color: red”>You can activate the interactive mode by adding a `keyboard` command in your script, then type in your visualization on the terminal.
To terminate debug mode and continue execution, use the dbcont command.
To terminate debug mode and exit the file without completing execution, use the dbquit command.</pre>



Similarly to the subboxing section, is interesting to see how the filtering strategy worked. So:

<pre style="color: red">Add a break point in the line where ckpt variable is declared.
Run coarse_gold_standard_alignment.m script section</pre>


For each stack of particles, the threshold determination and the fitted distributions in the graphs below.
 

[[File:gmm-population-clasification-sta.png|centre|thumb|950px| GMM fit to the CC score distribution of the particles from each tomogram and the automatically defined threshold for particle exclusion.]]




Following filtering, a new average per tomogram is computed. In the next step, the height of every particle is homogenised to a common center unit cell. To ensure that, independently of the tomogram, the z-height should be consistent across sub-volumes. To this end a synthetic reference is derived from a model:

<pre style="color: red">Compute a mask from a membrane object. </pre> 

The object's properties, thickness, radius, and side-length are already changed. In practice, you will need to visualise what you are doing. So <code>dview</code> is what you are looking for:

<pre style="color: red">dview theNameOfMyMask </pre> 


[[File:Synthetic-height-reference-sta.png|centre|thumb|550px|Synthetic reference for alignment to common height. Mask computed from membrane object X-View, left; Z-view, right. ]]


Both, template and reference, are band-passed and all the particles are aligned against it. Remember, axis orientation and x/y-shifts are already consistent due to the imposed C6 symmetry. Thus, the only permitted shift is along z-direction:

<pre style="color: red">Fill in the dalign allowing only 20 pixels of z-shift (flag: "lim").</pre> 


Then, the resulting transformation parameters are applied to the table via <code>dynamo_table_rigid</code>. Again, fixing coordinates that describe the same particle and saving the table to disk.

Finally, all subvolumes containing one centered unit cell are re-extracted one last time. The side length is of 192 pixels is preserved. This time, per tomogram, about 2,000 to 5,000 particles are expected. 


The point has been reached: the quality of the particles and their initial orientation are good enough to serve as the basis for a gold standard alignment. Prior to launching the project, as explained before:

<pre style="color: red">Merge everything to a single container
Get even and odd index
Prepare each half template.</pre> 


The configuration of the 2 independent alignment projects using the starting references that have been created before is displayed below. A total of 3 rounds with 3 iterations each are set up. The parameter search space is reduced after each round. The low-pass filter is increased to 32 pixels (8.1 Å), and the last 2 rounds are run on full-sized particles.

[[File:coarse-half-maps-alignment-configuration-sta.png|centre|thumb|600px|  Numerical parameters configuration for coarse half maps alignment project. ]]

The estimated times are, for registration, and re-cropping, 20  min (8 workers). The coarse even/odd project is around 48 hour per half-set (8 gpus). Therefore, if you have already completed the task filling in the missing parts:

<pre style="color: red">run coarse_gold_standard_alignment.m</pre> 


To check the resulting maps:
<pre style="color: red">aEPath = ddb([pr_EO{1} ':as']);
aOPath = ddb([pr_EO{2} ':as']);
dshow([dread(aEPath); dread(aOPath)])</pre>

'''To do the same but with the pre-computed results, do instead:'''

<pre style="color: blue">pr_EO = [root,'data/STA/projects/pr_EO'];
aEPath = [pr_EO,'/results/ite_0009/averages/average_symmetrized_ref_002_ite_0009.em'];;
aOPath = [pr_EO,'/results/ite_0009/averages/average_symmetrized_ref_001_ite_0009.em'];;
dshow([dread(aEPath); dread(aOPath)])</pre>


[[File:last-average-for-coarse-half-maps-sta.png|centre|thumb|900px| Last average from coarse gold standard alignment project. ]]

=== Extraction of refined half maps ===
The script  <code>refined_gold_standard_alignment.m</code> is responsible for refining previous results by pruning again the dataset and running a final alignment project with adapted parameters.

Open the corresponding script:

<pre style="color: blue">open refined_gold_standard_alignment.m </pre> 

This time, to remove particles of low quality, the dataset is reduced through a normalised, per tomogram, CC filtering. The reason for thresholding  separately is to avoid reducing the coverage of the zeros of the CTF. The risk emerges from different defocus values of the tomograms influencing the overall CC score. With the help of <code>dpksta.filters.byCC</code>, subvolumes with a score lower than 1 standard deviation below the mean are removed. This function performs normalisation by the angle of latitude. This is important, since particles at the poles of the VLPs generally have a lower CC score than particles on the equator due to missing wedge effects. More precisely, this function works by fitting the first 2 terms of a general Fourier series to the data. As usual, after the tables are pruned, a new average is computed that will serve as reference for the last alignment project.

For this 2 projects configuration the set up is the same parameters as the last round of the previous project but with a stricter low-pass filter up to 38 fourier pixels (~ 6.8 Å). Although, the intuition is that there is still room for improvement in terms of resolution. Therefore, the last bullet is to focus the refinement more onto the area of interest. 

[[File:refined-half-maps-alignment-configuration-sta.png|centre|thumb|600px|   Numerical parameters configuration for refined half maps alignment project. ]]


For this, the idea is to create a tighter mask. While one option is to draw the mask by hand, the choice is to exploit the already existing motives in dynamo (as last time, there is trial and error until the properties satisfy you biological target). For target the capsid:

<pre style="color: red">Create a custom tight mask as the intersection of a membrane (parameters given), and a cylinder ( radius : 87; size : 192; position : [97 97 97]).</pre> 



Be aware that the default mask option is still used. The new mask is overwritten via <code>dvput(theNameOfMyProject,'file_mask',theNameOfMyMask)</code>. 


[[File:adapted-mask-sta.png|centre|thumb|600px|   Custom mask for refinement gold standard project. ]]

The computation time of this for the refinement project is estimated to be 24 hours per half-set (8 gpus). Launch it via:



<pre style="color: red">run refined_gold_standard_alignment.m</pre> 


Et voilà! The latest results have been obtained. And with them, we have finished the subtomogram averaging projects.

[[File:final-sta-result-from-refined-gold-standard-sta.png|centre|thumb|600px| HIV-1 immature Gag structure, shown as cryoET final subtomogram averaging map though two orthogonal views via <code>dshow</code>, top left.]]

Finally, the resulting odd/even maps are reformatted and prepared. To this end, the even half-map is aligned to the odd half-map with the usual function <code>dalign</code>. The resulting transformation parameters are used to transform the table corresponding to the even half-map though <code>dynamo_table_rigid</code>. The table is then used to re-average the even particles using<code>daverage</code> imposing a '''c6''' symmetry, <code>dynamo_csym</code>.

The reason for this re-averaging is to avoid edge artifacts caused by the half-map alignment. Then, they are saved (<code>dwrite</code>) according to RELION convention so that they can be used for the gold standard Fourier shell correlation (FSC) computation. This is important, since unless the map's name starts with ‘half1_*.mrc’ and  ‘half2_*.mrc’. Else, RELION will not be able to parse them. 


A final density map is additionally computed by re-averaging all particles of the full dataset. This isosurface visualisation is displayed below. Note the necessity of re-averaging  to ensure a correct Fourier compensation of the final map. Therefore, the corresponding flag, '''fc''' is set to true (<code>daverage(nameOfTheStarFile,'t',dynamo_table_merge({theNameOfHalfMapEven, theNameOfHalfMapOdd}),'fc’,1)</code>). Then, simply:

<pre style="color: red">dynamo_chimera(theNameOfMyFinalMap)</pre> 

'''To do the same but with the pre-computed results, do instead:'''

 <pre style="color: blue"> theNameOfMyFinalMap = [root,'data/STA/projects/final_reaveraged_ctf_corrected_thick_tomos.mrc'];
 dynamo_chimera(theNameOfMyFinalMap)</pre> 

[[File:isosurface-final-structure-density-map-chimera-sta.png|centre|thumb|600px| Isosurface visualisation via Chimera of the final structure showing the immature CA-SP1 lattice viewed from outside the virus, left; from inside, right. ]]

== Outputs ==

- Final high-resolution dynamo average.<br>
- Two high-resolution half-maps.<br>
- Particle table of high-resolution dynamo average.<br>

= Fourier shell correlation =
== Initial comments ==

At the end of a gold-standard alignment, we can estimate the resolution of the final average by comparing the Fourier shell correlation (FSC) of the two half-maps.<br>

An FSC mask should be generated to limit the FSC computation to the areas where the protein densities are present. A better FSC mask can be used to deduce better resolutions, but a poor FSC mask can either deduce a worse resolution or even an artefactual resolution.<br>

So this analysis can be divided in two steps: the FSC mask creation and the actual FSC computation. Both of these steps can be done either in dynamo or in RELION. In this walktrough the alignment mask will be generated in dynamo and the FSC computation will be done in RELION.<br>

This choice was made because the mask creation in RELION (using the Mask creation job) can limit the user freedom. In fact, it will always be a Gaussian filtering and binarization of an input density map.<br>

On the other hand, the FSC computation in dynamo using the <code>dynamo_fsc</code> function will only compute the FSC curve (https://www.dynamo-em.org/w/index.php?title=Fourier_Shell_Correlation). RELION Post-processing job will calculate the phase randomized FSC in addition to the usual FSC.<br><br>


Here is an example of an FSC computation using RELION on a different dataset. Note the effect of the mask on the FSC curve.

[[File:Example_of_healty_FSC.png|centre|thumb|800px| FSC curves computation using RELION on another dataset for reference.]]

In general, an healthy FSC computation will have the following characteristics:<br>
- RELION assigns the final FSC close to the FSC of the masked half-maps (blue line, see legend).<br>
- The phase randomized FSC (red line) will go to zero quickly.<br>
- At high frequency all FSC curves will go to zero.<br><br><br>

We must take care to note the behavior of all of the FSC curves. Some typical causes of unhealty FSC curves can be: unmasked symmetrization artefacts, sharp edges caused by the FSC mask or duplicate particles in the average.

== Input ==

- two half-maps generated using the gold standard protocol.<br>

== Practical notes ==

The FSC computation will be done using RELION. Additionally, the script <code>FSC_mask_generator.m</code> will be used to generate the FSC mask.<br><br>

<pre style="color: blue">Remember to have the root variable as the path to the data folder</pre>

'''Change your MATLAB working folder to the FSC_computation/scripts subfolder.'''

The even and odd averages found through dynamo have been precomputed. You can load them:<br>

<pre style="color: blue">h1 = dread([root,'data/Dynamo_final_results/half1_final_reaveraged_unfil_ctf_corrected_thick_tomos.mrc']);

h2 = dread([root,'data/Dynamo_final_results/half2_final_reaveraged_unfil_ctf_corrected_thick_tomos.mrc']);

</pre>

You can inspect the differences between the two half-maps by:

<pre style="color: blue">dshow(abs(h1-h2));</pre>

You should see something like:
[[File:H1_d_h2.png|centre|thumb|600px| Difference between halfmaps]]

Now we can start to build the mask.

<pre style="color: blue">Inspect a half map and determine the dimensions and center of cylinder that would contain the middle "wineglass" structure.</pre>

Generate a cylindrical mask:

<pre style="color: blue">mask = dcylinder([35,35],192,[96,96,126]);</pre>

And compare it with the half-maps:

<pre style="color: blue">dshow(h1+mask)</pre>

You should see something like:
[[File:HIV_mask_s.png|centre|thumb|600px| Mask on template]]

As mentioned in the initial comments, the mask will be mollified to avoid sharp edges. Remove them and inspect the result:
<pre style="color: blue">mask_m = dynamo_mollify(mask,3);
dshow(h1+mask_m)</pre>

You should see something like:
[[File:HIV_mask_m.png|centre|thumb|600px| Smooth mask on template]]

The FSC mask will be created using the usual <code>dwrite</code> command followed by the <code>dpkio.tomo.mrc.setPixelSize</code> command to set the appropriate pixel size in Angstrom in the header of the mask .mrc file. Copy the half-maps too:

<pre style="color: blue">px_size = 1.35;

dwrite(mask_m,'FSC_mask.mrc');
dpkio.tomo.mrc.setPixelSize(px_size,'FSC_mask.mrc');

dwrite(h1,'half1.mrc');
dwrite(h2,'half2.mrc');

dpkio.tomo.mrc.setPixelSize(px_size,'half1.mrc');
dpkio.tomo.mrc.setPixelSize(px_size,'half2.mrc');</pre>

Now RELION will be used to estimate the resolution. The RELION job to do this analysis is the Post-processing job.

Some relevant RELION documentation can be found in:<br>
https://relion.readthedocs.io/en/release-4.0/STA_tutorial/Postprocess.html<br>
https://relion.readthedocs.io/en/release-4.0/SPA_tutorial/Mask.html<br>

In this case, we will use the the generated FSC mask and ask RELION to estimate the sharpening B-factor automatically.

Note that the two half-maps need to be in the same folder and follow the naming convention: 'half1_XXXX.mrc' and 'half2_XXXX.mrc'. It may be necessary to rename the two half-maps to be able to set them in RELION.<br><br>

<pre style="color: red">load RELION4 with relion --tomo</pre>

<pre style="color: red">Click to the Post-processing window</pre>

You will see the following GUI:
[[File:Post_process_IO2.png|centre|thumb|600px| RELION Post-process I/O field. Note the name of the half-map file. The second half map will be in the same folder and be named 'half2.mrc'.]]
[[File:Post_process_Sharpen2.png|centre|thumb|600px| RELION Post-process I/O sharpen with automatic B-factor estimation.]]
[[File:Post_process_Running2.png|centre|thumb|600px| RELION Post-process I/O running with forced mask.]]

<pre style="color: red">Fill in the GUI with half1.mrc as the map and FSC_mask.mrc as the FSC mask.

Set the pixel size at 1.35 A, force RELION to use your FSC mask and the press Run!

Inspect the result by using the RELION GUI Display option of the logfile.pdf
</pre>

RELION will compute the FSC curves of the two half map with and without the FSC mask and then will select the best curve to correct and select for the final resolution estimation. It is sometimes the case, especially for membrane protein, that the unmasked FSC will be completely artefactual and this can cause problem for the RELION final resolution estimation. In these cases, RELION will suggest to force the usage of the FSC mask by the '--force_mask' option.

The sharpened final average is saved as the postprocess.mrc file in the Postprocess job folder that RELION just created.

== Outputs ==

- FSC curve.<br>
- The corresponding FSC mask.<br>

= RELION conversion =
== Initial comments ==

RELION has some very interesting capabilities that can be used to complement dynamo. In particular, RELION 3d auto-refine job type is useful because it allows to run the equivalent of dynamo STA on locally back-projected particles. One advantage of this approach compared to cropping particles from the tomograms is that RELION can exploit explicitly the information that each particle is, in fact, a reconstruction based on information of multiple patches. However, a more immediately practical advantage is that RELION does not need to store in memory all of the full-scale tomograms.<br><br>
Some disavantages of RELION are that it needs relatively good initial particle positions to start a 3d auto-refine job correctly and that RELION allows for less control on the angular alignment search spaces that dynamo. In addition to this, RELION is not meant to receive reconstructed tomograms as inputs, but instead it requires the unaligned tilt series to allow for the local back-projection procedure. This means that RELION also needs to receive the alignment parameter of the tilt series as well as CTF and dose information, which can be confusing to implement in practice.<br><br>
However, these disadvantages can be removed by using the dynamo-RELION conversion tools. Two main tools are used in conjunction:<br>
- a tomogram.star generator: this script will use the results of the dynamo fiducial-based tilt series alignment to create a file for RELION that allows it to link together the tilt series alignment parameters, the CTF parameters, the dose information, the microscope information and the paths to the unaligned tilt series.<br>
- a particle.star generator: this script converts a dynamo particle table in a RELION-readable format. This allows RELION to start its 3d auto-refine with particles that of sufficient quality.<br><br><br>
By using these tools users of dynamo can exploit the functions available in RELION. In this case, RELION will be used to refine the average generated by dynamo. <br><br>
Other interesting ways to exploits the dynamo<->RELION converters can be:<br>
- reconstructing only binned tomograms with dynamo and generate an initial set of particles. Next, refining the unbinned particles in RELION. This approach allows for not needing to reconstruct all of the full-sized tomograms. Indeed, for large datasets and depending on the resources available, it might not be practically possible to store in memory all of the full-sized tomograms.<br>
- By converting the particles from RELION to dynamo, a user can exploit the dynamo functionalities that are not available in RELION, both in data analysis and visualization.

== Input ==

- the unaligned, not dose weighted, not CTF corrected tilt series.<br>
- the AWF folder paths.<br>
- the CTF information.<br>
- the dose information.<br>
- the microscope parameters.<br>
- a particle table.<br>

== Practical notes ==


The dynamo-RELION will be done using the script <code>convert_data.m</code>. Firstly, the script will recover the CTF information, the dose information, the microscope parameters and the alignment parameters previously computed and will generate a tomogram.star file. This file allows RELION to know how to locally reconstruct subtomograms directly from the unaligned, not dose weighted, not CTF corrected tilt series.

<pre style="color: blue">Change the root of convert_data.m in RELION_conversion/scripts</pre>

You can see that the first part of the script are the path to all the files required. RELION implicitly handles dose weighting, ctf correction, tilt series alignment and local reconstructions, so RELION needs access to this information. RELION has an implementation to automatically do this but only for IMOD output folders. A custom function to do the same with dynamo alignments is used in the script <code>convert_data.m</code>.<br><br>
The RELION also need the microscope parameters.
<pre style="color: blue">Update the parameters of convert_data.m</pre>
The second part of the script <code>convert_data.m</code> will convert dynamo particle table in a particle.star file that is meant to be compatible with the previously generated tomogram.star file. A particle.star file is paired with a tomogram.star file in similarly to how a dynamo particles table is paired with a dynamo catalogue of tomograms.<br><br>
Try the conversion:
<pre style="color: red">run convert_data.m</pre>
Two files have been generated: particles_converted.star and tomogram_converted.star.<br>
Inspect tomogram_converted.star:
<pre style="color: blue">open('tomogram_converted.star');</pre>
The script is divided in 6 parts: an header and 5 tomogram blocks. The header includes the information of the microscope parameters as well as the location of the raw tilt series. Each tomogram block is split in as many lines as micrographs present in that tilt series. Each lines encodes the alignment parameter, the ctf parameters and the dose.<br><br>
Inspect particles_converted.star:
<pre style="color: blue">open('particles_converted.star');</pre>
You will see the files has an header as well as a body with as many lines as particles. Each line includes the alignment parameters of a particle and the identification of which tilt series generates that particle.<br><br>
From here on everything can be done with RELION. For example, the particles can be shifted in Z so that the GAG layer will be centered. Do this with:
<pre style="color: red">system('relion_star_handler --i particles_converted.star --o particles_converted_shifted.star --center --center_Z 32');</pre>

== Outputs ==

- a tomogram.star.<br>
- a particle.star.<br>

= Final RELION refinement =
== Initial comments ==

As said in the previous section, RELION will be used to refine the average generated by dynamo. Since now the converted star file have been generated, all of the folowing steps are done in RELION.<br>
RELION will be used to locally backproject the particles, create an initial template and run a 3d auto-refine project. These are the RELION equivalent steps of the dynamo steps of particles cropping, particle averaging and STA.<br><br>
This step is strongly inspired by the RELION Subtomogram tutorial: https://relion.readthedocs.io/en/release-4.0/STA_tutorial/index.html. In particular the "Make pseudo-subtomograms", "Reconstruct particle", "High-resolution 3D refinement" and "Masks & Postprocessing" steps. Inspection of the RELION Subtomogram tutorial can be used to find how other RELION tools can be used to improve the resolution even further. However, in this workflow, only the basic RELION steps will be exploited.

== Input ==

- the tomogram.star.<br>
- the particle.star.<br>
- The usual FSC mask.<br>

== Practical notes ==

The Reconstruct particle job will be the first step. This will generate the template for the alignment.<br>
The inputs will be the tomogram.star and the particle.star generated previously and the template will be generated of size 192x192x192 unbinned pixels by local back-projection of 2D patches of dimension 512x512.<br><br>

'''Change your working folder to the Final_RELION_refinement/scripts subfolder.'''

The following steps indicate how to set up the RELION job.
<pre style="color: red">open the RELION4 tomogram GUI (relion --tomo)</pre>
<pre style="color: red">Click on the RELION4 tomogram GUI Tomo reconstruct particles</pre>

It will appear like this:
[[File:RELION_reconstruct_IO_v2.png|centre|thumb|600px| RELION Tomo reconstruct particle I/O field.]]
<pre style="color: red">Fill in the tomogram.star and the SHIFTED particle.star that was created in the previous step.

A version of these star files was precomputed and can be found in:
data/RELION_initial_starfiles/tomo.star
data/RELION_initial_starfiles/particle_zp32.star

Then go on the Average page
</pre>

It will appear like this:
[[File:RELION_reconstruct_average.png|centre|thumb|600px| RELION Tomo reconstruct particle average field.]]

<pre style="color: red">Set the parameters like in the previous figure

Go on the Running page

This is where you can set the computing resources available to RELION.
</pre>
[[File:RELION_reconstruct_Running.png|centre|thumb|600px| RELION Tomo reconstruct particle running field.]]
<pre style="color: red">
Once this is done, click the Run! button.
</pre>

The generated <code>merged.mrc</code> average should present the general structure of the GAG. It should look like the left-most average of the following figure:

[[File:range_of_relion_initial_positions.png|centre|thumb|600px| RELION Make pseudo-subtomos I/O field.]]

The creation of all of the particles is done with the Make pseudo-subtomos job. The same size as the template are used.<br>
<pre style="color: red">This Make pseudo-subtomos step is described here put will not be done because of time limitations</pre>
A RELION job can be chained by giving as input of a job the optimisation set of a previous completed job.


[[File:RELION_pseudo_IO.png|centre|thumb|600px| RELION Make pseudo-subtomos I/O field.]]


[[File:RELION_psuedo_reconstuct.png|centre|thumb|600px| RELION Make pseudo-subtomos reconstruct field.]]

[[File:RELION_pseudo_Running.png|centre|thumb|600px| RELION Make pseudo-subtomos running field.]]

Finally a 3D auto-refinement job can be started. The following parameters can be given:<br><br>
<pre style="color: red">This Refine3D step is described here but will not be done because of time limitations. This step took 40h on 8 GPUs and the computed results are provided in the folder:

data/RELION_final_results/refinement
data/RELION_final_results/postprocess</pre>
Here are some notes on the chosen parameters:<br>
- the reference can be band-passed at 12 A. A non-bandpassed reference could generate artefacts.<br>
- the reference will be C6.<br>
- a convenient Mask diameter is 230 A. This created a spherical mask that is used in addition to the alignment mask.<br>
- an initial angular sampling of 3.7 degrees and an initial offset range of 5 pixels with step of 1 pixel.<br>
- If possible, GPU computing should be used with the identifications of all GPUs<br>
- RELION requires a certain number of mpi processes to be defined in the Running tab. In the example, by setting 8 GPUs in the compute tab, we require 9 mpi processes: 4 for the even particles, 4 for the odd particles and one mpi master process.<br>
- By default RELION uses all of the GPU memory available.  GPU device memory (in Mb) to leave free after allocation can be set with the flag --free_gpu_memory.<br>

[[File:RELION_refine_IO.png|centre|thumb|600px| RELION 3D auto-refine I/O field.]]
[[File:RELION_refine_Reference.png|centre|thumb|600px| RELION 3D auto-refine reference field. Note the low-pass and the symmetry.]]
[[File:RELION_refine_CTF.png|centre|thumb|600px| RELION 3D auto-refine CTF field.]]
[[File:RELION_refine_optimisation.png|centre|thumb|600px| RELION 3D auto-refine optimisation field. Note the mask diameter.]]
[[File:RELION_refine_autosampling.png|centre|thumb|600px| RELION 3D auto-refine auto-sampling field.]]
[[File:RELION_refine_compute.png|centre|thumb|600px| RELION 3D auto-refine compute field. Note that the ids of the GPUs are separated by the colon symbol ":".]]
[[File:RELION_refine_Running.png|centre|thumb|600px| RELION 3D auto-refine running field. Note the number of mpi processes is one more than the number of GPUs]]

The final refined and postprocessed average can be found:

<pre style="color: blue">relion_refined_average = dread([root,'data/RELION_final_results/postprocess/postprocess.mrc']);</pre>

<pre style="color: blue">Open the [root,'data/RELION_final_results/postprocess/postprocess.mrc'] file with Chimera, set the step to 1 and use the tool/Hide Dust function to clear the view.</pre>

[[File:final_chimera_HIV1.png|centre|thumb|600px| Chimera screenshot with Hide dust tool.]]

[[File:final_FSC_curve.png|centre|thumb|600px| RELION FSC computed curve of the HIV1 average.]]

== Outputs ==

- two half-maps generated using the gold standard protocol in RELION.<br>
- the FSC curve.<br>
- the high-resolution average RELION.<br>
- the aligned particle positions in RELION.<br>