PyTom: Align subtomograms using spherical harmonic analysis
This approach achieves basically the same task as stated here, but much faster and more accurate. Some other functionalities are also included, such as resolution determination according to gold-standard FSC and CTF correction using Wiener filter.
This is by no means a comprehensive document about the mathematical details. For that, please refer to the paper: Fast and accurate reference-free alignment of subtomograms, Y. Chen et al., JSB 2013.
If you are a advanced user and for some reason you want to integrate this fast algorithm into your own software, it is also possible and quite straight forward. Check out the underlying library SH Alignment here (written in C and Python) and follow the instruction described below.
Here we assume you already have the subtomograms to be aligned, either by template matching or any other means. What matters here is the subtomograms on disk and the corresponding particle list file describing all the relevant information (please check here if you already have the subtomograms and want to generate the particle list).
Depending on the purpose, there are four scripts that could be used for the alignment, which are all contained in the module
pytom.frm. They can be classified into two categories: with-/without gold standard, with-/without Wiener filtering. You should choose one that fits to your problem. The ways to use all these scripts are similar and they are described below.
|Non gold standard||Gold standard|
|Non Wiener filter||FRMAlignment.py||GFRMAlignment.py|
This script is the most basic one, without gold standard FSC and CTF correction. If you just want to align the volumes and nothing more, this is the one you should use. The call of this script with MPI is:
mpirun --hostfile "pathToYourHostfile" -c "numberOfCPUs" pytom PathToPytom/frm/FRMAlignment.py -j job.xml -v
frm/createJob.pyor any text editor if you are very familar with the XML format. An example of the job file would be:
<FRMJob Destination='.' BandwidthRange='[4, 64]' Frequency='10' MaxIterations='10' PeakOffset='10' AdaptiveResolution='0.1' FSC='0.5'>
<Reference PreWedge="" File="YOUR_REFERENCE.em" Weighting="">
<Mask Filename="YOUR_MASK.em" Binning="1" isSphere="True"/>
<SampleInformation PixelSize="XX_ANGSTROM" ParticleDiameter="XX_ANGSTROM"/>
Note not all the fields in this file need to be set. There are quite some parameters are unused and exist only for historical reasons or compatibility to other parts in Pytom. If
createJob.py is used, you would be only required to enter the relevant ones.
Here all the fields can be or need to be set are explained here:
10 ∗ (1 + 0.1) = 11band in 6th iteration. This parameter is useful to alleviate the noise bias and still let the alignment "see" more information in each round so that the local minimum is avoided.
<AngularConstraint Type="Fixed Angle" Phi="0" Psi="0" Theta="0" Nearby="10" />
<AngularConstraint Type="Fixed Axis" X="0.267261241912" Y="0.534522483825" Z="0.801783725737" Nearby="10" />
<AngularConstraint Type="Adaptive Angle" Nearby="10" />
For 'reference-free' alignment, i.e., alignment without using a
This function generates an average from the particles using random orientations.
This is the version with gold standard FSC, which means the dataset is randomly split into two half sets and aligned independantly. The resoluton is then determined by FSC between the averages from the two half sets. Although in the single particle field the suggested FSC criterion is 0.143, For alignment in electron tomography it would be better to set to 0.3 or 0.5, because this can prevent the determined frequency to jump too further at a time. In the end, you can still use 0.143 criterion to determine the final resolution.
As you can see, nothing much is different except to tell the program to run two independant alignments on one dataset. So is the job description file:
<Reference PreWedge="" File="YOUR_REFERENCE1.em" Weighting="">
<Reference PreWedge="" File="YOUR_REFERENCE2.em" Weighting="">
You have to give two references in the job file instead of one. These two references could be the same, or not. The way to call this script is identical to the one above except for the script name:
mpirun --hostfile "pathToYourHostfile" -c "numberOfCPUs" pytom PathToPytom/frm/GFRMAlignment.py -j job.xml -v
This is the version with CTF correction using Wiener filter. For that you have to prepare two sets of particles: the phase flipped volumes and the CTF convoluted volumes. Currently PyTom does not have the functionality to prepare them for you. It has to be done yourself using whatever software you like. For the Wiener filter equation, please refer to the paper.
The job description file is still similar and looks like this:
<ParticleListPair PhaseFlippedParticleList='PARTICLE_LIST1.xml' CTFConvolutedParticleList='CONVOLUTED_DIR' CTFSquared='CTF_SQUARED_VOLUME.em' SNR='10'/>
<ParticleListPair PhaseFlippedParticleList='PARTICLE_LIST2.xml' CTFConvolutedParticleList='CONVOLUTED_DIR2' CTFSquared='CTF_SQUARED_VOLUME2.em' SNR='10'/>
As you can see, the difference is in the particle list part:
The call of this script is the same and omitted here.
This is the version with gold standard FSC and CTF correction. Its job file can be set analogously to the above one but with one more reference field.
If you want to integrate this fast alignment algorithm into your own software, you can check out the low level library written in C and Python at this page. But please do not forget to cite the paper!
It contains some documentations and examples to use to help you understand. It should be quite straight forward to comprehend and adapt to your purpose. Most importantly, pay attention to the python function
sh_alignment.frm module, and the swig wrapper of C functions:
sh_alignment.swig_frm.frm_fourier_corr. The former is for finding the best translation and rotation between two volumes. And the latter is for computing the cross correlation of two spherical functions (real or complex).
If you have further questions, please email to: email@example.com.