Functions

Tools for aligning spatial transcriptomics data

Functions:

LDDMM(xI, I, xJ, J[, pointsI, pointsJ, L, ...])

Run LDDMM between a pair of images.

LDDMM_3D_to_slice(xI, I, xJ, J[, pointsI, ...])

LDDMM for 3D to 2D slice mapping.

L_T_from_points(pointsI, pointsJ)

Compute an affine transformation from points.

analyze3Dalign(labelfile, xv, v, mat, xJ, ...)

Create dataframe with region annotations and 3D coordinates of target brain slice

build_transform(xv, v, A[, direction, XJ])

Create sample points to transform atlas to target from affine and velocity.

build_transform3D(xv, v, A[, direction, XJ])

Create sample points to transform atlas to target from affine and velocity.

clip(I)

clip an arrays values between 0 and 1.

download_aba_image_labels(imageurl, ...)

Create 3D altas image and region annotations

download_aba_ontology(url, file_name)

Create 3D altas ontology

extent_from_x(xJ)

Given a set of pixel locations, returns an extent 4-tuple for use with np.imshow.

interp(x, I, phii, **kwargs)

Interpolate the 2D image I, with regular grid positions stored in x (1d arrays), at the positions stored in phii (3D arrays with first channel storing component)

interp3D(x, I, phii, **kwargs)

Interpolate the 3D image I, with regular grid positions stored in x (1d arrays), at the positions stored in phii (4D arrays with first channel storing component)

make_scree(W, name[, p])

Create a scree plot for a given set of eigenvalues.

normalize(arr[, t_min, t_max])

Linearly normalizes an array between two specifed values.

plot_brain_regions(df)

Plot all brain regions in target brain slice with different color.

plot_subset_brain_regions(df, brain_regions)

Plot subset of brain regions in target brain slice with different color.

rasterize(x, y[, g, dx, blur, expand, draw, ...])

Rasterize a spatial transcriptomics dataset into a density image

rasterizePCA(x, y, G)

Rasterize a spatial transcriptomics dataset into a density image and perform PCA on it.

rasterize_with_signal(x, y[, s, dx, blur, ...])

Rasterize a spatial transcriptomics dataset into a density image

saveRasters(X, Y, W, V, Z, scree_fig, nrows, ...)

Create an RGB image principal components 1-3 and 4-6.

to_A(L, T)

Convert a linear transform matrix and a translation vector into an affine matrix.

to_A_3D(L, T)

Convert a linear transform matrix and a translation vector into an affine matrix.

transform_image_atlas_to_target(xv, v, A, xI, I)

Transform an image

transform_image_atlas_with_A(A, XI, I, XJ)

Transform an image with an affine matrix

transform_image_target_to_atlas(xv, v, A, xJ, J)

Transform an image

transform_points_atlas_to_target(xv, v, A, ...)

Transform points.

transform_points_target_to_atlas(xv, v, A, ...)

Transform points.

v_to_phii(xv, v)

Integrate a velocity field over time to return a position field (diffeomorphism).

v_to_phii_3D(xv, v)

Integrate a 3D velocity field over time to return a 3D position field (diffeomorphism).

STalign.STalign.LDDMM(xI, I, xJ, J, pointsI=None, pointsJ=None, L=None, T=None, A=None, v=None, xv=None, a=500.0, p=2.0, expand=2.0, nt=3, niter=5000, diffeo_start=0, epL=2e-08, epT=0.2, epV=2000.0, sigmaM=1.0, sigmaB=2.0, sigmaA=5.0, sigmaR=500000.0, sigmaP=20.0, device='cpu', dtype=torch.float64, muB=None, muA=None)

Run LDDMM between a pair of images.

This jointly estimates an affine transform A, and a diffeomorphism phi. The map is off the form x -> A phi x

Parameters

xIlist of torch tensor

Location of voxels in atlas image I

Itorch tensor

Atlas image I, with channels along first axis

xJlist of torch tensor

Location of voxels in target image J

Jtorch tensor

Target image J, with channels along first axis

Ltorch tensor

Initial guess for linear transform (2x2 torch tensor). Defaults to None (identity).

Ttorch tensor

Initial guess for translation (2 element torch tensor). Defaults to None (identity)

Atorch tensor

Initial guess for affine matrix. Either L and T can be specified, or A, but not both. Defaults to None (identity).

vtorch tensor

Initial guess for velocity field

xvtorch tensor

pixel locations for velocity field

afloat

Smoothness scale of velocity field (default 500.0)

pfloat

Power of Laplacian in velocity regularization (default 2.0)

expandfloat

