diff --git a/.github/workflows/deploy-docs.yml b/.github/workflows/deploy-docs.yml new file mode 100644 index 00000000..499d4806 --- /dev/null +++ b/.github/workflows/deploy-docs.yml @@ -0,0 +1,37 @@ +name: docs + +on: + push: + branches: + - main + +# This job installs dependencies, builds the docs, and pushes it to +# the `gh-pages` branch of the same repository. +jobs: + deploy-book: + if: github.repository == ‘SBC-Utrecht/pytom-match-pick’ + runs-on: ubuntu-latest + steps: + - name: Checkout repository + uses: actions/checkout@v4 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: 3.x + + - name: Install dependencies + run: | + pip install mkdocs mkdocs-material markdown-exec + + # Build the book + - name: Build the book + run: | + mkdocs build + + # Push the site to github-pages + - name: GitHub Pages action + uses: peaceiris/actions-gh-pages@v4 + with: + publish_dir: ./site + github_token: ${{ secrets.GITHUB_TOKEN }} diff --git a/README.md b/README.md index 3b213711..51ac9761 100644 --- a/README.md +++ b/README.md @@ -4,9 +4,16 @@ # pytom-match-pick: GPU template matching for cryo-ET -GPU template matching, originally developed in [PyTom](https://github.com/SBC-Utrecht/PyTom), as a standalone python package that can be run from the command line. +GPU template matching, originally developed in [PyTom](https://github. +com/SBC-Utrecht/PyTom), as a standalone python package that is run from the command +line. -![cover_image](images/tomo200528_100_illustration.png) +![cover_image](docs/images/tomo200528_100_illustration.png) + + +#--8<-- [start:docs] ## Requires @@ -54,7 +61,7 @@ The installation above also adds the optional dependencies `[matplotlib, seaborn python -m pip install pytom-match-pick ``` -## Cupy warning +### Cupy warning Having issues running the software? If cupy is not correctly installed, ```commandline python -c "import pytom_tm" @@ -72,10 +79,6 @@ specific build compatible with the installed cuda toolkit. ## Usage -Detailed usage instructions are available on the wiki: https://github.com/SBC-Utrecht/pytom-match-pick/wiki - -Also, a tutorial can be found on the same wiki: https://github.com/SBC-Utrecht/pytom-match-pick/wiki/Tutorial - The following scripts are available to run with `--help` to see parameters: - create a template from an mrc file containing a density map: `pytom_create_template.py --help` @@ -85,6 +88,8 @@ The following scripts are available to run with `--help` to see parameters: - estimate an ROC curve from a job file (.json): `pytom_estimate_roc.py --help` - merge multiple star files to a single starfile: `pytom_merge_stars.py --help` +Detailed usage instructions and a tutorial are available on [our site](https://SBC-Utrecht.github.io/pytom-match-pick). + ## Usage questions, ideas and solutions, engagement, etc Please use our [github discussions](https://github.com/SBC-Utrecht/pytom-match-pick/discussions) for: - Asking questions about bottlenecks. @@ -113,7 +118,13 @@ you make PRs: pre-commit install ``` -This uses Ruff to check and format whenever you make commits. +This uses Ruff to check and format whenever you make commits. + +If you update anything in the (documentation) `docs/` folder make sure to test build the website locally: + +```commandline +mkdocs serve +``` ## Tests @@ -148,3 +159,8 @@ For a reference on GPU accelerated template matching in tomograms please see the DOI = {10.3390/ijms241713375} } ``` + + +#--8<-- [end:docs] diff --git a/docs/Developers.md b/docs/Developers.md new file mode 100644 index 00000000..0d890519 --- /dev/null +++ b/docs/Developers.md @@ -0,0 +1,47 @@ +Star files written out by pytom-match-pick should be easily integratable with other software as it follows RELION star file conventions. The only difference are three column headers with extraction statistics, which are important to maintain for annotations. + +The exact header is: + +``` +# Created by the starfile Python package (version x.x.x) at xx:xx:xx on xx/xx/xxxx + + +data_ + +loop_ +_rlnCoordinateX #1 +_rlnCoordinateY #2 +_rlnCoordinateZ #3 +_rlnAngleRot #4 +_rlnAngleTilt #5 +_rlnAnglePsi #6 +_rlnLCCmax #7 +_rlnCutOff #8 +_rlnSearchStd #9 +_rlnDetectorPixelSize #10 +_rlnMicrographName #11 +``` + +With RELION5 compatibility mode in `pytom_extract_candidates.py` the star file +columns are this (correcting the position in to Angstrom and relative to the +tomogram center): + +``` +# Created by the starfile Python package (version x.x.x) at xx:xx:xx on xx/xx/xxxx + + +data_ + +loop_ +_rlnCenteredCoordinateXAngst #1 +_rlnCenteredCoordinateYAngst #2 +_rlnCenteredCoordinateZAngst #3 +_rlnAngleRot #4 +_rlnAngleTilt #5 +_rlnAnglePsi #6 +_rlnLCCmax #7 +_rlnCutOff #8 +_rlnSearchStd #9 +_rlnDetectorPixelSize #10 +_rlnTomoName #11 +``` \ No newline at end of file diff --git a/docs/Frequently-asked-questions-(FAQ).md b/docs/Frequently-asked-questions-(FAQ).md new file mode 100644 index 00000000..dd338f83 --- /dev/null +++ b/docs/Frequently-asked-questions-(FAQ).md @@ -0,0 +1,34 @@ +## How to deal with gold beads? + +Gold beads (and other artifacts) can often interfer in template matching due to their high electron scattering potential. The easiest way to deal with them is by removing them on the micrograph level prior to reconstruction. The following IMOD commands can do the trick: + +``` +imodfindbeads -size [GOLD_BEAD_DIAMETER_IN_PIXELS] -input TS_ID.st -output TS_ID.fid -spacing 0.8 +ccderaser --input TS_ID.st --output TS_ID_erased.st -model TS_ID.fid -order 0 -circle / -better [1.5 * GOLD_BEAD_RADIUS_IN_PIXELS] -merge 1 -exclude -skip 1 -expand 3 +``` + +## What template box size should I use? + +For the simple missing wedge model a box size that tightly fits the template and mask is easiest (and slightly faster). + +For the full per-tilt-weighting model its better to have some overhang, a rought estimate could be a box size of `4 * particle_diameter`. This aids in sampling the CTF function and the individual tilts in Fourier space. However, the mask radius should still snuggly fit the template. Although the larger box size slows down rotations slightly, the search benefits more from better sampling of the weighting function. + +## Is my particle handedness correct? + +If template matching is giving unexpectedly poor results for the particle of interest, it might that the template has the wrong handedness. In that case, the `pytom_create_template.py` has the option `--mirror` to produce a mirrored version of the template. We advice to create a mirrored and non-mirrored version of the template and to run template matching with both. After extracting ~1000 particle from both jobs, while setting the `--cut-off` to -1 (in `pytom_extract_candidates.py`). You can plot the results with the following python code: + +``` +import starfile +import matplotlib.pyplot as plt + +raw = starfile.read('[TOMO_ID]_particles.star') +mirror = starfile.read('[TOMO_ID]_mirror_particles.star') + +fig, ax = plt.subplots(figsize=(5,5)) +ax.plot(raw['ptmLCCmax'], label='raw') +ax.plot(mirror['ptmLCCmax'], label='mirror') +ax.set_xlabel('Particle index') +ax.set_ylabel('LCCmax') +ax.legend() +plt.show() +``` \ No newline at end of file diff --git a/docs/Timings.md b/docs/Timings.md new file mode 100644 index 00000000..67c60ec9 --- /dev/null +++ b/docs/Timings.md @@ -0,0 +1,14 @@ + +This page provides a quick overview of software timings. All were run on a tomogram with dimension (462x478x250) and a template box size of 34x34x34. + +``` +Single GPU (GTX 1080 Ti) +12.85 (deg) - 193 sec. (3.22 min.) +``` + +``` +4 GPU's (GTX 1080 Ti) +12.85 (deg) - 48 sec. (0.80 min.) +7.00 (deg) - 306 sec. (5.10 min.) +3.00 (deg) - 3801 sec. (63.35 min.) +``` \ No newline at end of file diff --git a/docs/Usage.md b/docs/Usage.md new file mode 100644 index 00000000..93c5c1d6 --- /dev/null +++ b/docs/Usage.md @@ -0,0 +1,216 @@ + +Features in a tomogram that resemble a structural 'template' can be localized in an automated fashion using 'template matching'. In this approach a 3D template is correlated with a given tomogram. In this procedure the different possible rotations and translations are sampled exhaustively using the algorithm described in [Förster et al. (2010)](http://dx.doi.org/10.1016/S0076-6879(10)83011-3). + +## Requirements + +For usage you need at least a set of reconstructed tomograms in the MRC format and a template structure in the MRC format. Tomograms in IMOD format (.rec) are also okay but need to be renamed (or softlinked!) to have the correct extension (.mrc). Tomograms are ideally binned 6x or 8x to prevent excessive runtimes. The template can be an EM reconstruction (from the EMDB) or a PDB that was coverted to a density (for example via Chimera molmap). + +## Template matching workflow + +Using template matching in this software consists of the following steps: + +1. Creating a template and mask +2. Matching the template in a tomogram +3. Extracting particles +4. Merging annotations for export to other software + + +## 1. Creating a template and mask + +**Important**: +- The template and mask need to have the same box size. +- The template needs to have the same contrast as the tomogram (e.g. the particles + are black in both the tomogram and template). Contrast can be adjusted with the + `--invert` option. + +### pytom_create_template.py + +Using an EM map as a reference structure generally leads to the best results. Alternatively a structure from the PDB can be converted in Chimera(X) using the molmap command to create an MRC file that models the electrostatic potential. + +```python exec="on" result="ansi" +import argparse +code = (""" +---8<--- "entry_points.py:create_template_usage" +""") +cleaned_lines = [line.lstrip() for line in code.splitlines() + if not any(keyword in line for keyword in ('action=', 'default=', 'type=')) +] +cleaned_code = '\n'.join(cleaned_lines) +exec(cleaned_code) +print(parser.format_help()) +``` + +### pytom_create_mask.py + +The mask around the template can be quite tight to remove as much noise as possible around the particles of interest. We recommend around 10%-20% overhang relative to the particle radius. You can also generate an ellipsoidal mask for particles that do not approximate well as a sphere. Though you will probably need to reorient this mask in chimera and resample to the grid of the template. Optionally you could also create a structured mask around the template in external software (via thresholding and dilation for example). Take into account that non-spherical masks roughly double the template matching computation time. + +```python exec="on" result="ansi" +import argparse +code = (""" +---8<--- "entry_points.py:create_mask_usage" +""") +cleaned_lines = [line.lstrip() for line in code.splitlines() + if not any(keyword in line for keyword in ('action=', 'default=', 'type=')) +] +cleaned_code = '\n'.join(cleaned_lines) +exec(cleaned_code) +print(parser.format_help()) +``` + +## 2. Matching the template in a tomogram + +### pytom_match_template.py + +This script requires at least a tomogram, a template, a mask, the min and max tilt angles (for missing wedge constraint), an angular search, and a GPU index to run. The search can be limited along any dimension with the `--search-x`, `--search-y`, and `--search-z` parameters; for example to skip some empty regions in the z-dimension where the ice layer is not yet present, or to remove some reconstruction artifact region along the x-dimension. With the `--volume-split` option, a tomogram can be split into chunks to allow them to fit in GPU memory (useful for large tomograms). Providing multiple GPU's will allow the program to split the angular search (or the subvolume search) over multiple cards to speed up the algorithm. + +The software automatically calculates the angular search based on the available +resolution and provided particle diameter. The required search is found from the +Crowther criterion $\Delta \alpha = \frac{180}{\pi r_{max} d}$. For the maximal +resolution the voxel size is used, unless a low-pass filter is specified as this +limits the available maximal resolution. You can exploit this to reduce the angular +search! For non-spherical particles we suggest choosing the particle diameter as the +longest axis of the macromolecule. + +In case the template matching is run with a non-spherical mask, it is essential to set the `--non-spherical-mask` flag. It requires a slight modification of the calculation that will roughly double the computation time, so only use non-spherical masks if absolutely necessary. + +#### Optimizing results: per tilt weighting with CTFs and dose accumulation + +Optimal results are obtained by also incorporating information for the 3D CTF. You +can pass the following files (and parameters): + +* Tilt angles: a `.rawtlt` or `.tlt` file to the `--tilt-angles` parameter with all the + tilt + angles used to reconstruct the tomogram. You should then also set the + `--per-tilt-weighting` flag. +* CTF data: a `.defocus` file from IMOD or `.txt` file to `--defocus-file`. The `. + txt` file + should specify the defocus of each tilt in $\mu m$. You can also give a + single defocus value (in $\mu m$). The CTF will also require input for + `--voltage`, `--amplitude-contrast`, and `--spherical-abberation`. +* Dose weighting: a `.txt` file to `--dose-accumulation` with the accumulated dose per + tilt (assuming the same ordering as `.tlt`). Each line contains a single float + specifying the accumulated dose in $e^{-}/\text{Å}^{2}$. Dose weighting only works in + combination with `--per-tilt-weighting`. + +_(As a side note, you can also only enable `--per-tilt-weighting` **without** dose accumulation and CTFs, or **with either** dose accumulation or CTFs.)_ + +When enabling the CTF model here (with the defocus file), it is important that the template is not multiplied with a CTF before passing it to this script. The template only needs to be scaled to the correct pixel size and the contrast should be adjusted to match the contrast in the tomograms. + +Secondly, if the tomogram was CTF corrected, for example by using IMODs +strip-based CTF correction or NovaCTF. Its important to add the parameter +`--tomogram-ctf-model phase-flip` which modifies the template CTF to match the +tomograms CTF correction. + +#### Background corrections + +The software contains two background correction methods that might improve results: +`--spectral-whitening` or `--random-phase-correction` (from STOPGAP). In our +experience the random phase correction is most reliable, while spectral whitening +never seemed to clearly improve results. + +```python exec="on" result="ansi" +import argparse +code = (""" +---8<--- "entry_points.py:match_template_usage" +""") +cleaned_lines = [line.lstrip() for line in code.splitlines() + if not any(keyword in line for keyword in ('action=', 'default=', 'type=')) +] +cleaned_code = '\n'.join(cleaned_lines) +exec(cleaned_code) +print(parser.format_help()) +``` + +## 3. Extracting particles + +Both scripts run on the job file created in `pytom_match_template.py` which contains details about correlation statistics and the output files. The job file always has the format `[TOMO_ID]_job.json`. + +**IMPORTANT** For both scripts the `[-r, --radius-px]` option needs to be considered carefully. The particle extraction will mask out spheres with this radius around each peak in the score volume and prevents selecting the same macromolecule twice. It is specified as an integer **number of pixels** (not Angstrom!) and ideally it should be the radius of the particle of interest. It can be found by dividing the particle radius by the pixel size, e.g. a ribosome (r = 290Å / 2) in a 15Å tomogram should gets a pixel radius of 9.6. As it needs to be an integer value and ribosomes are not perfect spheres, it is best to round it down to 9 pixels. + +### pytom_extract_candidates.py + +#### STAR file metadata + +Resulting STAR files from extraction have three colums with extraction statistics (`LCCmax`, `CutOff`, `SearchStd`). Dividing the `LCCmax` and the `CutOff` by the `SearchStd`, will express them as a number of $\sigma$ or (3D SNR; similar to [Rickgauer et al. (2017)](https://doi.org/10.7554/eLife.25648). + +STAR files written out by the template matching module will have RELION compliant column headers, i.e. `rlnCoordinateX` and `rlnAgleRot`, to simplify integration with other software. The Euler angles that are written out therefore also follow the same conventions as RELION and Warp, i.e. `rlnAngleRot`, `rlnAngleTilt`, `rlnAnglePsi` are intrinsic clockwise ZYZ Euler angles. Hence they can be directly used for subtomogram averaging in RELION. See here for more info: [https://www.ccpem.ac.uk/user_help/rotation_conventions.php](https://www.ccpem.ac.uk/user_help/rotation_conventions.php). + +Please see the [For developers](Developers.md) section for more details on the +metadata. + +#### Default true positive estimation + +The particle extraction has been updated to use the formula in [Rickgauer et al. (2017)](https://doi.org/10.7554/eLife.25648) for finding the extraction threshold based on the false alarm rate. This was not yet described in our [IJMS publication](https://www.mdpi.com/1422-0067/24/17/13375) but is essentially very similar to the Gaussian fit that we used. However, it is more reliable and also specific to the standard deviation $\sigma$ of the search in each tomogram. `pytom_match_template.py` keeps track of $\sigma$ and stores it in the job file. The user can specify a number of false positives to allow per tomogram with a minimum value of 1. It can be increased to make the extraction threshold more lenient which might increase the number of true positives at the expense of more false positives. The parameter should roughly correspond to the number of false positives that end up in the extracted particle list. + +Template matching has a huge search space $N_{voxels} * N_{rotations}$ which is mainly +false positives, and has in comparison a tiny fraction of true positives. If we have +a Gaussian for the background (with expected mean 0 and some standard deviation), +the false alarm rate can be calculated for a certain cut-off value, as it is +dependent on the size of the search space. For example, a false alarm rate of $(N_ +{voxels} * N_{rotations})^{-1}$, indicates it would expect 1 false positive in the +whole search. This can be calculated with the error function, + +$$N^{-1} = \text{erfc}( \theta / ( \sigma \sqrt{2} ) ) / (2 n_{\text{FP}})$$ + +, where theta is the cut-off, sigma the standard deviation of the Gaussian, and N the search space. $n_{\text{FP}}$ represents the scaling by the user of tolerated number of false positives. + +#### Tophat transform filter + +This option can be used to filter the score map for sharp peaks (steep local maxima) +which usually correspond to true positives. This will be described in a forthcoming +publication. For now, you can check out Marten's poster at CCPEM that shows some +preliminary results: [10.5281/zenodo.13165643](https://doi.org/10.5281/zenodo.13165643). + +```python exec="on" result="ansi" +import argparse +code = (""" +---8<--- "entry_points.py:extract_candidates_usage" +""") +cleaned_lines = [line.lstrip() for line in code.splitlines() + if not any(keyword in line for keyword in ('action=', 'default=', 'type=')) +] +cleaned_code = '\n'.join(cleaned_lines) +exec(cleaned_code) +print(parser.format_help()) +``` + +### pytom_estimate_roc.py + +This script runs the Gaussian fit as described in the [IJMS publication](https://www.mdpi.com/1422-0067/24/17/13375). It requires installation with plotting dependencies as it writes out or displays a figure showing the Gaussian fit and estimated ROC curve. The benefit is that it estimates some classification statistics (such as false discovery rate and sensitivity). You can use it to esimate an extraction threshold for a representative tomogram and then supply this threshold as the `[-c, --cut-off]` parameter for `pytom_extract_candidates.py`. + +```python exec="on" result="ansi" +import argparse +code = (""" +---8<--- "entry_points.py:estimate_roc_usage" +""") +cleaned_lines = [line.lstrip() for line in code.splitlines() + if not any(keyword in line for keyword in ('action=', 'default=', 'type=')) +] +cleaned_code = '\n'.join(cleaned_lines) +exec(cleaned_code) +print(parser.format_help()) +``` + +## 4. Merging annotations for export to other software + +After running template matching and candidate extraction on multiple tomograms, each tomogram will have an individual starfile with particle annotations. Each starfile will contain the `MicrographName` column which refers back to the tomogram name. Multiple starfiles can therefore be appended to results in a large list which can be used in other software (such as RELION, WarpM) to load annotations. These software will link the annotations to specific tilt-series using the `MicrographName` column. + +### pytom_merge_stars.py + +Without providing any parameters the script will try to merge all the starfiles in the current working directory and save them to a new file `particles.star`. + +```python exec="on" result="ansi" +import argparse +code = (""" +---8<--- "entry_points.py:merge_stars_usage" +""") +cleaned_lines = [line.lstrip() for line in code.splitlines() + if not any(keyword in line for keyword in ('action=', 'default=', 'type=')) +] +cleaned_code = '\n'.join(cleaned_lines) +exec(cleaned_code) +print(parser.format_help()) +``` + + + diff --git a/images/tomo200528_100_illustration.png b/docs/images/tomo200528_100_illustration.png similarity index 100% rename from images/tomo200528_100_illustration.png rename to docs/images/tomo200528_100_illustration.png diff --git a/docs/index.md b/docs/index.md new file mode 100644 index 00000000..7c7396a7 --- /dev/null +++ b/docs/index.md @@ -0,0 +1,9 @@ +# pytom-match-pick: GPU template matching for cryo-ET + +GPU template matching, originally developed in [PyTom](https://github. +com/SBC-Utrecht/PyTom), as a standalone python package that is run from the command +line. + +![cover_image](images/tomo200528_100_illustration.png) + +--8<-- "README.md:docs" diff --git a/docs/javascripts/katex.js b/docs/javascripts/katex.js new file mode 100644 index 00000000..8c642137 --- /dev/null +++ b/docs/javascripts/katex.js @@ -0,0 +1,10 @@ +document$.subscribe(({ body }) => { + renderMathInElement(body, { + delimiters: [ + { left: "$$", right: "$$", display: true }, + { left: "$", right: "$", display: false }, + { left: "\\(", right: "\\)", display: false }, + { left: "\\[", right: "\\]", display: true } + ], + }) +}) diff --git a/docs/tutorials/Tutorial.md b/docs/tutorials/Tutorial.md new file mode 100644 index 00000000..48f052fa --- /dev/null +++ b/docs/tutorials/Tutorial.md @@ -0,0 +1,255 @@ +In this tutorial we are going to use template matching to detect human 80S ribosomes in a sample of ER-derived microsomes. As the tutorial dataset imaged isolated ER-derived vesicles with a 200 keV microscope, and the ice layer thickness is on the order of ~180 nm, the contrast is quite good. So, to make it more challenging, we are also going to test matching with only the 60S subunit. This way we can demonstrate some optimizations that can be made to maximize correlation. Besides, using [baited reconstruction](https://elifesciences.org/reviewed-preprints/90486) can also be a good way to prevent bias in the template search. + +## Preparing the dataset and template + +The tutorial dataset can be downloaded from DataverseNL: . The raw tilt-series can also be found on DataverseNL: . Tilt-series were collected identically to [EMPIAR-11751](https://www.ebi.ac.uk/empiar/EMPIAR-11751/). A quick overview of the preprocessing workflow: (1) MotionCor2 was used to correct the motion in raw movie frames, (2) ctfplotter from IMOD was used to estimate defocus values, (3) ctfphaseflip from IMOD was used to do strip-based phaseflipping, and (4) AreTomo was used for reconstructing tomograms (weighted-back-projection) with 5x5 local patch alignment. Besides the tomograms (in MRC format), we provide for each tomogram a `.rawtlt` file with tilt angles, a `_dose.txt` file with accumulated dose per frame, and a `.defocus` file with per tilt defocus values. + +We will initially use a cryo-EM reconstruction of the human 80S ribosome, which you will need to download from the EMDB: . + +Secondly, we will use a PDB model of the human 80S ribosome, [6qzp](https://www.rcsb.org/structure/6QZP). To prepare the structure, first download the `.cif` file from the PDB. Then, the structure can be opened in ChimeraX and with the following command in the ChimeraX command line, a simple electrostatic potential model of the 60S large ribosomal subunit can be generated: + +```molmap #1/L? 5 gridSpacing 1.724```. + +This means we sample the map on a 1.724 Å grid (which corresponds to the tilt-series' pixel size of this dataset) and at a resolution of 5 Å. `#1/L?` selects only chains starting with L from loaded model 1, i.e. chains that belong to the large subunit. Afterwards save the generated potential as an .mrc file: `6qzp_60S.mrc`. (It can also be downloaded [here](data/6qzp_60S.mrc)) + +Finally, we assume the following folder structure for this tutorial, with the current working directory `tm_tutorial` : + +``` +tm_tutorial/ ++- dataset/ +¦ +- tomo200528_100.mrc +¦ +- tomo200528_100.rawtlt +¦ +- tomo200528_100.defocus +¦ +- tomo200528_100_dose.txt +¦ +- tomo200528_101.mrc +¦ +- ... ++- templates/ +¦ +- 6qzp_60S.mrc +¦ +- emd_2938.mrc ++- results_80S/ ++- results_60S/ +``` + +### Tomogram voxel size fix + +I made the mistake to not annotate the voxel size in the tomograms correctly (thanks to @jinxsfe for pointing this out!). Please run the following command to fix the voxel size in the MRC headers: + +``` bash +for x in dataset/.mrc; do python -c "import mrcfile; mrc = mrcfile.mmap('$x', 'r+'); mrc.voxel_size = 13.79"; done +``` + +## Part 1: Matching the 80S ribosome with a binary wedge and simple CTF + +This section corresponds to the base method where the template is prepared by convolution with a single CTF and is adjusted with a binary wedge during the template matching run. + +### Preparing the template and mask + +Let's generate a template for the template matching run: + +``` bash +pytom_create_template.py \ + -i templates/emd_2938.map \ + -o templates/80S.mrc \ + --input-voxel 1.1 \ + --output-voxel 13.79 \ + --center \ + --invert \ + -b 60 +``` + +The `--ctf-correction` option convolutes the template with a single 3D-CTF with the specified defocus to roughly match the appearance of the particles in the tomograms. Even though each tilt-series will have different defocus values, for full ribosomes this is usually sufficient. The `--flip-phase` option takes the `absolute()` of the CTF before applying it. This is because the tomograms have been CTF-corrected which is done through phase-flipping. Since the `absolute()` operation makes all the values positive and our input map is white, we get back a white template. However, the tomograms of the dataset have black contrast (veryify by opening a tomogram), so finally the `--invert` flag needs to be set to make the template black. + +Secondly, we need a mask with the same box size as the template and a radius that fully encompasses our ribosome structure. The diameter of a ribosome is roughly 300 Å, which means we need at least (300 Å / 13.79 Å =) 22 pixels to cover the ribosome. With some overhang we extend it to 24 pixels diameter, and therefore set the mask radius (`--radius`) to 12 pixels. The mask will also need a smooth fall-off to prevent aliasing artifacts, which we set with `--sigma`: + +``` bash +pytom_create_mask.py \ + -b 60 \ + -o templates/mask.mrc \ + --voxel-size 13.79 \ + --radius 12 \ + --sigma 1 +``` + +The `--voxel-size` is not really required for the mask, but in case you want to open the mask and template in ChimeraX (or napari) they are automatically scaled to the right size. If you have ChimeraX/napari/3dmod available its anyway a good idea to open the mask and template at this point to get a feel for the data. The mask should always be 1 in the center and 0 outside with a smooth fall-off. The template should have contrast matching your tomogram (i.e. black if the tomogram has black objects, and white vice versa). + +
+ ![template_and_mask](images/0_template_mask.png){ width="400" } +
Template and mask slice
+
+ +The large template and mask box sizes might seem surprising, but they are needed for proper weighting during template matching. Similar to large box size during downstream high-resolution averaging, the overhang here also provides better sampling of weighting functions in Fourier space. The larger size increases the number of sampling points for CTFs and tomographic point spread functions. + +### Running template matching + +The template matching for the 80S ribosome can be started as follows: + +``` bash +pytom_match_template.py \ + -t templates/80S.mrc \ + -m templates/mask.mrc \ + -v dataset/tomo200528_100.mrc \ + -d results_80S/ \ + --particle-diameter 300 \ + -a dataset/tomo200528_100.rawtlt \ + --low-pass 40 \ + --defocus 3 \ + --amplitude 0.08 \ + --spherical 2.7 \ + --voltage 200 \ + --tomogram-ctf-model phase-flip \ + -g 0 +``` + +A low-pass filter of 40 Å is set to prevent bias. Then, the angular search is calculated from the Crowther criterion using a diameter $d$ of 300 Å and a $r_{max}$ of 1/(40 Å) (due to the low-pass filter), which results in a roughly 7° increment ([see this wiki page for more info on the required sampling](https://github.com/SBC-Utrecht/pytom-template-matching-gpu/wiki/template-matching)). Here, the program was started on the system GPU with index 0. If you have more GPU's available for the tutorial you can simply add them and pytom-match-pick will automatically distribute the angular search, e.g. `-g 0 1 2 3` will run on the first 4 GPU's in the system. + +**Out of memory!** Although unlikely to happen during the tutorial, in later usage of the software your GPUs might run out of memory. This is especially likely for 4x or lower binned tomograms. If this happens, you can use the `--split` option to tile the tomogram into multiple sections. You need to specify the number of splits along each dimension, for example `--split 2 1 1` will split into 2 sections along the x-axis, `--split 2 2 1` will split once in x and once in y, resulting in 4 subvolumes! + +To visualize the results it can be very informative to open the tomogram and score map side-by-side in Napari. Which should look like this: + +
+ ![part1_score_map](images/1_tomo200528_100_slice69_comb.png){ width="800" } +
Slice 69 of dataset/tomo200528_100.mrc and results_80S/tomo200528_100_scores.mrc (in Napari)
+
+ +Then we calculate the receiver operator curve (ROC) to estimate true positives using: + +``` bash +pytom_estimate_roc.py \ + -j results_80S/tomo200528_100_job.json \ + -n 800 \ + -r 8 \ + --bins 16 \ + --crop-plot > results_80S/tomo200528_100_roc.log +``` + +This will automatically write a file to the folder where the job file is located (`tomo200528_100_roc.svg`). Additionally, here the terminal output is also written to the file `results_80S/tomo200528_100_roc.log` as it contains some info on the estimated number of particles in the sample. Catting the log file in the terminal (`cat results_80S/tomo200528_100_roc.log`), the estimated number of total true positives is approx. 180 ribosomes. + +
+ ![part1_roc](images/1_tomo200528_100_roc_80S.svg){ width="800" } +
ROC curve (results_80S/tomo200528_100_roc.svg)
+
+ +In this case the result is quite good: the Rectangle Under the Curve (RUC) is around 0.9 (1 is optimal). Also the fit of the gaussians seems very appropriate for the histogram, so we can use the estimated cut-off value (`0.194`) to extract particle coordinates and rotations into a STAR file: + +``` bash +pytom_extract_candidates.py -j results_80S/tomo200528_100_job.json -n 300 -r 8 -c 0.194 +``` + +This will extract ~160 particles and save their annotations to the file `results_80S/tomo200528_100_particles.star`. + +Alternatively, `pytom_extract_candidates.py` also has a default extraction cut-off estimation that can work well. Run the command without the `-c` parameter: + +``` bash +pytom_extract_candidates.py -j results_80S/tomo200528_100_job.json -n 300 -r 8 +``` + +It produces a very similar result: the script prints the estimated cut-off to the terminal (`0.197`) (very similar to the ROC estimated cut-off) and it also extracts ~160 particles. (Take care: the previous STAR file is overwritten.) + +## Part 2: Matching the 60S ribosomal subunit with a tilt-weighted PSF + +Now, let's apply the same base method using only the 60S subunit. + +### Preparing the template and mask + +The 60S template can be generated with this command: + +``` bash +pytom_create_template.py \ + -i templates/6qzp_60S.mrc \ + -o templates/60S.mrc \ + --input-voxel 1.724 \ + --output-voxel 13.79 \ + --center \ + --invert \ + -b 60 +``` + +For the mask the radius is slightly reduced to constrain it around the 60S subunit: + +``` bash +pytom_create_mask.py \ + -b 60 \ + -o templates/mask_60S.mrc \ + --voxel-size 13.79 \ + --radius 10 \ + --sigma 1 +``` + +### Running template matching + +Template matching is run for the 60S subunit, where the low-pass filter is now set +to 35 Å to include more high-resolution information in the template matching. This +results in an angular sampling ~6.5°. We +could also entirely remove it, but this would increase the sampling to 5.3° +which is less feasible for the tutorial. This potentially introduces more bias, +however: with only the large subunit as a reference we can always test for detection +of the small subunit in subtomogram averages. + +The tilt-weighted PSF is also used here with the flag `--per-tilt-weighting` and by +providing an IMOD-style defocus file (.defocus) and a text file dose accumulation. +Furthermore, we run with `--random-phase-correction` to flatten the background noise +using the method introduced by STOPGAP. + +``` bash +pytom_match_template.py \ + -t templates/60S.mrc \ + -m templates/mask_60S.mrc \ + -v dataset/tomo200528_100.mrc \ + -d results_60S/ \ + --particle-diameter 300 \ + -a dataset/tomo200528_100.rawtlt \ + --per-tilt-weighting \ + --low-pass 35 \ + --defocus dataset/tomo200528_100.defocus \ + --amplitude 0.08 \ + --spherical 2.7 \ + --voltage 200 \ + --tomogram-ctf-model phase-flip \ + --dose-accumulation dataset/tomo200528_100_dose.txt \ + --random-phase \ + -g 0 +``` + +Again, the slice and score map can be visualized in Napari: + +
+ ![part2_score_map](images/2_tomo200528_100_slice69_simple.png){ width="800" } +
Slice 69 of dataset/tomo200528_100.mrc and results_60S/tomo200528_100_scores.mrc (in Napari)
+
+ +Matching only with the 60S subunit greatly reduces the peak height at the actual ribosome locations. Instead gold markers now start correlating heavily with the template and get higher values than the particles of interest. The reason is two-fold: (1) the size of the template has become smaller and (2) the switch from an EM-map to a PDB based model usually reduces correlation. The second point (2) can be caused by the electrostatic potential modelling with the simplistic molmap, which does not consider atom number and solvent embedding, but also because EM-maps might contain some additional heterogeneous densities. + +The ROC curve is now estimated with: + +``` bash +pytom_estimate_roc.py \ + -j results_60S/tomo200528_100_job.json \ + -n 800 \ + -r 8 \ + --bins 16 \ + --crop-plot > results_60S/tomo200528_100_roc.log +``` + +
+ ![part2_roc](images/2_tomo200528_100_roc_simple.svg){ width="800" } +
ROC curve (results_60S/tomo200528_100_roc.svg)
+
+ +The loss of quality is mainly visible from the gaussian fit (left plot). Even though the reported quality is still quite good (RUC), it is clear that the curve fit is not proper. This is due to the three mixed distributions that now appear: the left-most correctly fitted background noise (blue curve), the middle peak is from 60S subunits and the right-most peak from gold markers. The particle of interest and gold markers are both fitted as true positives (orange curve). + +Sadly, gold markers (but also carbon film, and ice) can quite often interfere with +template matching because of their high scattering potential leading to edges with very high signal-to-noise ratio (SNR). One way of dealing with this, is gold marker removal for which their are tools in IMOD, and also deep-learning based tools (e.g. [fidder](https://teamtomo.org/fidder/)), that remove gold markers on the tilt-image level before tomographic reconstruction. + +## Increasing the angular sampling + +To improve the results, the angular sampling needs to be increased. For +part 2 of this tutorial, we set the low-pass filter to 35 Å. This means template +matching is now including resolution up to 1/(35 Å), +which changes the required angular sampling to 6.5°. Running the template matching +command from part 2 without a low-pass filter increases the angular increment to 5.3°. +This roughly doubles the runtime. For me, it takes 1h20m on one RTX3060. + +## Conclusion + +We hope you enjoyed our tutorial on GPU-acculated template matching in +cryo-ET! If you ran into any issues during the tutorial please let us know on our [github issues page](https://github.com/SBC-Utrecht/pytom-match-pick/issues). On this repository's [Discussions page](https://github.com/SBC-Utrecht/pytom-match-pick/discussions) you can also reach out with specific questions about running the software for your dataset, and find previous discussions. diff --git a/docs/tutorials/data/6qzp_60S.mrc b/docs/tutorials/data/6qzp_60S.mrc new file mode 100644 index 00000000..14749f1c Binary files /dev/null and b/docs/tutorials/data/6qzp_60S.mrc differ diff --git a/docs/tutorials/images/0_template_mask.png b/docs/tutorials/images/0_template_mask.png new file mode 100644 index 00000000..e6e30d79 Binary files /dev/null and b/docs/tutorials/images/0_template_mask.png differ diff --git a/docs/tutorials/images/1_tomo200528_100_roc_80S.svg b/docs/tutorials/images/1_tomo200528_100_roc_80S.svg new file mode 100644 index 00000000..385a4a96 --- /dev/null +++ b/docs/tutorials/images/1_tomo200528_100_roc_80S.svg @@ -0,0 +1,1579 @@ + + + + + + + + 2023-12-11T12:18:19.475489 + image/svg+xml + + + Matplotlib v3.8.0, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/tutorials/images/1_tomo200528_100_slice69_comb.png b/docs/tutorials/images/1_tomo200528_100_slice69_comb.png new file mode 100644 index 00000000..7b25860f Binary files /dev/null and b/docs/tutorials/images/1_tomo200528_100_slice69_comb.png differ diff --git a/docs/tutorials/images/2_tomo200528_100_roc_simple.svg b/docs/tutorials/images/2_tomo200528_100_roc_simple.svg new file mode 100644 index 00000000..814503a8 --- /dev/null +++ b/docs/tutorials/images/2_tomo200528_100_roc_simple.svg @@ -0,0 +1,1593 @@ + + + + + + + + 2023-12-11T16:00:53.661043 + image/svg+xml + + + Matplotlib v3.8.0, https://matplotlib.org/ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/tutorials/images/2_tomo200528_100_slice69_simple.png b/docs/tutorials/images/2_tomo200528_100_slice69_simple.png new file mode 100644 index 00000000..278fa33c Binary files /dev/null and b/docs/tutorials/images/2_tomo200528_100_slice69_simple.png differ diff --git a/docs/tutorials/tests/test_tutorial.py b/docs/tutorials/tests/test_tutorial.py new file mode 100644 index 00000000..6d84e634 --- /dev/null +++ b/docs/tutorials/tests/test_tutorial.py @@ -0,0 +1,10 @@ +# this tests the tutorial.md +import subprocess +from mdextractor import extract_md_blocks + +lines = "".join(open("Tutorial.md").readlines()) +blocks = extract_md_blocks(lines) +for block in blocks: + if block.split()[0].endswith(".py"): + print(f"Running: {block}") + subprocess.run(block) diff --git a/mkdocs.yml b/mkdocs.yml new file mode 100644 index 00000000..03c3d875 --- /dev/null +++ b/mkdocs.yml @@ -0,0 +1,53 @@ +site_name: pytom-match-pick +site_url: https://SBC-Utrecht.github.io/pytom-match-pick +site_author: Marten Chaillet +repo_name: pytom-match-pick +repo_url: https://github.com/SBC-Utrecht/pytom-match-pick +copyright: Copyright © 2024 + +nav: + - Home: + - Overview: index.md + - Usage: Usage.md + - Timings: Timings.md + - FAQ: Frequently-asked-questions-(FAQ).md + - For developers: Developers.md + - Tutorial: + - Ribosomes on ER microsomes: tutorials/Tutorial.md + - Discussions: https://github.com/SBC-Utrecht/pytom-match-pick/discussions + - Report a bug: https://github.com/SBC-Utrecht/pytom-match-pick/issues + +theme: + name: material + features: + - navigation.instant + # - navigation.instant.progress # possible progress bar, but site is too small to justify + - navigation.tracking + - navigation.tabs + - navigation.tabs.sticky + - navigation.expand + - toc.follow + - content.code.copy + +markdown_extensions: + # - include + - pymdownx.snippets: + base_path: ["README.md", "src/pytom_tm/entry_points.py"] + check_paths: true + - pymdownx.superfences + - pymdownx.arithmatex: + generic: true + - attr_list + - md_in_html + +extra_javascript: + - javascripts/katex.js + - https://unpkg.com/katex@0/dist/katex.min.js + - https://unpkg.com/katex@0/dist/contrib/auto-render.min.js + +extra_css: + - https://unpkg.com/katex@0/dist/katex.min.css + +plugins: + - search + - markdown-exec diff --git a/pyproject.toml b/pyproject.toml index bceae9e9..bfc48181 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -39,7 +39,7 @@ dependencies = [ [project.optional-dependencies] plotting = ["matplotlib", "seaborn"] -dev = ["ruff", "pre-commit"] +dev = ["ruff", "pre-commit", "mkdocs", "mkdocs-material", "markdown-exec"] all = ["pytom-match-pick[plotting]", "pytom-match-pick[dev]"] diff --git a/src/pytom_tm/entry_points.py b/src/pytom_tm/entry_points.py index ec5eab35..fc349fd1 100644 --- a/src/pytom_tm/entry_points.py +++ b/src/pytom_tm/entry_points.py @@ -34,9 +34,12 @@ def pytom_create_mask(argv=None): argv = _parse_argv(argv) + # ---8<--- [start:create_mask_usage] + parser = argparse.ArgumentParser( + prog="pytom_create_mask.py", description="Create a mask for template matching. " - "-- Marten Chaillet (@McHaillet)" + "-- Marten Chaillet (@McHaillet)", ) parser.add_argument( "-b", @@ -98,6 +101,9 @@ def pytom_create_mask(argv=None): "Values in the range from 0.5-1.0 are usually sufficient for tomograms with " "20A-10A voxel sizes.", ) + + # ---8<--- [end:create_mask_usage] + argv = _parse_argv(argv) args = parser.parse_args(argv) @@ -126,9 +132,13 @@ def pytom_create_template(argv=None): from pytom_tm.template import generate_template_from_map argv = _parse_argv(argv) + + # ---8<--- [start:create_template_usage] + parser = argparse.ArgumentParser( + prog="pytom_create_template.py", description="Generate template from MRC density. " - "-- Marten Chaillet (@McHaillet)" + "-- Marten Chaillet (@McHaillet)", ) parser.add_argument( "-i", @@ -190,15 +200,13 @@ def pytom_create_template(argv=None): help="Specify a desired size for the output box of the template. " "Only works if it is larger than the downsampled box size of the input.", ) - ( - parser.add_argument( - "--invert", - action="store_true", - default=False, - required=False, - help="Multiply template by -1. " - "WARNING: not needed if ctf with defocus is already applied!", - ), + parser.add_argument( + "--invert", + action="store_true", + default=False, + required=False, + help="Multiply template by -1. " + "WARNING: not needed if ctf with defocus is already applied!", ) parser.add_argument( "-m", @@ -216,6 +224,9 @@ def pytom_create_template(argv=None): action=ParseLogging, help="Can be set to `info` or `debug`", ) + + # ---8<--- [end:create_template_usage] + args = parser.parse_args(argv) logging.basicConfig(level=args.log, force=True) @@ -273,9 +284,12 @@ def estimate_roc(argv=None): argv = _parse_argv(argv) from pytom_tm.plotting import plist_quality_gaussian_fit + # ---8<--- [start:estimate_roc_usage] + parser = argparse.ArgumentParser( + prog="pytom_estimate_roc.py", description="Estimate ROC curve from TMJob file. " - "-- Marten Chaillet (@McHaillet)" + "-- Marten Chaillet (@McHaillet)", ) parser.add_argument( "-j", @@ -350,6 +364,9 @@ def estimate_roc(argv=None): action=ParseLogging, help="Can be set to `info` or `debug`", ) + + # ---8<--- [end:estimate_roc_usage] + args = parser.parse_args(argv) logging.basicConfig(level=args.log, force=True) @@ -389,8 +406,12 @@ def estimate_roc(argv=None): def extract_candidates(argv=None): argv = _parse_argv(argv) + + # ---8<--- [start:extract_candidates_usage] + parser = argparse.ArgumentParser( - description="Run candidate extraction. -- Marten Chaillet (@McHaillet)" + prog="pytom_extract_candidates.py", + description="Run candidate extraction. -- Marten Chaillet (@McHaillet)", ) parser.add_argument( "-j", @@ -486,6 +507,9 @@ def extract_candidates(argv=None): action=ParseLogging, help="Can be set to `info` or `debug`", ) + + # ---8<--- [end:extract_candidates_usage] + args = parser.parse_args(argv) logging.basicConfig(level=args.log, force=True) @@ -514,8 +538,12 @@ def match_template(argv=None): from pytom_tm.parallel import run_job_parallel argv = _parse_argv(argv) + + # ---8<--- [start:match_template_usage] + parser = argparse.ArgumentParser( - description="Run template matching. -- Marten Chaillet (@McHaillet)" + prog="pytom_match_template.py", + description="Run template matching. -- Marten Chaillet (@McHaillet)", ) io_group = parser.add_argument_group("Template, search volume, and output") io_group.add_argument( @@ -581,7 +609,7 @@ def match_template(argv=None): "diameter. If given a float it will generate an angle list with healpix " "for Z1 and X1 and linear search for Z2. The provided angle will be used " "as the maximum for the " - "linear search and for the mean angle difference from healpix.\n" + "linear search and for the mean angle difference from healpix." "Alternatively, a .txt file can be provided with three Euler angles " "(in radians) per line that define the angular search. " "Angle format is ZXZ anti-clockwise (see: " @@ -811,6 +839,9 @@ def match_template(argv=None): action=ParseLogging, help="Can be set to `info` or `debug`", ) + + # ---8<--- [end:match_template_usage] + args = parser.parse_args(argv) logging.basicConfig(level=args.log, force=True) @@ -903,11 +934,14 @@ def match_template(argv=None): def merge_stars(argv=None): import pandas as pd + # ---8<--- [start:merge_stars_usage] + parser = argparse.ArgumentParser( + prog="pytom_merge_stars.py", description=( "Merge multiple star files in the same directory. " "-- Marten Chaillet (@McHaillet)" - ) + ), ) parser.add_argument( "-i", @@ -937,6 +971,9 @@ def merge_stars(argv=None): action=ParseLogging, help="Can be set to `info` or `debug`", ) + + # ---8<--- [end:merge_stars_usage] + args = parser.parse_args(argv) logging.basicConfig(level=args.log, force=True)