Single-cell transcriptomes of the three cell lineages of the Human blastocyst

Jialei Duan

2022-02-07

Blakeley, P., Fogarty, N.M.E., del Valle, I., Wamaitha, S.E., Hu, T.X., Elder, K., Snell, P., Christie, L., Robson, P., and Niakan, K.K. (2015). Defining the three cell lineages of the human blastocyst by single-cell RNA-seq. Development 142, 3151–3165.

BioProject Accession: PRJNA277181
GEO Accession: GSE66507

GitHub repository: jlduan/Replica@GitHub

Rendered at 2022-02-07 02:05:19 America/Chicago.

Load required packages.

library(tidyverse)
library(Matrix)
library(patchwork)
library(extrafont)

Sys.time()

[1] "2022-02-07 02:05:20 CST"

Sys.timezone()

[1] "America/Chicago"

Data preparation

Functions loading

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

Data loading

PROJECT_DIR <- file.path(
    "/Users/jialei/Dropbox/Data/Projects/UTSW/Peri-implantation/"
)

ad <- reticulate::import(module = "anndata", convert = TRUE)
print(ad$`__version__`)

[1] "0.7.6"

adata_files <- purrr::map(c("PRJNA277181"), \(x) {
    file.path(
        PROJECT_DIR,
        "raw",
        "public",
        x,
        "matrix",
        "adata.h5ad"
    )
})
purrr::map_lgl(adata_files, file.exists)

[1] TRUE

BACKED <- NULL
matrix_readcount_use <- purrr::map(adata_files, function(x) {
    ad$read_h5ad(
        filename = x, backed = BACKED
    ) |>
        convert_adata()
}) |>
    purrr::reduce(cbind)

Metadata

BACKED <- "r"
cell_metadata <- purrr::map(adata_files, function(x) {
    ad$read_h5ad(
        filename = x, backed = BACKED
    )$obs |>
        tibble::rownames_to_column(var = "cell") |>
        dplyr::select(cell, everything())
}) |>
    dplyr::bind_rows()

cell_metadata |>
    head() |>
    knitr::kable()

cell	batch	num_umis	num_features	mt_percentage
GSM1624222	PRJNA277181	16601969	8180	0.0551772
GSM1624223	PRJNA277181	10942129	7742	0.0780155
GSM1624224	PRJNA277181	5953198	6561	0.0700953
GSM1624225	PRJNA277181	16240393	7977	0.0462827
GSM1624226	PRJNA277181	7168882	6764	0.0710903
GSM1624227	PRJNA277181	2558268	7219	0.1186041

cell_metadata_PRJNA277181 <- vroom::vroom(
    file = file.path(
        PROJECT_DIR,
        "raw",
        "public",
        "PRJNA277181",
        "matrix",
        "cell_metadata.csv"
    )
)

cell_metadata_PRJNA277181 |>
    dplyr::count(lineage)

# A tibble: 3 × 2
  lineage                n
  <chr>              <int>
1 Epiblast               8
2 Primitive Endoderm     6
3 Trophectoderm         16

cell_metadata_PRJNA277181 |>
    dplyr::count(developmental_stage)

# A tibble: 1 × 2
  developmental_stage     n
  <chr>               <int>
1 Blastocyst             30

Clustering of human blastocysts

Clustering

SEED <- 20210719
N_COMPONENTS <- 30
N_COMPONENTS_SELECTED <- 20

MINIMAL_NUM_CELLS_REQUIRED_FOR_GENE <- 1
MINIMAL_NUM_COUNTS_REQUIRED_FOR_GENE <- 60
RESOLUTION <- 0.8

N_FEATURES <- 3000
REDUCTION <- "pca"
PERPLEXITY <- 5

features_selected <- rownames(matrix_readcount_use)[
    (Matrix::rowSums(matrix_readcount_use > 0) >=
        MINIMAL_NUM_CELLS_REQUIRED_FOR_GENE) & (
        Matrix::rowSums(matrix_readcount_use) >=
            MINIMAL_NUM_COUNTS_REQUIRED_FOR_GENE)
]
features_selected |> length()

[1] 14671

