Walkthrough on localized reconstruction

In this walkthrough, we use localized reconstruction to reconstruct a set of subtomograms at full pixel resolution without the requirement of creating a full reconstruction first. In this approach, we:

• create a reconstruction at low resolution
• mark the coordinates of interest in the low resolution reconstruction
• reconstruct separately the different coordinates at full resolution

For the sake of simplicity, in this walkthrough, the landmarks of interest will be simply the gold beads, which are trivially recognizable both in tomograms and in projections in the tilt series.

Example data set

The data set is an aligned tilt series. No filtering has been performed on it.

Source

The tilt series was used in the publication Cryo-EM structure of the extended type VI secretion system sheath-tube complex (J Wang et. al - Nature microbiology, 2017).

Download tilt series file

The aligned tilt series can be download with the command:

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

in a linux shell, or through:

curl -O  https://wiki.dynamo.biozentrum.unibas.ch/w/doc/data/t6ss/vc17_align.mrc

under MacOS.

This will create a file called vc17_align.ali in your current working directory. The tilt angles go from -60 to 60 in intervals of 3 degrees.

The file can be brought into memory through the command

ts = dread('vc17_align.ali');

creating the variable ts (name has been chosen arbitrarily) which contains the tilt series as a numerical matrix.

Inspecting the tilt series

The tilt series at full pixel resolution can be inspected by:

dtmshow(ts)

tmshow of raw tilt series

Prebinning during inspection

You will notice that inspecting the full series in dtmshow is bound to rather slow transitions. It is more convenient to bin the images to a smaller bin level. This can be done by switching on the binning radiobutton. Thus, the accessed slices will be binned on the fl when required.

Note that the coordinates of the axis are kept to the dimension of the unbinned tilt series.

Also, it is possible to create a prebinned version of the stack for different binning levels, so that the slice is accessed already in its prebinned dimension. This can be done after dtmshow has been initialized

or already when starting it (at the cost of an slower initialization)

dtmshow(ts,'bz',3)

tmshow of 3x binned tilt series

Geometric conventions

In Dynamo, we call stack 3d center to a point (x_sc, y_sc,z_sc) univocally defined by a stack of micrographs:

• x_sc and z_sc and determined by the rotation axis implicitly defined by the alignment.
• y_sc is the center of the stack along direction y, defined as floor(Ny/2)+0.5; for a tilt series with Ny pixels along y.

Creation of a binned reconstruction

Binned tilt series

We bin the original tilt series:

sbin = dpktilt.bin(ts,3;)

This command applies a binning of level two to the micorgraphs in the matrix ts separately. Note that the dimensions of the new matrix are floor( [sx,sy]/2^2<2>).

Filtered tilt series

In this walkthrough, we use Weighted Backprojection as reconstruction algorithm. We apply here the simplest possible filter onto the original projections: a ramp filter.

 sf = dpktomo.rec.applyRampFilter(sbin);

Reconstruction

As we have binned the tiltseries, a reconstruction at the full extent of the tomogram will fit easily in memory. We need to chose a size of the output tomogram, expressed in pixels of the stack By default, the center of the reconstructed volume will be located at the stack 3d center.

tiltAngles = -60:3:60;  

 outputBP = dpktilt.math.backproject([400,300,200],sf,tiltAngles);

The reconstruction can be normally accessed through the field rec in the output.

dtmshow(outputBP.rec);

tmshow on the reconstruction, z and x' views

In order to store this tomogram into a file, type:

dwrite(outputBP.rec,'temp.mrc');

Reconstruction options

The command dpktilt.math.backproject allows several modifications. The most used ones are:

• 'mw' to use several cores working in parallel.
• 'maxGb' to state the maximum size that can be kept at once in memory. If the reconstruction needs more space, it will be computed in chunks along z. Each chunk will be computed separately and immediately written into disk.

Reconstruction output object

Other than the final reconstruction, the output of the dpktilt.math.backproject command is an object that stores different information, including a description of the geometrical relationship between the reconstucted volume and the tilt series. We will use this information later in this walkthrough.

In actual practice, this reconstruction object can be saved into disk (as a binary file with extension .obp), so that the information on how a reconstruction relates to a given tilt series can be recalled in a later session. The command

outputBP.save('myInfo');

will create a file called myInfo.obp in disk, whose contents can be brought later into memory through

oRead = dread('myInfo.obp');

CTF considerations

For simplicity, we are skipping the CTF correction step in this walkthrough, where we are only interested in the geometrical manipulations. In a real case, you would just use a version of the tilt series where each micrograph has been phase-flipped according to a defocus estimation computed on the full micrograph.

Annotation on binned reconstruction

