mini_visium

Date:

2022-09-16

library(Giotto)

Install Python modules

To run this vignette you need to install all the necessary Python modules.

1. This can be done manually, see https://rubd.github.io/Giotto_site/articles/installation_issues.html#python-manual-installation

2. This can be done within R using our installation tools (installGiottoEnvironment), see https://rubd.github.io/Giotto_site/articles/tut0_giotto_environment.html for more information.

(optional) set giotto instructions

# to automatically save figures in save_dir set save_plot to TRUE
temp_dir = '~/Temp/'
myinstructions = createGiottoInstructions(save_dir = temp_dir,
                                          save_plot = FALSE,
                                          show_plot = F)

1. A. Create a Giotto object

minimum requirements:

- matrix with expression information (or path to)

- x,y(,z) coordinates for cells or spots (or path to)

# giotto object
expr_path = system.file("extdata", "visium_DG_expr.txt.gz", package = 'Giotto')
loc_path = system.file("extdata", "visium_DG_locs.txt", package = 'Giotto')
mini_visium <- createGiottoObject(raw_exprs = expr_path,
                                   spatial_locs = loc_path,
                                   instructions = myinstructions)

How to work with Giotto instructions that are part of your Giotto object:

- show the instructions associated with your Giotto object with showGiottoInstructions

- change one or more instructions with changeGiottoInstructions

- replace all instructions at once with replaceGiottoInstructions

- read or get a specific giotto instruction with readGiottoInstructions

Of note, the python path can only be set once in an R session. See the reticulate package for more information.

# show instructions associated with giotto object (mini_visium)
showGiottoInstructions(mini_visium)

1. B. Add image

## 1. read image
png_path = system.file("extdata", "deg_image.png", package = 'Giotto')
mg_img = magick::image_read(png_path)

## 2. test and modify image alignment
mypl = spatPlot(mini_visium, return_plot = T, point_alpha = 0.8)
orig_png = createGiottoImage(gobject = mini_visium, mg_object = mg_img, name = 'image',
                             xmax_adj = 450, xmin_adj = 550,
                             ymax_adj = 200, ymin_adj = 200)
mypl_image = addGiottoImageToSpatPlot(mypl, orig_png)
mypl_image

## 3. add images to Giotto object ##
image_list = list(orig_png)
mini_visium = addGiottoImage(gobject = mini_visium,
                             images = image_list)
showGiottoImageNames(mini_visium)

2. processing steps

• filter genes and cells based on detection frequencies

• normalize expression matrix (log transformation, scaling factor and/or z-scores)

• add cell and gene statistics (optional)

• adjust expression matrix for technical covariates or batches (optional). These results will be stored in the custom slot.

# explore gene and cell distribution
filterDistributions(mini_visium, detection = 'feats')
filterDistributions(mini_visium, detection = 'cells')
filterCombinations(mini_visium,
                   expression_thresholds = c(1),
                   gene_det_in_min_cells = c(20, 20, 50, 50),
                   min_det_genes_per_cell = c(100, 200, 100, 200))

# filter and normalize
mini_visium <- filterGiotto(gobject = mini_visium,
                            expression_threshold = 1,
                            gene_det_in_min_cells = 50,
                            min_det_genes_per_cell = 100,
                            expression_values = c('raw'),
                            verbose = T)
mini_visium <- normalizeGiotto(gobject = mini_visium, scalefactor = 6000, verbose = T)
mini_visium <- addStatistics(gobject = mini_visium)

3. dimension reduction

• identify highly variable genes (HVG)

• perform PCA

• identify number of significant prinicipal components (PCs)

• run UMAP and/or TSNE on PCs (or directly on matrix)

mini_visium <- calculateHVG(gobject = mini_visium)

mini_visium <- runPCA(gobject = mini_visium)
screePlot(mini_visium, ncp = 30)
plotPCA(gobject = mini_visium)

mini_visium <- runUMAP(mini_visium, dimensions_to_use = 1:10)
plotUMAP(gobject = mini_visium)
mini_visium <- runtSNE(mini_visium, dimensions_to_use = 1:10)
plotTSNE(gobject = mini_visium)

4. clustering

• create a shared (default) nearest network in PCA space (or directly on matrix)

• cluster on nearest network with Leiden or Louvan (kmeans and hclust are alternatives)

mini_visium <- createNearestNetwork(gobject = mini_visium, dimensions_to_use = 1:10, k = 15)
mini_visium <- doLeidenCluster(gobject = mini_visium, resolution = 0.4, n_iterations = 1000)

pDataDT(mini_visium)

# visualize UMAP cluster results
plotUMAP(gobject = mini_visium, cell_color = 'leiden_clus', show_NN_network = T, point_size = 2.5)

# visualize UMAP and spatial results
spatDimPlot(gobject = mini_visium, cell_color = 'leiden_clus', spat_point_shape = 'voronoi')