# begin
seurat_obj <- Seurat::CreateSeuratObject(
    counts = matrix_readcount_use[features_selected, ],
    project = "Duan",
    assay = "RNA",
    names.field = 1,
    names.delim = "_",
    meta.data = data.frame(
        cell_metadata,
        row.names = cell_metadata$cell
    )[colnames(matrix_readcount_use), ],
    #
    min.cells = 0,
    min.features = 0
) |>
    Seurat::NormalizeData(
        normalization.method = "LogNormalize",
        scale.factor = colSums(matrix_readcount_use[features_selected, ]) |>
            median(),
        margin = 1,
        verbose = TRUE
    ) |>
    Seurat::FindVariableFeatures(
        assay = NULL,
        selection.method = "vst",
        loess.span = 0.3,
        clip.max = "auto",
        # mean.function = FastExpMean,
        # dispersion.function = FastLogVMR,
        num.bin = 20,
        binning.method = "equal_width",
        nfeatures = N_FEATURES,
        mean.cutoff = c(0.1, 8),
        dispersion.cutoff = c(1, Inf),
        verbose = TRUE
    ) |>
    Seurat::ScaleData(
        # features = features_selected,
        vars.to.regress = NULL,
        model.use = "linear",
        use.umi = FALSE,
        do.scale = TRUE,
        do.center = TRUE,
        scale.max = Inf,
        block.size = 500,
        min.cells.to.block = 3000,
        verbose = TRUE
    ) |>
    Seurat::RunPCA(
        # features = features_selected,
        npcs = N_COMPONENTS,
        rev.pca = FALSE,
        weight.by.var = TRUE,
        verbose = TRUE,
        ndims.print = 1:2,
        nfeatures.print = 5,
        reduction.name = "pca",
        reduction.key = "PC_",
        seed.use = SEED,
        approx = FALSE
    ) |>
    Seurat::RunUMAP(
        dims = seq_len(N_COMPONENTS_SELECTED),
        reduction = REDUCTION,
        features = NULL,
        graph = NULL,
        assay = NULL,
        nn.name = NULL,
        slot = "data",
        umap.method = "uwot",
        reduction.model = NULL,
        return.model = FALSE,
        n.neighbors = 30L,
        n.components = 2L,
        metric = "cosine",
        n.epochs = NULL,
        learning.rate = 1,
        min.dist = 0.1,
        spread = 1,
        set.op.mix.ratio = 1,
        local.connectivity = 1L,
        repulsion.strength = 1,
        negative.sample.rate = 5L,
        a = NULL,
        b = NULL,
        uwot.sgd = FALSE,
        seed.use = SEED,
        metric.kwds = NULL,
        angular.rp.forest = FALSE,
        verbose = TRUE,
        reduction.name = "umap",
        reduction.key = "UMAP_"
    ) |>
    Seurat::RunTSNE(
        reduction = REDUCTION,
        cells = NULL,
        dims = seq_len(N_COMPONENTS_SELECTED),
        features = NULL,
        seed.use = SEED,
        tsne.method = "Rtsne",
        dim.embed = 2,
        distance.matrix = NULL,
        reduction.name = "tsne",
        reduction.key = "tSNE_",
        #
        perplexity = PERPLEXITY,
        max_iter = 3000,
        verbose = TRUE,
        #
        check_duplicates = TRUE
    ) |>
    Seurat::FindNeighbors(
        reduction = REDUCTION,
        dims = seq_len(N_COMPONENTS_SELECTED),
        assay = NULL,
        features = NULL,
        k.param = 20,
        return.neighbor = FALSE,
        # compute.SNN = !return.neighbor,
        prune.SNN = 1 / 15,
        nn.method = "annoy",
        n.trees = 50,
        annoy.metric = "euclidean",
        nn.eps = 0,
        verbose = TRUE,
        force.recalc = FALSE,
        do.plot = FALSE,
        graph.name = NULL,
        l2.norm = FALSE,
        cache.index = FALSE
    ) |>
    Seurat::FindClusters(
        graph.name = NULL,
        modularity.fxn = 1,
        initial.membership = NULL,
        node.sizes = NULL,
        resolution = RESOLUTION,
        method = "igraph",
        algorithm = 4,
        n.start = 100,
        n.iter = 100,
        random.seed = SEED,
        group.singletons = TRUE,
        temp.file.location = NULL,
        edge.file.name = NULL,
        verbose = TRUE
    )

Registered S3 method overwritten by 'spatstat.geom':
  method     from
  print.boxx cli

Warning: Feature names cannot have underscores ('_'), replacing with dashes
('-')

Centering and scaling data matrix

PC_ 1 
Positive:  ENSG00000241186-TDGF1, ENSG00000077782-FGFR1, ENSG00000137285-TUBB2B, ENSG00000157601-MX1, ENSG00000187193-MT1X 
Negative:  ENSG00000111371-SLC38A1, ENSG00000163382-NAXE, ENSG00000083307-GRHL2, ENSG00000102243-VGLL1, ENSG00000003989-SLC7A2 
PC_ 2 
Positive:  ENSG00000112419-PHACTR2, ENSG00000108883-EFTUD2, ENSG00000118137-APOA1, ENSG00000109654-TRIM2, ENSG00000115414-FN1 
Negative:  ENSG00000130055-GDPD2, ENSG00000168907-PLA2G4F, ENSG00000112759-SLC29A1, ENSG00000162757-C1orf74, ENSG00000124766-SOX4

Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session

02:05:23 UMAP embedding parameters a = 1.577 b = 0.8951

02:05:23 Read 30 rows and found 20 numeric columns

02:05:23 Using Annoy for neighbor search, n_neighbors = 30

02:05:23 Building Annoy index with metric = cosine, n_trees = 50

02:05:23 Writing NN index file to temp file /var/folders/yr/vrcclyl91kd_dg97yrqtwqz40000gn/T//RtmpyEZmal/file9e38364121a8

02:05:23 Searching Annoy index using 1 thread, search_k = 3000

02:05:23 Annoy recall = 100%

02:05:23 Commencing smooth kNN distance calibration using 1 thread

02:05:23 Initializing from normalized Laplacian + noise (using irlba)

02:05:23 Commencing optimization for 500 epochs, with 870 positive edges

02:05:24 Optimization finished