The reconstruction in bin level 3 fits in memory and can be explored easily with dtmslice

We open it creating a catalogue on the fly to contain our annotations:

dtmslice('temp.rec','c','gbcat');

Boxes-type models are the most suitable for picking in 3D

The flag 'c' orders dtmslice to prepare a catalogue arbitrarily called gbcat

Annotation on 3D

The most trivial strategy for picking particles is to use the 3D viewer. Use C to input a 3d coordinate into the model. Further controls are explained here

Select a few gold beads.

Note that you don't need to be very precise to catch the center of the gold bead.

Gold beads can be selected approximately

Annotation on projections

An alternative way to quickly annotate the gold beads is by concurrent projection. In this approach we launch a projection view of the full tomogram [Project > project full shown fragment alongz z].

Location of a single gold bead in the total projection viewed along z. A right click (two-finger tap on the Mac) marks the x-y position

Two concurrent viewers will pop up showing cuts along x and y, i.e, YZ and XZ planes where the z coordinate can be clicked (in either of them).

A left click marks the z

The (x,y,z) point will update in the tomoslice window. Also, the windows that show the cuts long x and y will update upon new right-clicks on the z projection.

Saving annotation as model

On either way you chose to perform your annotations (by 3d or by concurrent 2d projections and cut), you can save your model in the catalogue with the disk icon or under the active model menu.

Location of 3d annotations on 2d tilt series

Now we want to localize the 3d coordinates of the gold beads in the binned tilt serie. We first bring the coordinates to a memory variable that we can manipulate. 'If your dtmslice is still open, you can put the points into the workspace memory by Active Model > Export to workspace > Points, which will create in the memory a variable called temp_points with Nx3 entries. In this variable, each row are the xyz coordinates of a single gold bead.

If you already closed the dtmslice session, you need to read the model that you saved into catalogue:

dcmodels gbcat -ws oum
m = dread(oum.files{1}); 

The points are stored as Nx3 matrix in the points property of the model m

 gb = m.points;

So, we have defined in memory a matrix Nx3 called gb. Each row is the 3d position of the center on one of the manually picked gold beads, expressed in the coordinates of the 3x binned reconstruction. You can check the geometry with low level depicting commands as dslice3d

figure; dslice3d('temp.mrc','r',gb(1,:));

which plots a single slice of the tomogram cutting a point of interest (by default, in direction z and covering the full extent of the volume), followed by a depiction of the 3d point of interest. For instance, for the first gold bead (in row 1)

hold on;h=dpkgeom.plotCloud(gb(1,:),'or');

dslice3d

A different way to check it is to crop a subtomogram around the 3d coordinates of interest. For instance, for the second gold bead:

gb3d = dcrop('temp.rec',40,gb(2,:));

and show it as montage of slices with dview:

dview(gb3d);

dview of an area cropped with dview

Computation

To find the 2d coordinates of all the goldbeads at once in the micrograph corresponding to a tilt angle of -60, we would write

tilt = 1;
r2d = dpktilt.point3dInMicrograph(gb,tiltAngles(tilt),'sizeStack',size(sf),'sc',[200.5,150.5,100.5]);

Note that you need to express the geometric context of r3d.

• 'sc' stands for 'stack3dcenter'
• 'ss' stands for 'stackSize'

In other words, gb are the coordinates of the points of interest in the reconstructed volume. You need to express how the reconstructed volume related to the original stack dpktilt.point3dInMicrograph needs this information explicitly. Alternatively, you can just use the method 'point3dInAlignedMicrograph' in the obp object, assuming it is still in memory. As this object carries the geometric information, the syntax is simpler.

 r2d = outputBP.point3dInAlignedMicrograph(gb,-60);

or with the notation {} to pass position in the tilt series instead of the angle value.

 r2d = outputBP.point3dInAlignedMicrograph(r3d,{1});

Visualization

Basic Matlab commands

So, did it work? You can check it by just plotting the image of the micrograph

figure(); 
hi = dshow(sbin(:,:,tilt)); 

and then plotting on it the positions of the markers.

 hold on;
h = plot(r2d(:,1),r2d(:,2),'o');

further tunning of the graphical output is possible through access to the properties of the handle h

h.LineWidth=1;
h.MarkerSize=10;
h.Color='r'; 

Projections of tilts 1 and 21 (-60 and 0 degrees)

Note that the default convention in dshow is to plot the rows of the matrix along the 'x' axis of the plot, with the origin of coordinates in the left bottom corner of the plot. The intensity of matrix entry [i,j] is represented in the pixel centered at (i-0.5,j-0.5).

Interactive GUI

