Difference between revisions of "Walkthrough for template matching"

From Dynamo
Jump to navigation Jump to search
Line 110: Line 110:
 
[[  File:TemplateMatchingControlThickness.png |thumb|center|500px | Thickness controls  ]]
 
[[  File:TemplateMatchingControlThickness.png |thumb|center|500px | Thickness controls  ]]
  
Then, you go through each particle and indicate a point that you considere that could lay in the central axis of the particle by clicking on <tt>N</yy>. The particle will rotate on screen to match this direction.  The center of the particle can also be set by clicking on <tt>N</tt>
+
Then, you go through each particle and indicate a point that you considere that could lay in the central axis of the particle by clicking on <tt>N</tt>. The particle will rotate on screen to match this direction.  The center of the particle can also be set by clicking on <tt>N</tt>
  
 
[[  File:TemplateMatchingClickNorth.png |thumb|center|500px| Select the north direction with <tt>N</tt>  ]]
 
[[  File:TemplateMatchingClickNorth.png |thumb|center|500px| Select the north direction with <tt>N</tt>  ]]
Line 125: Line 125:
  
 
===== Measure distances =====
 
===== Measure distances =====
 +
 +
 
[[  File:TemplateMatchingMarkZaxis.png |thumb|center|500px  ]]
 
[[  File:TemplateMatchingMarkZaxis.png |thumb|center|500px  ]]
  
Line 137: Line 139:
 
[[  File:TemplateMatchingMaskSave.png |thumb|center|500px  ]]
 
[[  File:TemplateMatchingMaskSave.png |thumb|center|500px  ]]
  
 +
====Cropping the borders of a template====
  
====Cropping the borders of a template====
+
<tt>shifted =circshift(dsym(o.average,'c17'),[0,0,-1]);dslices(final,'y','overlay','maskTight.em','overlay_as','mask'); </tt>
  
 
[[  File:TemplateMatchingCircshiftTemplate.png |thumb|center|500px  ]]
 
[[  File:TemplateMatchingCircshiftTemplate.png |thumb|center|500px  ]]
 +
 +
and we can actually extract the part of the template volume where we actually have signal.
 +
 +
<tt>cropped = dynamo_crop(shifted ,32);</tt>
 
[[ File:TemplateMatchingTemplateInterior.png |thumb|center|500px| Template cropped to 32 pixels ]]
 
[[ File:TemplateMatchingTemplateInterior.png |thumb|center|500px| Template cropped to 32 pixels ]]
  
Line 146: Line 153:
 
Alternatively, you can use <tt>dynamo_mask</tt> or <tt>dynamo_tube</tt> to create a synthetic model.
 
Alternatively, you can use <tt>dynamo_mask</tt> or <tt>dynamo_tube</tt> to create a synthetic model.
  
== Creating a cross correlation process==
+
== Creating a template matching process==
 +
The main function is <tt>dynamo_match</tt>
  
  <nowiki>ptm = dynamo_match(fileFull,'average48.em','mask','maskData48.em',....
+
  <nowiki>pts = dynamo_match(fileFull,'average32.em','mask','maskData32.em',....
 
     'outputFolder','cs30',...
 
     'outputFolder','cs30',...
     'ytilt',[-39,36],'sc',[1000,1000,300],'cr',360,'cs',30,'bin',1);</nowiki>
+
     'ytilt',[-39,36],'sc',[256,256,256],'cr',360,'cs',30,'bin',1);</nowiki>
 +
 
 +
* 'bin' <br/>  allows binning on the fly. Template, mask and tomograms are input in original size.
 +
* 'sc' : size of chunk <br/> This is the maximum extent of a block in the unbinned tomogram that will be kept in memory by each one of the processors working in parallel (only one in this experiment). We need to be careful not to crowd the available memory.
 +
* 'cr' : cone range <br/> orientations will be looked for inside a cone. In this context, the most usual values are 360 (sample the full sphere) or 0 (do not scan orientations)
 +
 
 +
 
 +