Factor to expand size of velocity field around image boundaries (default 2.0)

ntint

Number of timesteps for integrating velocity field (default 3). Ignored if you input v.

pointsItorch tensor

N x 2 set of corresponding points for matching in atlas image. Default None (no points).

pointsJtorch tensor

N x 2 set of corresponding points for matching in target image. Default None (no points).

niterint

Number of iterations of gradient descent optimization

diffeo_startint

Number of iterations of gradient descent optimization for affine only, before nonlinear deformation.

epLfloat

Gradient descent step size for linear part of affine.

epTfloat

Gradient descent step size of translation part of affine.

epVfloat

Gradient descent step size for velocity field.

sigmaMfloat

Standard deviation of image matching term for Gaussian mixture modeling in cost function. This term generally controls matching accuracy with smaller corresponding to more accurate. As an common example (rule of thumb), you could chose this parameter to be the variance of the pixels in your target image.

sigmaBfloat

Standard deviation of backtround term for Gaussian mixture modeling in cost function. If there is missing tissue in target, we may label some pixels in target as background, and not enforce matching here.

sigmaAfloat

Standard deviation of artifact term for Gaussian mixture modeling in cost function. If there are artifacts in target or other lack of corresponding between template and target, we may label some pixels in target as artifact, and not enforce matching here.

sigmaR: float

Standard deviation for regularization. Smaller sigmaR means a smoother resulting transformation. Regularization is of the form: 0.5/sigmaR^2 int_0^1 int_X |Lv|^2 dx dt.

sigmaP: float

Standard deviation for matching of points. Cost is of the form 0.5/sigmaP^2 sum_i (atlas_point_i - target_point_i)^2

device: str

torch device. defaults to ‘cpu’. Can also be ‘cuda:0’ for example.

dtype: torch dtype

torch data type. defaults to torch.float64

muA: torch tensor whose dimension is the same as the target image

Defaults to None, which means we estimate this. If you provide a value, we will not estimate it. If the target is a RGB image, this should be a tensor of size 3. If the target is a grayscale image, this should be a tensor of size 1.

muB: torch tensor whose dimension is the same as the target image

Defaults to None, which means we estimate this. If you provide a value, we will not estimate it.

Returns a dictionary

{ ‘A’: torch tensor

Affine transform

‘v’: torch tensor

Velocity field

‘xv’: list of torch tensor

Pixel locations in v

‘WM’: torch tensor

Resulting weight 2D array (matching)

‘WB’: torch tensor

Resulting weight 2D array (background)

‘WA’: torch tensor

Resulting weight 2D array (artifact)

}

STalign.STalign.LDDMM_3D_to_slice(xI, I, xJ, J, pointsI=None, pointsJ=None, L=None, T=None, A=None, v=None, xv=None, a=500.0, p=2.0, expand=1.25, nt=3, niter=5000, diffeo_start=0, epL=1e-06, epT=10.0, epV=1000.0, sigmaM=1.0, sigmaB=2.0, sigmaA=5.0, sigmaR=100000000.0, sigmaP=20.0, device='cpu', dtype=torch.float64, muA=None, muB=None)

LDDMM for 3D to 2D slice mapping.

muA: torch tensor whose dimension is the same as the target image

Defaults to None, which means we estimate this. If you provide a value, we will not estimate it. If the target is a RGB image, this should be a tensor of size 3. If the target is a grayscale image, this should be a tensor of size 1.

muB: torch tensor whose dimension is the same as the target image

Defaults to None, which means we estimate this. If you provide a value, we will not estimate it.

STalign.STalign.L_T_from_points(pointsI, pointsJ)

Compute an affine transformation from points.

Note for an affine transformation (6dof) we need 3 points.

Outputs, L,T should be rconstructed blockwize like [L,T;0,0,1]

Parameters

pointsIarray

An Nx2 array of floating point numbers describing atlas points in ROW COL order (not xy)

pointsJarray

An Nx2 array of floating point numbers describing target points in ROW COL order (not xy)

Returns

Larray

A 2x2 linear transform array.

Tarray

A 2 element translation vector

STalign.STalign.analyze3Dalign(labelfile, xv, v, mat, xJ, dx, scale_x, scale_y, x, y, X_, Y_, namesdict, device='cpu')

Create dataframe with region annotations and 3D coordinates of target brain slice

labelfilearray

File name of 3D atlas with region annotations.

xvlist of array

Sample points for velocity

varray

time dependent velocity field

Aarray

Affine transformation matrix

directionchar

‘f’ for forward and ‘b’ for backward. ‘b’ is default and is used for transforming images. ‘f’ is used for transforming points.

