Human Heart

Spatial Gene Expression Dataset by Space Ranger 1.1.0

The library (T1T2-G12) was prepared following the Visium Spatial Gene Expression Reagent Kits User Guide (CG000239) and was sequenced on an Illumina NovaSeq 6000. Detailed info of this dataset from 10x can be found at here.


Load required packages.

library(data.table)
library(tidyverse)
library(magrittr)
library(patchwork)
library(extrafont)
library(Matrix)
Sys.time()
## [1] "2021-02-03 20:43:51 CST"


Data preparation


Functions

source(
    file = file.path(
        SCRIPT_DIR,
        "utilities.R"
    )
)

Matrix

PROJECT_DIR <- "/Users/jialei/Dropbox/Data/Projects/UTSW/Spatial_sequencing/datasets/Human_heart"
matrix_readcount_use <- read_10x_matrix_h5(
    h5_file = file.path(
        PROJECT_DIR,
        "raw",
        "V1_Human_Heart_filtered_feature_bc_matrix.h5"
    ),
    cell_id_prefix = NULL,
    feature_type = "Gene Expression"
)

Embedding

Cluster embedding

embedding <- readr::read_csv(
    file = file.path(
        PROJECT_DIR,
        "clustering",
        "exploring",
        "embedding_ncomponents20_ccc1_seed20200416.csv.gz"
    )
)
## 
## ── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## cols(
##   cell = col_character(),
##   batch = col_double(),
##   louvain = col_double(),
##   leiden = col_double(),
##   x_tsne = col_double(),
##   y_tsne = col_double(),
##   x_umap = col_double(),
##   y_umap = col_double(),
##   x_fitsne = col_double(),
##   y_fitsne = col_double(),
##   x_phate = col_double(),
##   y_phate = col_double(),
##   `x_umap_min_dist=0.1` = col_double(),
##   `y_umap_min_dist=0.1` = col_double(),
##   x_multicoretsne = col_double(),
##   y_multicoretsne = col_double()
## )

Visium embedding

tissue_positions <- readr::read_csv(
    file = file.path(
        PROJECT_DIR,
        "raw/spatial/tissue_positions_list.csv"
    ),
    col_names = c(
        "barcode", "in_tissue", "array_row",
        "array_col", "pxl_col_in_fullres", "pxl_row_in_fullres"
    )
)
## 
## ── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## cols(
##   barcode = col_character(),
##   in_tissue = col_double(),
##   array_row = col_double(),
##   array_col = col_double(),
##   pxl_col_in_fullres = col_double(),
##   pxl_row_in_fullres = col_double()
## )


Visualization of Visium


Barebone

geom_point_size <- 1.2

# UMIs
values <- colSums(matrix_readcount_use)
values <- tissue_positions %>%
    mutate(
        value = case_when(
            barcode %in% colnames(matrix_readcount_use) ~ as.character(values[barcode]),
            TRUE ~ "NA"
        ),
        value = log10(as.integer(value))
    ) %>%
    pull(value)

p_spatial_UMI <- plot_embedding(
    embedding = tissue_positions %>%
        select(pxl_row_in_fullres, pxl_col_in_fullres),
    color_values = values,
    show_color_legend = TRUE,
    geom_point_size = geom_point_size
) +
    scale_y_continuous(trans = scales::reverse_trans()) +
    coord_fixed() +
    theme_customized_spatial() +
    ggplot2::guides(
        color = ggplot2::guide_colourbar(
            title = expression(paste("Log"[10], " (UMIs)"))
        )
    )

# MT
mt_ratio <- (colSums(matrix_readcount_use[rownames(matrix_readcount_use) %>%
    stringr::str_remove(pattern = "^E.+_") %>%
    stringr::str_detect(pattern = "MT-"), ])) / colSums(matrix_readcount_use)

values <- tissue_positions %>%
    mutate(
        value = case_when(
            barcode %in% colnames(matrix_readcount_use) ~ as.character(mt_ratio[barcode]),
            TRUE ~ "NA"
        ),
        value = as.numeric(value)
    ) %>%
    pull(value)

