Functions
Tools for aligning spatial transcriptomics data
Functions:

Run LDDMM between a pair of images. 

LDDMM for 3D to 2D slice mapping. 

Compute an affine transformation from points. 

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

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

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

clip an arrays values between 0 and 1. 

Create 3D altas image and region annotations 

Create 3D altas ontology 

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

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) 

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) 

Create a scree plot for a given set of eigenvalues. 

Linearly normalizes an array between two specifed values. 
Plot all brain regions in target brain slice with different color. 


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

Rasterize a spatial transcriptomics dataset into a density image 

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

Rasterize a spatial transcriptomics dataset into a density image 

Create an RGB image principal components 13 and 46. 

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

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

Transform an image 

Transform an image with an affine matrix 

Transform an image 

Transform points. 

Transform points. 

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

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=2e08, 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=1e06, 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 4tuple 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’, ‘DGsg’])
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 1D 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 13 and 46. 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 1D 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