You should get the following response on the screen (will take around 5 minutes to complete)
  
 
  <nowiki>------------------------------------------------------------
 
  <nowiki>------------------------------------------------------------
Line 161: Line 176:
 
   - Mb per tile (operation) : 663.84
 
   - Mb per tile (operation) : 663.84
 
  ... initializing output elements
 
  ... initializing output elements
Preparing to run on 1 blocks
+
Preparing to run on 25 blocks
 
Running on single processor
 
Running on single processor
 
  Computing tile 1
 
  Computing tile 1
  -Range (original block): x[1:958] y[1:926] z[1:214]
+
 
  -Range (binned block)  : x[1:479] y[1:463] z[1:107]
+
  <Information for each block>
 +
 
 
   ... tile 1 finished in 387 seconds (8 for setup;  7.73 per triplet)
 
   ... tile 1 finished in 387 seconds (8 for setup;  7.73 per triplet)
 
   [ok] ... template matching process completed upto creation of cross correlation
 
   [ok] ... template matching process completed upto creation of cross correlation
 
       you can proceed to peak location and particle extraction.
 
       you can proceed to peak location and particle extraction.
 
------------------------------------------------------------</nowiki>
 
------------------------------------------------------------</nowiki>
 +
 +
The main result of <tt>dynamo_match</tt> is  a volume (''cross correlation'' volume) that assigns a score to each pixel. This tomogram is written in the file <tt>cs30/CC.mrc</tt>. Some auxiliary outputs are other volumes that store the angles [tdrot,tilt,narot] that maximize the cross-correlation for each pixel in the volume.
  
 
=== The <tt>Process</tt> object ===
 
=== The <tt>Process</tt> object ===
The return of the function <notwiki>ptm</nowiki> is an object of class XXX. In short, this is an object created
+
The return of the function <notwiki>ptm</nowiki> is an object that keeps track of the location of all the created output, and is used in the subsequent steps that extract the actual location of the peaks
 
 
 
 
=== Considerations when creating a process===
 
  
 
== Locating cross correlation peaks  ==
 
== Locating cross correlation peaks  ==
Line 181: Line 196:
 
=== Looking at the cross correlation map ===
 
=== Looking at the cross correlation map ===
  
[[  File:TemplateMatchingCCRaw.png |thumb|center|500px  ]]
+
<tt>pts.showCC</tt>
 +
 
 +
[[  File:TemplateMatchingCCRaw.png |thumb|center|500px| <tt>tmshow</tt> on the CC volume ]]
  
 
=== Looking at the cross correlation profile ===
 
=== Looking at the cross correlation profile ===
  
[[  File:TemplateMatchingPeakGood.png |thumb|center|500px  ]]
+
We can get a plot of the cross-correlation value found on the local maxima of the cc volume wiith the order
[[  File:TemplateMatchingPeakGoodVolume.png |thumb|center|500px  ]]
+
 
[[  File:TemplateMatchingPeakKink.png |thumb|center|500px  ]]
+
<tt>pts.peaks.plotCCPeaks('sidelength',32)</tt>
[[  File:TemplateMatchingPeakBad.png |thumb|center|500px  ]]
+
 
 +
The cross-correlation values of the peaks appear in ascending order. We can check the quality of the peaks by auxiliary clicking on the curve to select one particle and then selecting some visualization option.
 +
 
 +
[[  File:TemplateMatchingPeakGood.png |thumb|center|500px| auxiliary click on the area of good peaks ]]
 +
 
 +
The <tt>sidelength</tt> that we pass will be used for this visualization options to crop the particle around the peak.
 +
 
 +
[[  File:TemplateMatchingPeakGoodVolume.png |thumb|center|500px  |  particle cropped around the selected correlation peak]]
 +
 
 +
We click on a couple of particles in the area of kink in the cross-correlation to roughly estimate the cross-correlation threshold. In out case, it seems to be around 0.37 (being rather conservative)
 +
 
 +
[[  File:TemplateMatchingPeakKink.png |thumb|center|500px | Looking for an area of transition between good and bad peaks]]
 +
 
 +