Read the 30 x 20 data matrix successfully!
Using no_dims = 2, perplexity = 5.000000, and theta = 0.500000
Computing input similarities...
Building tree...
Done in 0.00 seconds (sparsity = 0.720000)!
Learning embedding...
Iteration 50: error is 74.069343 (50 iterations in 0.00 seconds)
Iteration 100: error is 65.065538 (50 iterations in 0.00 seconds)
Iteration 150: error is 71.986840 (50 iterations in 0.00 seconds)
Iteration 200: error is 62.424255 (50 iterations in 0.00 seconds)
Iteration 250: error is 61.934033 (50 iterations in 0.00 seconds)
Iteration 300: error is 1.723239 (50 iterations in 0.00 seconds)
Iteration 350: error is 1.235360 (50 iterations in 0.00 seconds)
Iteration 400: error is 1.098426 (50 iterations in 0.00 seconds)
Iteration 450: error is 1.061904 (50 iterations in 0.00 seconds)
Iteration 500: error is 0.992590 (50 iterations in 0.00 seconds)
Iteration 550: error is 0.904049 (50 iterations in 0.00 seconds)
Iteration 600: error is 0.882830 (50 iterations in 0.00 seconds)
Iteration 650: error is 0.870322 (50 iterations in 0.00 seconds)
Iteration 700: error is 0.854095 (50 iterations in 0.00 seconds)
Iteration 750: error is 0.727305 (50 iterations in 0.00 seconds)
Iteration 800: error is 0.727670 (50 iterations in 0.00 seconds)
Iteration 850: error is 0.705883 (50 iterations in 0.00 seconds)
Iteration 900: error is 0.699640 (50 iterations in 0.00 seconds)
Iteration 950: error is 0.703992 (50 iterations in 0.00 seconds)
Iteration 1000: error is 0.699231 (50 iterations in 0.00 seconds)
Iteration 1050: error is 0.707053 (50 iterations in 0.00 seconds)
Iteration 1100: error is 0.705576 (50 iterations in 0.00 seconds)
Iteration 1150: error is 0.702926 (50 iterations in 0.00 seconds)
Iteration 1200: error is 0.704696 (50 iterations in 0.00 seconds)
Iteration 1250: error is 0.703848 (50 iterations in 0.00 seconds)
Iteration 1300: error is 0.703590 (50 iterations in 0.00 seconds)
Iteration 1350: error is 0.699929 (50 iterations in 0.00 seconds)
Iteration 1400: error is 0.696757 (50 iterations in 0.00 seconds)
Iteration 1450: error is 0.702982 (50 iterations in 0.00 seconds)
Iteration 1500: error is 0.700752 (50 iterations in 0.00 seconds)
Iteration 1550: error is 0.701990 (50 iterations in 0.00 seconds)
Iteration 1600: error is 0.699863 (50 iterations in 0.00 seconds)
Iteration 1650: error is 0.690780 (50 iterations in 0.00 seconds)
Iteration 1700: error is 0.696472 (50 iterations in 0.00 seconds)
Iteration 1750: error is 0.700205 (50 iterations in 0.00 seconds)
Iteration 1800: error is 0.697540 (50 iterations in 0.00 seconds)
Iteration 1850: error is 0.703494 (50 iterations in 0.00 seconds)
Iteration 1900: error is 0.694863 (50 iterations in 0.00 seconds)
Iteration 1950: error is 0.701220 (50 iterations in 0.00 seconds)
Iteration 2000: error is 0.698841 (50 iterations in 0.00 seconds)
Iteration 2050: error is 0.699521 (50 iterations in 0.00 seconds)
Iteration 2100: error is 0.700054 (50 iterations in 0.00 seconds)
Iteration 2150: error is 0.693202 (50 iterations in 0.00 seconds)
Iteration 2200: error is 0.699307 (50 iterations in 0.00 seconds)
Iteration 2250: error is 0.699259 (50 iterations in 0.00 seconds)
Iteration 2300: error is 0.697938 (50 iterations in 0.00 seconds)
Iteration 2350: error is 0.698317 (50 iterations in 0.00 seconds)
Iteration 2400: error is 0.693987 (50 iterations in 0.00 seconds)
Iteration 2450: error is 0.696271 (50 iterations in 0.00 seconds)
Iteration 2500: error is 0.695814 (50 iterations in 0.00 seconds)
Iteration 2550: error is 0.696823 (50 iterations in 0.00 seconds)
Iteration 2600: error is 0.698417 (50 iterations in 0.00 seconds)
Iteration 2650: error is 0.698130 (50 iterations in 0.00 seconds)
Iteration 2700: error is 0.695174 (50 iterations in 0.00 seconds)
Iteration 2750: error is 0.695155 (50 iterations in 0.00 seconds)
Iteration 2800: error is 0.694876 (50 iterations in 0.00 seconds)
Iteration 2850: error is 0.695021 (50 iterations in 0.00 seconds)
Iteration 2900: error is 0.695160 (50 iterations in 0.00 seconds)
Iteration 2950: error is 0.699060 (50 iterations in 0.00 seconds)
Iteration 3000: error is 0.697582 (50 iterations in 0.00 seconds)
Fitting performed in 0.03 seconds.

Computing nearest neighbor graph

Computing SNN

