Human CyCIF PDAC#

library(Giotto)

Install Python Modules#

Warning

This tutorial was written with Giotto version 0.3.6.9046, your version is 1.0.3. This is a more recent version and results should be reproducible.

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

Important

Python module installation can be done either automatically via our installation tool (from within R) (see step 2.2A) or manually (see step 2.2B).

See Part 2.2 Giotto-Specific Python Packages of our Giotto Installation section for step-by-step instructions.

Dataset Explanation#

The CyCIF data to run this tutorial can be found here. Alternatively you can use the getSpatialDataset to automatically download this dataset like we do in this example.

We will re-analyze the PDAC data generated by tissue-CyCIF as explained in the paper by Lin et al.

ADD IMAGE

Dataset Download#

The merFISH data to run this tutorial can be found here. Alternatively you can use the getSpatialDataset to automatically download this dataset like we do in this example.

# download data to working directory ####
# if wget is installed, set method = 'wget'
# if you run into authentication issues with wget, then add " extra = '--no-check-certificate' "
getSpatialDataset(dataset = 'cycif_PDAC', directory = results_folder, method = 'wget')

1. Giotto Global Instructions and Preparations#

1.1 Optional: Set Giotto Instructions#

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

1.2 Giotto Object#

# 2. create giotto object from provided paths ####
expr_path = paste0(results_folder, "cyCIF_PDAC_expression.txt.gz")
loc_path = paste0(results_folder, "cyCIF_PDAC_coord.txt")
meta_path = paste0(results_folder, "cyCIF_PDAC_annot.txt")
png_path = paste0(results_folder, "canvas.png")

2. Create Giotto Object and Process The Data#

# read in data information

# expression info
pdac_expression = readExprMatrix(expr_path, transpose = T)
# cell coordinate info
pdac_locations = data.table::fread(loc_path)
# metadata
pdac_metadata = data.table::fread(meta_path)
pdac_test <- createGiottoObject(raw_exprs = pdac_expression,
            spatial_locs = pdac_locations,
            instructions = instrs,
            cell_metadata = pdac_metadata)

## normalize & adjust
pdac_test <- normalizeGiotto(gobject = pdac_test, scalefactor = 10000, verbose = T)
pdac_test <- addStatistics(gobject = pdac_test)

## visualize original annotations ##
spatPlot(gobject = pdac_test, point_size = 0.5, coord_fix_ratio = 1, cell_color = 'frame',
 background_color = 'black', legend_symbol_size = 3,
 save_param = list(save_name = '2_a_spatPlot'))
2_a_spatPlot.png
spatPlot(gobject = pdac_test, point_size = 0.3, coord_fix_ratio = 1,
    cell_color = 'COL', background_color = 'black', legend_symbol_size = 3,
    save_param = list(save_name = '2_b_spatPlot_column'))
2_b_spatPlot_column.png
spatPlot(gobject = pdac_test, point_size = 0.3, coord_fix_ratio = 1,
    cell_color = 'ROW', background_color = 'black', legend_symbol_size = 3,
    save_param = list(save_name = '2_c_spatPlot_row'))
2_c_spatPlot_row.png
## add external histology information
pdac_metadata = pDataDT(pdac_test)
pancreas_frames = c(1:6, 27:31, 15:19, 40:44)
PDAC_frames = c(23:26, 35:37, 51:52, 64:65, 77)
small_intestines_frames = c(49:50, 63, 75:76, 88:89, 100:103, 112:116, 125:129, 137:140)

# detailed histology
hist_info = ifelse(pdac_metadata$frame %in% pancreas_frames, 'pancr',
       ifelse(pdac_metadata$frame %in% PDAC_frames, 'PDAC',
          ifelse(pdac_metadata$frame %in% small_intestines_frames, 'small_intest', 'other')))
pdac_test = addCellMetadata(pdac_test, new_metadata = hist_info)

spatPlot(gobject = pdac_test, point_size = 0.3, coord_fix_ratio = 1, cell_color = 'hist_info',
    background_color = 'black', legend_symbol_size = 3,
    save_param = list(save_name = '2_d_spatPlot_hist'))
2_d_spatPlot_hist.png
# coarse histology
hist_info2 = ifelse(pdac_metadata$frame %in% pancreas_frames, 'pancr',
        ifelse(pdac_metadata$frame %in% small_intestines_frames, 'small_intest','PDAC'))
pdac_test = addCellMetadata(pdac_test, new_metadata = hist_info2)