In the gaussian distribution below this ara , the peaks do not correspond to particles
  
 +
[[  File:TemplateMatchingPeakBad.png |thumb|center|500px | false positives ]]
  
 
=== Extracting a table ===
 
=== Extracting a table ===
Line 199: Line 230:
  
 
=== Looking at  table positions===
 
=== Looking at  table positions===
 +
 
[[  File:TemplateMatchingTable3d.png |thumb|center|500px  ]]
 
[[  File:TemplateMatchingTable3d.png |thumb|center|500px  ]]
 +
 
[[  File:TemplateMatchingTable3dSide.png |thumb|center|500px  ]]
 
[[  File:TemplateMatchingTable3dSide.png |thumb|center|500px  ]]
  
Line 205: Line 238:
 
We can check how the individual particles look like on a gallery modus:
 
We can check how the individual particles look like on a gallery modus:
  
<tt>pff.peaks.browse();</tt>
+
<tt>pr.peaks.browse();</tt>
  
 
This order just opens <tt>ddbrowse</tt>. We are using here the support object <tt>peaks</tt>, but this command is equivalent to just invoking <tt>ddbrowse</tt>
 
This order just opens <tt>ddbrowse</tt>. We are using here the support object <tt>peaks</tt>, but this command is equivalent to just invoking <tt>ddbrowse</tt>

Revision as of 21:23, 21 August 2017


Dynamo includes a set of tools for location of particles inside tomograms. The most simple one is template matching.

Template matching

In this technique, a template representing a molecule of interest is systematically cross-correlated against a tomogram, producing a cross-correlation map of the tomogram. Each pixel in this map represents a score assigned the corresponding pixel in the tomogram map. This score measures the similarity of the neighbourhood of the tomogram pixel to the used template. This similarity is measured exclusively inside a mask.

Data set

Tomogram description

The tomogram contains a buffer with T20S proteasome on a holey carbon grid collected on a Krios + K2. Original pixelsize was 1.76 angstroms. The tomogram provided here has bin binned twice (yielding thus 7.04 ang), defocus is 4.4 microns, no CTF correction, no energy filter.

Acknowledgements

The tomogram has been kindly provided by Alex Noble, from the New York Structural Biology Center. Data collection was performed using Leginon and Appion-Protomo at the Simons Electron Microscopy Center and National Resource for Automated Molecular Microscopy located at the New York Structural Biology Center, supported by grants from the Simons Foundation (349247), NYSTAR, and the NIH National Institute of General Medical Sciences (GM103310) with additional support from Agouron Institute [Grant Number: F00316].

Getting the tomogram

In Linux, you can write:

wget  https://wiki.dynamo.biozentrum.unibas.ch/w/doc/data/t20s/t20s.mrc

the corresponding command is curl -O in the Mac.

Visualizing the tomogram

We can get a first glance on how the tomogram looks like:

dtmshow -otf t20s.mrc  

to use the on-the-fly access to the slice shown at each given moment, or

v = dread('t20s.mrc');dtmshow(v); 

to preload the full tomogram into a memory variable (arbitrarily called v). In either option, you will see, the proteasomes are densely packed in an layer. The layer is slightly oblique, what can be seen browsing through z or y.

dtmshow on the pt20s tomogram, View along z

Navigating on-the-fly you'll see that transitions in y are slower than transitions in z, because all the pixels of the same slice are stored sequentially in the disk.

dtmshow on the pt20s tomogram, View along y


Estimation of the missing wedge

Estimation of the missing wedge'

Creating a template

There are different strategies to create the first template. In the case where the general shape of each protein is roughly recognisable by eye, it is not difficult to just crop and align manually some of the particles. When this is not possible, you have the option of using a density map that mimics the general topology of your protein.

Through manual alignment

We will use a model inside tmslice to pick a set of ~10 particles. We will then crop them, and align them manually using dgallery

Manual selection of some particles

We use for this our tool dtmslice. As the tomogram is provided is fairly small you can probably just open it without any further binning.

dtmslice t20s.mrc -cc ct20

dtmslice on tomogram

Use the model pool menu to add a new model of type General