XJarray

Sample points for target (meshgrid with ij index style). Defaults to None

to keep sampling on the xv.

dxint

Step of rasterized image.

scale_xdouble

Value that x positions of atlas were scaled by

scale_ydouble

Value that y positions of atlas were scaled by

xarray

X positions of target brain slice

yarray

Y positions of target brain slice

X_array
Y_array
namesdictDictionary

Dictionary of brain regions corresponding to region numbers

df : Dataframe containing each cell in original target brain slice, region annotations and 3D coordinates in terms of the altas.

STalign.STalign.build_transform(xv, v, A, direction='b', XJ=None)

Create sample points to transform atlas to target from affine and velocity.

Parameters

xvlist of array

Sample points for velocity

varray

time dependent velocity field

Aarray

Affine transformation matrix

directionchar

‘f’ for forward and ‘b’ for backward. ‘b’ is default and is used for transforming images. ‘f’ is used for transforming points.

XJarray

Sample points for target (meshgrid with ij index style). Defaults to None to keep sampling on the xv.

Returns

Xsarray

Sample points in mehsgrid format.

STalign.STalign.build_transform3D(xv, v, A, direction='b', XJ=None)

Create sample points to transform atlas to target from affine and velocity.

Parameters

xvlist of array

Sample points for velocity

varray

time dependent velocity field

Aarray

Affine transformation matrix

directionchar

‘f’ for forward and ‘b’ for backward. ‘b’ is default and is used for transforming images. ‘f’ is used for transforming points.

XJarray

Sample points for target (meshgrid with ij index style). Defaults to None to keep sampling on the xv.

Returns

Xsarray

Sample points in mehsgrid format.

STalign.STalign.clip(I)

clip an arrays values between 0 and 1. Useful for visualizing

Parameters

Itorch tensor

A torch tensor to clip.

Returns

Ictorch tensor

Clipped torch tensor.

STalign.STalign.download_aba_image_labels(imageurl, labelurl, imagefile, labelfile)

Create 3D altas image and region annotations

Parameters

imageurlarray

Link containing atlas cell stained image.

labelurlarray

Link containing atlas altas images for each voxel in imageurl.

imagefilearray

File location to save imageurl information.

labelfilearray

File location to save labelurl information.

Returns

imagefilearray

File location to save imageurl information.

labelfilearray

File location to save labelurl information.

STalign.STalign.download_aba_ontology(url, file_name)

Create 3D altas ontology

Parameters

url : Link to url contain atlas ontology file_name : File name to save atlas ontology

Returns

ontology_name : File name to store ontologies namesdict : Dictionary of all brain rgeion names

STalign.STalign.extent_from_x(xJ)

Given a set of pixel locations, returns an extent 4-tuple for use with np.imshow.

Note inputs are locations of pixels along each axis, i.e. row column not xy.

Parameters

xJlist of torch tensors

Location of pixels along each axis

Returns

extenttuple

(xmin, xmax, ymin, ymax) tuple

Examples

>>> extent_from_x(xJ)
>>> fig,ax = plt.subplots()
>>> ax.imshow(J,extent=extentJ)
STalign.STalign.interp(x, I, phii, **kwargs)

Interpolate the 2D image I, with regular grid positions stored in x (1d arrays), at the positions stored in phii (3D arrays with first channel storing component)

Parameters

xlist of arrays

List of arrays storing the pixel locations along each image axis. convention is row column order not xy.

Iarray

Image array. First axis should contain different channels.

phiiarray

Sampling array. First axis should contain sample locations corresponding to each axis.

**kwargsdict

Other arguments fed into the torch interpolation function torch.nn.grid_sample

Returns

outtorch tensor

The image I resampled on the points defined in phii.

Notes

Convention is to use align_corners=True.

This uses the torch library.

STalign.STalign.interp3D(x, I, phii, **kwargs)

Interpolate the 3D image I, with regular grid positions stored in x (1d arrays), at the positions stored in phii (4D arrays with first channel storing component)

Parameters

xlist of arrays

List of arrays storing the pixel locations along each image axis. convention is row column order not xy.

Iarray

Image array. First axis should contain different channels.

phiiarray

Sampling array. First axis should contain sample locations corresponding to each axis.

**kwargsdict

Other arguments fed into the torch interpolation function torch.nn.grid_sample

Returns

outtorch tensor

The image I resampled on the points defined in phii.

Notes

Convention is to use align_corners=True.

This uses the torch library.

STalign.STalign.make_scree(W, name, p=6)

Create a scree plot for a given set of eigenvalues.

Parameters

Wnumpy array

Eigenvalues in descending order