spatPlot(gobject = pdac_test, point_size = 0.3, coord_fix_ratio = 1, cell_color = 'hist_info2',
 background_color = 'black', legend_symbol_size = 3, point_border_stroke = 0.001,
 save_param = list(save_name = '2_e_spatPlot_hist2'))
2_e_spatPlot_hist2.png

2.1 Add and Align Image#

2.1.1 Read Image with magick#

# read
mg_img = magick::image_read(png_path)

2.1.2 Optional: Modify Image#

Examples: flip axis, negate, change background, etc.

# flip/flop (convert x and y axes)
mg_img = magick::image_flip(mg_img)
mg_img = magick::image_flop(mg_img)

# negate image
mg_img2 = magick::image_negate(mg_img)

2.1.3 Test Image#

Check to see if it is aligned

## align image ##
# 1. create spatplot
mypl = spatPlot(gobject = pdac_test, point_size = 0.3, coord_fix_ratio = NULL, cell_color = 'hist_info2',
    legend_symbol_size = 3, point_border_stroke = 0.001,
    save_plot = F, return_plot = T)

# 2.create giotto image and make adjustments (xmax_adj, xmin_adj, ...)
hist_png = createGiottoImage(gobject = pdac_test, mg_object = mg_img2, name = 'image_hist',
             xmax_adj = 5000, xmin_adj = 2500, ymax_adj = 1500, ymin_adj = 1500)

# 3. add giotto image to spatplot to check alignment
mypl_image = addGiottoImageToSpatPlot(mypl, hist_png)
mypl_image

2.1.4 Add Giotto Image(s) to Object(s)#

## add images to Giotto object ##
image_list = list(hist_png)
pdac_test = addGiottoImage(gobject = pdac_test,
               images = image_list)
showGiottoImageNames(pdac_test)

3. Dimension Reduction#

# PCA
pdac_test <- runPCA(gobject = pdac_test, expression_values = 'normalized',
        scale_unit = T, center = F, method = 'factominer')
signPCA(pdac_test, scale_unit = T, scree_ylim = c(0, 3),
        save_param = list(save_name = '3_a_signPCA'))
3_a_signPCA.png
plotPCA(gobject = pdac_test, point_shape = 'no_border', point_size = 0.2,
        save_param = list(save_name = '3_b_PCAplot'))
3_b_PCAplot.png
# UMAP
pdac_test <- runUMAP(pdac_test, dimensions_to_use = 1:14, n_components = 2, n_threads = 12)
plotUMAP(gobject = pdac_test, point_shape = 'no_border', point_size = 0.2,
    save_param = list(save_name = '3_c_UMAP'))
3_c_UMAP.png

4. Clustering#

## sNN network (default)
pdac_test <- createNearestNetwork(gobject = pdac_test, dimensions_to_use = 1:14, k = 20)

## 0.2 resolution
pdac_test <- doLeidenCluster(gobject = pdac_test, resolution = 0.2, n_iterations = 100, name = 'leiden')

# create customized color palette for leiden clustering results
pdac_metadata = pDataDT(pdac_test)
leiden_colors = Giotto:::getDistinctColors(length(unique(pdac_metadata$leiden)))
names(leiden_colors) = unique(pdac_metadata$leiden)
color_3 = leiden_colors['3'];color_10 = leiden_colors['10']
leiden_colors['3'] = color_10; leiden_colors['10'] = color_3

plotUMAP(gobject = pdac_test, cell_color = 'leiden', point_shape = 'no_border',
    point_size = 0.2, cell_color_code = leiden_colors,
    save_param = list(save_name = '4_a_UMAP'))
4_a_UMAP.png
plotUMAP(gobject = pdac_test, cell_color = 'hist_info',point_shape = 'no_border',   point_size = 0.2,
    save_param = list(save_name = '4_b_UMAP'))
4_b_UMAP.png
spatPlot(gobject = pdac_test, cell_color = 'leiden', point_shape = 'no_border', point_size = 0.2,
    cell_color_code = leiden_colors, coord_fix_ratio = 1,
    save_param = list(save_name = '4_c_spatplot'))
4_c_spatplot.png

4.1 Add Background Image#

showGiottoImageNames(pdac_test)
spatPlot(gobject = pdac_test, show_image = T, image_name = 'image_hist',
    cell_color = 'leiden',
    point_shape = 'no_border', point_size = 0.2, point_alpha = 0.7,
    cell_color_code = leiden_colors, coord_fix_ratio = 1,
    save_param = list(save_name = '4_d_spatPlot'))
4_d_spatPlot.png

5. Visualize the Spatial and Expression Space#