embedding <- seurat_obj@reductions$umap@cell.embeddings |>
    as.data.frame() |>
    tibble::rownames_to_column(var = "cell") |>
    dplyr::left_join(
        seurat_obj@meta.data,
        by = c("cell" = "cell")
    ) |>
    dplyr::select(
        -c("orig.ident", "nCount_RNA", "nFeature_RNA")
    ) |>
    dplyr::mutate(
        leiden = Seurat::Idents(seurat_obj)
    ) |>
    dplyr::left_join(
        seurat_obj@reductions$pca@cell.embeddings |>
            as.data.frame() |>
            tibble::rownames_to_column(var = "cell") |>
            dplyr::select(cell, x_pca = PC_1, y_pca = PC_2),
        by = "cell"
    ) |>
    dplyr::left_join(
        seurat_obj@reductions$tsne@cell.embeddings |>
            as.data.frame() |>
            tibble::rownames_to_column(var = "cell") |>
            dplyr::select(cell, x_tsne = tSNE_1, y_tsne = tSNE_2),
        by = "cell"
    ) |>
    dplyr::select(
        cell,
        batch,
        leiden,
        x_pca,
        y_pca,
        x_umap = UMAP_1,
        y_umap = UMAP_2,
        x_tsne,
        y_tsne,
        everything()
    )

Visualization

GEOM_POINT_SIZE <- 0.6
RASTERISED <- TRUE
x_column <- "x_umap"
y_column <- "y_umap"

GEOM_POINT_SIZE <- 2
RASTERISED <- FALSE

embedding <- embedding |>
    dplyr::left_join(
        cell_metadata_PRJNA277181 |>
            dplyr::select(
                cell = sample_name,
                lineage,
                embryo_number
            ),
        by = "cell"
    )

Code

embedding |>
    dplyr::left_join(
        cell_metadata_PRJNA277181
    ) |>
    dplyr::group_by(embryo_number) |>
    dplyr::summarise(
        num_cells = n(),
        median_umis = median(num_umis)
    ) |>
    gt::gt() |>
    gt::data_color(
        columns = c(median_umis),
        colors = scales::col_numeric(
            palette = c(
                "green", "orange", "red"
            ),
            domain = NULL
        )
    ) |>
    gt::fmt_number(
        columns = c(median_umis),
        decimals = 0,
        use_seps = TRUE,
        suffixing = FALSE
    ) |>
    gt::summary_rows(
        columns = c(embryo_number),
        fns = list(
            Count = ~ n()
        ),
        decimals = 0
    ) |>
    gt::summary_rows(
        columns = c(median_umis),
        fns = list(
            Mean = ~ mean(.)
        ),
        decimals = 0
    ) |>
    gt::summary_rows(
        columns = c(num_cells),
        fns = list(
            Sum = ~ sum(.)
        ),
        decimals = 0
    )

	embryo_number	num_cells	median_umis
	Embryo 2	1	11,917,699
	Embryo 3	37	7,168,882
	Embryo 4	11	17,894,715
	Embryo 5	20	8,267,444
	Embryo 6	1	9,956,208
	Embryo 8	65	12,239,684
	Embryo 9	1	3,152,485
Count	7	—	—
Mean	—	—	10,085,302
Sum	—	136	—

Code

embedding |>
    dplyr::left_join(
        cell_metadata_PRJNA277181
    ) |>
    dplyr::group_by(lineage) |>
    dplyr::summarise(
        num_cells = n(),
        median_umis = median(num_umis)
    ) |>
    gt::gt() |>
    gt::data_color(
        columns = c(median_umis),
        colors = scales::col_numeric(
            palette = c(
                "green", "orange", "red"
            ),
            domain = NULL
        )
    ) |>
    gt::fmt_number(
        columns = c(median_umis),
        decimals = 0,
        use_seps = TRUE,
        suffixing = FALSE
    ) |>
    gt::summary_rows(
        columns = c(num_cells),
        fns = list(
            Count = ~ n()
        ),
        decimals = 0
    ) |>
    gt::summary_rows(
        columns = c(median_umis),
        fns = list(
            Mean = ~ mean(.)
        ),
        decimals = 0
    ) |>
    gt::summary_rows(
        columns = c(num_cells),
        fns = list(
            Sum = ~ sum(.)
        ),
        decimals = 0
    )

	lineage	num_cells	median_umis
	Epiblast	16	12,280,532
	Primitive Endoderm	18	9,539,092
	Trophectoderm	102	10,456,983
Count	—	3	—
Mean	—	—	10,758,869
Sum	—	136	—

Code

embedding |>
    dplyr::left_join(
        cell_metadata_PRJNA277181
    ) |>
    dplyr::group_by(lineage, embryo_number) |>
    dplyr::summarise(
        num_cells = n(),
        median_umis = median(num_umis)
    ) |>
    gt::gt() |>
    gt::data_color(
        columns = c(median_umis),
        colors = scales::col_numeric(
            palette = c(
                "green", "orange", "red"
            ),
            domain = NULL
        )
    ) |>
    gt::data_color(
        columns = c(embryo_number),
        colors = scales::col_factor(
            palette = paletteer::paletteer_d(
                n = 6, palette = "colorblindr::OkabeIto"
            ) |> as.character(),
            domain = NULL
        )
    ) |>
    gt::fmt_number(
        columns = c(median_umis),
        decimals = 0,
        use_seps = TRUE,
        suffixing = FALSE
    ) |>
    gt::summary_rows(
        columns = c(num_cells),
        fns = list(
            Sum = ~ sum(.)
        ),
        decimals = 0
    )

	embryo_number	num_cells	median_umis
