In jussi-kupari/kumar-et-al-2022: Rakesh NtsCre project

# Load libraries and some functions
source(here::here("R", "libraries.R"))
R.utils::sourceDirectory(here("R", "functions"))

# Load SS2 data and annotation
FosTRAP <-
  read.delim(here("data", "raw", "SS420_counts.txt"), sep = " ") %>%
  as.matrix()

FosTRAP_annotation <-
  read.delim(
    here("data", "reference", "SS420_cells_annot_for_Yizhou.txt"),
    sep = " "
  ) %>%
  rename(cell_type = x) %>%
  rownames_to_column("cell")

# Create Seurat object
FosTRAP %<>%
  CreateSeuratObject(project = "FosTRAP", min.features = 500) %>%
  NormalizeData() %>%
  ScaleData(features = rownames(FosTRAP)) %>%
  FindVariableFeatures() %>%
  RunPCA()

# Plot some QC metrics
VlnPlot(FosTRAP, c("nFeature_RNA", "nCount_RNA"))

ElbowPlot(FosTRAP, ndims = 50)

# Cluster cells
FosTRAP %<>%
  RunUMAP(reduction = "pca", dims = 1:10, n.neighbors = 50) %>%
  FindNeighbors(reduction = "pca", dims = 1:10) %>%
  FindClusters(resolution = 0.5) %>%
  BuildClusterTree(
    features = NULL, dims = 1:10, graph = NULL,
    reorder = TRUE, reorder.numeric = TRUE, verbose = TRUE
  )

annot <-
  rownames([e-mail]) %>%
  as_tibble() %>%
  left_join(FosTRAP_annotation, by = c("value" = "cell")) %>%
  rename(cell = value)

[e-mail] %<>%
  rownames_to_column("cell") %>%
  left_join(FosTRAP_annotation, by = "cell") %>%
  column_to_rownames("cell")

FeaturePlot(FosTRAP, "Rbfox3", label = TRUE) +
  VlnPlot(FosTRAP, c("Rbfox3"))+ NoLegend()

# Remove nonneuronal cells
FosTRAP %<>%
  subset(idents = 3:5) %>%
  NormalizeData() %>%
  ScaleData(features = rownames(FosTRAP)) %>%
  FindVariableFeatures() %>%
  RunPCA() %>%
  RunUMAP(reduction = "pca", dims = 1:10, n.neighbors = 50) %>%
  FindNeighbors(reduction = "pca", dims = 1:10) %>%
  FindClusters(resolution = 0.5) %>%
  BuildClusterTree(
    features = NULL, dims = 1:10, graph = NULL,
    reorder = TRUE, reorder.numeric = TRUE, verbose = TRUE
  )

FeaturePlot(FosTRAP, "Rbfox3", label = TRUE) +
  VlnPlot(FosTRAP, "Rbfox3") + NoLegend()

# Load and prepare Häring et al. data
haring <- readRDS(here("data", "reference", "haring.rds"))

haring %>%
  set_idents("standard_names") %>%
  NormalizeData() %>%
  ScaleData(features = rownames(haring)) %>%
  FindVariableFeatures() %>%
  RunPCA()

# Create classifier with scPred using MDA as model
haring %<>%
  getFeatureSpace("standard_names") %>%
  trainModel(model = "mda")

# Check model performance
get_scpred(haring)

# Predict labels for FosTRAP neurons
FosTRAP %<>% scPredict(haring, threshold = 0.75)

# Max score distributions, red vertical line (0.55) is the assignment boundary
[e-mail] %>%
  ggplot(aes(scpred_max)) +
  geom_histogram(
    bins = 50,
    color = "white",
    fill = "darkblue",
    alpha = 0.7
  ) +
  geom_vline(xintercept = 0.75, color = "red", size = 1) +
  cowplot::theme_cowplot() +
  scale_y_continuous(expand = expansion(mult = c(0, .1))) +
  labs(
    title = "",
    x = "Max prediction score",
    y = "Number of cells"
  )