p_spatial_MT <- plot_embedding(
    embedding = tissue_positions %>%
        select(pxl_row_in_fullres, pxl_col_in_fullres),
    color_values = values,
    show_color_legend = TRUE,
    geom_point_size = geom_point_size
) +
    scale_y_continuous(trans = scales::reverse_trans()) +
    coord_fixed() +
    theme_customized_spatial() +
    ggplot2::guides(color = ggplot2::guide_colourbar(title = "MT %"))

# MYL2
selected_features <- "ENSG00000111245_MYL2"

values <- matrix_readcount_use[selected_features, ]
values <- log10(values * 1000000 / colSums(matrix_readcount_use) + 1)

values <- tissue_positions %>%
    mutate(
        value = case_when(
            barcode %in% colnames(matrix_readcount_use) ~ as.character(values[barcode]),
            TRUE ~ "NA"
        ),
        value = as.numeric(value)
    ) %>%
    pull(value)

p_spatial_MYL2 <- plot_embedding(
    embedding = tissue_positions %>%
        select(pxl_row_in_fullres, pxl_col_in_fullres),
    color_values = values,
    show_color_legend = TRUE,
    geom_point_size = geom_point_size
) +
    scale_y_continuous(trans = scales::reverse_trans()) +
    coord_fixed() +
    theme_customized_spatial() +
    ggplot2::guides(
        color = ggplot2::guide_colourbar(
            # title = expression(paste("Log"[10], " (CPM + 1)"))
            title = stringr::str_remove(
                string = selected_features,
                pattern = "^E.+_"
            )
        )
    )

# leiden
values <- tissue_positions %>%
    left_join(
        embedding %>%
            dplyr::select(barcode = cell, leiden)
    ) %>%
    pull(leiden)

p_spatial_leiden <- plot_embedding(
    embedding = tissue_positions %>%
        select(pxl_row_in_fullres, pxl_col_in_fullres),
    color_values = values %>% as.factor(),
    show_color_legend = TRUE,
    geom_point_size = geom_point_size
) +
    scale_y_continuous(trans = scales::reverse_trans()) +
    coord_fixed() +
    theme_customized_spatial() +
    guides(colour = guide_legend(
        title = "Leiden", override.aes = list(size = 2), ncol = 1
    ))
list(
    p_spatial_UMI,
    p_spatial_MT,
    p_spatial_leiden,
    p_spatial_MYL2
) %>% purrr::reduce(`+`) +
    patchwork::plot_layout(
        nrow = 2
    ) +
    patchwork::plot_annotation(
        theme = theme(plot.margin = margin())
    )

mt_ratio %>% fivenum() %>% unname()
## [1] 0.09500247 0.34412550 0.38047412 0.41770273 0.67411658

w/ image

img_info <- summarize_img(
    data_directory = file.path(
        PROJECT_DIR,
        "raw/spatial"
    ),
    sample_id = "Human_heart"
)
img_info
p_spatial_img <- tissue_positions %>%
    filter(in_tissue != 0) %>%
    mutate(
        pxl_row_in_fullres = pxl_row_in_fullres * img_info$scalef,
        pxl_col_in_fullres = pxl_col_in_fullres * img_info$scalef
    ) %>%
    ggplot(
        aes(
            x = pxl_row_in_fullres,
            y = pxl_col_in_fullres
        )
    ) +
    geom_spatial(
        data = img_info[1, ], aes(grob = grob), x = 0.5, y = 0.5
    ) +
    theme_customized()

p_spatial_img

img_x_start <- 71
img_x_end <- 570

img_y_start <- 71
img_y_end <- 570

img_info <- summarize_img(
    data_directory = file.path(
        PROJECT_DIR,
        "raw/spatial"
    ),
    sample_id = "Human_heart",
    image_xlim = c(img_x_start, img_x_end),
    image_ylim = c(img_y_start, img_y_end)
)

