## page was renamed from EMAN2/e2gmm
= Structural variability analysis with Gaussian mixture model =

Available in EMAN release after 10/01/2020. Still under development. For details of the method, refer to [[https://arxiv.org/abs/2101.10356| Chen2021]].
Updates:
 * After 03/22/2021, a new GUI program, '''e2gmm_analysis''' is added to visualize the motion trajectory of the Gaussian model after training.
 * After 05/05/2021, a new '''--conv''' option is added that uses a convolutional neural network structure and improves convergence. 
 * If you are interested in how this works, there are now Jupyter notebooks for the two tutorial datasets in Github [[https://github.com/cryoem/eman2-jupyter/tree/master/tutorials|here]]. Notebooks of more advanced usage of the method will be uploaded later.

== Dataset ==

First, for the compositional heterogeneity analysis, we use the dataset of L17-Depleted 50S Ribosomal Intermediates from [[https://www.ebi.ac.uk/pdbe/emdb/empiar/entry/10076| EMPIAR-10076]]. In this tutorial, we use a subset of only 3000 particles of box size 128. The small dataset can be downloaded [[https://drive.google.com/file/d/19jP_XmmTtlGPznHS5LtnOauQ6tkPY8Ww |here]]. The particles are artificially selected from the full dataset to keep the tutorial is as simple as possible. Inside the zip file is a minimum single particle refinement project from the new e2spa pipeline. This includes a '''particles''' folder that contains the particles, and a '''r3d_00''' folder with particle alignment information, mask, and the 3D structure. There is also a '''strucfac.txt''' file that is used for optimal sharpening of the density maps.

 {{attachment:particle_input.png| Particle input |width=400}}

If you have previously run e2refine_easy on this data, you can extract the necessary files from a refine_XX folder using e2evalrefine.py with the --extractorientptcl option, and then just use the corresponding threed_XX file as a starting volume for the analysis. 

The single model refinement can also be performed with the new e2spa_refine.py pipeline, which directly outputs the aligned list file and average density map that can be used as input.

To convert a refinement from Relion, a script called e2convertrelion.py can be found in the examples folder (after 01/2021) that produces phase flipped particles and the aligned list file.


== Build the Gaussian model ==

First, create a folder to contain related files.
{{{
mkdir gmm_00
}}}

To limit memory consumption, the model is trained with 2D projections of the model instead of the 3D map. So here we make projections of the averaged map first. If you start from a e2refine_easy run, the projections can also be found in the refinement folder.
{{{
e2project3d.py r3d_00/threed_01.hdf --outfile gmm_00/projections.hdf --orientgen=eman:delta=5 --parallel thread:12
}}}

Now we begin training the Gaussian model from scratch. We start from a tiny box size and a small amount of Gaussian functions so it converges faster. Then run with a Fourier box size of 48 so the Gaussian model matches the density map at ~15 Å. When increasing the number of Gaussian functions, the new Gaussian will be seeded near the old ones, so the overall shape of the Gaussian model is preserved at each run of the refinement. The mask will exclude the Gaussian centers outside the volume. The last command can be repeated if the Gaussian model does not match well with the density map. 

{{{
e2gmm_refine.py --projs gmm_00/projections.hdf --npt 512 --maxboxsz 24 --modelout gmm_00/model_gmm.txt --niter 20 --mask r3d_00/mask.hdf
e2gmm_refine.py --projs gmm_00/projections.hdf --npt 1024 --maxboxsz 24 --modelout gmm_00/model_gmm.txt --niter 20 --mask r3d_00/mask.hdf --model gmm_00/model_gmm.txt
e2gmm_refine.py --projs gmm_00/projections.hdf --npt 1600 --maxboxsz 48 --modelout gmm_00/model_gmm.txt --niter 20 --mask r3d_00/mask.hdf --model gmm_00/model_gmm.txt --evalmodel gmm_00/model_projs.hdf
}}}

The '''--evalmodel''' option of e2gmm_refine.py will write the projection of the model to an image stack so we can take a look at the 3D reconstruction and check its FSC with the input density map.
{{{
e2make3dpar.py --input gmm_00/model_projs.hdf --output gmm_00/model_avg.hdf --pad 128 --mode trilinear --keep 1 --threads 12 --setsf strucfac.txt
e2proc3d.py gmm_00/model_avg.hdf gmm_00/fsc_map_model.txt --calcfsc r3d_00/threed_01.hdf
}}}

 {{attachment:avgmodel.png| Gaussian model of neutral structure |width=600}}


== Heterogeneity analysis ==

Now we have a Gaussian model of the neutral structure, we can look into the heterogeneity of the particles. Run the e2gmm_refine.py command with particle input and --heter option to let the program embed them to a conformation space. The --gradout option writes the gradient with respect to the neutral structure to a file so the step can be skipped next time (read the saved file using --gradin). The middle layer output is saved in mid_00.txt. The --pas option controls the adjustment of position, amplitude, and sigma of the Gaussian functions. Setting it to 010 means changing amplitude only. Here use the new '''--conv''' mode to speed up convergence. For this dataset the improvement of the convolutional mode is minimal, but there are cases that it makes a significant difference.
{{{
e2gmm_refine.py --model gmm_00/model_gmm.txt --ptclsin r3d_00/ptcls_01.lst --heter --gradout gmm_00/allgrads.hdf --maxboxsz 48 --midout gmm_00/mid_00.txt --pas 010 --conv
}}}

Finally, the e2gmm_eval.py script runs PCA on the middle layer output, classifies the particles, and reconstructs 3D maps of each class. 
{{{
e2gmm_eval.py --pts gmm_00/mid_00.txt --pcaout gmm_00/mid_pca.txt --ptclsin r3d_00/ptcls_01.lst --ptclsout gmm_00/ptcl_cls.lst --setsf strucfac.txt --ncls 3
}}}

 {{attachment:classifyresult.png| Classification result |width=600}}

== Continuous heterogeneity ==

We use another dataset of pre-catalytic spliceosome from [[https://www.ebi.ac.uk/pdbe/emdb/empiar/entry/10180| EMPIAR-10180]] as an example of continuous conformational change. The sample dataset includes 20,000 particles with alignment information in their headers. The temporary download link can be found [[https://drive.google.com/file/d/1th4OP4McoZFWYQcxDzo4SkahawqvNNTq | here]]. Similar to the ribosome dataset, here we intentionally selected a subset of particles that follows a certain motion trajectory. 

First, create a new folder and generate projections of the given 3D structure.
{{{
mkdir gmm_00
e2project3d.py r3d_00/threed_01.hdf --outfile gmm_00/projections.hdf --orientgen=eman:delta=5 --parallel thread:12
}}}

Now build the Gaussian model. We start from lower resolution and fewer Gaussian functions, then increase both to have better convergence. 
{{{
e2gmm_refine.py --projs gmm_00/projections.hdf --npt 512 --maxboxsz 24 --modelout gmm_00/model_gmm.txt --niter 20 --mask r3d_00/mask.hdf
e2gmm_refine.py --projs gmm_00/projections.hdf --npt 1024 --maxboxsz 64 --model gmm_00/model_gmm.txt --modelout gmm_00/model_gmm.txt --niter 20 --mask r3d_00/mask.hdf --evalmodel gmm_00/model_projs.hdf
}}}

Check the Gaussian model and make sure it is consistent with the density map. 
{{{
e2make3dpar.py --input gmm_00/model_projs.hdf --output gmm_00/model_avg.hdf --pad 128 --mode trilinear --keep 1 --threads 12 --setsf strucfac.txt
e2proc3d.py gmm_00/model_avg.hdf gmm_00/fsc_map_model.txt --calcfsc r3d_00/threed_01.hdf
}}}

Now perform the heterogeneity analysis. Here we look at both movement trajectory and amplitude change, so the '''--pas''' option can be omitted. 
{{{
e2gmm_refine.py --model gmm_00/model_gmm.txt --ptclsin r3d_00/ptcls_01.lst --heter --maxboxsz 64 --midout gmm_00/mid_00.txt --conv
}}}

Finally, classify particles based on their distribution in the conformational space. Since this is a continuous motion, we use '''--mode=regress''' to generate 3D density maps along the first eigenvector at a constant spacing. 
{{{
e2gmm_eval.py --pts gmm_00/mid_00.txt --pcaout gmm_00/mid_pca.txt --ptclsin r3d_00/ptcls_01.lst --ptclsout gmm_00/ptcls_cls.lst --setsf strucfac.txt --mode regress --ncls 9
}}}


 {{attachment:regressout0.gif| Regression result |width=300}}  {{attachment:regressout1.gif| Regression result |width=300}}