Opening a model

Try to pick particles on different orientations, both in top and side views.

Pick particles with the C key

You can move the slice with the z slidebar in the control panel, or dragging it while the mouse button is pressed. Top views are more abundant below the crowded accumulation of side views in the central slice.

Use the wheel to go up and down in the tomogram to locate top views

When you have clicked around ten, we will prepare to crop them out of the tomogram. For this, we need to get first an idea of how big should the box. With the keys [1] and [2] you can place anchor points in the scene. Auxiliary click in the bar between them shows the distance (in pixels of the unbinned tomogram).

Measure distances with the [1] and [2] anchors

Now we can open the dtcrop utility for this model, by visiting the Active model" menu.

TemplateMatchingMenuForDtcrop.png

In the GUI that pops up we can fix the sidelength.

TemplateMatchingDtcrop.png

After completion of the cropping process, you will get a data folder called t20s_Particles, which will also contain a table file called t20s_Particles/crop.tbl

Manual alignment of some particles

Now we want to align the particles. We open them inside the dgallery browser:

dgallery -d t20s_Particles -t t20s_particles/crop.tbl 

This GUI has several functionalities. We will focus in using it to quickly get an alignment of the particles. For this, we first load the particles expressed in the table into the internal memory of the gallery by clicking into Load.

Preload all the particles into memory

This seems not to have any effect in the depiction: you need to fix the range of the depicted particles moving the slider bar. When you do so, each particle will appear in the scene as a single red box.

Adjust the scene to show all the particles

You can play with different settings for viewing, like the direction (x, y or z)

viewing direction controls

or the number of slices that are projected into the depicted image for each particle.

Thickness controls

Then, you go through each particle and indicate a point that you considere that could lay in the central axis of the particle by clicking on N. The particle will rotate on screen to match this direction. The center of the particle can also be set by clicking on N

Select the north direction with N
Select the north direction with N</tt

You may need to complete two rounds of operating on the particles, changing the viewing direction (first in x, then in y)

Averaging manually aligned particles

TemplateMatchingAverageManual.png

Creating a tight mask

Measure distances
TemplateMatchingMarkZaxis.png
TemplateMatchingMarkXYPlane.png
TemplateMatchingClickNorthTopViewX.png
Create mask
TemplateMatchingMaskOptions.png


TemplateMatchingMaskSave.png

Cropping the borders of a template

shifted =circshift(dsym(o.average,'c17'),[0,0,-1]);dslices(final,'y','overlay','maskTight.em','overlay_as','mask'); 
TemplateMatchingCircshiftTemplate.png
and we can actually extract the part of the template volume where we actually have signal.
cropped = dynamo_crop(shifted ,32);
Template cropped to 32 pixels

Through geometrical shapes

Alternatively, you can use dynamo_mask or dynamo_tube to create a synthetic model.

Creating a template matching process

The main function is dynamo_match

pts = dynamo_match(fileFull,'average32.em','mask','maskData32.em',....
    'outputFolder','cs30',...
    'ytilt',[-39,36],'sc',[256,256,256],'cr',360,'cs',30,'bin',1);
  • 'bin'
    allows binning on the fly. Template, mask and tomograms are input in original size.
  • 'sc' : size of chunk
    This is the maximum extent of a block in the unbinned tomogram that will be kept in memory by each one of the processors working in parallel (only one in this experiment). We need to be careful not to crowd the available memory.
  • 'cr' : cone range
    orientations will be looked for inside a cone. In this context, the most usual values are 360 (sample the full sphere) or 0 (do not scan orientations)


You should get the following response on the screen (will take around 5 minutes to complete)

------------------------------------------------------------
 
 Template matching process.
 computing in CPU
 Output folder: cs30.TM
A total of 1 tiles have been created
  - Mb per tile (reading)   : 724.19
  - Mb per tile (operation) : 663.84
 ... initializing output elements
Preparing to run on 25 blocks
Running on single processor
 Computing tile 1

  <Information for each block>

  ... tile 1 finished in 387 seconds (8 for setup;  7.73 per triplet)
  [ok] ... template matching process completed upto creation of cross correlation
       you can proceed to peak location and particle extraction.