p_spatial_img <- tissue_positions %>%
    filter(in_tissue != 0) %>%
    mutate(
        pxl_row_in_fullres = pxl_row_in_fullres * img_info$scalef,
        pxl_col_in_fullres = pxl_col_in_fullres * img_info$scalef
    ) %>%
    ggplot(
        aes(
            x = pxl_row_in_fullres,
            y = pxl_col_in_fullres
        )
    ) +
    geom_spatial(
        data = img_info[1, ], aes(grob = grob), x = 0.5, y = 0.5
    ) +
    theme_customized()

geom_point_size <- 1

# UMI
p_spatial_img_UMI <- tissue_positions %>%
    filter(in_tissue != 0) %>%
    mutate(
        value = Matrix::colSums(matrix_readcount_use[, barcode]),
        pxl_row_in_fullres = pxl_row_in_fullres * img_info$scalef,
        pxl_col_in_fullres = pxl_col_in_fullres * img_info$scalef
    ) %>%
    ggplot(
        aes(
            x = pxl_row_in_fullres,
            y = pxl_col_in_fullres,
            color = value %>% log10()
        )
    ) +
    geom_spatial(
        data = img_info[1, ], aes(grob = grob), x = 0.5, y = 0.5
    ) +
    geom_point(stroke = 0, shape = 16, size = geom_point_size) +
    scale_y_continuous(trans = scales::reverse_trans()) +
    coord_fixed(
        xlim = c(img_x_start, img_x_end),
        ylim = c(img_y_end, img_y_start),
        expand = FALSE
    ) +
    theme_customized_spatial_img(x = 0.06, y = 0.22) +
    ggplot2::guides(
        color = ggplot2::guide_colourbar(
            title = expression(paste("Log"[10], " (UMIs)")),
            ticks.colour = "black"
        )
    )

# MT
p_spatial_img_MT <- tissue_positions %>%
    filter(in_tissue != 0) %>%
    mutate(
        value = mt_ratio[barcode],
        pxl_row_in_fullres = pxl_row_in_fullres * img_info$scalef,
        pxl_col_in_fullres = pxl_col_in_fullres * img_info$scalef
    ) %>%
    ggplot(
        aes(
            x = pxl_row_in_fullres,
            y = pxl_col_in_fullres,
            color = value
        )
    ) +
    geom_spatial(
        data = img_info[1, ], aes(grob = grob), x = 0.5, y = 0.5
    ) +
    geom_point(stroke = 0, shape = 16, size = 1) +
    scale_y_continuous(trans = scales::reverse_trans()) +
    coord_fixed(
        xlim = c(img_x_start, img_x_end),
        ylim = c(img_y_end, img_y_start),
        expand = FALSE
    ) +
    theme_customized_spatial_img(x = 0.06, y = 0.22) +
    ggplot2::guides(
        color = ggplot2::guide_colourbar(
            title = "MT %", ticks.colour = "black"
        )
    )

# Leiden
p_spatial_img_leiden <- tissue_positions %>%
    left_join(
        embedding %>% dplyr::select(cell, leiden),
        by = c("barcode" = "cell")
    ) %>%
    # filter(in_tissue != 0) %>%
    filter(!is.na(leiden)) %>%
    mutate(
        pxl_row_in_fullres = pxl_row_in_fullres * img_info$scalef,
        pxl_col_in_fullres = pxl_col_in_fullres * img_info$scalef
    ) %>%
    ggplot(
        aes(
            x = pxl_row_in_fullres,
            y = pxl_col_in_fullres,
            color = leiden %>% as.factor()
        )
    ) +
    geom_spatial(
        data = img_info[1, ], aes(grob = grob), x = 0.5, y = 0.5
    ) +
    geom_point(stroke = 0, shape = 16, size = 1) +
    scale_y_continuous(trans = scales::reverse_trans()) +
    coord_fixed(
        xlim = c(img_x_start, img_x_end),
        ylim = c(img_y_end, img_y_start),
        expand = FALSE
    ) +
    theme_customized_spatial_img(x = 0.06) +
    scale_color_manual(
        values = ggthemes::tableau_color_pal("Tableau 10")(10)
    ) +
    guides(
        colour = guide_legend(
            title = "Leiden", override.aes = list(size = 2), ncol = 1
        )
    )
