all_times <- list()  # store the time for each chunk
knitr::knit_hooks$set(time_it = local({
  now <- NULL
  function(before, options) {
    if (before) {
      now <<- Sys.time()
    } else {
      res <- difftime(Sys.time(), now, units = "secs")
      all_times[[options$label]] <<- res
    }
  }
}))
knitr::opts_chunk$set(
  tidy = TRUE,
  tidy.opts = list(width.cutoff = 95),
  message = FALSE,
  warning = FALSE,
  time_it = TRUE
)

Marker Identification & Cluster Annotation

scCustomize has several helper functions to assist with identification of marker genes and annotation of clusters. Load packages & Data

# Load Packages
library(ggplot2)
library(dplyr)
library(magrittr)
library(Seurat)
library(scCustomize)
library(qs)

# Load example dataset for tutorial
pbmc <- pbmc3k.SeuratData::pbmc3k.final

Helper functions available are:

We will use pbmc3K dataset from SeuratData package as example for this vignette.

# Update pbmc check
pbmc <- UpdateSeuratObject(pbmc)

Percent Difference in Expression

If we run FindAllMarkers() like normal, this is the output.

# Run basic FindAllMarkers DE test
all_markers <- FindAllMarkers(object = pbmc, verbose = F)

all_markers <- all_markers %>%
  group_by(cluster) %>%  
  arrange(desc(avg_log2FC)) %>%   
  arrange(cluster)
head(all_markers, 5) %>%
  kableExtra::kbl(row.names = TRUE) %>%
  kableExtra::kable_styling(bootstrap_options = c("bordered", "condensed", "responsive", "striped")) 

While Seurat::FindAllMarkers()returns the percent of cells in identity 1 (pct.1) and identity 2 (pct.2) that express a marker it can be helpful to view the difference in these two measures in addition to the values alone.

scCustomize contains helper function: Add_Pct_Diff() to add the percent difference between two clusters. Add_Pct_Diff can be used with any output from Seurat, SeuratWrappers, or other DE tests and can be piped with the function itself.

Use Add_Pct_Diff piped with the marker function.

# run piped with Seurat command
all_markers_pct <- FindAllMarkers(object = pbmc, verbose = F) %>% 
  Add_Pct_Diff()

all_markers_pct <- all_markers_pct %>% 
  group_by(cluster) %>%
  arrange(desc(avg_log2FC)) %>% 
  arrange(cluster)
head(all_markers_pct, 5) %>%
  kableExtra::kbl(row.names = TRUE) %>%
  kableExtra::kable_styling(bootstrap_options = c("bordered", "condensed", "responsive", "striped")) 

Use Add_Pct_Diff post-hoc

Add_Pct_Diff can also simply be run after generating marker data.frame:

all_markers_pct <- Add_Pct_Diff(marker_dataframe = all_markers)

Use with non-Seurat data.frames

Add_Pct_Diff can be used with any data.frame as long as the column names that specify percent expressed in the 2 groups are provided

all_markers_pct <- Add_Pct_Diff(marker_dataframe = all_markers, pct.1_name = "percent1", pct.2_name = "percent2")

Extract the top N marker genes per cluster, identity, treatment etc

scCustomize contains function Extract_Top_Markers to get the top N genes from each cluster (or other group of choice). It has the option to return either a vector of genes or a filtered data.frame.

Returns a Named Vector (Default)

# By default it returns a named vector of genes (each gene named for cluster it was associated with)
top_5 <- Extract_Top_Markers(marker_dataframe = all_markers, num_genes = 5, rank_by = "avg_log2FC")
head(top_5, 10) 

Return data.frame (Optional)

top_5_df <- Extract_Top_Markers(marker_dataframe = all_markers, num_genes = 5, data_frame = TRUE, rank_by = "avg_log2FC")

head(top_5_df, 10) 
top_5_df <- Extract_Top_Markers(marker_dataframe = all_markers, num_genes = 5, data_frame = TRUE, rank_by = "avg_log2FC")

head(top_5_df, 10) %>%
  kableExtra::kbl(row.names = TRUE) %>%
  kableExtra::kable_styling(bootstrap_options = c("bordered", "condensed", "responsive", "striped")) 

Use with non-Seurat data.frames

Extract_Top_Markers uses the data.frame output of Seurat::FindAllMarkers() by default but can be used with any data.frame by supplying few extra parameters. Let's use the function on the output from presto package and order the results by auc instead of log2FC.

# Read in pre-computed data.frame results
presto_markers <- qread("assets/presto_markers.qs")

# Extract top markers
top_5_presto <- Extract_Top_Markers(marker_dataframe = presto_markers, num_genes = 5, group_by = "cluster", rank_by = "auc", gene_column = "gene")
head(top_5_presto, 10)

Unique Vectors

When returning an unnamed vector you may want to return a unique list without any duplicate features. This is often true when supplying a list of features to plotting function like DoHeatmap() or Cluster_DotPlot(). You can do this by setting parameter make_unique = TRUE.

# Default parameters
top_5 <- Extract_Top_Markers(marker_dataframe = all_markers, num_genes = 5, rank_by = "avg_log2FC", named_vector = FALSE)

any(duplicated(x = top_5))

# Set `make_unique = TRUE`
top_5_unique <- Extract_Top_Markers(marker_dataframe = all_markers, num_genes = 5, rank_by = "avg_log2FC", named_vector = FALSE, make_unique = TRUE)

any(duplicated(x = top_5_unique))

Pairing Extract_Top_Markers with Iterate_FeaturePlot_scCustom()

The advantage of exporting a named vector from Extract_Top_Markers is that Iterate_FeaturePlot_scCustom can detect the names and either add those to the title of the plot (if single_pdf = TRUE) or the file names if FALSE.

Iterate_FeaturePlot_scCustom(seurat_object = pbmc, gene_list = top_5_presto, single_pdf = T)

Example plot showing the gene name and additional name of the plot.

Create simple annotation files

Now for a simple helper we can use Create_Cluster_Annotation_File to create simple csv file to use for annotation.

Create_Cluster_Annotation_File(file_path = "data/", file_name = "ExpXX_Round01_Annotation")

Here is an example that has already been filled out (using Excel, Numbers, or other editor program).

annotation_example <- read.csv(file = "assets/ExpXX_Round01_Annotation_Completed.csv", stringsAsFactors = FALSE)
head(annotation_example, 5) %>%
  kableExtra::kbl(row.names = TRUE) %>%
  kableExtra::kable_styling(bootstrap_options = c("bordered", "condensed", "responsive", "striped")) 

Extract annotation information

There is also complementary function to pull information back out of the this (or other) csv files or data.frames regarding clustering annotation.

annotation_info <- Pull_Cluster_Annotation(annotation = "assets/ExpXX_Round01_Annotation_Completed.csv")

By default the function expects data.frame or file in same format as Create_Cluster_Annotation_File but has optional parameters that can be adjusted to file/data.frame with different column headers (see cluster_name_col and cell_type_col).

Pull_Cluster_Annotation results

The result of Pull_Cluster_Annotation is a list that contains:

Example output of `Pull_Cluster_Annotation

Renaming Clusters

Finally can use Rename_Clusters to easily rename and set active.ident in Seurat Object. This is simple wrapper around Seurat commands but adds additional checks/warnings and is simple one-liner.

obj_renamed <- Rename_Clusters(seurat_object = obj, new_idents = annotation_info$new_cluster_idents)

Rename_Clusters can take input from either Pull_Cluster_Annotation "new_cluster_idents" or any correctly ordered vector of new idents.



samuel-marsh/scCustomize documentation built on Dec. 20, 2024, 7:41 a.m.