## ---- echo=FALSE, include=FALSE----------------------------------------------- library(knitr) knitr::opts_chunk$set( cache=TRUE, warning=FALSE, message=FALSE, cache.lazy=FALSE ) ## ---- eval=FALSE-------------------------------------------------------------- # if (!requireNamespace("BiocManager", quietly=TRUE)) # install.packages("BiocManager") # # BiocManager::install("tidySingleCellExperiment") ## ----message=FALSE------------------------------------------------------------ # Bioconductor single-cell packages library(scater) library(scran) library(SingleR) library(SingleCellSignalR) # Tidyverse-compatible packages library(ggplot2) library(purrr) library(tidyHeatmap) # Both library(tidySingleCellExperiment) ## ----------------------------------------------------------------------------- pbmc_small_tidy <- tidySingleCellExperiment::pbmc_small ## ----------------------------------------------------------------------------- pbmc_small_tidy ## ----------------------------------------------------------------------------- assay(pbmc_small_tidy) ## ----------------------------------------------------------------------------- pbmc_small_tidy$file[1:5] ## ----------------------------------------------------------------------------- # Create sample column pbmc_small_polished <- pbmc_small_tidy %>% extract(file, "sample", "../data/([a-z0-9]+)/outs.+", remove=FALSE) # Reorder to have sample column up front pbmc_small_polished %>% select(sample, everything()) ## ----------------------------------------------------------------------------- # Use colourblind-friendly colours friendly_cols <- dittoSeq::dittoColors() # Set theme custom_theme <- list( scale_fill_manual(values=friendly_cols), scale_color_manual(values=friendly_cols), theme_bw() + theme( panel.border=element_blank(), axis.line=element_line(), panel.grid.major=element_line(size=0.2), panel.grid.minor=element_line(size=0.1), text=element_text(size=12), legend.position="bottom", aspect.ratio=1, strip.background=element_blank(), axis.title.x=element_text(margin=margin(t=10, r=10, b=10, l=10)), axis.title.y=element_text(margin=margin(t=10, r=10, b=10, l=10)) ) ) ## ----plot1-------------------------------------------------------------------- pbmc_small_polished %>% tidySingleCellExperiment::ggplot(aes(nFeature_RNA, fill=groups)) + geom_histogram() + custom_theme ## ----plot2-------------------------------------------------------------------- pbmc_small_polished %>% tidySingleCellExperiment::ggplot(aes(groups, nCount_RNA, fill=groups)) + geom_boxplot(outlier.shape=NA) + geom_jitter(width=0.1) + custom_theme ## ----------------------------------------------------------------------------- pbmc_small_polished %>% join_features(features=c("HLA-DRA", "LYZ")) %>% ggplot(aes(groups, .abundance_counts + 1, fill=groups)) + geom_boxplot(outlier.shape=NA) + geom_jitter(aes(size=nCount_RNA), alpha=0.5, width=0.2) + scale_y_log10() + custom_theme ## ----preprocess--------------------------------------------------------------- # Identify variable genes with scran variable_genes <- pbmc_small_polished %>% modelGeneVar() %>% getTopHVGs(prop=0.1) # Perform PCA with scater pbmc_small_pca <- pbmc_small_polished %>% runPCA(subset_row=variable_genes) pbmc_small_pca ## ----pc_plot------------------------------------------------------------------ # Create pairs plot with GGally pbmc_small_pca %>% as_tibble() %>% select(contains("PC"), everything()) %>% GGally::ggpairs(columns=1:5, ggplot2::aes(colour=groups)) + custom_theme ## ----cluster------------------------------------------------------------------ pbmc_small_cluster <- pbmc_small_pca # Assign clusters to the 'colLabels' of the SingleCellExperiment object colLabels(pbmc_small_cluster) <- pbmc_small_pca %>% buildSNNGraph(use.dimred="PCA") %>% igraph::cluster_walktrap() %$% membership %>% as.factor() # Reorder columns pbmc_small_cluster %>% select(label, everything()) ## ----cluster count------------------------------------------------------------ # Count number of cells for each cluster per group pbmc_small_cluster %>% tidySingleCellExperiment::count(groups, label) ## ----------------------------------------------------------------------------- # Identify top 10 markers per cluster marker_genes <- pbmc_small_cluster %>% findMarkers(groups=pbmc_small_cluster$label) %>% as.list() %>% map(~ .x %>% head(10) %>% rownames()) %>% unlist() # Plot heatmap pbmc_small_cluster %>% join_features(features=marker_genes) %>% group_by(label) %>% heatmap(.feature, .cell, .abundance_counts, .scale="column") ## ----umap--------------------------------------------------------------------- pbmc_small_UMAP <- pbmc_small_cluster %>% runUMAP(ncomponents=3) ## ----umap plot, eval=FALSE---------------------------------------------------- # pbmc_small_UMAP %>% # plot_ly( # x=~`UMAP1`, # y=~`UMAP2`, # z=~`UMAP3`, # color=~label, # colors=friendly_cols[1:4] # ) ## ----eval=FALSE--------------------------------------------------------------- # # Get cell type reference data # blueprint <- celldex::BlueprintEncodeData() # # # Infer cell identities # cell_type_df <- # # assays(pbmc_small_UMAP)$logcounts %>% # Matrix::Matrix(sparse = TRUE) %>% # SingleR::SingleR( # ref = blueprint, # labels = blueprint$label.main, # method = "single" # ) %>% # as.data.frame() %>% # as_tibble(rownames="cell") %>% # select(cell, first.labels) ## ----------------------------------------------------------------------------- # Join UMAP and cell type info pbmc_small_cell_type <- pbmc_small_UMAP %>% left_join(cell_type_df, by="cell") # Reorder columns pbmc_small_cell_type %>% tidySingleCellExperiment::select(cell, first.labels, everything()) ## ----------------------------------------------------------------------------- # Count number of cells for each cell type per cluster pbmc_small_cell_type %>% count(label, first.labels) ## ----------------------------------------------------------------------------- pbmc_small_cell_type %>% # Reshape and add classifier column pivot_longer( cols=c(label, first.labels), names_to="classifier", values_to="label" ) %>% # UMAP plots for cell type and cluster ggplot(aes(UMAP1, UMAP2, color=label)) + geom_point() + facet_wrap(~classifier) + custom_theme ## ----------------------------------------------------------------------------- pbmc_small_cell_type %>% # Add some mitochondrial abundance values mutate(mitochondrial=rnorm(dplyr::n())) %>% # Plot correlation join_features(features=c("CST3", "LYZ"), shape="wide") %>% ggplot(aes(CST3 + 1, LYZ + 1, color=groups, size=mitochondrial)) + geom_point() + facet_wrap(~first.labels, scales="free") + scale_x_log10() + scale_y_log10() + custom_theme ## ----------------------------------------------------------------------------- pbmc_small_nested <- pbmc_small_cell_type %>% filter(first.labels != "Erythrocytes") %>% mutate(cell_class=dplyr::if_else(`first.labels` %in% c("Macrophages", "Monocytes"), "myeloid", "lymphoid")) %>% nest(data=-cell_class) pbmc_small_nested ## ----warning=FALSE------------------------------------------------------------ pbmc_small_nested_reanalysed <- pbmc_small_nested %>% mutate(data=map( data, ~ { .x <- runPCA(.x, subset_row=variable_genes) variable_genes <- .x %>% modelGeneVar() %>% getTopHVGs(prop=0.3) colLabels(.x) <- .x %>% buildSNNGraph(use.dimred="PCA") %>% igraph::cluster_walktrap() %$% membership %>% as.factor() .x %>% runUMAP(ncomponents=3) } )) pbmc_small_nested_reanalysed ## ----------------------------------------------------------------------------- pbmc_small_nested_reanalysed %>% # Convert to tibble otherwise SingleCellExperiment drops reduced dimensions when unifying data sets. mutate(data=map(data, ~ .x %>% as_tibble())) %>% unnest(data) %>% # Define unique clusters unite("cluster", c(cell_class, label), remove=FALSE) %>% # Plotting ggplot(aes(UMAP1, UMAP2, color=cluster)) + geom_point() + facet_wrap(~cell_class) + custom_theme ## ---- eval=FALSE-------------------------------------------------------------- # pbmc_small_nested_interactions <- # pbmc_small_nested_reanalysed %>% # # # Unnest based on cell category # unnest(data) %>% # # # Create unambiguous clusters # mutate(integrated_clusters=first.labels %>% as.factor() %>% as.integer()) %>% # # # Nest based on sample # tidySingleCellExperiment::nest(data=-sample) %>% # tidySingleCellExperiment::mutate(interactions=map(data, ~ { # # # Produce variables. Yuck! # cluster <- colData(.x)$integrated_clusters # data <- data.frame(assays(.x) %>% as.list() %>% .[[1]] %>% as.matrix()) # # # Ligand/Receptor analysis using SingleCellSignalR # data %>% # cell_signaling(genes=rownames(data), cluster=cluster) %>% # inter_network(data=data, signal=., genes=rownames(data), cluster=cluster) %$% # `individual-networks` %>% # map_dfr(~ bind_rows(as_tibble(.x))) # })) # # pbmc_small_nested_interactions %>% # select(-data) %>% # unnest(interactions) ## ----------------------------------------------------------------------------- tidySingleCellExperiment::pbmc_small_nested_interactions ## ----------------------------------------------------------------------------- sessionInfo()