list(
    p_spatial_img,
    p_spatial_img_UMI +
        scale_color_gradientn(
            colors = colorRampPalette(
                colors = rev(
                    x = RColorBrewer::brewer.pal(n = 11, name = "Spectral")
                )
            )(n = 100)
        ),
    p_spatial_img_MT + scale_color_gradientn(
        colors = colorRampPalette(
            colors = rev(
                x = RColorBrewer::brewer.pal(n = 11, name = "Spectral")
            )
        )(n = 100)
    ),
    p_spatial_img_leiden
) %>% purrr::reduce(`+`) +
    patchwork::plot_layout(
        nrow = 2
    ) +
    patchwork::plot_annotation(
        theme = theme(plot.margin = margin())
    )


Extended Data Fig. 5

FEATURES_SELECTED <- c(
    "ENSG00000179776_CDH5",
    "ENSG00000010319_SEMA3G",
    "ENSG00000265107_GJA5",
    "ENSG00000213088_ACKR1",
    "ENSG00000130300_PLVAP",
    "ENSG00000133392_MYH11",
    "ENSG00000107796_ACTA2",
    "ENSG00000101384_JAG1",
    "ENSG00000134250_NOTCH2"
)

purrr::map(FEATURES_SELECTED, function(x) {
    selected_features <- x

    tissue_positions %>%
        filter(in_tissue != 0) %>%
        mutate(
            value = matrix_readcount_use[selected_features, barcode],
            pxl_row_in_fullres = pxl_row_in_fullres * img_info$scalef,
            pxl_col_in_fullres = pxl_col_in_fullres * img_info$scalef,
            value = log10(value * 1000000 / colSums(matrix_readcount_use) + 1)
        ) %>%
        ggplot(
            aes(
                x = pxl_row_in_fullres,
                y = pxl_col_in_fullres,
                color = value
            )
        ) +
        geom_spatial(
            data = img_info[1, ], aes(grob = grob), x = 0.5, y = 0.5
        ) +
        geom_point(stroke = 0, shape = 16, size = geom_point_size) +
        scale_y_continuous(trans = scales::reverse_trans()) +
        coord_fixed(
            xlim = c(img_x_start, img_x_end),
            ylim = c(img_y_end, img_y_start),
            expand = FALSE
        ) +
        theme_customized_spatial_img(x = 0.06, y = 0.22) +
        ggplot2::guides(
            color = ggplot2::guide_colourbar(
                title = expression(paste("Log"[10], " (UMIs)")),
                ticks.colour = "black"
            )
        ) +
        scale_color_gradientn(
            colors = colorRampPalette(
                colors = rev(
                    x = RColorBrewer::brewer.pal(n = 11, name = "Spectral")
                )
            )(n = 100)
        ) +
        ggtitle(selected_features)
}) %>%
    purrr::reduce(`+`) +
    patchwork::plot_layout(
        ncol = 2
    ) +
    patchwork::plot_annotation(
        theme = theme(plot.margin = margin())
    )



R session info

devtools::session_info()$platform
##  setting  value                       
##  version  R version 4.0.3 (2020-10-10)
##  os       macOS Big Sur 10.16         
##  system   x86_64, darwin20.2.0        
##  ui       unknown                     
##  language (EN)                        
##  collate  en_US.UTF-8                 
##  ctype    en_US.UTF-8                 
##  tz       America/Chicago             
##  date     2021-02-03
devtools::session_info()$pack %>%
    as_tibble() %>%
    dplyr::select(
        package,
        loadedversion,
        date,
        `source`
    ) %>%
    # print(n = Inf)
    gt::gt() %>%
    gt::tab_options(table.font.size = "median")