namestr

Name of the dataset the eigenvalue originate from Will be included in plot title

pint

Number of eigenvalue from W to plot. Defaults to 6

Returns

figmatplotlib Figure

Figure handle for scree plot

Raises

ValueError

If p is larger than the number to eigenvalues in W.

STalign.STalign.normalize(arr, t_min=0, t_max=1)

Linearly normalizes an array between two specifed values.

Parameters

arrnumpy array

array to be normalized

t_minint or float

Lower bound of normalization range

t_maxint or float

Upper bound of normalization range

Returns

norm_arrnumpy array

1D array with normalized arr values

STalign.STalign.plot_brain_regions(df)

Plot all brain regions in target brain slice with different color.

Parameters

dfDataFrame

Dataframe containing each cell in original target brain slice, region annotations and 3D coordinates in terms of the altas.

Returns

None

STalign.STalign.plot_subset_brain_regions(df, brain_regions)

Plot subset of brain regions in target brain slice with different color.

Parameters

dfDataFrame

Dataframe containing each cell in original target brain slice, region annotations and 3D coordinates in terms of the altas.

brain_regionsarray

Subset of brain regions of interest (i.e. [‘CA1’, ‘CP’, ‘DG-sg’])

Returns

None

STalign.STalign.rasterize(x, y, g=array([1.]), dx=30.0, blur=1.0, expand=1.1, draw=10000, wavelet_magnitude=False, use_windowing=True)

Rasterize a spatial transcriptomics dataset into a density image

Paramters

xnumpy array of length N

x location of cells

ynumpy array of length N

y location of cells

gnumpy array of length N

RNA count of cells If not given, density image is created

dxfloat

Pixel size to rasterize data (default 30.0, in same units as x and y)

blurfloat or list of floats

Standard deviation of Gaussian interpolation kernel. Units are in number of pixels. Can be aUse a list to do multi scale.

expandfloat

Factor to expand sampled area beyond cells. Defaults to 1.1.

drawint

If True, draw a figure every draw points return its handle. Defaults to False (0).

wavelet_magnitudebool

If True, take the absolute value of difference between scales for raster images. When using this option blur should be sorted from greatest to least.

Returns

Xnumpy array

Locations of pixels along the x axis

Ynumpy array

Locations of pixels along the y axis

Mnumpy array

A rasterized image with len(blur) channels along the first axis

figmatplotlib figure handle

If draw=True, returns a figure handle to the drawn figure.

Raises

Exception

If wavelet_magnitude is set to true but blur is not sorted from greatest to least.

Examples

Rasterize a dataset at 30 micron pixel size, with three kernels.

>>> X,Y,M,fig = tools.rasterize(x,y,dx=30.0,blur=[2.0,1.0,0.5],draw=10000)

Rasterize a dataset at 30 micron pixel size, with three kernels, using difference between scales.

>>> X,Y,M,fig = tools.rasterize(x,y,dx=30.0,blur=[2.0,1.0,0.5],draw=10000, wavelet_magnitude=True)
STalign.STalign.rasterizePCA(x, y, G)

Rasterize a spatial transcriptomics dataset into a density image and perform PCA on it.

Parameters

xnumpy array of length N

x location of cells

ynumpy array of length N

y location of cells

Gpandas Dataframe of shape (N,M)

gene expression level of cells for M genes

Returns

Xnumpy array

A rasterized image for each gene with the channel along the first axis

Ynumpy array

A rasterized image unraveled to 1-D for each gene; data is centered

Wnumpy array

The eigenvalues for the principal components in descending order

Vnumpy array

The normalized eigenvectors for the principal components V[:, i] corresponds to eigenvalue at W[i]

Znumpy array

X rotated to align with the principal component axes

nrowsint

Row dimension of rasterized image

ncolsint

Column dimension of rasterized image

Notes

Each value/row at the same index in x, y, and G should all correspond to the same cell. x[i] <-> y[i] <-> G[i,:]

STalign.STalign.rasterize_with_signal(x, y, s=None, dx=30.0, blur=1.0, expand=1.1, draw=0, wavelet_magnitude=False, use_windowing=True)

Rasterize a spatial transcriptomics dataset into a density image

Paramters

xnumpy array of length N

x location of cells

ynumpy array of length N

y location of cells

snumpy array of length N by M for M signals

signal value should be length NxM

dxfloat

Pixel size to rasterize data (default 30.0, in same units as x and y)

blurfloat or list of floats

Standard deviation of Gaussian interpolation kernel. Units are in number of pixels. Can be aUse a list to do multi scale.

expandfloat