Epiblast
	Embryo 2	1	11,917,699
	Embryo 3	1	831,175
	Embryo 4	9	17,894,715
	Embryo 5	4	9,819,580
	Embryo 6	1	9,956,208
Primitive Endoderm
	Embryo 4	1	18,879,216
	Embryo 5	16	8,100,936
	Embryo 8	1	12,239,684
Trophectoderm
	Embryo 3	36	9,055,506
	Embryo 4	1	16,678,709
	Embryo 8	64	11,224,321
	Embryo 9	1	3,152,485
Sum	—	136	—

PCA

total_variance <- sum(
    matrixStats::rowVars(
        Seurat::GetAssayData(seurat_obj, assay = "RNA", slot = "scale.data")
    )
)

pc_variance_explained <- (seurat_obj@reductions$pca@stdev)^2 / total_variance

x <- "pca"
purrr::map(c("lineage", "embryo_number"), \(y) {
    plot_embedding(
        data = embedding[, paste0(c("x_", "y_"), x)],
        color = embedding[[y]] |> as.factor(),
        label = glue::glue("{toupper(x)}; {stringr::str_to_title(y)}"),
        color_labels = FALSE,
        color_legend = TRUE,
        sort_values = FALSE,
        shuffle_values = FALSE,
        rasterise = RASTERISED,
        geom_point_size = GEOM_POINT_SIZE
    ) +
        theme_customized_embedding(
            axis_text = TRUE,
            axis_title = TRUE
        ) +
        labs(
            x = paste("PC1", scales::percent(pc_variance_explained[1], accuracy = 0.01)),
            y = paste("PC2", scales::percent(pc_variance_explained[1], accuracy = 0.01))
        )
}) |>
    purrr::reduce(`+`) +
    patchwork::plot_layout(ncol = 2, byrow = TRUE) +
    patchwork::plot_annotation(
        theme = ggplot2::theme(plot.margin = ggplot2::margin())
    )

Non-linear

purrr::map(c("umap", "tsne"), \(x) {
    purrr::map(c("lineage", "embryo_number"), \(y) {
        plot_embedding(
            data = embedding[, paste0(c("x_", "y_"), x)],
            color = embedding[[y]] |> as.factor(),
            label = glue::glue("{toupper(x)}; {stringr::str_to_title(y)}"),
            color_labels = FALSE,
            color_legend = TRUE,
            sort_values = FALSE,
            shuffle_values = FALSE,
            rasterise = RASTERISED,
            geom_point_size = GEOM_POINT_SIZE
        ) +
            theme_customized_embedding()
    })
}) |>
    unlist(recursive = F) |>
    purrr::reduce(`+`) +
    patchwork::plot_layout(ncol = 2, byrow = TRUE) +
    patchwork::plot_annotation(
        theme = ggplot2::theme(plot.margin = ggplot2::margin())
    )

Expression

features_selected <- c(
    "ENSG00000204531_POU5F1",
    "ENSG00000203909_DPPA5",
    "ENSG00000125798_FOXA2",
    "ENSG00000134853_PDGFRA",
    "ENSG00000107485_GATA3",
    "ENSG00000118777_ABCG2"
)

groups_selected_violin <- c(
    "Epiblast", "Primitive Endoderm", "Trophectoderm"
)

# blastoid
cells_violin <- purrr::map(groups_selected_violin, \(x) {
    embedding |>
        dplyr::filter(
            lineage == x
        ) |>
        dplyr::pull(cell)
})
names(cells_violin) <- groups_selected_violin

plot_violin(
    cells = cells_violin,
    features = features_selected,
    matrix_cpm = calc_cpm(matrix_readcount_use),
    x_range_breaks = NULL
) +
    ggplot2::scale_color_manual(
        name = NULL,
        values = scales::hue_pal()(n = 3) |> rev()
    ) +
    ggplot2::scale_fill_manual(
        name = NULL,
        values = scales::hue_pal()(n = 3) |> rev()
    ) +
    theme_customized_violin() +
    ggplot2::theme(
        plot.title = ggplot2::element_text(
            family = "Arial",
            size = 8,
            margin = ggplot2::margin(
                t = 0, r = 0, b = 1, l = 0,
                unit = "mm"
            ),
            hjust = 0.5
        )
    )

R session info

devtools::session_info()