# Discard unassigned neurons and cell types with less than 3 cells
FosTRAP %<>% set_idents("scpred_prediction")

FosTRAP %<>%
  subset(idents = as.character(unique(Idents(FosTRAP))) %>%
    keep(~ .x %not_in% "unassigned"))

clusters_to_keep <-
  [e-mail] %>% 
  count(scpred_prediction, sort = TRUE) %>% 
  mutate(keep = n >= 3) %>% 
  select(scpred_prediction, keep)

[e-mail] %<>%
  rownames_to_column("cell") %>% 
  left_join(clusters_to_keep, by = "scpred_prediction") %>% 
  column_to_rownames("cell") 

FosTRAP %<>% 
  set_idents("keep") %>% 
  subset(idents = TRUE) %>% 
  set_idents("scpred_prediction")

# scPred score heatmap
[e-mail] %>%
  select(
    starts_with("scpred_"),
    -c(scpred_max, scpred_prediction, scpred_no_rejection)
  ) %>%
  rename_with(~ str_remove(.x, "scpred_")) %>%
  t() %>%
  pheatmap::pheatmap(
    cluster_rows = FALSE,
    cluster_cols = TRUE,
    treeheight_col = 0,
    show_colnames = FALSE,
    color = viridisLite::viridis(10),
    fontsize = 12,
    title = ""
  )

# Plot UMAP with predicted identites
DimPlot(FosTRAP, label = TRUE, repel = TRUE) 

# Plot neurotransmitter type (Glut/Gaba)
(FeaturePlot(FosTRAP, 
             c("Slc17a6", "Slc32a1"), 
             label = TRUE, 
             pt.size = 1, 
             repel = TRUE)) /
  (VlnPlot(FosTRAP, c("Slc17a6", "Slc32a1")))

# Plot number of neurons per celltype
[e-mail] %>%
  count(scpred_prediction, sort = TRUE) %>%
  mutate(prop = n / sum(n)) %>%
  ggplot(aes(fct_inorder(scpred_prediction), prop)) +
  geom_col(fill = "darkblue", alpha = 0.7) +
  cowplot::theme_cowplot() +
  scale_y_continuous(expand = expansion(mult = c(0, .1))) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
  labs(x = "Predicted cell type", y = "Proportion of cells")

# Tabulate cell counts per predicted identity
[e-mail] %>% 
  count(scpred_prediction, sort = TRUE) %>% 
  knitr::kable()

# Get markers for Häring et al. data
haring_markers <-
  haring %>%
  set_idents("standard_names") %>%
  FindAllMarkers(only.pos = TRUE)

# Get top genes for Häring et al. cell types
top_genes_haring <-
  haring_markers %>%
  filter(cluster %in% unique([e-mail]$scpred_prediction)) %>%
  group_by(cluster) %>%
  slice_min(p_val_adj, n = 5, with_ties = FALSE)

# Re-order FosTRAP data (with assigned IDs) for plotting
order <-
  c(
    "Glut1",
    "Glut2",
    "Glut3",
    "Glut4",
    "Glut6",
    "Glut7",
    "Glut8",
    "Gaba10",
    "Gaba12",
    "Gaba13"
  )

Idents(FosTRAP) <- factor(x = Idents(FosTRAP), levels = order)

# Dotplot of top Häring markers in FosTRAP neurons
FosTRAP %>%
  DotPlot(
    features = unique(c("Slc32a1", "Slc17a6", "Prkcg", top_genes_haring$gene))
  ) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  coord_flip()

# Plot Nts and Prkcq in FosTRAP neurons
VlnPlot(FosTRAP, c("Nts", "Prkcg"), pt.size = 0.1)

sessionInfo()

jussi-kupari/kumar-et-al-2022 documentation built on June 15, 2022, 9:37 p.m.