diff --git a/Tutorial/Euler_Angle_List b/Tutorial/Euler_Angle_List index dea9c265aa0a290e0cba40ec430bdd9bb924d5b6..89248ed319bc3587fd77a13edf455d1834472fb6 100644 --- a/Tutorial/Euler_Angle_List +++ b/Tutorial/Euler_Angle_List @@ -1,2 +1,3 @@ +2 0.836000 1.369000 2.640000 -1.375000 0.400000 3.018000 diff --git a/Tutorial_BioEM/BioEM_Manual.pdf b/Tutorial_BioEM/BioEM_Manual.pdf index f67b1aa94709cb0eb665d871c52a547fcc8e903f..db8d019259bfe86dfd954fc35f0897eef281c089 100644 Binary files a/Tutorial_BioEM/BioEM_Manual.pdf and b/Tutorial_BioEM/BioEM_Manual.pdf differ diff --git a/Tutorial_BioEM/BioEM_Manual.tex b/Tutorial_BioEM/BioEM_Manual.tex new file mode 100644 index 0000000000000000000000000000000000000000..a269fcca234ce1aa9ec7ce1be0324ff9401a4222 --- /dev/null +++ b/Tutorial_BioEM/BioEM_Manual.tex @@ -0,0 +1,1573 @@ +\documentclass[12pt, psamsfonts]{book} +\usepackage[english]{babel} +\usepackage[latin1]{inputenc} +\usepackage[reqno]{amsmath} + +\usepackage{threeparttable,rotating,fancybox} +\usepackage[colorlinks=true]{hyperref} %Insert hyperlinks in latex files + +\usepackage{graphicx} +\usepackage{amsfonts} +\usepackage{latexsym} +\usepackage[dvips]{color} +%% \usepackage{amssymb, amsmath, amsfonts} +%% \usepackage{latin1} +%% \usepackage{color} +%% \usepackage[dvips]{graphics} + +\usepackage{fancyhdr} +\setlength{\headheight}{15.2pt} +\pagestyle{fancy} + +%\usepackage[utf8]{inputenc} + +\pagenumbering{roman} + + + +\newcommand*{\titleGM}{\begingroup % Create the command for including the title page in the document +\hbox{ % Horizontal box +\hspace*{0.1\textwidth} % Whitespace to the left of the title page +\rule{3pt}{\textheight} % Vertical line +\hspace*{0.05\textwidth} % Whitespace between the vertical line and title page text +\parbox[b]{0.75\textwidth}{ % Paragraph box which restricts text to less than the width of the page + +{\noindent\Huge BioEM Manual}\\[2\baselineskip] % Title +{\large \textit{A software for bayesian analysis of EM images}}\\[4\baselineskip] % Tagline or further description +{\Large \textsc{ }}%Pilar Cossio \\ David Rohr \\ Volker Linderstruth \\ Gerhard Hummer}} % Author name + +\vspace{0.5\textheight} % Whitespace between the title block and the publisher +{\noindent Max Planck Institute of Biophysics }\\[\baselineskip] % Publisher and logo +}} +\endgroup} + + +\title{BioEM Manual} + +\author{P.Cossio, D. Rohr, V. Lindestruth, G. Hummer } +\date{March 2015} +\begin{document} +\pagestyle{plain} + +\titleGM + +\tableofcontents + +%\pagestyle{fancy} + + +\chapter{The BioEM software} +\pagenumbering{arabic} + +\section{Introduction} + +The BioEM method calculates the posterior probability +of a model to multiple experimental EM images, using Bayesian analysis. +The key idea is not to modify the raw images but to create a calculated image, from the original model, as similar as possible to the observed experimental image. +The calculated image takes into account the relevant factors in the experiment and image formation, such as the molecule orientation, interference effects, uncertainties in the particle center, offset in the intensities, etc., and, importantly, noise. +Technically, the model is first rotated to a given orientation, then projected along the $z$-axis, +then it is convoluted with a point spread function in Fourier space to cope with imaging artefacts, +next it is shifted by a certain number of pixels to account for the center displacement, +and finally this modified projection is compared to a reference particle-image. +The similarity between the calculated projection and the experimental image, for a given parameter set, is assessed through a likelihood function. +The posterior probability of the model is obtained by integrating the likelihood function, and prior probabilities, over a all parameter ranges. +The BioEM software is used to perform this integration numerically. A detailed description of the BioEM method is found in Refs. \cite{CossioHummer_2013,BioEM_server}. + +\section{Installation} + + +Before installation, there are several programs or libraries that +should be installed on the compute node. In the following, +we will give a brief explanation of the mandatory, and the optional prerequisite programs. + + +\subsection{Prerequisite programs and libraries} + +\begin{itemize} + +\item {\it FFTW library:} is a subroutine library for computing the discrete Fourier transform. +It is specifically used in BioEM, to calculate the convolution of the ideal image with the point spread function (PSF), and +the cross-correlation of the calculated image to the experimental image. FFTW can be downloaded from the webpage www.fftw.org. + +\item {\it BOOST library:} provides support for tasks and structures such as linear algebra, pseudorandom number generation, multithreading, +image processing, and unit testing. In particular, they are used in the code to access and organize input-data. +BOOST can be downloaded from www.boost.org. + +\item {\it OpenMP:} is a programming interface that supports multi-platform shared memory parallel programming. +It is normally, included in the standard GNU or Intel c++ compliers (so no downloading should be necessary). For more information see http://openmp.org/. + +\end{itemize} + +The optional (but +{\bf encouraged} to use) programs for an easy compilation, and optimal performance, are described bellow. + +\begin{itemize} + +\item {\it CMake:} is a cross-platform software for managing the build process of software using a compiler-independent method ({\it i.e.} creating a Makefile). +CMake can be downloaded from www.cmake.org. + +\item {\it CUDA:} is a parallel computing platform implemented by the graphics processing units (GPUs) that NVIDIA\cite{} produce. +Thus, NVIDIA graphics cards are necessary for running BioEM with the CUDA implementation. For more information see +www.nvidia.com. + +\item {\it MPI:} Message Passing Interface is a standardized and portable message-passing system designed to function on a wide variety of parallel computers, with and without shared-memory. +{\bf PC:} Difference between openMPI and MPICH (?). Is one recommend over the other? + +\end{itemize} + + +After these programs are successfully installed in your compute node, +it will be possible to install BioEM. + + + +{\it Note:} It is recommended that the same complier that is used for the +libraries, is also used BioEM. + +\subsection{Download} + +A compressed directory of the BioEM software can be downloaded from [mpi biophys]. +After downloading the {\it tar} file, uncompress by executing +\vspace{0.5cm} + +\fbox{% +\parbox{12cm}{ +{\footnotesize \texttt{tar -zxvf BioEM.tar.gz}}}} + +\vspace{0.5cm} + +In the uncompressed {\it BioEM} directory, there are +\begin{itemize} +\item[--]the source code {\it cpp} files, and {\it include} directory with corresponding header files. +\item[--]the copyright license, and README file. +\item[--]the {\it CMakeLists.txt} file that is necessary for CMake (see section bellow). +\item[--]the {\it Tutorial} directory that includes the example files used in the tutorial (see chapter \ref{tutorial}). Inside this directory, +there is also a directory called {\it MODEL\_COMPARISON}. +\item[--]the {\it Quaternions} directory that includes list files of quaternions that sample uniformly +the rotational group (see section \ref{intor}). +\end{itemize} + + +\subsection{Installing with CMake} + +The easiest installation is done with the CMake program +that generates automatically a Makefile, according to the +specific CPU/GPU architecture, and desired features. CMake +uses the {\it CMakeLists.txt} file, that contains all the instructions to generate the Makefile. +This file is provided in the uncompressed BioEM directory. +At the beginning of the {\it CMakeLists.txt} the modifiable options are provided. +These options should be enabled/disabled ({\bf ON}/{\bf OFF}) according to the desired functionalities. + +\begin{itemize} + +\item To enable or disable {\it OpenMP} + +\fbox{% +\parbox{12cm}{ +{\footnotesize \texttt{option (INCLUDE\_OPENMP "Build BioEM with OpenMP support" ON/OFF)}}}} + +\item To enable or disable {\it MPI} + +\fbox{% +\parbox{12cm}{ +{\footnotesize \texttt{option (INCLUDE\_MPI "Build BioEM with MPI support" ON/OFF)}}}} + +\item To enable or disable {\it CUDA} + +\fbox{% +\parbox{12cm}{ +{\footnotesize \texttt{option (INCLUDE\_CUDA "Build BioEM with CUDA support" ON/OFF)}}}} + + +\item To print out the CMake variables + +\fbox{% +\parbox{12cm}{ +{\footnotesize \texttt{option (PRINT\_CMAKE\_VARIABLES "List all CMAKE Variables" ON/OFF)}}}} + +\item To use {\it CUDA} with the INTEL complier ({\it icc, mpiicc}) +one needs to turn {\bf ON} the \texttt{CUDA\_FORCE\_GCC} variable + +\fbox{% +\parbox{12cm}{ +{\footnotesize \texttt{option (CUDA\_FORCE\_GCC "Force GCC as host compiler for CUDA part (If +standard host compiler is incompatible with CUDA)" ON/OFF)}}}} + +\end{itemize} + +{\it Note:} For certain architectures, additional files ({\it e.g.} FindFFTW.cmake) could also be needed for successful CMake +execution. For more information, on specific CMake features ({\it e.g.} changing compiler) see www.cmake.org. + +\subsubsection{Steps for basic installation} + + +\begin{itemize} + +\item[--] Setup the desired features in the CMakeLists.txt file. +\item[--] Generate a build directory in the main BioEM directory + +\fbox{% +\parbox{10cm}{ +{\footnotesize \texttt{mkdir build}}}} + +\item[--]Access the build directory, and run CMake with the {\it CMakeLists.txt} file + +\fbox{% +\parbox{10cm}{ +{\footnotesize \texttt{cd build \\ + cmake ../CMakeLists.txt}}}} +\item[--] If this process is successful, a Makefile and CMakeFiles directory should be generated in the build directory. +If this is not the case, turn on the \texttt{PRINT\_CMAKE\_VARIABLES} option in the CMakeLists.txt file, and re-run +CMake with verbosity. +\item[--] After generating the Makefile, run it in the build directory + +\fbox{% +\parbox{10cm}{ +{\footnotesize \texttt{make}}}} +\item[--] If this process is successful a bioEM executable should be generated. + +\end{itemize} + +For a simple test, run +\vspace{0.5cm} + +\fbox{% +\parbox{12cm}{ +{\footnotesize \texttt{./bioEM}}}} + +\vspace{0.5cm} + +The output on the screen should be +\vspace{0.5cm} + +\fbox{% +\parbox{12cm}{ +{\footnotesize \texttt{ +Command line inputs:\\ +--Modelfile arg (Mandatory) Name of model file\\ +--Particlesfile arg if BioEM (Mandatory) Name of particles file\\ +--Inputfile arg if BioEM (Mandatory) Name of input parameter file\\ +--PrintBestCalMap arg (Optional) Only print best calculated map (file nec.). + NO BioEM (!)\\ +--ReadOrientation arg (Optional) Read orientation list instead of uniform + grid (file nec.)\\ +--ReadPDB (Optional) If reading model file in PDB format\\ +--ReadMRC (Optional) If reading particle file in MRC format \\ +--ReadMultipleMRC (Optional) If reading Multiple MRCs \\ +--DumpMaps (Optional) Dump maps after they were red from maps file\\ +--LoadMapDump (Optional) Read Maps from dump instead of maps file\\ +--OutputFile arg (Optional) For changing the outputfile name\\ +--help (Optional) Produce help message\\ +}}}} + +\vspace{1cm} + + + +\chapter{BioEM input} + + +The BioEM software's main input is from the command line. Here the filenames +of the model, particle-images and parameter file should be provided. +We will now give a detailed description of each input item. + +\section{Model file} + + +The model is represented as spheres (or points) in 3-dimensional space, with corresponding radius and density. +%either in PDB or +% txt format. +Through the command line one has to provide the name of the file +that contains these items: + +\vspace{0.5cm} + +\fbox{% +\parbox{12cm}{ +{\footnotesize \texttt{ --Modelfile arg}}}} + +\vspace{0.5cm} + +Where \texttt{arg} is the filename. There are two types of model file formats that are read by BioEM: + +\begin{itemize} +\item[--] {\it *.txt *.dat file:} +Useful for all atom representation or 3D voxel representation of density maps. +With format "\%f \%f \%f \%f \%f", + the first three columns are the coordinates of atoms or + voxels, the fourth column is the radius or voxel side length ($\AA$) and the + last column is the corresponding electron density +(Format: {\footnotesize \texttt{ x --- y --- z --- radius --- density }}). + +\item[--] {\it pdb file:} BioEM reads only the C$_\alpha$ atom positions from standard {\it pdb} files, with +their corresponding residue type. The residues are modeled as a sphere, centered at the + C$_\alpha$ with corresponding van-der-Waals radii, and electron density (as in ref. + \cite{CossioHummer_2013}). To read pdb files the following commandline keyword is needed: + +\fbox{% +\parbox{12cm}{ +{\footnotesize \texttt{--ReadPDB}}}} + +{\it Note:} the .pdb extension is not mandatory. + +\end{itemize} + +{\it Additional Feature:} It is possible to model the elements simply as points +(instead of spheres) by projecting only the density of each element. For this, add the keyword \texttt{"NO\_PROJECT\_RADIUS"}, +in the input parameter file (see section \ref{Physparm}). + + + +\section{Particle-image file} + +The name of the experimental image file is needed, as a commandline input in BioEM: +\vspace{0.5cm} + +\fbox{% +\parbox{12cm}{ +{\footnotesize \texttt{ --Particlesfile arg}}}} + +\vspace{0.5cm} + +Where \texttt{arg} is the name of the file containing the experimental particle images. +Two format options are allowed for the these files: + +\begin{itemize} + \item[--] {\it *.txt or .dat file:} Data are formatted as + "\%8d\%8d\%16.8f" where the first two columns are + the pixel indexes, and the third column is the image intensity. + Multiple particles are read in the same file with the + separator \texttt{"PARTICLE" \& particle number}. + Pixel indexes should start at 0, and all pixels should be + included. It is recommended that particles + are normalized to zero average and unit standard deviation. + + \item[--] {\it .mrc file:} BioEM also reads standard {\it .mrc} particle-image + files. To do so, the additional commandline keyword is needed: + +\fbox{% +\parbox{12cm}{ +{\footnotesize \texttt{ --ReadMRC }}}} + + If reading multiple {\it mrc} files, the name of the file containing the {\it list} + of all the {\it mrc} files should be provided. The additional + command is also required + +\fbox{% +\parbox{12cm}{ +{\footnotesize \texttt{--ReadMultipleMRC}}}} + +{\it Example:} + +\fbox{% +\parbox{12cm}{ +{\footnotesize \texttt{ --Particlesfile LIST --ReadMRC --ReadMultipleMRC}}}} + +\texttt{LIST} is the name of the file containing the list of the multiple {\it mrc} files. + + +{\it Notes:} The {\it .mrc} extension is not mandatory. +By default, when reading {\it mrc} particles, the intensities are normalized to +zero average and unit standard deviation. Use keyword \texttt{NO\_MAP\_NORM} to unset this default. + +\end{itemize} + + +{\it Additional Features:} These are some useful commandline keywords for executing multiple times BioEM +with a large image set: + +\vspace{0.5cm} + +\fbox{% +\parbox{12cm}{ +{\footnotesize \texttt{--DumpMaps}}}} + +\vspace{0.5cm} + +writes out a file {\it maps.dump} containing the particle maps in binary format for a faster re-reading. + + To read the dumped maps, use + +\vspace{0.5cm} + +\fbox{% +\parbox{12cm}{ +{\footnotesize \texttt{--LoadMapDump }}}} + + +\vspace{0.5cm} + + See the chapter \ref{tutorial} for a detailed description. + + + + + + + +\section{Input-parameter file} + +\label{Physparm} + +BioEM has two sets of parameters. +One set describes the physical problem, like the number of pixels. +Another set describes the runtime configuration, which involves how to parallelize, whether to +use a GPU, and some other algorithmic settings. These parameters do not (or only slightly) +change the output, but have a large influence on the compute performance. They are passed to +BioEM via environment variables. +The two sets are treated differently, because the first set is related to the actual problem, while +the second set belongs to the compute node where the problem is processed. +In this section, we will only describe the input parameters related to the physical +problem. For a detailed description of the performance variables see chapter \ref{perfparm}. + +Additionally from the model and experimental images input, BioEM needs +an input file describing the physical constraints, and integration +limits of the algorithm: + +\vspace{0.5cm} +\fbox{% +\parbox{12cm}{ +{\footnotesize \texttt{ --Inputfile arg}}}} + +\vspace{0.5cm} +where \texttt{arg} is the filename of the input file. This file contains specific +keywords that will describe the physical conditions, and constraints for the algorithm. +Bellow we will describe, in detail, each keyword for the input parameter file. + +\subsection{Micrograph parameters} +Mandatory inputs for the description of the experimental map are + +\begin{itemize} + +\item \texttt{PIXEL\_SIZE (float)} + + Pixel size in $\AA$ of the experimental micrograph. + +\item \texttt{NUMBER\_PIXELS (int)} + +We assume a square particle-image. Here, \texttt{(int)} is the number of pixels in each + dimension, {\it e.g.} particle-image of 220 x 220 pixles, then \texttt{(int)= 220}. +\end{itemize} + + +In a standard BioEM calculation, the integration over the model orientations, +PSF parameters, and center displacement are performed numerically. +To do so, one needs to define the integration ranges, +and grid spacing for each parameter. These quantities will +depend on the experimental conditions (such as defocus ranges, and preferred particle orientations), +thus should be specified by the user. + +\subsection{Integration of orientations} +\label{intor} +In BioEM, there are two manners of describing the orientations of the model +in 3D space, with the Euler angles or with quaternions. +There are several possibilities for sampling the space of rotations, and +performing the integration. + + +\begin{itemize} +\item {\it Grid-sampling of the Euler Angles ($\alpha,\beta,\gamma$):} +Sampling of the full Euler angle space with uniform distribution of $\alpha \in [-\pi,\pi]$, $\mathrm cos(\beta) \in [-1,1]$ +and $\gamma \in [-\pi,\pi]$. For this case one needs only to provide the number of grid points +in $\alpha$, and $\beta$ (because those in $\gamma$ will be the same as in $\alpha$). +The keywords in the parameter file are + +\texttt{GRIDPOINTS\_ALPHA (int)}\\ +\texttt{GRIDPOINTS\_BETA (int)} + +{\it Note:} To sample uniformly it is recommended that +\texttt{GRIDPOINTS\_ALPHA$\sim$ 2*GRIDPOINTS\_BETA}. + +\item {\it Grid-sampling of quaternions:} +With BioEM it is also possible to generate a grid in quaternion space for +the integration over the orientations. For this option, one should provide the +keyword: + +\texttt{USE\_QUATERNIONS} + +Additionally, to sample the quaternions in a hypercubic-grid, one should include + +\texttt{GRIDPOINTS\_QUATERNION (int)} + +where \texttt{(int)} is the grid spacing in each quaternion dimension $\in [-1,1]$. + + + +\item {\it Read orientations from file:} With this feature the orientational grid points are not +calculated directly in the code, but are instead read from a file. This provides more flexibility in the sampling, and +integration. For this feature, an extra keyword in the command line is necessary: + +\fbox{% +\parbox{12cm}{ +{\footnotesize \texttt{--ReadOrientation arg}}}} + +Where \texttt{arg} is the name of the file containing the orientations. +The first row of the file should have \texttt{(int)} equal to the total number of orientations. + +The orientations can be described with the standard Euler angles, +or with quaternions. +\begin{itemize} +\item[--] The format for the Euler angle file should be "\%12.6f\%12.6f\%12.6f", ordered +as $\alpha,\beta,\gamma$, respectively. +\item[--] The format for the file containing the quaternions should be: +"\%12.6f\%12.6f\%12.6f\%12.6f". +One should recall that to use quaternions the extra keyword \texttt{USE\_QUATERNIONS} in the input parameter file is necessary. + +\item[--] {\it Prior for orientations:} For this input, it is possible to +assign prior probabilities for each orientation. To do so, one should add an extra +column (of format "\%12.6f") that indicates the value of the prior probability for each orientation. + +\end{itemize} + + +\end{itemize} +{\bf Important Remark:} We note that not all of these possibilities sample uniformly the group +of rotations in 3D space ({\it SO3}). +For uniform sampling of {\it SO3}, we recommend +the successive orthonormal images method from Ref. \cite{quaternionsMitchell}, with quaternions. +In the directory {\it Quaternions}, we provide lists of quaternions that sample uniformly {\it SO3} +using this method. Several lists with different number of total orientations are provided. +We highly recommend to use these files instead of using trivial grid-sampling of the Euler angles or quaternions. + + + + +\subsection{Integration of the CTF or PSF} + +\subsubsection{Integration in Fourier space using the CTF and Envelop:} +To take into account the interference effects in the +experiment, we multiply the idea image $I_0$ from a model with the {\bf C}ontrast {\bf T}ransfer {\bf F}unction (CTF) in Fourier space, + $\mathrm{CTF}(s)\mathcal{F}(I_{0})$, where $s$ is the radial spatial frequency, and $\mathcal{F}$ indicates the Fourier transform. +An approximate expression for the CTF is $-A\cos(as^2/2)-\sqrt{1-A^2}\sin(as^2/2)$, where the parameter $a=2\pi \lambda \Delta f$ and $\lambda$ is the +electron wavelength, and $\Delta f$ is the defocus. $A \in [0,1]$ is a parameter that establishes the contributions of the cosine and sine components. +In addition, the CTF is normally modulated by an envelop function (see Eqs. 2 and 3 in Ref. \cite{bioEM}) with $\mathrm{Env}(s)=e^{-bs^2/2}$, +where the parameter $b=B/2$ is half of the B-factor\cite{PenzekXXX}. %We note that all this analysis is done in Fourier space. +To calculate the posterior probability, one must integrate numerically the three parameter $a,b$ and $A$, or equivalently, the defocus $\Delta f$, B-factor $B$, +and amplitude $A$. +To do so, one should include in the input parameter file the start and end limits, and number of grid points in each integration: + +\textit{Parameter -- (start) -- (end) -- (gridpoints)}\\ +\indent \texttt{CTF\_DEFOCUS (float) (float) (int)}\\ +\indent \texttt{CTF\_B\_FACTOR (float) (float) (int)}\\ +\indent \texttt{CTF\_AMPLITUDE (float) (float) (int)}\\ + +The defocus units should be in micro-meters, and that of the B-factor should be in \AA$^2$. + +\subsubsection{Integration in real space using the PSF:} + +The {\bf P}oint {\bf S}pread {\bf F}unction (PSF), that is the real-space equivalent of the contrast transfer function, +is defined in the Supplementary Information of Ref. \cite{BioEM_server}. +Here, an ideal calculated image is convoluted with the PSF to mimic the interference effects in the imaging experiment. +Similarly as the CTF, the PSF has three real-space variables: amplitude $A_R$, envelop $\chi$ and phase $\theta$ (following the same +notation as in Ref. \cite{BioEM_server}). +The amplitude defines the contribution of the PSF real-space sine or cosine parts, and is within $[0,1]$. +The envelop is the real space equivalent of the B-factor, and should be given in +units of $\AA^{-2}$. The phase is also in real space, and should be given in units of $\AA^{-2}$. +The specific ranges of these parameters will depend on the imaging conditions of each experiment +(defocus, astigmatism etc). + +To use the PSF (CTF analysis is default), one should include the keyword: + +\texttt{USE\_PSF} + +The information of the integration limits and grid points should be included in the input file as: + +\textit{Parameter -- (start) -- (end) -- (gridpoints)}\\ +\indent \texttt{PSF\_ENVELOPE (float) (float) (int)}\\ +\indent \texttt{PSF\_PHASE (float) (float) (int)}\\ +\indent \texttt{PSF\_AMPLITUDE (float) (float) (int)}\\ + + +{\it Additional feature:} One can print out the corresponding CTF parameters that maximize the posterior, + with the keyword in the parameter file + \texttt{WRITE\_CTF\_PARAM}\\ + +{\it Note:} It is important to remark that one of these two procedures (either the CTF or PSF) is mandatory in the BioEM analysis. + + +\subsection{Integration of center displacement} + +The integration of the particle center is done over a square and uniform grid. +The particle, along both directions, is translated from its center up to a maximum distance ({\it max displ.}). +Users should provide this maximum displacement and the integration grid spacing in units of pixels. + +The keywords in parameter file are:\\ +\textit{Parameter - (max displ.) - (grid-space)}\\ + \texttt{DISPLACE\_CENTER (int) (int)}\\ + + \vspace{0.2cm} +{\it Example:} If \texttt{[DISPLACE\_CENTER 10 2]}, the integration will be done +between $[x_c-10,x_c+10]$ along $x$ (where $c$ denotes center), and $[y_c-10,y_c+10]$ along $y$, with sampling every 2 pixels. + +\vspace{0.5cm} +{\bf Important Remark:} +The integration of the {\it normalization}, {\it offset} and {\it noise} parameters is carried out analytically. See +Supplementary Information of Ref. \cite{CossioHummer_2013}. + + +\section{Additional features} +\label{crosscor} + +Besides calculating the BioEM probability, one can use the code to extract other features +that can be useful in the analysis of the data. + +These extra features are: + +\begin{itemize} + +\item {\it Model prior probability:} It is possible to +include the model's prior using the keyword in the parameter file: \\ +\texttt{PRIOR\_MODEL (float)}\\ +where \texttt{(float)} is the numerical value of the model's prior. + + +\item {\it Posterior probability as a function of orientations:} +One can write out the log-posterior as a function of each orientation. +In this case, the integration is performed over PSF, center, normalization, offset and noise, +but not over the orientations. The keyword in parameter file is\\ + \texttt{WRITE\_PROB\_ANGLES}\\ +With this feature there is an additional output file \texttt{ANG\_PROB} where the results are written. + +\item {\it Cross Correlation calculation:} +It is also possible to write out the best cross-correlation as a function +of the pixels. In this case, the standard integration over the specified parameters is performed, and the cross correlation is stored, +using Bayesian analysis. The optimal cross-correlation at pixel $x$,$y$ is +\begin{equation} + CC^{opt}(x,y)=-\ln \int e^{-\left[I^{\mathrm{obs}}(x,y)-I(x,y|m,\boldsymbol\theta)\right]^{2}} d\boldsymbol\theta +\end{equation} +where $I^{\mathrm{obs}}(x,y)$ is the experimental image, and $I(x,y|m,\boldsymbol\theta)$ is the calculated image +generated from model $m$ with parameters $\boldsymbol\theta$. +The integral is over all parameters apart from the center displacement. + +%To optimize speed, a grid can be used ({\it i.e.} instead of storing $CC^{opt}(x,y)$ at every pixel, +%on can store it every \texttt{[X,Y]} pixels). + The necessary keywords in the input parameter file is \\ + \texttt{ WRITE\_CROSSCOR}\\ + Additional keywords in the parameter file are: + \begin{itemize} + \item[] \texttt{FLIPPED}\\ + {\it Important:} If the micrograph intensities are inverted (particles are white instead of black).\\ + \item[] \texttt{CROSSCOR\_NOTBAYESIAN}\\ + If only the maximum of $CC^{opt}(x,y)$ is stored instead of the sum over the parameters. + \end{itemize} +There is an additional outputfile called \texttt{CROSS\_CORRELATION} where the cross-correlations are written +(see chapter \ref{tutorial} for details on the output file). +This feature can be used for particle identification in entire micrographs, so it it recommend +to have few input particle images, and few integration grid points. +The cross-correlation is written every 10 pixels, in each dimension, by default. + +\item {\it Print projection map:} +This option is {\bf completely independent} of the BioEM posterior probability calculation. +The idea is to print a synthetic image from a model given a specific set of parameters ({\it e.g.} those +that maximize the posterior). The commandline keyword is + +\vspace{0.5cm} + +\fbox{% +\parbox{12cm}{ +{\texttt{--PrintBestCalMap arg}}}} + +\vspace{0.5cm} + +where \texttt{arg} is the name of the inputfile containing the specific parameter set (note that +this input file is different from that used in the standard BioEM calculations). In this case the only file +\texttt{BESTMAP} with the printed map will be generated. See chapter \ref{tutorial} for examples. + +\end{itemize} + +{\it Note:} The previous additional options are mutually exclusive. + +\section{List of keywords for input-file} + +\subsection{BioEM posterior probability calculation:} +\begin{itemize} + + +\item \texttt{PIXEL\_SIZE (float)}: Indicate the micrograph pixel size in \AA. +\item \texttt{NUMBER\_PIXELS (int)}: Assuming a square particle-image, one should indicate the number of pixels along an axis. +This should coincide with the number of pixels read form the micrograph. + +\item \texttt{GRIDPOINTS\_ALPHA (int)}: (Integration of Orientations) Number of grid points used in the integration +over Euler angle $\alpha$. For this keyword of uniform-grid Euler angle sampling, Euler angle $\gamma$ is integrated in the same manner. +\item \texttt{GRIDPOINTS\_BETA (int)}: (Integration of Orientations) Number of grid points used in the integration +over $\cos(\beta)$. + +\item \texttt{USE\_QUATERNIONS}: (Integration of Orientations) If using quaternions to the describe the orientations. +\item \texttt{GRIDPOINTS\_QUATERNION (int)}: (Integration of Orientations) For a cubic-grid for quaternion sampling ($\in [-1,1]$). +\texttt{(int)} is the number of grid points per dimension. + +\item \texttt{CTF\_DEFOCUS (float) (float) (int)}: (CTF Integration) Grid sampling of CTF defocus, $\Delta f$. Units of micro-meters. +\texttt{(float) (float)} are the starting and ending limits (e.g., 1.5 6.5 [$\mu$m]), and \texttt{(int)} is the number of grid points. + +\item \texttt{CTF\_B\_FACTOR (float) (float) (int)}: (CTF Integration) Grid sampling of CTF B-factor, $B$. Units of \AA$^2$. +\texttt{(float) (float)} are the starting and ending limits (e.g., 100.4 300.5 [\AA$^2$]), and \texttt{(int)} is the number of grid points. + +\item \texttt{CTF\_AMPLITUDE (float) (float) (int)}: (CTF Integration) Grid sampling of the CTF amplitude, $A$ (adimensional $\in [0,1]$). +\texttt{(float) (float)} are the starting and ending limits (e.g., 0.1 0.7), and \texttt{(int)} is the number of grid points. + +\item \texttt{USE\_PSF}: (PSF Integration) If using the point spread function to describe the interference effects. + +\item \indent \texttt{PSF\_PHASE (float) (float) (int)}: (PSF Integration) Grid sampling of PSF phase, $\theta$. Units of \AA$^{-2}$. +\texttt{(float) (float)} are the starting and ending limits (e.g., 0.002 0.008 [\AA$^{-2}$]), and \texttt{(int)} is the number of grid points. + +\item \indent \texttt{PSF\_ENVELOPE (float) (float) (int)}: (PSF Integration) Grid sampling of PSF Envelop, $\chi$. Units of \AA$^{-2}$. +\texttt{(float) (float)} are the starting and ending limits (e.g., 0.0002 0.0008 [\AA$^{-2}$]), and \texttt{(int)} is the number of grid points. + +\item \indent \texttt{PSF\_AMPLITUDE (float) (float) (int)}: (PSF Integration) Grid sampling of the PSF amplitude, $A$ (adimensional $\in [0,1]$). +\texttt{(float) (float)} are the starting and ending limits (e.g., 0.1 0.7), and \texttt{(int)} is the number of grid points. + +\item \texttt{DISPLACE\_CENTER (int) (int)}: (Integration of particle center displacement) Sampling within a square grid. Units of pixels. + \texttt{(int) (int)} are the maximum displacement from the center in both directions, and the grid spacing, respectively. + + +\end{itemize} + +\subsection{Optional:} + +\begin{itemize} +\item \texttt{NO\_PROJECT\_RADIUS}: If atomic positions instead of coarse-grained or Van der Waals radii are used. +\item \texttt{ELECTRON\_WAVELENGTH (float)}: To change the default value of the electron wavelength \texttt{(float)} used in the CTF phase (with defocus) calculation. +Default 0.0220 $\AA$. +\item \texttt{PRIOR\_MODEL (float)}: Prior probability of model. Default 1. +\item \texttt{PRIOR\_ANGLES}: To read the prior of each orientation in the input file of orientations (\texttt{--ReadOrientation arg} from command-line). +\item \texttt{NO\_MAP\_NORM}: Condition to not normalize to zero mean and unit variance the input maps. + + +\end{itemize} + + +\subsection{Additional features:} + +\begin{itemize} +\item \texttt{WRITE\_CTF\_PARAM}: To writeout the CTF parameters that maximize the posterior. + +\item \texttt{WRITE\_PROB\_ANGLES}: To writeout the posterior as a function of the orientation. + +\item \texttt{WRITE\_CROSSCOR}: To writeout the posterior as a function of the pixels in the micrograph. + +\item \texttt{CROSSCOR\_NOTBAYESIAN}: To be used with \texttt{WRITE\_CROSSCOR} for calculating the cross-correlation +as a function of pixels in a none Bayesian manner (see above). + +\item \texttt{FLIPPED}: To be used with \texttt{WRITE\_CROSSCOR} if the micrograph intensities are flipped (particles are white instead of black). + +\item \texttt{IGNORE\_CROSSCORR\_OFFSET}: To be used with \texttt{WRITE\_CROSSCOR} if the intensity offsets should be disregarded. + + + +\end{itemize} + +\subsection{Input for printing a 2D map from model (no BioEM):} + +This option is {\bf completely independent} of the BioEM posterior probability calculation. +The idea is to print a synthetic image from a model given a specific set of parameters. +\vspace{0.5cm} +\fbox{% +\parbox{12cm}{ +{\texttt{--PrintBestCalMap arg}}}} +\vspace{0.5cm} + +where \texttt{arg} is the name of the inputfile containing the specific parameters to construct the map: + +\begin{itemize} +\item \texttt{NUMBER\_PIXELS (int)}: Same as above. +\item \texttt{PIXEL\_SIZE (float)}: Same as above. +\item \texttt{BEST\_DX (int)}: Displacement in direction $x$. +\item \texttt{BEST\_DY (int)}: Displacement in direction $y$. +\item \texttt{BEST\_CTF\_DEFOCUS (float)}: CTF defocus. +\item \texttt{BEST\_CTF\_B\_FACTOR (float)}: CTF B-factor. +\item \texttt{BEST\_CTF\_AMP (float)}: CTF amplitude. +\item \texttt{BEST\_NORM (float)}: Map normalization. +\item \texttt{BEST\_OFFSET (float)}: Map offset. +\item \texttt{BEST\_ALPHA (float)}: Euler angle $\alpha$. +\item \texttt{BEST\_GAMMA (float)}: Euler angle $\gamma$. +\item \texttt{BEST\_BETA (float)}: Euler angle $\beta$. +\item \texttt{WITHNOISE (float)}: (Optional) if noise is added. \texttt{(float)} is the signal-to-noise-ratio. +\item \texttt{NO\_PROJECT\_RADIUS}: (Optional) if atomic positions instead of coarse-grained or Van der Waals radii are used. +\item \texttt{PRINT\_ROTATED\_MODELS}: (Optional) to print 3D coordinates of rotated model. +\item \texttt{USE\_PSF}: (Optional) To use the PSF. +\item \texttt{BEST\_PSF\_PHASE (float)} PSF phase. +\item \texttt{BEST\_PSF\_ENVELOPE (float)} PSF envelop. +\item \texttt{BEST\_PSF\_AMP (float)}: PSF amplitude. + +\end{itemize} + +\label{perfparm} + +These variables enhance or modify the code's computational performance according to the specific +node characteristics, without modifying the numerical results. +They should be tuned for every specific problem where BioEM is executed. + +\chapter{Performance} +\section{Ways of parallelization} +\label{wayparallel} + +BioEM compares various projections of a model to a set of reference particle-images. +Technically, the model is first projected along a given angular orientation, +then it is convoluted with a point spread function in Fourier space to cope with imaging artefacts, +next it is shifted by a certain number of pixels to account for the center displacement, +and finally this modified projection is compared to a reference particle-image. +From a computational complexity perspective, the comparison of a given projection and convolution +with all the reference map is by far the most time consuming part. Calculating the projection and +performing the convolution both make up for only 0-5\% of the overall execution time (depending on +the number of maps and convolution kernels). +BioEM facilitates this task using an all model-projection to an all-particle image comparison through a nested loop. +As shown in Fig. 2 (BioEM software paper), +the outermost loop is over the orientations. Then there is a loop over the different PSF kernels, +finally the innermost loop iterates over all maps and center displacements. +This offers, in principle, many dimensions for parallelization. + +In the following , we will explain the types of parallelization used within BioEM, list all relevant +environment variables, and provide some suggestions for runtime configurations in different +situations. + +\begin{itemize} +\item {\it MPI:} BioEM uses MPI to parallelize over the orientations in the outermost loop. In this case +the probabilities for all reference maps / PSF kernels / center displacements are calculated for a +certain subset of orientations by each MPI process. Afterward, the probabilities computed by +every MPI process are reduced to the final probabilities. If started via \texttt{'mpirun'}, BioEM will +automatically distribute the orientations evenly among all MPI processes. +\item {\it OpenMP:} BioEM can use OpenMP to parallelize over the particle images in the innermost loop. As processing of +these maps is totally independent, there is no synchronization required at all. BioEM will +automatically multithread over the maps. The number of employed threads can be controlled with +the standard + +\fbox{% +\parbox{12cm}{ +{\texttt{export OMP\_NUM\_THREADS=[x]}}}} + +environment variable for OpenMP, where \texttt{[x]} is the number of OpenMP threads. +\item {\it Graphics Processing Unites (GPUs):} BioEM can use GPUs to speed up the processing. In this case, +the innermost loop over all particle-images, and with all center displacements, is processed by the GPU. +The projections and the PSF convolutions are still processed by the +CPU. This process is pipelined such that the CPU prepares the next projections, and PSF convolutions +while the GPU calculates the probabilities to all particle-images for the previous calculated images. Hence, +this is a horizontal parallelization layer among the particle images with an additional +vertical layer through the pipeline. Usage of GPUs must be enabled with the + +\fbox{% +\parbox{12cm}{ +{\texttt{export GPU=1 }}}} + +environment variable. One BioEM process will always only use one GPU, by default the fastest one. A GPU +device can be explicitly configured with the environment variable: + +\fbox{% +\parbox{12cm}{ +{\texttt{export GPUDEVICE=[x]}}}} + +Multiple GPUs can be used through MPI. In this case, every GPU will process all particle-images but calculate the +probabilities only for a subset of the orientations (see description of MPI above). Selection +of GPU devices for each MPI processed must be carried out by: +\begin{itemize} +\item Either the scheduler must mask all GPUs but one, such that the MPI processes sees only 1 +GPU device. +\item Or the \texttt{GPUDEVICE=[x]} environment variable must be set differently for each MPI process. +\item Or it is possible to set \texttt{GPUDEVICE=-1}. In this case the MPI process with rank N on a system +with G GPUs will take the GPU with ID (N \% G). +\end{itemize} +\item {\it GPU / CPU combined processing:} Besides the pipeline approach described in the previous point, +which employs the CPU for creating the calculated image, and the GPU for calculating +the likelihood to all particle-images, there is also the possibility to split the +set of particle-images among the CPU and the GPU. This is facilitated by the + +\fbox{% +\parbox{12cm}{ +{\texttt{export GPUWORKLOAD=[x]}}}} + +environment +variable, which sets the percentage of the number of particle-images processed by the GPU to x\%. The +pipeline is still used, {\it i.e.} in an optimal situation the CPU will: +\begin{itemize} +\item Issue a GPU kernel call such that the GPU calculates the probabilities for x\% of the particle-images +for the current orientation and convolution. +\item Process its own fraction of (100-x)\% of the particle-images for the current orientation and PSF convolution +in parallel to the GPU. +\item Afterward, finish the preparation of the next orientation and convolution before the GPU has +finished calculating the probabilities for the current orientation and PSF convolution. +\end{itemize} + +\item {\it Multiple Projections at once via OpenMP:} BioEM can prepare the projection of multiple +orientations at once using OpenMP. The benefit compared to the pure OpenMP parallelization over +the particle images, however, is tiny. This is relevant if MPI is not used, OpenMP is used, GPU is used, +and if the number of reference maps is small. The number of projections at once +is determined by the environment variable + +\fbox{% +\parbox{12cm}{ +{\texttt{export BIOEM\_PROJECTIONS\_AT\_ONCE=[x]}}}} + +Where \texttt{[x]} is the number of projections that will be calculated simultaneously. + +\item {\it Fourier-algorithm to process all center displacements in parallel:} +BioEM can operate using two algorithms, +a Fourier-algorithm and a non-Fourier-algorithm. Besides numerical effects, both produce +identical results assuming certain boundary conditions. The Fourier-algorithm automatically +takes all displacements into account without having to loop over them. Hence, its runtime is +almost independent from the number of center displacements. The non-Fourier-algorithms is faster for very +few shifts. Benchmarks have shown, that the Fourier-algorithm is faster if the number of shifts +per direction is four or larger, which is the case for almost every relevant problem. +Hence, the Fourier-algorithm should almost always be used. It can be activated / deactivated +through the +\fbox{% +\parbox{12cm}{ +{\texttt{export FFTALGO=[x]}}}} + +environment variable, which defaults to 1. With the Fourier-algorithm +used normally, the number of nested loops is reduced to three, which makes the loop over the +particle-images the innermost loop. + +\end{itemize} + +{\it Note:} For parallelization over the CPU cores of one node, there are in principle two options: +\begin{itemize} +\item[--] One can simply use OpenMP to parallelize over the particle images (optionally using the +\texttt{BIOEM\_PROJECTIONS\_AT\_ONCE=[x]}). +\item[--] One can use MPI with as many MPI processes as there are CPU cores, and with \texttt{OMP\_NUM\_THREAD=1}. +In this case, the parallelization is done over the orientations only. +\end{itemize} +In general, the first option is better for many maps while the second option is faster for +few maps. + +\vspace{0.5cm} +Naturally, different methods of parallelization can be combined: +\begin{itemize} +\item[--] The Fourier-algorithm can be combined (and is by default combined) with any chosen option +automatically. +\item[--] One can combine MPI with the GPU algorithm to use multiple GPUs at once (as described in +the GPU section). +\item[--] As stated in the previous paragraph, on one node OpenMP works better for many maps while MPI +is better with few maps. For a medium number of particle-images, both methods can be combined. For +instance, \texttt{OMP\_NUM\_THREADS=[x]} can be set to \texttt{x = 1/4 N}, where \texttt{N} is the number of CPU cores on the system, +and BioEM can be called with 'mpirun', and 4 MPI processes. In that case always four +orientations, are processed in parallel, and \texttt{x} particle-images are processed in parallel. +\item[--] On multiple nodes, one can either use MPI exclusively, using as many MPI processes as there are +CPU cores in total. Alternatively, one can use one MPI process per node and do the intra-node +parallelization with OpenMP. Naturally this can also be combined with the previous option, to +mix MPI and OpenMP inside a single node. +\item[--] One can use GPUs and CPU cores jointly to calculate the probabilities for all particle-images. For more +than one GPU, MPI must be employed. In this case, the number of MPI processes must match the +number of GPUs. So the above method to combine MPI, and OpenMP inside one node must be employed, +in order to use all CPU cores. +\end{itemize} + +\section{List of environment variables} + +\begin{itemize} + +\item \texttt{FFTALGO=[x]} (Default: 1) +Set to 1 to enable the Fourier-algorithm (default) or to 0 to use the non-Fourier-algorithm. + +\item \texttt{GPU=[x]} (Default: 0) +Set to 1 to enable GPU usage, set to 0 to use only the CPU. +GPU will be used to calculate the probabilities for all particle-images. The preparation of projections +and PSF convolutions will be processed by the CPU. This is arranged in a pipeline to ensure +continuous GPU utilization. + +\item \texttt{GPUALGO=[x]} (Default: 2) +This option is only relevant if \texttt{GPU=1} and \texttt{FFTALGO=0}. +Hence, it is commonly not used, since \texttt{FFTALGO} defaults to 1. +For the non-Fourier-algorithm there are three \texttt{GPUALGO} implementations: +\begin{itemize} +\item \texttt{x=2}: This will parallelize over the particle-images, and over the center displacements. The approach requires less + memory bandwidth than \texttt{GPUALGO=0} or \texttt{GPUALGO=1}. However, it poses several constraints + on the problem configuration: {\it i)} The number of center displacements per dimension must be a power of 2, and + {\it ii)} must be a factor of the number of CUDA threads per block. +\item \texttt{x=1}: This will parallelize over the particle-images, and then loop over the center displacements on the GPU. It is + usually slower than \texttt{GPUALGO=2} but there are no constraints on the problem + configuration. The particle-images are not processed all at once but in chunks. +\item \texttt{x=0}: As \texttt{GPUALGO=1}, all particle images are processed at once. It is always slower than \texttt{GPUALGO=1} + and should not be used anymore. +\end{itemize} + +\item \texttt{GPUDEVICE=[x]} (Default: fastest) +Only relevant if \texttt{GPU=1}. +\begin{itemize} +\item If this is not set, BioEM will autodetect the fastest GPU. +\item If \texttt{x >= 0}, BioEM will use GPU number \texttt{x}. +\item If \texttt{x = -1}, BioEM runs with \texttt{N} MPI threads, and the system has \texttt{G} GPUs, then BioEM will use GPU +with number (\texttt{N \% G}). The idea is that one can place multiple MPI processes on one node, and +each will use a different GPU. For a multi-node configuration, one must make sure that +consecutive MPI ranks are placed on the same node, {\it i.e.} four processes on two nodes (node0 and +node1) must be placed as (node0, node0, node1, node1) and not as (node0, node1, node0, node1), +because in the latter case only 1 GPU per node will be used (by two MPI processes each). +\end{itemize} + +\item \texttt{GPUWORKLOAD=[x]} (Default: 100) +Only relevant if \texttt{GPU=1}. +This defines the fraction of the workload in percent. To be precise: the fraction of the number +of particle-images processed by the GPU. The remaining maps will be processed by the CPU. Preparation +of projection and convolution will be processed by the CPU in any case. + +\item \texttt{GPUASYNC=[x]} (Default: 1) +Only relevant if \texttt{GPU=1}. +This uses a pipeline to overlap the processing on the GPU, the preparation of projections and +convolutions on the CPU, and the DMA transfer. There is no reason to disable this except for +debugging purposes. + +\item \texttt{GPUDUALSTREAM=[x]} (Default: 1) +Only relevant if \texttt{GPU=1}. +If this is set to 1, the GPU will use two streams in parallel. This can help to improve the GPU +utilization. Benchmarks have shown that there is a very little positive effect by this setting. + +\item \texttt{OMP\_NUM\_THREADS=[x]} (Default: Number of CPU cores) +This is the standard OpenMP environment variable to define the number of OpenMP threads. +It can be used for profiling purposes to analyze the scaling. +If one choses to use MPI for intra-node parallelization (see last two paragraphs of section \ref{wayparallel}). +this can be set to \texttt{x=1} (to use MPI exclusively) or to other values for a mixed MPI / OpenMP +configuration. + +\item \texttt{BIOEM\_PROJECTIONS\_AT\_ONCE=[x]} (Default: 1) +This defines the number of projections prepared at once. Benchmarks have shown that the effect is +negligible. OpenMP is used to prepare these projections in parallel. It is mostly relevant, if +OpenMP is used, no GPU is used, and/or the number of reference maps is very small. + +\item \texttt{BIOEM\_DEBUG\_BREAK=[x]} (Default: deactivated) +This is a debugging option. It will reduce the number of projection and PSF convolutions to a maximum +of \texttt{x} both. It can be used for profiling to analyze scaling, and for fast sanity tests. + +\item \texttt{BIOEM\_DEBUG\_NMAPS=[x]} (Default: deactivated) +As \texttt{BIOEM\_DEBUG\_BREAK}, with the difference that this limits the number of reference particle-images to a +maximum of \texttt{x}. + +\item \texttt{BIOEM\_DEBUG\_OUTPUT=[x]} (Default: 2) +Change the verbosity of the output. Higher means more output, lower means less output. +\begin{itemize} +\item 0: Stands for no debug output. +\item 1: Limited timing output. +\item 2: Standard timing output showing durations of projection, convolution, and comparison (creation +of probabilities. +\end{itemize} +Values above 1 add successively more extensive output. + +\end{itemize} + +\section{Suggestions for runtime configurations} + +{\it Default Settings:} It is recommended that the following settings should be left at theirs defaults: +\texttt{FFTALGO} (Default 1), \texttt{GPUALGO} (Default 2), \texttt{GPUASYNC} (Default 1), \texttt{GPUDUALSTREAM} (Default 1). + +{\it Profiling:} one should limit the number of orientations and projections. +Good settings are \texttt{BIOEM\_DEBUG\_BREAK=4-64} +\texttt{BIOEM\_DEBUG\_OUTPUT=2} is a good choice to get the timing. For a few number of particle-images it might +make sense to switch to \texttt{BIOEM\_DEBUG\_OUTPUT=1}. + +{\it Production environment:} + +\begin{itemize} +\item \texttt{BIOEM\_DEBUG\_OUTPUT=0} can reduce the size of the text output. + +\item \texttt{BIOEM\_PROJECTIONS\_AT\_ONCE=[x]} has usually a tiny or none, but never a negative effect. +The memory footprint increases with \texttt{x}, so it should not be set unnecessarily large. +For best performance, choose a multiple of the number of OpenMP threads. + +\item \texttt{GPU=1} should usually be used if a GPU is available. +Performance wise, one Titan GPU corresponds roughly to 20 cores at 3 GHz. +If the CPU has significant compute capabilities, one should use \texttt{GPUWORKLOAD=[x]} for combined +CPU / GPU processing. + +\item If one uses combined CPU / GPU processing, a good value for \texttt{GPUWORKLOAD=[x]} must be determined. +A good starting point is to measure CPU and GPU individually with + +\texttt{BIOEM\_DEBUG\_OUTPUT=2; BIOEM\_DEBUG\_BREAK=4; GPU=0/1} + +and compare the time of CPU to that of the GPU. Assuming the CPU takes \texttt{c} seconds for the comparison, and the +GPU takes \texttt{g} second, a good starting point for \texttt{GPUWORKLOAD=[x]} is \texttt{x = 100 * g / (c + g)}. + +\item On a single node, one should use OpenMP parallelization for many maps (\texttt{> 1000}) and MPI +parallelization for few maps (\texttt{< 100}). Assume a system with \texttt{N} CPU cores, the command +for many particle-images would be + +\texttt{BIOEM\_PROJECTIONS\_AT\_ONCE=[4*N] OMP\_NUM\_THREADS=[N]} + +and for few + +\texttt{OMP\_NUM\_THREADS=1 ; mpirun -n [N] } + +\item For a medium number of maps, a combined MPI / OpenMP configuration can be better. + + +{\it Example:} +Assume 20 CPU cores, possible options would be (among others): + +\begin{itemize} +\item[--] 20 MPI processes with 1 OMP thread each: + +\texttt{OMP\_NUM\_THREADS=1 mpirun -n 20} +\item[--] 10 MPI processes with 2 OMP threads each: + +\texttt{OMP\_NUM\_THREADS=2 mpirun -n 10} +\item[--] 5 MPI processes with 4 OMP threads each: + +\texttt{OMP\_NUM\_THREADS=5 mpirun -n 4 } +\item[--] 2 MPI processes with 10 OMP threads each: + +\texttt{OMP\_NUM\_THREADS=10 mpirun -n 2} +\end{itemize} +In any case, you should make sure that the number of MPI processes times the number of OMP threads +per process equals the number of (virtual) CPU cores. +In general, the more MPI processes, the better for few particle-images, the more OMP threads, the better for +many maps. +\end{itemize} + +If a system offers multiple GPUs, all GPUs should be used. This must be accomplished via MPI. +In this case, the number of MPI processes per node must match the number of GPUs per node. +There are different ways to make sure every MPI process uses a different GPU (as discussed in +the GPU paragraph of section \ref{wayparallel}). Assuming the MPI processes are placed such, that consecutive MPI +ranks are placed on one node, one can use the \texttt{GPUDEVICE=-1} setting. This is assumed here. +Let us assume an example of \texttt{N} nodes with \texttt{C} CPU cores each and \texttt{G} GPUs each. +The following command will use all GPUs, and ignore the CPUs: + +\texttt{OMP\_NUM\_THREADS=1 GPU=1 GPUDEVICE=-1 mpirun -n [N*G]} + +With the following command, you can use all the CPU cores as well. A combined MPI / OpenMP setting +as discussed in the previous paragraph must be used, under the constraint that the number of MPI +processes matches the number of GPUs: + +\texttt{BIOEM\_PROJECTIONS\_AT\_ONCE=[4*C/G] OMP\_NUM\_THREADS=[C/G] GPU=1} + +\texttt{GPUDEVICE=-1 GPUWORKLOAD=[x] mpirun -n [N*G]} + +Here, \texttt{GPUWORKLOAD} must be tuned as described before. + + +\chapter{Tutorial} +\label{tutorial} + +\section{Posterior probability using BioEM} +\label{commandline} +We now show examples of the different input formats, features and execution commands +for calculating the posterior probability of model given a particle-image set, using BioEM. +All files are provided in the Tutorial directory that comes with the BioEM package. +Here, we will concentrate only on the commandline and physical parameter setup. +Examples of the environmental variable setup for optimizing the performance +are shown in chapter \ref{perfparm}. + +\begin{itemize} + +\item {\it Text Model - Text Image:} +For calculating the probability of a model, in text format, given particle images also in text form. The example parameter inputfile +has few grid points in the Euler angle representation of the orientation. +This should be changed for a larger number in a production run ($\sim 10000$ total orientations). + +{\bf Files:} +\begin{itemize} +\item[] Model file: {\it Model\_Text} +\item[] Parameter input file: {\it Param\_Input} +\item[] Particle-image file: {\it Text\_Image\_Form} +\end{itemize} + +{\bf Commandline execution:} + +\fbox{% +\parbox{12cm}{ +{\texttt{./bioEM --Inputfile Param\_Input --Modelfile Model\_Text --Particlesfile Text\_Image\_Form}}}} + +{\bf Outputfile:} {\it Output\_Probabilities}. + +{\it Note:} Check coordinates in the output \texttt{COORDREAD} file to verify that the model is correct. + +\item {\it PDB Model - Text Image:} +To perform the BioEM calculation with a model in {\it pdb} format. + +{\bf New Command:}\texttt{ --ReadPDB} + +{\bf Files:} +\begin{itemize} +\item[] Model file: {\it Model.pdb} +\item[] Parameter file: {\it Param\_Input} +\item[] Particle-image file: {\it Text\_Image\_Form} +\end{itemize} + +{\bf Commandline execution:} + +\fbox{% +\parbox{12cm}{ +{\texttt{./bioEM --Inputfile Param\_Input --Modelfile Model.pdb --ReadPDB --Particlesfile Text\_Image\_Form}}}} + +{\bf Outputfile:} {\it Output\_Probabilities}. + +\item {\it PDB Model - One MRC Image:} +To perform the BioEM calculation, when reading one single {\it mrc} particle image file. + +{\bf New Command:}\texttt{ --ReadMRC} + +{\bf Files:} +\begin{itemize} +\item[] Model file: {\it Model.pdb} +\item[] Parameter file: {\it Param\_Input} +\item[] Particle-image file: {\it OneImage.mrc} +\end{itemize} + +{\bf Commandline execution:} + +\fbox{% +\parbox{12cm}{ +{\texttt{./bioEM --Inputfile Param\_Input --Modelfile Model.pdb --ReadPDB --Particlesfile OneImage.mrc --ReadMRC}}}} + +{\bf Outputfile:} {\it Output\_Probabilities}. +\item {\it PDB Model - Multiple MRCs:} To perform the BioEM calculation, when multiple {\it mrc} files are read. +In this case, the file name containing the list of all {\it mrc} images should be provided. + +{\bf New Command:}\texttt{ --ReadMultipleMRC} + +{\bf Files:} +\begin{itemize} +\item[] Model file: {\it Model.pdb} +\item[] Parameter file: {\it Param\_Input} +\item[] FILE with names of MRC files : {\it ListMRC} +\end{itemize} + +{\it Note:} The file {\it ListMRC} contains the names of files {\it OneImage.mrc} and +{\it TwoImages.mrc} that are provided in the Tutorial directory. + + +{\bf Commandline execution:} + +\fbox{% +\parbox{12cm}{ +{\texttt{./bioEM --Inputfile Param\_Input --Modelfile Model.pdb --ReadPDB --Particlesfile ListMRC --ReadMRC --ReadMultipleMRC --OutputFile Output\_Probabilities\_MultipleMRC}}}} + +{\bf Example outputfile:} {\it Output\_Probabilities\_MultipleMRC} + +{\it Note:} Both commands \texttt{--ReadMRC --ReadMultipleMRC} are required. + +\end{itemize} + +\subsection{Additional commandline options} + +Several additional features using the command line input are available with BioEM: + + +\begin{itemize} + +\item {\it Dump particle-images:} +This extra feature writes out the particle-images in binary ({\bf PC:?}) format. +This makes them faster to read in a further BioEM execution. + +{\bf New Command:}\texttt{ --DumpMaps} + +{\bf Files:} +\begin{itemize} +\item[] Model file: {\it Model.pdb} +\item[] Parameter file: {\it Param\_Input} +\item[] FILE with names of MRC files : {\it ListMRC} +\end{itemize} + +{\bf Commandline execution:} + +\fbox{% +\parbox{12cm}{ +{\texttt{./bioEM --Inputfile Param\_Input --Modelfile Model.pdb --ReadPDB --Particlesfile ListMRC --ReadMRC --ReadMultipleMRC --DumpMaps}}}} + +{\bf Outputfiles: }{\it Output\_Probabilities} and {\it maps.dump} (that will be useful in a further execution). + +\item {\it Load particle-images:} +This feature reads in the particle images in binary format from file {\it maps.dump} +that was described previously. In this case, no particle file name is necessary, +but the {\it maps.dump} file should be in the current directory. + +{\bf New Command:}\texttt{ --LoadMapDump} + +{\bf Files:} +\begin{itemize} +\item[] Model file: {\it Model.pdb} +\item[] Parameter file: {\it Param\_Input} +\item[] Dumped Mapfile: {\it maps.dump} +\end{itemize} + +{\bf Commandline execution:} + +\fbox{% +\parbox{12cm}{ +{\texttt{./bioEM --Inputfile Param\_Input --Modelfile Model.pdb --ReadPDB --LoadMapDump }}}} + +{\bf Outputfile:} {\it Output\_Probabilities} + +\item {\it Read orientations from file:} +Related to section \ref{intor}. +With this feature the orientations are read from an inputfile. +The first line of the file should contain the total number of orientations. +The format should of the file be "\%12.6f\%12.6f\%12.6f" for $\alpha$, $\beta$ and $\gamma$ +Euler angles, respectively, and "\%12.6f\%12.6f\%12.6f\%12.6f" for the quaternions. + +{\bf New Command:}\texttt{ --ReadOrientation arg} + +If quaternions are used, additionally, in the input parameter file one has to add +the keyword + +\texttt{QUATERNIONS} + +An example, for reading an Euler angle list file is + + +{\bf Files:} +\begin{itemize} +\item[] Model file: {\it Model.pdb} +\item[] Parameter file: {\it Param\_Input} +\item[] Particle image file: {\it Text\_Image\_Form} +\item[] EulerAngle File: {\it Euler\_Angle\_List} +\end{itemize} +{\bf Commandline execution:} + +\fbox{% +\parbox{12cm}{ +{\texttt{./bioEM --Inputfile Param\_Input\_ReadEulerAng --Modelfile Model.pdb --ReadPDB --Particlesfile Text\_Image\_Form --ReadOrientation Euler\_Angle\_List}}}} + +{\bf Outputfile:} {\it Output\_Probabilities} + +In the directory Quaternions, there are several quaternion list input files that +sample uniformly the rotational group in 3D space. These are files +recommended to use. If using quaternions, remember to include the keyword +\texttt{QUATERNIONS} in the input parameter file, for an example see {\it Param\_Input\_Quaternions} in +the tutorial directory. + +\item {\it Print map from Model:} This option is {\bf completely independent } of the BioEM +posterior probability calculation. +It simply prints out a synthetic image +from a model given a set of parameters. This is useful for post-analysis +characterization, for example in comparing the best calculated map +to the experimental image (as in Fig.1 Ref.[1]). Or for creating synthetic +image sets with artificial noise. + + +{\bf New Command:}\texttt{--PrintBestCalMap arg} + +{\bf Files:} +\begin{itemize} +\item[] Model file: {\it Model.pdb} +\item[] Print Map Parameter file: {\it Param\_Print\_MAP} +\end{itemize} + +{\it Important:} Print Map input parameters are different from those of standard BioEM parameter input. +For this option, (for the moment) one can only use the Euler angles. +An example file, {\it Param\_Print\_MAP}, is in the Tutorial directory. +Specifically, the keyword \texttt{BEST\_[X]} indicates the specific value of parameter \texttt{[X]} used for the printed map. +To add random Gaussian noise, use the keyword \texttt{WITHNOISE (float)} where \texttt{(float)} is +the signal-to-noise ratio. If this keyword is omitted no noise is added. + +{\bf Commandline execution:} + +\fbox{% +\parbox{12cm}{ +{\texttt{./bioEM --Modelfile Model.pdb --ReadPDB --PrintBestCalMap Param\_Print\_MAP}}}} + +{\bf Outputfile:} \texttt{BESTMAP} with format ( +\texttt{"MAP" --- MAP number --- Pixel X --- Pixel Y --- Intensity}). +%After map row there is an extra line for easily plotting with pm3d map in Gnuplot +%(see http://gnuplot.info/). + + + +{\it Note:} BioEM Output probability files are not generated for this case. + +\end{itemize} + + + +\section{Optional calculations} + +Some extra calculations are available in with the BioEM program. Additionally, +from calculating the BioEM posterior probability, the options +assess further information of the Bayesian analysis. +The keywords to enable these calculations should be included in the input parameter file. + +%, allow assessing further information +%of the analysis. + + +\begin{itemize} +\item {\it Posterior as a function of orientations:} This option prints out the +posterior probabilities of the model as a function of the orientation. +In this case, no integration over the orientations is performed, allowing the assessment +of the probability distribution for each orientation. +The keyword in the parameter file needed is + +\fbox{% +\parbox{12cm}{ +{\texttt{WRITE\_PROB\_ANGLES}}}} + +an additional outputfile \texttt{ANG\_PROB} is generated. An example of the parameter +input is provide in the {\it Param\_Input\_WritePAng} file. + +\item {\it Including prior probabilities}: + +To include the prior probabilities both for the model see the {\it Param\_Input\_Priors} file. The prior probabilities for the orientations +should be included in an additional file (see {\it Euler\_Angle\_List\_Prior}). An example of this execution is: + + +\fbox{% +\parbox{12cm}{ +{\texttt{./bioEM --Modelfile Model\_Text --Particlesfile Text\_Image\_Form --Inputfile Param\_Input\_Priors --ReadOrientation Euler\_Angle\_List\_Prior}}}} + + +\item{\it Cross Correlation Calculation:} This option prints the optimal cross correlation +of the model at as a function of the pixels in the micrograph (see section \ref{crosscor}). +This can be useful in the preliminary steps of particle picking (identification). +The keywords in the parameter file are + +\fbox{% +\parbox{12cm}{ +{\texttt{WRITE\_CROSSCOR \\ +FLIPPED (optional)\\ +CROSSCOR\_NOTBAYESIAN (optional)}}}} + +An additional file, \texttt{CROSS\_CORRELATION}, shows the cross-correlation as a function of the pixels. +By default the cross-correlation is written every 10 pixels. This feature is still in exploratory phase. + + + +\end{itemize} + + + +{\it Note:} These options can be included in with the standard BioEM posterior probability setup +executions. However, the print map option does not work with these keywords. + +\section{BioEM output} +The main output file is called \texttt{Output\_Probabilities} by default. +To change the name of the output file include the following keyword +in the commandline +\fbox{% +\parbox{12cm}{ +{\texttt{--OutputFile arg }}}} + +where \texttt{arg} is the desired name of the output file. +This file contains the logarithm of the posterior probability of +the model to each individual experimental image. It also +includes the parameters that maximize the posterior. + +\fbox{% +\parbox{12cm}{ +{\texttt{RefMap [number Particle Map] LogProb [$\ln (P)$]\\ + RefMap [number Particle Map] Maximizing Param: [Orientation] [CTF or PSF parameters] [center displacement]}}}} + +{\it Important remark:} The posterior probability is not normalized. Thus, it is always recommended to compare +$\ln (P)$ of different models or relative to noise as in \cite{CossioHummer_2013}. See section \ref{modcom}. + +Before executing a production run, it is recommended to check that the values of the log-posterior are finite, +and that the maximizing parameters are in a reasonable range ({\it e.g.} not at the borders of the integration limits). + +The output file \texttt{COORDREAD} is always generated. It is to check that the model coordinates, radius and density +are read correctly. + +\subsection{Optional outputs} + +The optional output files for BioEM are: + +\begin{itemize} + +\item[--] \texttt{ANG\_PROB}: This file has the posterior probabilities for each orientation. +The keyword \texttt{WRITE\_PROB\_ANGLES} should be included in the parameter inputfile. +For the Euler angles, the output file format is + +\fbox{% +\parbox{12cm}{ +{\texttt{[Map number -- alpha -- beta -- gamma -- log Probability]} }}} + +For the quaternions, its format is + +\fbox{% +\parbox{12cm}{ +{\texttt{[Map number -- q1 -- q2 -- q3 -- log Probability]}}}} + +where the fourth quaternion is obtained by $q_4=\sqrt{1-q_1^2-q_2^2-q_3^2}$. + +\item[--] \texttt{CROSS\_CORRELATION}: This file contains the cross-correlation as a function of the pixels. +The keyword \texttt{WRITE\_CROSSCOR} should be included in the parameter inputfile. +The output file format is + +\fbox{% +\parbox{12cm}{ +{ \texttt{[Map number -- Pixel x -- Pixel y -- Cross-Correlation] }}}} + +\item[--] \texttt{BESTMAP}: This file contains a model projection image from a specific set of parameters (this excludes a BioEM posterior probability calculation). +The commandline keyword is \texttt{--PrintBestCalMap}. +The output file format is + +\fbox{% +\parbox{12cm}{ +{\texttt{[Map number -- Pixel X -- Pixel Y -- Intensity]}}}} +\end{itemize} + +\section{Model comparison using BioEM} +\label{modcom} + +BioEM is standardly used for model comparison and ranking. +Here, we provide a complete example of how to analyse the +output files of BioEM for discriminating between structural models. +The relevant files are found in the {\it MODEL\_COMPARISON} directory +that is inside the {\it Tutorial} directory. There you will find:\\ + +\fbox{% +\parbox{12cm}{ +{ - {\it MODEL\_1}: First model in text format.\\ +- {\it MODEL\_2}: Second model in text format.\\ +- {\it Param\_Input\_ModelComparision}: example of parameter input.\\ +- {\it Quaternion\_List}: Here, we are using quaternions instead of Euler angles, and they are read +them from this additional file.\\ +- {\it 20\_ParticleImages}: Image file in text format. \\ +- {\it subtract\_LogP.sh}: Bash script to extract the difference in log posterior. +}}} + +For the model comparison, one need to run the BioEM program for each model. +We note that for these models, the radius of the atoms +is less than the pixel size. If this is the case for your system, +you should include the keyword \texttt{NO\_PROJECT\_RADIUS} in your input parameter file +(in this regard, if the radius is less than the pixel size, this option is mandatory). + +\begin{itemize} +\item[] For Model 1: + +\fbox{% +\parbox{12cm}{ +{\texttt{BIOEM\_DEBUG\_OUTPUT=0 ./bioEM --Modelfile MODEL\_1 --Particlesfile 20\_ParticleImages --Inputfile Param\_Input\_ModelComparision +--ReadOrientation Quaternion\_List --OutputFile Output\_MODEL\_1}}}} + +\item[] For Model 2: + +\fbox{% +\parbox{12cm}{ +{\texttt{BIOEM\_DEBUG\_OUTPUT=0 ./bioEM --Modelfile MODEL\_2 --Particlesfile 20\_ParticleImages --Inputfile Param\_Input\_ModelComparision +--ReadOrientation Quaternion\_List --OutputFile Output\_MODEL\_2}}}} +\end{itemize} + +After executing these two commands, the two output files containing the posterior probabilities of each model ( +\texttt{Output\_MODEL\_1} and \texttt{Output\_MODEL\_2}) are generated. +Since the input parameter and particle-image files are the same, then the output files +should only differ in the specific numeric results. + +To extract the difference in log-posterior of model 1 with respect to model 2, one +can simply run in terminal the bash script {\it subtract\_LogP.sh}:\\ +%\vspace{0.5cm} + +\fbox{% +\parbox{12cm}{ +{\texttt{./subtract\_LogP.sh Output\_MODEL\_1 Output\_MODEL\_2}}}} +\vspace{0.5cm} + +This prints out the particle number, log-posterior of model 1, log-posterior of model 2, +difference in log-posteriors (Model 1- Model 2), and cumulative difference of the log-posterior. +The first line of the output on the screen should be: \\ + + +\fbox{% +\parbox{12cm}{ +{\texttt{{\footnotesize +MapNum: 0 Mod1: -70902.70000 Mod2: -70900.40000 Dif: -2.30000 Cum: -2.30000 + }}}}} + +\vspace{0.5cm} + +One can modify, and treat the output of the posterior probabilities in many different manners, as done in +\cite{CossioHummer_2013}. The data analysis, and interpretation is left to +the criteria of each individaul user. + + + +\bibliographystyle{plain} +\bibliography{Biblio_BioEM.bib} + + + + +\end{document} diff --git a/Tutorial_BioEM/MODEL_COMPARISION/Param_Input_ModelComparision b/Tutorial_BioEM/MODEL_COMPARISION/Param_Input_ModelComparision index 82f4afd673a5d3f7d97ff75c0c72421a889b62d3..4b1f4fe376b5ef93f32118c3bdd203b90af39e24 100644 --- a/Tutorial_BioEM/MODEL_COMPARISION/Param_Input_ModelComparision +++ b/Tutorial_BioEM/MODEL_COMPARISION/Param_Input_ModelComparision @@ -1,15 +1,17 @@ +##### Micrograph Parameters ####### NUMBER_PIXELS 224 PIXEL_SIZE 1.32 -GRIDPOINTS_ENVELOPE 3 -START_ENVELOPE 0.000 -END_ENVELOPE 0.001 -GRIDPOINTS_PSF_PHASE 3 -START_PSF_PHASE 0.004 -END_PSF_PHASE 0.024 -GRIDPOINTS_PSF_AMP 1 -START_PSF_AMP 1. -END_PSF_AMP 1. -MAX_D_CENTER 20 -PIXEL_GRID_CENTER 2 + +###### Sincs its a full atom structure ##### NO_PROJECT_RADIUS -QUATERNIONS + +##### Quaterion grid points: ####### +USE_QUATERNIONS + +##### Constrast transfer integration: ####### +CTF_B_FACTOR 100.0 300.5 8 +CTF_DEFOCUS 1.0 6.0 5 +CTF_AMPLITUDE 0.1 0.1 1 + +##### Center displacement: ####### +DISPLACE_CENTER 10 2 diff --git a/Tutorial_BioEM/Param_Print_MAP b/Tutorial_BioEM/Param_Print_MAP index 748438c154b9eb94296f0f2854c57eaa2c1e7bd2..537607e4db6c7f3aa2ea463304fbd4904754be5f 100644 --- a/Tutorial_BioEM/Param_Print_MAP +++ b/Tutorial_BioEM/Param_Print_MAP @@ -1,14 +1,34 @@ +##### Micrograph Parameters ####### + NUMBER_PIXELS 220 PIXEL_SIZE 1.77 -BEST_AMP 1 + +######## Center displacement ################ BEST_DX -10 BEST_DY 5 + +######## Normalization and Offset ################ BEST_NORM 1 BEST_OFFSET 0 -BEST_PHASE 0.00068 -BEST_ENVELOPE 0.03 -WITHNOISE 5.476 + +######## Euler Angles ################ BEST_ALPHA -0.464793 BEST_GAMMA -0.409391 BEST_BETA 1.97283 -#PRINT_ROTATED_MODELS + +######## Using PSF ################ +#USE_PSF +BEST_PSF_ENVELOPE 0.0002 +BEST_PSF_PHASE 0.004 +BEST_PSF_AMP 1 + +##### Using CTF ####### +BEST_CTF_B_FACTOR 200.5 +BEST_CTF_DEFOCUS 3.5 +BEST_CTF_AMP 0.3 + +######## ADDING NOISE ########### +WITHNOISE 5.476 + +########## Printing the 3d Model rotation ########### +PRINT_ROTATED_MODELS