─ Session info ───────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.1.2 (2021-11-01)
 os       macOS Monterey 12.2
 system   aarch64, darwin20.6.0
 ui       unknown
 language (EN)
 collate  en_US.UTF-8
 ctype    en_US.UTF-8
 tz       America/Chicago
 date     2022-02-07
 pandoc   2.17.1.1 @ /opt/homebrew/bin/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────
 package         * version     date (UTC) lib source
 abind             1.4-5       2016-07-21 [1] CRAN (R 4.1.1)
 assertthat        0.2.1       2019-03-21 [1] CRAN (R 4.1.1)
 backports         1.4.1       2021-12-13 [1] CRAN (R 4.1.2)
 bit               4.0.4       2020-08-04 [1] CRAN (R 4.1.1)
 bit64             4.0.5       2020-08-30 [1] CRAN (R 4.1.1)
 brio              1.1.3       2021-11-30 [1] CRAN (R 4.1.2)
 broom             0.7.12      2022-01-28 [1] CRAN (R 4.1.2)
 cachem            1.0.6       2021-08-19 [1] CRAN (R 4.1.1)
 callr             3.7.0       2021-04-20 [1] CRAN (R 4.1.1)
 cellranger        1.1.0       2016-07-27 [1] CRAN (R 4.1.1)
 checkmate         2.0.0       2020-02-06 [1] CRAN (R 4.1.1)
 cli               3.1.1       2022-01-20 [1] CRAN (R 4.1.2)
 cluster           2.1.2       2021-04-17 [2] CRAN (R 4.1.2)
 codetools         0.2-18      2020-11-04 [2] CRAN (R 4.1.2)
 colorspace        2.0-2       2021-06-24 [1] CRAN (R 4.1.1)
 cowplot           1.1.1       2020-12-30 [1] CRAN (R 4.1.1)
 crayon            1.4.2       2021-10-29 [1] CRAN (R 4.1.1)
 data.table        1.14.2      2021-09-27 [1] CRAN (R 4.1.1)
 DBI               1.1.2       2021-12-20 [1] CRAN (R 4.1.2)
 dbplyr            2.1.1       2021-04-06 [1] CRAN (R 4.1.1)
 deldir            1.0-6       2021-10-23 [1] CRAN (R 4.1.1)
 desc              1.4.0       2021-09-28 [1] CRAN (R 4.1.1)
 devtools          2.4.3.9000  2022-01-27 [1] Github (r-lib/devtools@41280ac)
 digest            0.6.29      2021-12-01 [1] CRAN (R 4.1.2)
 dplyr           * 1.0.7.9000  2022-02-04 [1] Github (tidyverse/dplyr@e4cd8c2)
 dtplyr            1.2.1.9000  2022-02-06 [1] Github (tidyverse/dtplyr@aea63d9)
 ellipsis          0.3.2       2021-04-29 [1] CRAN (R 4.1.1)
 evaluate          0.14        2019-05-28 [1] CRAN (R 4.1.1)
 extrafont       * 0.17        2014-12-08 [1] CRAN (R 4.1.1)
 extrafontdb       1.0         2012-06-11 [1] CRAN (R 4.1.1)
 fansi             1.0.2       2022-01-14 [1] CRAN (R 4.1.2)
 farver            2.1.0       2021-02-28 [1] CRAN (R 4.1.1)
 fastmap           1.1.0       2021-01-25 [1] CRAN (R 4.1.1)
 fitdistrplus      1.1-6       2021-09-28 [1] CRAN (R 4.1.1)
 forcats         * 0.5.1.9000  2022-01-27 [1] Github (tidyverse/forcats@b4dade0)
 fs                1.5.2.9000  2022-01-27 [1] Github (r-lib/fs@6d1182f)
 future            1.23.0      2021-10-31 [1] CRAN (R 4.1.1)
 future.apply      1.8.1       2021-08-10 [1] CRAN (R 4.1.1)
 gargle            1.2.0       2021-07-02 [1] CRAN (R 4.1.1)
 generics          0.1.2       2022-01-31 [1] CRAN (R 4.1.2)
 ggplot2         * 3.3.5.9000  2022-01-27 [1] Github (tidyverse/ggplot2@c89c265)
 ggrepel           0.9.1       2021-01-15 [1] CRAN (R 4.1.1)
 ggridges          0.5.3       2021-11-24 [1] Github (wilkelab/ggridges@01db7fc)
 globals           0.14.0      2020-11-22 [1] CRAN (R 4.1.1)
 glue              1.6.1.9000  2022-01-27 [1] Github (tidyverse/glue@3da70df)
 goftest           1.2-3       2021-10-07 [1] CRAN (R 4.1.1)
 googledrive       2.0.0       2021-07-08 [1] CRAN (R 4.1.1)
 googlesheets4     1.0.0       2021-07-21 [1] CRAN (R 4.1.1)
 gridExtra         2.3         2017-09-09 [1] CRAN (R 4.1.1)
 gt                0.3.1       2021-08-07 [1] CRAN (R 4.1.2)
 gtable            0.3.0.9000  2022-01-27 [1] Github (r-lib/gtable@a0bd272)
 haven             2.4.3       2021-08-04 [1] CRAN (R 4.1.1)
 highr             0.9         2021-04-16 [1] CRAN (R 4.1.1)
 hms               1.1.1       2021-09-26 [1] CRAN (R 4.1.1)
 htmltools         0.5.2       2021-08-25 [1] CRAN (R 4.1.1)
 htmlwidgets       1.5.4       2021-09-08 [1] CRAN (R 4.1.1)
 httpuv            1.6.5       2022-01-05 [1] CRAN (R 4.1.2)
 httr              1.4.2       2020-07-20 [1] CRAN (R 4.1.1)
 ica               1.0-2       2018-05-24 [1] CRAN (R 4.1.1)
 igraph            1.2.11.9123 2022-01-27 [1] Github (igraph/rigraph@cfef324)
 irlba             2.3.5       2021-12-06 [1] CRAN (R 4.1.2)
 jsonlite          1.7.3       2022-01-17 [1] CRAN (R 4.1.2)
 KernSmooth        2.23-20     2021-05-03 [2] CRAN (R 4.1.2)
 knitr             1.37        2021-12-16 [1] CRAN (R 4.1.2)
 labeling          0.4.2       2020-10-20 [1] CRAN (R 4.1.1)
 later             1.3.0       2021-08-18 [1] CRAN (R 4.1.1)
 lattice           0.20-45     2021-09-22 [2] CRAN (R 4.1.2)
 lazyeval          0.2.2       2019-03-15 [1] CRAN (R 4.1.1)
 leiden            0.3.9       2021-07-27 [1] CRAN (R 4.1.1)
 lifecycle         1.0.1       2021-09-24 [1] CRAN (R 4.1.1)
 listenv           0.8.0       2019-12-05 [1] CRAN (R 4.1.1)
 lmtest            0.9-39      2021-11-07 [1] CRAN (R 4.1.2)
 lubridate         1.8.0       2022-01-27 [1] Github (tidyverse/lubridate@566590f)
 magrittr          2.0.2       2022-01-26 [1] CRAN (R 4.1.2)
 MASS              7.3-55      2022-01-13 [2] CRAN (R 4.1.2)
 Matrix          * 1.4-0       2021-12-08 [2] CRAN (R 4.1.2)
 matrixStats       0.61.0      2021-09-17 [1] CRAN (R 4.1.1)
 memoise           2.0.1       2021-11-26 [1] CRAN (R 4.1.2)
 mgcv              1.8-38      2021-10-06 [2] CRAN (R 4.1.2)
 mime              0.12        2021-09-28 [1] CRAN (R 4.1.1)
 miniUI            0.1.1.1     2018-05-18 [1] CRAN (R 4.1.1)
 modelr            0.1.8.9000  2022-01-27 [1] Github (tidyverse/modelr@16168e0)
 munsell           0.5.0       2018-06-12 [1] CRAN (R 4.1.1)
 nlme              3.1-155     2022-01-13 [2] CRAN (R 4.1.2)
 paletteer         1.4.0.9000  2022-02-01 [1] Github (EmilHvitfeldt/paletteer@c814e0d)
 parallelly        1.30.0      2021-12-17 [1] CRAN (R 4.1.2)
 patchwork       * 1.1.0.9000  2022-01-27 [1] Github (thomasp85/patchwork@79223d3)
 pbapply           1.5-0       2021-09-16 [1] CRAN (R 4.1.1)
 pillar            1.7.0       2022-02-01 [1] CRAN (R 4.1.2)
 pkgbuild          1.3.1       2021-12-20 [1] CRAN (R 4.1.2)
 pkgconfig         2.0.3       2019-09-22 [1] CRAN (R 4.1.1)
 pkgload           1.2.4       2021-11-30 [1] CRAN (R 4.1.2)
 plotly            4.10.0      2021-10-09 [1] CRAN (R 4.1.1)
 plyr              1.8.6       2020-03-03 [1] CRAN (R 4.1.1)
 png               0.1-7       2013-12-03 [1] CRAN (R 4.1.1)
 polyclip          1.10-0      2019-03-14 [1] CRAN (R 4.1.1)
 prettyunits       1.1.1.9000  2022-01-27 [1] Github (r-lib/prettyunits@8706d89)
 prismatic         1.1.0       2021-10-17 [1] CRAN (R 4.1.2)
 processx          3.5.2       2021-04-30 [1] CRAN (R 4.1.1)
 promises          1.2.0.1     2021-02-11 [1] CRAN (R 4.1.1)
 ps                1.6.0       2021-02-28 [1] CRAN (R 4.1.1)
 purrr           * 0.3.4.9000  2022-01-27 [1] Github (tidyverse/purrr@5aca9df)
 R.cache           0.15.0      2021-04-30 [1] CRAN (R 4.1.1)
 R.methodsS3       1.8.1       2020-08-26 [1] CRAN (R 4.1.1)
 R.oo              1.24.0      2020-08-26 [1] CRAN (R 4.1.1)
 R.utils           2.11.0      2021-09-26 [1] CRAN (R 4.1.1)
 R6                2.5.1.9000  2022-01-27 [1] Github (r-lib/R6@1b05b89)
 ragg              1.2.1.9000  2022-01-27 [1] Github (r-lib/ragg@c68c666)
 RANN              2.6.1       2019-01-08 [1] CRAN (R 4.1.1)
 RColorBrewer      1.1-2       2014-12-07 [1] CRAN (R 4.1.1)
 Rcpp              1.0.8       2022-01-13 [1] CRAN (R 4.1.2)
 RcppAnnoy         0.0.19      2021-07-30 [1] CRAN (R 4.1.1)
 readr           * 2.1.2.9000  2022-02-04 [1] Github (tidyverse/readr@8296ae5)
 readxl            1.3.1.9000  2022-01-27 [1] Github (tidyverse/readxl@2ccb82c)
 rematch2          2.1.2       2020-05-01 [1] CRAN (R 4.1.1)
 remotes           2.4.2       2022-01-27 [1] Github (r-lib/remotes@9355549)
 reprex            2.0.1       2021-08-05 [1] CRAN (R 4.1.1)
 reshape2          1.4.4       2020-04-09 [1] CRAN (R 4.1.1)
 reticulate        1.24        2022-01-26 [1] CRAN (R 4.1.2)
 rlang             1.0.1.9000  2022-02-03 [1] Github (r-lib/rlang@f4c6d14)
 rmarkdown         2.11.15     2022-02-04 [1] Github (rstudio/rmarkdown@ede2d7c)
 ROCR              1.0-11      2020-05-02 [1] CRAN (R 4.1.1)
 rpart             4.1.16      2022-01-24 [2] CRAN (R 4.1.2)
 rprojroot         2.0.2       2020-11-15 [1] CRAN (R 4.1.1)
 rstudioapi        0.13.0-9000 2022-02-04 [1] Github (rstudio/rstudioapi@942cebe)
 Rtsne             0.16        2022-01-27 [1] Github (jkrijthe/Rtsne@52a07c2)
 Rttf2pt1          1.3.9       2021-07-22 [1] CRAN (R 4.1.1)
 rvest             1.0.2       2021-10-16 [1] CRAN (R 4.1.1)
 sass              0.4.0       2021-05-12 [1] CRAN (R 4.1.1)
 scales            1.1.1.9000  2022-01-29 [1] Github (r-lib/scales@13b01f0)
 scattermore       0.7         2020-11-24 [1] CRAN (R 4.1.1)
 sctransform       0.3.3       2022-01-13 [1] Github (satijalab/sctransform@e9e52a4)
 sessioninfo       1.2.2       2021-12-06 [1] CRAN (R 4.1.2)
 Seurat            4.1.0       2022-01-15 [1] Github (satijalab/seurat@eaaf675)
 SeuratObject      4.0.4       2021-11-23 [1] CRAN (R 4.1.2)
 shiny             1.7.1       2021-10-02 [1] CRAN (R 4.1.1)
 spatstat.core     2.3-2       2021-11-26 [1] CRAN (R 4.1.2)
 spatstat.data     2.1-2       2021-12-17 [1] CRAN (R 4.1.2)
 spatstat.geom     2.3-1       2021-12-10 [1] CRAN (R 4.1.2)
 spatstat.sparse   2.1-0       2021-12-17 [1] CRAN (R 4.1.2)
 spatstat.utils    2.3-0       2021-12-12 [1] CRAN (R 4.1.2)
 stringi           1.7.6       2021-11-29 [1] CRAN (R 4.1.2)
 stringr         * 1.4.0.9000  2022-01-27 [1] Github (tidyverse/stringr@85f6140)
 styler          * 1.6.2.9000  2022-02-06 [1] Github (r-lib/styler@c43c99d)
 survival          3.2-13      2021-08-24 [2] CRAN (R 4.1.2)
 systemfonts       1.0.3       2021-10-13 [1] CRAN (R 4.1.2)
 tensor            1.5         2012-05-05 [1] CRAN (R 4.1.1)
 testthat          3.1.2.9000  2022-01-27 [1] Github (r-lib/testthat@54b9db2)
 textshaping       0.3.6       2021-10-13 [1] CRAN (R 4.1.1)
 tibble          * 3.1.6.9001  2022-02-04 [1] Github (tidyverse/tibble@281c74d)
 tidyr           * 1.2.0.9000  2022-02-03 [1] Github (tidyverse/tidyr@5a1ea09)
 tidyselect        1.1.1       2021-04-30 [1] CRAN (R 4.1.1)
 tidyverse       * 1.3.1.9000  2022-01-27 [1] Github (tidyverse/tidyverse@6186fbf)
 tzdb              0.2.0       2021-10-27 [1] CRAN (R 4.1.1)
 usethis           2.1.5.9000  2022-01-27 [1] Github (r-lib/usethis@57b109a)
 utf8              1.2.2       2021-07-24 [1] CRAN (R 4.1.1)
 uwot              0.1.11.9000 2022-02-06 [1] Github (jlmelville/uwot@2be098a)
 vctrs             0.3.8       2021-04-29 [1] CRAN (R 4.1.1)
 viridisLite       0.4.0       2021-04-13 [1] CRAN (R 4.1.1)
 vroom             1.5.7.9000  2022-01-30 [1] Github (r-lib/vroom@d91f2aa)
 withr             2.4.3       2021-11-30 [1] CRAN (R 4.1.2)
 xfun              0.29        2021-12-14 [1] CRAN (R 4.1.2)
 xml2              1.3.3       2021-11-30 [1] CRAN (R 4.1.2)
 xtable            1.8-4       2019-04-21 [1] CRAN (R 4.1.1)
 yaml              2.2.2       2022-01-25 [1] CRAN (R 4.1.2)
 zoo               1.8-9       2021-03-09 [1] CRAN (R 4.1.1)

 [1] /opt/homebrew/lib/R/4.1/site-library
 [2] /opt/homebrew/Cellar/r/4.1.2/lib/R/library

─ Python configuration ───────────────────────────────────────────────────────
 python:         /Users/jialei/.pyenv/shims/python
 libpython:      /Users/jialei/.pyenv/versions/miniforge3-4.10.1-5/lib/libpython3.9.dylib
 pythonhome:     /Users/jialei/.pyenv/versions/miniforge3-4.10.1-5:/Users/jialei/.pyenv/versions/miniforge3-4.10.1-5
 version:        3.9.5 | packaged by conda-forge | (default, Jun 19 2021, 00:24:55)  [Clang 11.1.0 ]
 numpy:          /Users/jialei/.pyenv/versions/miniforge3-4.10.1-5/lib/python3.9/site-packages/numpy
 numpy_version:  1.20.3
 anndata:        /Users/jialei/.pyenv/versions/miniforge3-4.10.1-5/lib/python3.9/site-packages/anndata
 
 NOTE: Python version was forced by RETICULATE_PYTHON

──────────────────────────────────────────────────────────────────────────────