# heatmap and dendrogram
showClusterHeatmap(gobject = mini_visium, cluster_column = 'leiden_clus')
showClusterDendrogram(mini_visium, h = 0.5, rotate = T, cluster_column = 'leiden_clus')

5. differential expression

scran_markers = findMarkers_one_vs_all(gobject = mini_visium,
                                       method = 'scran',
                                       expression_values = 'normalized',
                                       cluster_column = 'leiden_clus')
# violinplot
topgenes_scran = scran_markers[, head(.SD, 2), by = 'cluster']$genes
violinPlot(mini_visium, genes = topgenes_scran, cluster_column = 'leiden_clus',
           strip_text = 10, strip_position = 'right')

# metadata heatmap
topgenes_scran = scran_markers[, head(.SD, 6), by = 'cluster']$genes
plotMetaDataHeatmap(mini_visium, selected_genes = topgenes_scran,
                    metadata_cols = c('leiden_clus'))

6. A. simple cell type annotation

clusters_cell_types = c('Ambiguous_1', 'Tcf7l2_cells', 'Ambiguous_2', 'Ipcef1_cells',
                        'Prdm8_cells', 'Spink8_cells', 'Ccn2_cells', 'Cdc153_cells')
names(clusters_cell_types) = 1:8
mini_visium = annotateGiotto(gobject = mini_visium, annotation_vector = clusters_cell_types,
                             cluster_column = 'leiden_clus', name = 'cell_types')

# check new cell metadata
pDataDT(mini_visium)

# visualize annotations
spatDimPlot(gobject = mini_visium, cell_color = 'cell_types', spat_point_size = 3, dim_point_size = 3)

6. B. simple cell type enrichment

7.2 spatial enrichment: cell type distribution

Here we will use known markers for different mouse brain cell types to identify which cell types are enriched in the individual spots or identified clusters.

Paper: eisel, A. et al. Molecular Architecture of the Mouse Nervous System. Cell 174, 999-1014.e22 (2018).

## cell type signatures ##
## combination of all marker genes identified in Zeisel et al
sign_matrix_path = system.file("extdata", "sig_matrix.txt", package = 'Giotto')
brain_sc_markers = data.table::fread(sign_matrix_path) # file don't exist in data folder
sig_matrix = as.matrix(brain_sc_markers[,-1]); rownames(sig_matrix) = brain_sc_markers$Event

## enrichment tests
mini_visium = runSpatialEnrich(mini_visium,
                               sign_matrix = sig_matrix,
                               enrich_method = 'PAGE') #default = 'PAGE'

## heatmap of enrichment versus annotation (e.g. clustering result)
cell_types = colnames(sig_matrix)
plotMetaDataCellsHeatmap(gobject = mini_visium,
                         metadata_cols = 'leiden_clus',
                         value_cols = cell_types,
                         spat_enr_names = 'PAGE',
                         x_text_size = 8, y_text_size = 8)


enrichment_results = mini_visium@spatial_enrichment$PAGE
enrich_cell_types = colnames(enrichment_results)
enrich_cell_types = enrich_cell_types[enrich_cell_types != 'cell_ID']

## spatplot
spatCellPlot(gobject = mini_visium, spat_enr_names = 'PAGE',
             cell_annotation_values = enrich_cell_types,
             cow_n_col = 3,coord_fix_ratio = NULL, point_size = 1)

7. spatial grid

Create a grid based on defined stepsizes in the x,y(,z) axes.

mini_visium <- createSpatialGrid(gobject = mini_visium,
                                 sdimx_stepsize = 300,
                                 sdimy_stepsize = 300,
                                 minimum_padding = 50)
showGrids(mini_visium)
spatPlot(gobject = mini_visium, show_grid = T, point_size = 1.5)

# extract grid and associated metadata spots
annotated_grid = annotateSpatialGrid(mini_visium)
annotated_grid_metadata = annotateSpatialGrid(mini_visium,
                                              cluster_columns = c('leiden_clus', 'cell_types', 'nr_genes'))

8. spatial network

• visualize information about the default Delaunay network

• create a spatial Delaunay network (default)

• create a spatial kNN network

plotStatDelaunayNetwork(gobject = mini_visium, maximum_distance = 300)
mini_visium = createSpatialNetwork(gobject = mini_visium, minimum_k = 2, maximum_distance_delaunay = 400)
mini_visium = createSpatialNetwork(gobject = mini_visium, minimum_k = 2, method = 'kNN', k = 10)
showNetworks(mini_visium)

# visualize the two different spatial networks
spatPlot(gobject = mini_visium, show_network = T,
         network_color = 'blue', spatial_network_name = 'Delaunay_network',
         point_size = 2.5, cell_color = 'leiden_clus')

spatPlot(gobject = mini_visium, show_network = T,
         network_color = 'blue', spatial_network_name = 'kNN_network',
         point_size = 2.5, cell_color = 'leiden_clus')

9. spatial genes

Identify spatial genes with 3 different methods:

- binSpect with kmeans binarization (default)