spatDimPlot2D(gobject = pdac_test, cell_color = 'leiden',
      spat_point_shape = 'no_border', spat_point_size = 0.2,
      dim_point_shape = 'no_border', dim_point_size = 0.2,
      cell_color_code = leiden_colors,
      save_param = list(save_name = '5_a_spatdimplot'))
5_a_spatdimplot.png
spatDimPlot2D(gobject = pdac_test, cell_color = 'leiden',
      spat_point_shape = 'border',
      spat_point_size = 0.2, spat_point_border_stroke = 0.01,
      dim_point_shape = 'border', dim_point_size = 0.2,
      dim_point_border_stroke = 0.01, cell_color_code = leiden_colors,
      save_param = list(save_name = '5_b_spatdimplot'))
5_b_spatdimplot.png
spatDimPlot2D(gobject = pdac_test, cell_color = 'hist_info2',
      spat_point_shape = 'border', spat_point_size = 0.2,
      spat_point_border_stroke = 0.01, dim_point_shape = 'border',
      dim_point_size = 0.2, dim_point_border_stroke = 0.01,
      save_param = list(save_name = '5_c_spatdimplot'))
5_c_spatdimplot.png

6. Cell-Type Marker Gene Detection#

# resolution 0.2
cluster_column = 'leiden'

# gini
markers_gini = findMarkers_one_vs_all(gobject = pdac_test,
                  method = "gini",
                  expression_values = "scaled",
                  cluster_column = cluster_column,
                  min_genes = 5)
markergenes_gini = unique(markers_gini[, head(.SD, 5), by = "cluster"][["genes"]])

plotMetaDataHeatmap(pdac_test, expression_values = "norm",
        metadata_cols = c(cluster_column),
        selected_genes = markergenes_gini,
        custom_cluster_order = c(1, 10, 3, 12, 8, 2, 9, 6, 11, 13, 4, 5, 7),
        y_text_size = 8, show_values = 'zscores_rescaled',
        save_param = list(save_name = '6_a_metaheatmap'))
6_a_metaheatmap.png
topgenes_gini = markers_gini[, head(.SD, 1), by = 'cluster']$genes
violinPlot(pdac_test, genes = unique(topgenes_gini), cluster_column = cluster_column,
    strip_text = 8, strip_position = 'right',
    save_param = c(save_name = '6_b_violinplot_gini', base_width = 5, base_height = 10))
6_b_violinplot_gini.png

7. Cell-Type Annotation#

7.1 Metadata Heatmap#

Inspect the potential cell type markers

## all genes heatmap
plotMetaDataHeatmap(pdac_test, expression_values = "norm", metadata_cols = 'leiden',
        custom_cluster_order = c(1, 10, 3, 12, 8, 2, 9, 6, 11, 13, 4, 5, 7),
        y_text_size = 8, show_values = 'zscores_rescaled',
        save_param = list(save_name = '7_a_metaheatmap'))
7_a_metaheatmap.png

Inspect the potential cell type markers stratified by tissue location

plotMetaDataHeatmap(pdac_test, expression_values = "norm", metadata_cols = c('leiden','hist_info2'),
        first_meta_col = 'leiden', second_meta_col = 'hist_info2',
        y_text_size = 8, show_values = 'zscores_rescaled',
        save_param = list(save_name = '7_b_metaheatmap'))
7_b_metaheatmap.png

7.2 Spatial Subsets#

Inspect subsets of the data based on tissue location

spatPlot(pdac_test, cell_color = 'leiden', cell_color_code = leiden_colors,
    point_shape = 'no_border', point_size = 0.75, group_by = 'hist_info2',
    save_param = list(save_name = '7_c_spatplot'))
7_c_spatplot.png
spatPlot(pdac_test, cell_color = 'leiden', cell_color_code = leiden_colors,
    point_shape = 'no_border', point_size = 0.3,
    group_by = 'hist_info2', group_by_subset = c('pancr'), cow_n_col = 1,
    save_param = list(save_name = '7_d_spatplot'))
7_d_spatplot.png
spatPlot(pdac_test, cell_color = 'leiden', cell_color_code = leiden_colors,
    point_shape = 'no_border', point_size = 0.3,
    group_by = 'hist_info2', group_by_subset = c('PDAC'), cow_n_col = 1,
    save_param = list(save_name = '7_e_spatplot'))
7_e_spatplot.png
spatPlot(pdac_test, cell_color = 'leiden', cell_color_code = leiden_colors,
    point_shape = 'no_border', point_size = 0.3,
    group_by = 'hist_info2', group_by_subset = c('small_intest'), cow_n_col = 1,
    save_param = list(save_name = '7_f_spatplot'))