package loadedversion date source
assertthat 0.2.1 2019-03-21 CRAN (R 4.0.0)
backports 1.2.1 2020-12-09 CRAN (R 4.0.3)
broom 0.7.4.9000 2021-02-02 Github (tidymodels/broom@b750e4f)
bslib 0.2.4.9000 2021-02-02 Github (rstudio/bslib@b3cd7a9)
cachem 1.0.1 2021-01-22 Github (r-lib/cachem@27c8d89)
callr 3.5.1.9000 2021-01-06 Github (r-lib/callr@743069f)
cellranger 1.1.0 2016-07-27 CRAN (R 4.0.0)
checkmate 2.0.0 2020-02-06 CRAN (R 4.0.0)
cli 2.3.0 2021-01-31 CRAN (R 4.0.3)
colorspace 2.0-0 2020-11-10 R-Forge (R 4.0.3)
crayon 1.4.0 2021-01-30 CRAN (R 4.0.3)
data.table 1.13.7 2021-01-26 local
DBI 1.1.1 2021-01-15 CRAN (R 4.0.3)
dbplyr 2.0.0 2020-11-03 CRAN (R 4.0.3)
desc 1.2.0 2018-05-01 CRAN (R 4.0.0)
devtools 2.3.1.9000 2021-01-21 Github (r-lib/devtools@ef962e4)
digest 0.6.27 2020-10-24 CRAN (R 4.0.3)
dplyr 1.0.4.9000 2021-02-02 Github (tidyverse/dplyr@2455c77)
ellipsis 0.3.1 2020-05-15 CRAN (R 4.0.3)
evaluate 0.14 2019-05-28 CRAN (R 4.0.0)
extrafont 0.17 2014-12-08 CRAN (R 4.0.2)
extrafontdb 1.0 2012-06-11 CRAN (R 4.0.0)
fansi 0.4.2 2021-01-15 CRAN (R 4.0.3)
farver 2.0.3 2020-01-16 CRAN (R 4.0.0)
fastmap 1.1.0 2021-01-25 CRAN (R 4.0.3)
forcats 0.5.1.9000 2021-01-28 Github (tidyverse/forcats@b5fce89)
fs 1.5.0 2020-07-31 CRAN (R 4.0.3)
generics 0.1.0 2020-10-31 CRAN (R 4.0.3)
ggplot2 3.3.3.9000 2021-01-30 Github (tidyverse/ggplot2@dbd7d79)
ggthemes 4.2.4 2021-01-20 Github (jrnold/ggthemes@4b5e80e)
glue 1.4.1.9000 2021-01-06 Github (tidyverse/glue@f0a7b2a)
gt 0.2.2 2020-11-25 Github (rstudio/gt@bae32f4)
gtable 0.3.0 2019-03-25 CRAN (R 4.0.0)
haven 2.3.1 2020-06-01 CRAN (R 4.0.0)
highr 0.8 2019-03-20 CRAN (R 4.0.0)
hms 1.0.0 2021-01-13 CRAN (R 4.0.3)
htmltools 0.5.1.9000 2021-01-23 Github (rstudio/htmltools@e7f0393)
httr 1.4.2 2020-07-20 CRAN (R 4.0.2)
jquerylib 0.1.3 2020-12-17 Github (rstudio/jquerylib@8f8e639)
jsonlite 1.7.2 2020-12-09 CRAN (R 4.0.3)
knitr 1.31.4 2021-01-29 Github (yihui/knitr@d83e8de)
labeling 0.4.2 2020-10-20 CRAN (R 4.0.3)
lattice 0.20-41 2020-04-02 CRAN (R 4.0.3)
lifecycle 0.2.0 2020-03-06 CRAN (R 4.0.0)
lubridate 1.7.9.2 2021-01-04 Github (tidyverse/lubridate@aab2e30)
magrittr 2.0.1.9000 2020-12-14 Github (tidyverse/magrittr@bb1c86a)
Matrix 1.3-2 2021-01-06 CRAN (R 4.0.3)
memoise 2.0.0 2021-01-26 CRAN (R 4.0.3)
modelr 0.1.8.9000 2021-01-23 Github (tidyverse/modelr@16168e0)
munsell 0.5.0 2018-06-12 CRAN (R 4.0.0)
patchwork 1.1.1 2020-12-17 CRAN (R 4.0.3)
pillar 1.4.99.9006 2021-02-02 Github (r-lib/pillar@c6f8311)
pkgbuild 1.2.0 2020-12-15 CRAN (R 4.0.3)
pkgconfig 2.0.3 2019-09-22 CRAN (R 4.0.0)
pkgload 1.1.0 2020-05-29 CRAN (R 4.0.0)
png 0.1-7 2013-12-03 CRAN (R 4.0.0)
prettyunits 1.1.1.9000 2020-11-23 Github (r-lib/prettyunits@b1cdad8)
processx 3.4.5 2020-11-30 CRAN (R 4.0.3)
ps 1.5.0 2020-12-05 CRAN (R 4.0.3)
purrr 0.3.4.9000 2020-11-23 Github (tidyverse/purrr@af06d45)
R6 2.5.0 2020-11-02 Github (r-lib/R6@6cf7d4e)
ragg 1.0.0.9000 2021-01-15 Github (r-lib/ragg@aebed7f)
RColorBrewer 1.1-2 2014-12-07 CRAN (R 4.0.0)
Rcpp 1.0.6 2021-01-15 CRAN (R 4.0.3)
readr 1.4.0.9000 2021-01-23 Github (tidyverse/readr@483ca6c)
readxl 1.3.1.9000 2021-01-23 Github (tidyverse/readxl@9f85fa5)
remotes 2.2.0.9000 2021-01-25 Github (r-lib/remotes@cf2b4a9)
reprex 1.0.0 2021-01-27 CRAN (R 4.0.3)
rhdf5 2.34.0 2020-10-27 Bioconductor
rhdf5filters 1.2.0 2020-10-27 Bioconductor
Rhdf5lib 1.12.1 2021-01-26 Bioconductor
rjson 0.2.20 2018-06-08 CRAN (R 4.0.0)
rlang 0.4.10.9000 2021-02-02 Github (r-lib/rlang@d15299e)
rmarkdown 2.6.6 2021-02-02 Github (rstudio/rmarkdown@d8e7a15)
rprojroot 2.0.2 2020-11-15 CRAN (R 4.0.3)
rstudioapi 0.13.0-9000 2021-01-02 Github (rstudio/rstudioapi@4baeb39)
Rttf2pt1 1.3.8 2020-01-10 CRAN (R 4.0.0)
rvest 0.3.6 2020-07-25 CRAN (R 4.0.2)
sass 0.3.1 2021-01-24 CRAN (R 4.0.3)
scales 1.1.1 2020-05-11 CRAN (R 4.0.3)
sessioninfo 1.1.1 2018-11-05 CRAN (R 4.0.3)
stringi 1.5.3 2020-09-09 CRAN (R 4.0.2)
stringr 1.4.0.9000 2021-01-23 Github (tidyverse/stringr@1f03eb0)
styler 1.3.2.9000 2021-02-02 Github (r-lib/styler@3560b39)
systemfonts 1.0.0.9000 2021-02-02 Github (r-lib/systemfonts@d85abe2)
testthat 3.0.1 2021-01-29 Github (r-lib/testthat@b19b5ac)
textshaping 0.2.1.9000 2021-01-15 Github (r-lib/textshaping@f6f2697)
tibble 3.0.6.9000 2021-01-31 Github (tidyverse/tibble@eb99cb6)
tidyr 1.1.2.9000 2021-01-23 Github (tidyverse/tidyr@c338aa9)
tidyselect 1.1.0 2020-05-11 CRAN (R 4.0.3)
tidyverse 1.3.0.9000 2020-11-23 Github (hadley/tidyverse@8a0bb99)
usethis 2.0.0.9000 2021-01-31 Github (r-lib/usethis@716b703)
utf8 1.1.4 2018-05-24 CRAN (R 4.0.0)
vctrs 0.3.6 2020-12-17 CRAN (R 4.0.3)
viridisLite 0.3.0 2018-02-01 CRAN (R 4.0.0)
withr 2.4.1 2021-01-26 CRAN (R 4.0.3)
xfun 0.20 2021-01-06 CRAN (R 4.0.3)
xml2 1.3.2 2020-04-23 CRAN (R 4.0.0)
yaml 2.2.1 2020-02-01 CRAN (R 4.0.0)