- binSpect with rank binarization

- silhouetteRank

Visualize top 4 genes per method.

km_spatialgenes = binSpect(mini_visium)
spatGenePlot(mini_visium, expression_values = 'scaled',
             genes = km_spatialgenes[1:4]$feats,
             point_shape = 'border', point_border_stroke = 0.1,
             show_network = F, network_color = 'lightgrey', point_size = 2.5,
             cow_n_col = 2)

rank_spatialgenes = binSpect(mini_visium, bin_method = 'rank')
spatGenePlot(mini_visium, expression_values = 'scaled',
             genes = rank_spatialgenes[1:4]$feats,
             point_shape = 'border', point_border_stroke = 0.1,
             show_network = F, network_color = 'lightgrey', point_size = 2.5,
             cow_n_col = 2)

silh_spatialgenes = silhouetteRank(gobject = mini_visium) # TODO: suppress print output
spatGenePlot(mini_visium, expression_values = 'scaled',
             genes = silh_spatialgenes[1:4]$genes,
             point_shape = 'border', point_border_stroke = 0.1,
             show_network = F, network_color = 'lightgrey', point_size = 2.5,
             cow_n_col = 2)

10. spatial co-expression patterns

Identify robust spatial co-expression patterns using the spatial network or grid and a subset of individual spatial genes.

1. calculate spatial correlation scores

2. cluster correlation scores

# 1. calculate spatial correlation scores
ext_spatial_genes = km_spatialgenes[1:100]$feats
spat_cor_netw_DT = detectSpatialCorGenes(mini_visium,
                                         method = 'network', spatial_network_name = 'Delaunay_network',
                                         subset_genes = ext_spatial_genes)

# 2. cluster correlation scores
spat_cor_netw_DT = clusterSpatialCorGenes(spat_cor_netw_DT, name = 'spat_netw_clus', k = 8)
heatmSpatialCorGenes(mini_visium, spatCorObject = spat_cor_netw_DT, use_clus_name = 'spat_netw_clus')


netw_ranks = rankSpatialCorGroups(mini_visium, spatCorObject = spat_cor_netw_DT, use_clus_name = 'spat_netw_clus')
top_netw_spat_cluster = showSpatialCorGenes(spat_cor_netw_DT, use_clus_name = 'spat_netw_clus',
                                            selected_clusters = 6, show_top_genes = 1)

cluster_genes_DT = showSpatialCorGenes(spat_cor_netw_DT, use_clus_name = 'spat_netw_clus', show_top_genes = 1)
cluster_genes = cluster_genes_DT$clus; names(cluster_genes) = cluster_genes_DT$feat_ID

mini_visium = createMetagenes(mini_visium, gene_clusters = cluster_genes, name = 'cluster_metagene')
spatCellPlot(mini_visium,
             spat_enr_names = 'cluster_metagene',
             cell_annotation_values = netw_ranks$clusters,
             point_size = 1.5, cow_n_col = 3)

11. spatial HMRF domains

hmrf_folder = paste0(temp_dir,'/','11_HMRF/')
if(!file.exists(hmrf_folder)) dir.create(hmrf_folder, recursive = T)

# perform hmrf
my_spatial_genes = km_spatialgenes[1:100]$genes
HMRF_spatial_genes = doHMRF(gobject = mini_visium,
                            expression_values = 'scaled',
                            spatial_genes = my_spatial_genes,
                            spatial_network_name = 'Delaunay_network',
                            k = 8,
                            betas = c(28,2,2),
                            output_folder = paste0(hmrf_folder, '/', 'Spatial_genes_brain/SG_top100_k8_scaled'))

# check and select hmrf
for(i in seq(28, 30, by = 2)) {
  viewHMRFresults2D(gobject = mini_visium,
                    HMRFoutput = HMRF_spatial_genes,
                    k = 8, betas_to_view = i,
                    point_size = 2)
}

mini_visium = addHMRF(gobject = mini_visium,
                      HMRFoutput = HMRF_spatial_genes,
                      k = 8, betas_to_add = c(28),
                      hmrf_name = 'HMRF')

giotto_colors = getDistinctColors(8)
names(giotto_colors) = 1:8
spatPlot(gobject = mini_visium, cell_color = 'HMRF_k8_b.28',
         point_size = 3, coord_fix_ratio = 1, cell_color_code = giotto_colors)

12. export Giotto Analyzer to Viewer

viewer_folder = paste0(temp_dir, '/', 'Mouse_cortex_viewer')

# select annotations, reductions and expression values to view in Giotto Viewer
exportGiottoViewer(gobject = mini_visium, output_directory = viewer_folder,
                   factor_annotations = c('cell_types',
                                          'leiden_clus',
                                          'HMRF_k8_b.28'),
                   numeric_annotations = 'total_expr',
                   dim_reductions = c('umap'),
                   dim_reduction_names = c('umap'),
                   expression_values = 'scaled',
                   expression_rounding = 3,
                   overwrite_dir = T)