7_f_spatplot.png

7.3 Spatial Distribution of Clusters#

Visually inspect the spatial distribution of different clusters

# spatial enrichment of groups
for(group in unique(pDataDT(pdac_test)$leiden)) {
    spatPlot(pdac_test, cell_color = 'leiden', point_shape = 'no_border',
        point_size = 0.3, other_point_size = 0.1,
        select_cell_groups = group, cell_color_code = 'red',
        save_param = list(save_name = paste0('7_g_spatplot_', group)))
}
7_g_spatplot_1.png 7_g_spatplot_2.png

7.4 Annotate Clusters by Position and Expression#

Annotate clusters based on spatial position and dominant expression patterns

cell_metadata =  pDataDT(pdac_test)
cluster_data = cell_metadata[, .N, by = c('leiden', 'hist_info2')]
cluster_data[, fraction:= round(N/sum(N), 2), by = c('leiden')]
data.table::setorder(cluster_data, leiden, hist_info2, fraction)

# final annotation
names = 1:13
location = c('pancr', 'intest', 'general', 'intest', 'pancr',
     'intest', 'pancr', 'canc', 'general', 'pancr',
     'general', 'pancr', 'intest')
feats = c('epithelial_I', 'fibroblast_VEGFR+', 'stroma_HER2+_pERK+', 'epithelial_lining_p21+', 'epithelial_keratin',
  'epithelial_prolif', 'epithelial_actin++', 'immune_PD-L1+', 'stromal_actin-', 'epithelial_tx_active',
  'epithelial_MET+_EGFR+', 'immune_CD45+', 'epithelial_pAKT')

annot_dt = data.table::data.table('names' = names, 'location' = location, 'feats' = feats)
annot_dt[, annotname := paste0(location,'_',feats)]
cell_annot = annot_dt$annotname;names(cell_annot) = annot_dt$names
pdac_test = annotateGiotto(pdac_test, annotation_vector = cell_annot, cluster_column = 'leiden')


# specify colors
leiden_colors
leiden_names = annot_dt$annotname; names(leiden_names) = annot_dt$names

cell_annot_colors = leiden_colors
names(cell_annot_colors) = leiden_names[names(leiden_colors)]

# covisual
spatDimPlot(gobject = pdac_test, cell_color = 'cell_types', cell_color_code = cell_annot_colors,
    spat_point_shape = 'border', spat_point_size = 0.2, spat_point_border_stroke = 0.01,
    dim_point_shape = 'border', dim_point_size = 0.2, dim_point_border_stroke = 0.01,
    dim_show_center_label = F, spat_show_legend = T, dim_show_legend = T, legend_symbol_size = 3,
    save_param = list(save_name = '7_h_spatdimplot'))
7_h_spatdimplot.png
# spatial only
spatPlot(gobject = pdac_test, cell_color = 'cell_types', point_shape = 'no_border', point_size = 0.2,
    coord_fix_ratio = 1, show_legend = T, cell_color_code = cell_annot_colors, background_color = 'black',
    save_param = list(save_name = '7_i_spatplot'))
7_i_spatplot.png
# dimension only
plotUMAP(gobject = pdac_test, cell_color = 'cell_types', point_shape = 'no_border', point_size = 0.2,
    show_legend = T, cell_color_code = cell_annot_colors, show_center_label = F, background_color = 'black',
    save_param = list(save_name = '7_j_umap'))
7_j_umap.png

1. Spatial Grid#

pdac_test <- createSpatialGrid(gobject = pdac_test,
               sdimx_stepsize = 150,
               sdimy_stepsize = 150,
               minimum_padding = 0)
spatPlot(pdac_test,
    cell_color = 'leiden',
    show_grid = T, point_size = 0.75, point_shape = 'no_border',
    grid_color = 'red', spatial_grid_name = 'spatial_grid',
    save_param = list(save_name = '8_a_spatplot'))
8_a_spatplot.png

9. Spatial Network#

pdac_test <- createSpatialNetwork(gobject = pdac_test, minimum_k = 2)

10. Cell-Cell Preferential Proximity#

## calculate frequently seen proximities
cell_proximities = cellProximityEnrichment(gobject = pdac_test,
                   cluster_column = 'cell_types',
                   spatial_network_name = 'Delaunay_network',
                   number_of_simulations = 200)
## barplot
cellProximityBarplot(gobject = pdac_test, CPscore = cell_proximities, min_orig_ints = 5, min_sim_ints = 5,
         save_param = list(save_name = '12_a_barplot'))