Factor to expand sampled area beyond cells. Defaults to 1.1.

drawint

If True, draw a figure every draw points return its handle. Defaults to False (0).

wavelet_magnitudebool

If True, take the absolute value of difference between scales for raster images. When using this option blur should be sorted from greatest to least.

Returns

Xnumpy array

Locations of pixels along the x axis

Ynumpy array

Locations of pixels along the y axis

Mnumpy array

A rasterized image with len(blur) channels along the last axis

figmatplotlib figure handle

If draw=True, returns a figure handle to the drawn figure.

Raises

Exception

If wavelet_magnitude is set to true but blur is not sorted from greatest to least.

Examples

Rasterize a dataset at 30 micron pixel size, with three kernels.

>>> X,Y,M,fig = tools.rasterize(x,y,dx=30.0,blur=[2.0,1.0,0.5],draw=10000)

Rasterize a dataset at 30 micron pixel size, with three kernels, using difference between scales.

>>> X,Y,M,fig = tools.rasterize(x,y,dx=30.0,blur=[2.0,1.0,0.5],draw=10000, wavelet_magnitude=True)
STalign.STalign.saveRasters(X, Y, W, V, Z, scree_fig, nrows, ncols, name, path)

Create an RGB image principal components 1-3 and 4-6. Save numpy arrays necessary for rasterizing and performing PCA. Save scree plot and RGB image as jpeg, and their figure handles.

Parameters

Xnumpy array

A rasterized image for each gene with the channel along the first axis

Ynumpy array

A rasterized image unraveled to 1-D for each gene; data is centered

Wnumpy array

The eigenvalues for the principal components in descending order

Vnumpy array

The normalized eigenvectors for the principal components V[:, i] corresponds to eigenvalue at W[i]

Znumpy array

X rotated to align with the principal component axes

scree_figmatplotlib Figure

Figure handle for scree plot

nrowsint

Row dimension of rasterized image

ncolsint

Column dimension of rasterized image

namestr

Name of the dataset the arrays originate from

pathstr

Absolute path to folder where files should be saved

Returns

None

STalign.STalign.to_A(L, T)

Convert a linear transform matrix and a translation vector into an affine matrix.

Parameters

Ltorch tensor

2x2 linear transform matrix

Ttorch tensor

2 element translation vector (note NOT 2x1)

Returns

Atorch tensor

Affine transform matrix

STalign.STalign.to_A_3D(L, T)

Convert a linear transform matrix and a translation vector into an affine matrix.

Parameters

Ltorch tensor

3x3 linear transform matrix

Ttorch tensor

3 element translation vector (note NOT 2x1)

Returns

Atorch tensor

Affine transform matrix

STalign.STalign.transform_image_atlas_to_target(xv, v, A, xI, I, XJ=None)

Transform an image

STalign.STalign.transform_image_atlas_with_A(A, XI, I, XJ)

Transform an image with an affine matrix

Parameters

Atorch tensor

Affine transform matrix

XIlist of numpy arrays

List of arrays storing the pixel location in image I along each image axis. convention is row column order not xy. i.e, locations of pixels along the y axis (rows) followed by locations of pixels along the x axis (columns)

Inumpy array

A rasterized image with len(blur) channels along the first axis

XJlist of numpy arrays

List of arrays storing the pixel location in image I along each image axis. convention is row column order not xy. i.e, locations of pixels along the y axis (rows) followed by locations of pixels along the x axis (columns)

Returns

AItorch tensor

image I after affine transformation A, with channels along first axis

STalign.STalign.transform_image_target_to_atlas(xv, v, A, xJ, J, XI=None)

Transform an image

STalign.STalign.transform_points_atlas_to_target(xv, v, A, pointsI)

Transform points. Note points are in row column order, not xy.

STalign.STalign.transform_points_target_to_atlas(xv, v, A, pointsI)

Transform points. Note points are in row column order, not xy.

STalign.STalign.v_to_phii(xv, v)

Integrate a velocity field over time to return a position field (diffeomorphism).

Parameters

xvlist of torch tensor

List of 1D tensors describing locations of sample points in v

vtorch tensor

5D (nt,2,v0,v1) velocity field

Returns

phii: torch tensor

Inverse map (position field) computed by method of characteristics

STalign.STalign.v_to_phii_3D(xv, v)

Integrate a 3D velocity field over time to return a 3D position field (diffeomorphism).

Parameters

xvlist of torch tensor

List of 1D tensors describing locations of sample points in v

vtorch tensor

5D (nt,3,v0,v1,v2) velocity field

Returns

phii: torch tensor

Inverse map (position field) computed by method of characteristics