------------------------------------------------------------

The main result of dynamo_match is a volume (cross correlation volume) that assigns a score to each pixel. This tomogram is written in the file cs30/CC.mrc. Some auxiliary outputs are other volumes that store the angles [tdrot,tilt,narot] that maximize the cross-correlation for each pixel in the volume.

The Process object

The return of the function <notwiki>ptm</nowiki> is an object that keeps track of the location of all the created output, and is used in the subsequent steps that extract the actual location of the peaks

Locating cross correlation peaks

Looking at the cross correlation map

pts.showCC
tmshow on the CC volume

Looking at the cross correlation profile

We can get a plot of the cross-correlation value found on the local maxima of the cc volume wiith the order

pts.peaks.plotCCPeaks('sidelength',32)

The cross-correlation values of the peaks appear in ascending order. We can check the quality of the peaks by auxiliary clicking on the curve to select one particle and then selecting some visualization option.

auxiliary click on the area of good peaks

The sidelength that we pass will be used for this visualization options to crop the particle around the peak.

particle cropped around the selected correlation peak

We click on a couple of particles in the area of kink in the cross-correlation to roughly estimate the cross-correlation threshold. In out case, it seems to be around 0.37 (being rather conservative)

Looking for an area of transition between good and bad peaks

In the gaussian distribution below this ara , the peaks do not correspond to particles

false positives

Extracting a table

A table can be extracted through:

myTable = pff.peaks.computeTable('mcc',0.378);

Visual evaluation of results

Looking at table positions

TemplateMatchingTable3d.png
TemplateMatchingTable3dSide.png

Looking at cropped particles

We can check how the individual particles look like on a gallery modus:

pr.peaks.browse();

This order just opens ddbrowse. We are using here the support object peaks, but this command is equivalent to just invoking ddbrowse

ddbrowse('d','cr30.TM','t',myTable);

ddbrowse with a tomogram as a data source, needing the input of a sidelength in pixels
Result launched by ddbrowse without a connected table
The table must be connected to align the particles
z-view of aligned particles
x-view of aligned particles
y-view of aligned particles

Looking at averages

o = pff.peaks.average();

TemplateMatchingAverageRaw.png
TemplateMatchingAverageRandomized.png

Getting the tables back into the catalogue

You don't need to operate with the catalogue: you could just use dtcrop on the just computed table to extract particles from the original tomogram into a data folder. However, cataloguing has to advantages:

  • Eases keeping track of all the steps you've performed, specially if you are going to use several tomograms, and..
  • allows you to visualize the peaks on the tomogram interactively, so that you can delete false positives or add peaks that were not located by the template matching.

Rescaling the table

Remember that the table that we are working with has the scale of the auxiliary binned volume that has been created along with the CC matrix. In the catalogue, we want to work with the particles in their original scale.

tableOriginalScale = dynamo_table_rescale(myTable,'factor',2);

In the syntax of dynamo_table_rescale, the factor is expressed in terms of how many times is the apix in the original table bigger than in the target table to be computed. In our case, the target table was computed with an apix of 14 Angstrom (one time binned in relation to the tomogram t20s.rec, with an apix of 7.04A). The factor is thus 2. In case of doubt, it is convenient to just run

dtinfo(tableOriginalScale);

to check the extent of the entries in the columns 24, 25 and 26, which are positions (in pixels) of the particles indexed by the table. Now, if we write the upscaled table into a file

dwrite(tableOriginalScale,'peaks.tbl');

we can just put this table back into the catalogue through:

dmimport -t peaks.tbl -c ct20s -i 1 -mn ccpeaks

which will create a model called ccpeaks and assign it to the first volume (-i 1) in the catalogue ct20s. You can check it by writing:

dcm -c ct20s -i 1 -l m

which asks for a listing of models (-l m) in the first volume of the catalogue.


TemplateMatchingIntoCatalogue.png
TemplateMatchingIntoCatalogueZoom.png