12_a_barplot.png
## network
cellProximityNetwork(gobject = pdac_test, CPscore = cell_proximities,
         remove_self_edges = T, only_show_enrichment_edges = F,
         save_param = list(save_name = '12_b_network'))
12_b_network.png
## visualization
spec_interaction = "1--5"
cellProximitySpatPlot2D(gobject = pdac_test, point_select_border_stroke = 0,
        interaction_name = spec_interaction,
        cluster_column = 'leiden', show_network = T,
        cell_color = 'leiden', coord_fix_ratio = NULL,
        point_size_select = 0.3, point_size_other = 0.1,
        save_param = list(save_name = '12_c_proxspatplot'))
12_c_proxspatplot.png

XX. Analyses for Paper#

XX.1 Heatmap for Cell Type Annotation#

cell_type_order_pdac = c("pancr_epithelial_actin++", "pancr_epithelial_I",
         "intest_epithelial_lining_p21+", "pancr_epithelial_keratin",
         "intest_epithelial_prolif" ,"general_epithelial_MET+_EGFR+",
         "intest_epithelial_pAKT", "pancr_epithelial_tx_active",
         "canc_immune_PD-L1+","general_stromal_actin-",
         "pancr_immune_CD45+", "intest_fibroblast_VEGFR+",
         "general_stroma_HER2+_pERK+")

plotMetaDataHeatmap(pdac_test, expression_values = "scaled", metadata_cols = c('cell_types'),
        custom_cluster_order = cell_type_order_pdac,
        y_text_size = 8, show_values = 'zscores_rescaled',
        save_param = list(save_name = 'xx_a_metaheatmap'))
xx_a_metaheatmap.png

XX.2 Pancreas#

Highlight region in pancreas:

## pancreas region ##
my_pancreas_Ids = pdac_metadata[frame == 17][['cell_ID']]
my_pancreas_giotto = subsetGiotto(pdac_test, cell_ids = my_pancreas_Ids)

spatPlot(my_pancreas_giotto, cell_color = 'leiden', point_shape = 'no_border',
    point_size = 1, cell_color_code = leiden_colors,
    save_param = list(save_name = 'xx_b_spatplot'))
xx_b_spatplot.png
spatGenePlot(my_pancreas_giotto, expression_values = 'scaled', point_border_stroke = 0.01,
     genes = c('E_Cadherin','Catenin', 'Vimentin',  'VEGFR'), point_size = 1,
     save_param = list(save_name = 'xx_c_spatgeneplot'))
xx_c_spatgeneplot.png
plotUMAP(my_pancreas_giotto, cell_color = 'leiden', point_shape = 'no_border', point_size = 0.5,
    cell_color_code = leiden_colors, show_center_label = F,
    save_param = list(save_name = 'xx_d_umap'))
xx_d_umap.png
dimGenePlot(my_pancreas_giotto, expression_values = 'scaled', point_border_stroke = 0.01,
    genes = c('E_Cadherin','Catenin', 'Vimentin',  'VEGFR'), point_size = 1,
    save_param = list(save_name = 'xx_e_dimGeneplot'))
xx_e_dimGeneplot.png

XX.3 Small Intestines#

Highlight region in small intestines (same as in original paper):

## intestine region ##
my_intest_Ids = pdac_metadata[frame == 115][['cell_ID']]
my_intest_giotto = subsetGiotto(pdac_test, cell_ids = my_intest_Ids)

spatPlot(my_intest_giotto, cell_color = 'leiden', point_shape = 'no_border',
    point_size = 1, cell_color_code = leiden_colors,
    save_param = list(save_name = 'xx_f_spatplot'))
xx_f_spatplot.png
spatGenePlot(my_intest_giotto, expression_values = 'scaled', point_border_stroke = 0.01,
     genes = c('PCNA','Catenin', 'Ki67',  'pERK'), point_size = 1,
     save_param = list(save_name = 'xx_g_spatGeneplot'))
xx_g_spatGeneplot.png
plotUMAP(my_intest_giotto, cell_color = 'leiden', point_shape = 'no_border', point_size = 0.5,
    cell_color_code = leiden_colors, show_center_label = F,
    save_param = list(save_name = 'xx_h_umap'))
xx_h_umap.png
dimGenePlot(my_intest_giotto, expression_values = 'scaled', point_border_stroke = 0.01,
    genes = c('PCNA','Catenin', 'Ki67',  'pERK'), point_size = 1,
    save_param = list(save_name = 'xx_i_dimGeneplot'))
xx_i_dimGeneplot.png