The basic tool for visualizing sets of 2d markers defined on a tilt series is dmarkers. This allows to navigate at once through all the annotations on all the images. We first compute the positions in 2d of all the 3d positions in gb:

r2dAll = outputBP.point3dInAlignedMicrograph(gb,'*');

and then create a shape set object with this information. This step is necessary to provide dmarkers with the right format.

ss     = dpktilt.aux.markgold.MarkerSet();
ss.setFrameCellArray(r2dAll); 

Now we can open the dmarkers GUI on the tilt series ts3 for this shape set:

og = dmarkers(sbin,'lm',ss);

dmarkers allow to slide along the projections. On the right, a zoom view (with the mouse wheel)

Remember that the gold beads had been picked loosely. If you zoom in, you'll see that the red spots are close to but not on the gold bead projections.

Localized reconstruction

Now we want to reconstruct locally each of the gold beads, i.e.

Reconstruct a single area

ogb = outputBP.localReconstruction(gb(1,:),[20,20,20],sf);
ogb.view(); 

Single area, binned

Reconstruction into a data folder

The command dtrec is the counterpart of dtcrop for tilt series. It loops on the positions inside a table. We will use a different approach.

outputBP.reconstructIntoDataFolder('binnedGoldBeads',20,sf,'points',gb,'fo',1);
figure;dslices binnedGoldBeads z -j * 

All gold beads, binned. Each red box represents a z projection.

Full resolution reconstruction

Now we can proceed to transport the locations annotated in the 3x binned geometry into the full size information.

fbp = outputBP.rescale('up',3);

We can apply the weight on the original tilt series at full resolution ts:

sWeightedFull = dpktomo.rec.applyRampFilter(ts);

It will take around a minute. The function accepts the mw flag to increase the speed.

fgb = fbp.localReconstruction(gb(1,:)*2^3,[60,60,60], sWeightedFull); 

The localized reconstruction can be shown through:

fgb.view();

Single gold bead, full resolution

Note that we don't really need to use the object fgb as support. Its only role is to work out the geometry, computing the position of the 3d stack center with reference to the subtomogram that we want to reconstruct locally. We can work out the maths and compute the vector: [1035.431 551.378 318.140] . We could use this vector directly in dpktilt.math.backproject and get the same result.

fbpDirect  = dpktilt.math.backproject([60,60,60], sWeightedFull,tiltAngles,'sc',[1035.431 551.378 318.140]);

We can now reconstruct all the gold beads at full resolution.

fbp.reconstructIntoDataFolder('fullGoldBeads',80, sWeightedFull,'points',gb*8,'fo',1); 

and create a depiction of the data set.

figure;dslices fullGoldBeads z -j * 

All gold beads, full resolution, compared to the reconstructions computed with the binned tilt series

Refining the full resolution reconstruction

Now, we could find the actual positions of the gold beads with an alignment project. We create a spherical template and mask:

dwrite(-dsphere(10,80),'template.mrc');
dwrite(dsphere(16,80),'mask.mrc');
dwrite(dynamo_table_blank('fullGoldBeads'),'blankFull.tbl'); 

ddelete goldAligner 
dcp.new('goldAligner','d','fullGoldBeads','t','blankFull.tbl','template','template.mrc','mask','mask.mrc','show',0);
dvput goldAligner -cr 0;
dvput goldAligner -ir 0;
dvput goldAligner -rf 0;
dvput goldAligner -lims 40;
dvput goldAligner -ite 1;
dvput goldAligner -dim 80;
check = dvcheck('goldAligner');

If the check doesn't report any problem, we can then just unfold the project:

 dvunfold goldAligner
 goldAligner 

Then we can check the result.

ddb goldAligner:a -v 

The resulting gold bead is now centered

Averaged gold bead (centered)

More importantly, the individual gold beads are now centered:

Averaged gold bead (centered)

The alignment parameters are stored in the refinement table, which we can read into a workspace variable (we arbitrarily call it t)

ddb goldAligner:rt -r t 

The position of gold beads in 3d is the addition of the computed shifts (columns 4 to 6 in the table) on the coarse positions of the gold beads:

gf3d = t(:,4:6)+gb*2^3;

r2coarse = outputBP.point3dInAlignedMicrograph(gf3d,'*');
sc = dpktilt.aux.markgold.MarkerSet();
sc.setFrameCellArray(r2coarse);
ogc = dmarkers(sbin,'lm',sc); 

New set of markers in 2d (red)

We can show in this session of dmarkers the same shape set of positions that we computed before the alignment project. We will paint it i blue for comparison with the current positions.

ss.name ='previous';
ogc.bindGenericShapeSet(ss,'pointColor','b');

Raw set of markers (blue)