Load libraries

suppressPackageStartupMessages({
  library(BiocStyle)
  library(tidyverse)
  library(here)
  library(glue)
  library(patchwork)
  library(paletteer)
  library(readxl)
  library(ggupset)
  library(stringr)
})

Load Data

files <- list.files(here("data",
                         paste0("C133_Neeland_batch", 0:6),
                         "data", 
                         "SCEs"), 
                      pattern = "doublets_called",
                      full.names = TRUE)
  
sceLst <- sapply(files, function(fn){
    readRDS(file = fn)
})

sceLst

$`/Users/maksimovicjovana/Work/Projects/MCRI/melanie.neeland/paed-inflammation-CITEseq/data/C133_Neeland_batch0/data/SCEs/C133_Neeland_batch0.doublets_called.SCE.rds`
class: SingleCellExperiment 
dim: 33538 26939 
metadata(0):
assays(1): counts
rownames(33538): ENSG00000243485 ENSG00000237613 ... ENSG00000277475
  ENSG00000268674
rowData names(20): ID Symbol ... is_mito is_pseudogene
colnames(26939): 1_AAACCCAAGCTAGTTC-1 1_AAACCCACAGTCGCTG-1 ...
  4_TTTGTTGTCTAGTACG-1 4_TTTGTTGTCTCGAACA-1
colData names(25): Barcode Capture ... scDblFinder.cxds_score batch
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):

$`/Users/maksimovicjovana/Work/Projects/MCRI/melanie.neeland/paed-inflammation-CITEseq/data/C133_Neeland_batch1/data/SCEs/C133_Neeland_batch1.doublets_called.SCE.rds`
class: SingleCellExperiment 
dim: 36601 21970 
metadata(0):
assays(1): counts
rownames(36601): ENSG00000243485 ENSG00000237613 ... ENSG00000278817
  ENSG00000277196
rowData names(20): ID Symbol ... is_mito is_pseudogene
colnames(21970): 1_AAACCCACACTTCCTG-1 1_AAACCCACAGACAAAT-1 ...
  2_TTTGTTGTCATTGGTG-1 2_TTTGTTGTCGATGGAG-1
colData names(32): Barcode Capture ... scDblFinder.cxds_score batch
reducedDimNames(0):
mainExpName: NULL
altExpNames(2): HTO ADT

$`/Users/maksimovicjovana/Work/Projects/MCRI/melanie.neeland/paed-inflammation-CITEseq/data/C133_Neeland_batch2/data/SCEs/C133_Neeland_batch2.doublets_called.SCE.rds`
class: SingleCellExperiment 
dim: 36601 45711 
metadata(0):
assays(1): counts
rownames(36601): ENSG00000243485 ENSG00000237613 ... ENSG00000278817
  ENSG00000277196
rowData names(20): ID Symbol ... is_mito is_pseudogene
colnames(45711): 1_AAACCCAAGACCTGGA-1 1_AAACCCAAGACTGTTC-1 ...
  2_TTTGTTGTCTCATGGA-1 2_TTTGTTGTCTCCAAGA-1
colData names(32): Barcode Capture ... scDblFinder.cxds_score batch
reducedDimNames(0):
mainExpName: NULL
altExpNames(2): HTO ADT

$`/Users/maksimovicjovana/Work/Projects/MCRI/melanie.neeland/paed-inflammation-CITEseq/data/C133_Neeland_batch3/data/SCEs/C133_Neeland_batch3.doublets_called.SCE.rds`
class: SingleCellExperiment 
dim: 36601 58038 
metadata(0):
assays(1): counts
rownames(36601): ENSG00000243485 ENSG00000237613 ... ENSG00000278817
  ENSG00000277196
rowData names(20): ID Symbol ... is_mito is_pseudogene
colnames(58038): 1_AAACCCAAGCAGCACA-1 1_AAACCCAAGCATCTTG-1 ...
  2_TTTGTTGTCTAGGCCG-1 2_TTTGTTGTCTCGGCTT-1
colData names(32): Barcode Capture ... scDblFinder.cxds_score batch
reducedDimNames(0):
mainExpName: NULL
altExpNames(2): HTO ADT

$`/Users/maksimovicjovana/Work/Projects/MCRI/melanie.neeland/paed-inflammation-CITEseq/data/C133_Neeland_batch4/data/SCEs/C133_Neeland_batch4.doublets_called.SCE.rds`
class: SingleCellExperiment 
dim: 36601 45760 
metadata(0):
assays(1): counts
rownames(36601): ENSG00000243485 ENSG00000237613 ... ENSG00000278817
  ENSG00000277196
rowData names(20): ID Symbol ... is_mito is_pseudogene
colnames(45760): 1_AAACCCAAGCGTTAGG-1 1_AAACCCAAGGATTTGA-1 ...
  2_TTTGTTGTCGACGATT-1 2_TTTGTTGTCTAGGCCG-1
colData names(32): Barcode Capture ... scDblFinder.cxds_score batch
reducedDimNames(0):
mainExpName: NULL
altExpNames(2): HTO ADT

$`/Users/maksimovicjovana/Work/Projects/MCRI/melanie.neeland/paed-inflammation-CITEseq/data/C133_Neeland_batch5/data/SCEs/C133_Neeland_batch5.doublets_called.SCE.rds`
class: SingleCellExperiment 
dim: 36601 44715 
metadata(0):
assays(1): counts
rownames(36601): ENSG00000243485 ENSG00000237613 ... ENSG00000278817
  ENSG00000277196
rowData names(20): ID Symbol ... is_mito is_pseudogene
colnames(44715): 1_AAACCCAAGAAGATCT-1 1_AAACCCAAGATGCAGC-1 ...
  2_TTTGTTGTCGGATTAC-1 2_TTTGTTGTCTGAGAGG-1
colData names(32): Barcode Capture ... scDblFinder.cxds_score batch
reducedDimNames(0):
mainExpName: NULL
altExpNames(2): HTO ADT

$`/Users/maksimovicjovana/Work/Projects/MCRI/melanie.neeland/paed-inflammation-CITEseq/data/C133_Neeland_batch6/data/SCEs/C133_Neeland_batch6.doublets_called.SCE.rds`
class: SingleCellExperiment 
dim: 36601 46485 
metadata(0):
assays(1): counts
rownames(36601): ENSG00000243485 ENSG00000237613 ... ENSG00000278817
  ENSG00000277196
rowData names(20): ID Symbol ... is_mito is_pseudogene
colnames(46485): 1_AAACCCAAGAAGCGCT-1 1_AAACCCAAGACTCATC-1 ...
  2_TTTGTTGTCGAGAATA-1 2_TTTGTTGTCTACTGAG-1
colData names(32): Barcode Capture ... scDblFinder.cxds_score batch
reducedDimNames(0):
mainExpName: NULL
altExpNames(2): HTO ADT

Load sample information

ext_metadata_df <- read_excel(
  here("data/CZI_samples_design_with_micro.xlsx"),
  col_types = 
    c("text", "numeric", "text", "text", "text", "numeric", "numeric", "text",
      "numeric", "text", "numeric", "numeric", "numeric", "numeric",
      "text", "text", "text")) %>%
  mutate(HASHTAG = paste0("Human_HTO_", HASHTAG)) %>%
  mutate(sample.id = ave(Donor_ID, Donor_ID, 
                         FUN = function(i) paste0("sample_", i, '.', seq_along(i)))) %>%
  dplyr::rename(Batch = BATCH, 
                Donor = Donor_ID, 
                Hashtag = HASHTAG,
                arest.id = AREST_ID) %>%
  mutate(Donor = as.factor(Donor))
knitr::kable(ext_metadata_df)

Assign sample IDs to cells

For genetically distinct samples, cells are assigned to a sample using the genetic assignment; the sample information is then assigned using the corresponding HTO. For longitudinal samples only, both cells and sample information are assigned to the individual samples using HTOs.

sceLst <- sapply(1:length(sceLst), function(i){
  message(i)
  sce <- sceLst[[i]]
  
  # convert genetic and HTO labels to character
  sce$GeneticDonor <- as.character(sce$GeneticDonor)
  sce$dmmHTO <- ifelse(!stringr::str_detect(sce$dmmHTO, "Human_HTO_") & 
                         !(stringr::str_detect(sce$dmmHTO, "Doublet|Negative")),
                       paste0("Human_HTO_", sce$dmmHTO),
                       as.character(sce$dmmHTO))
  
  # match genetics and HTOs
  tab <- table(sce$GeneticDonor, sce$dmmHTO)
  tab <- tab[!rownames(tab) %in% c("Doublet","Unknown"),
             !colnames(tab) %in% c("Doublet","Negative")]
  # extract all genetically unique samples in this batch
  sub <- tab[rowSums(tab/rowSums(tab) > 0.2) == 1, , drop = FALSE]
  # create dictionary to translate genetics to HTO
  dict <- apply(sub, 1, function(x){
    colnames(tab)[which(x == max(x))]
  })
  dict <- c(dict, "Doublet" = "Doublet", "Unknown" = "Negative")
 
  sce$ProcessedHTO <- NA
  # assign HTO based on genetics
  genetics <- sce$GeneticDonor %in% names(dict)
  sce$ProcessedHTO[genetics] <- unname(dict[sce$GeneticDonor[genetics]])
  
  # if there are longitudinal samples also assign HTO based on hashtags
  if(dim(sub)[1] != dim(sub)[2]){
    HTO <- !sce$GeneticDonor %in% names(dict)
    sce$ProcessedHTO[HTO] <- sce$dmmHTO[HTO]
    
  }
  
  # add sample metadata based on processed HTO
  colData(sce) %>%
    data.frame() %>%
    rownames_to_column("cell") %>%
    left_join(ext_metadata_df %>% 
                dplyr::filter(Batch == (i-1)) %>%
                data.frame, 
              by = c("ProcessedHTO" = "Hashtag")) %>%
    mutate(sample.id = ifelse(is.na(sample.id), 
                              ProcessedHTO,
                              sample.id)) %>%
    column_to_rownames("cell") %>% 
    DataFrame() -> colData(sce)
  
  sce
})

p <- lapply(1:length(sceLst), function(i){
  sce <- sceLst[[i]]
  colData(sce) %>%
    data.frame() %>%
    ggplot(aes(x = sample.id, fill = Donor)) +
    geom_bar(position = position_dodge()) + 
    geom_text(aes(label = ..count..),
              stat = "count", vjust = 1.5, 
              colour = "black", size = 2) +
    theme(legend.position = "none") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
    ggtitle(glue("batch {i-1}"))
})

wrap_plots(p, ncol = 2)

p <- lapply(1:length(sceLst), function(i){
  sce <- sceLst[[i]]
  colData(sce) %>%
    data.frame() %>%
    mutate(scds = ifelse(hybrid_call, "Doublet", "Singlet"),
           scdf = ifelse(scDblFinder.class == "doublet", "Doublet", "Singlet")) %>%
    dplyr::select(ProcessedHTO, scds, scdf, Capture) %>% 
    rownames_to_column(var = "cell") %>%
    mutate(vireo_dmm_dbl = (ProcessedHTO == "Doublet"),
           scds_dbl = (scds == "Doublet"),
           scdf_dbl = (scdf == "Doublet")) %>%
    dplyr::select(cell, vireo_dmm_dbl, scds_dbl, scdf_dbl, Capture) %>% 
    pivot_longer(cols = c(vireo_dmm_dbl, scds_dbl, scdf_dbl), names_to = "method") %>%
    dplyr::filter(value == TRUE) %>%
    group_by(cell) %>%
    summarise(data = list(method), Capture) %>%
    rowwise() %>%
    ggplot(aes(x = data)) +
    geom_bar() +
    scale_x_upset(n_intersections = 20) +
    geom_text(stat = 'count', aes(label = after_stat(count)), 
              vjust = -0.5, size = 2) +
    ggtitle(glue("batch {i-1}"))
})

wrap_plots(p, ncol = 2)

Doublet filtering

We HAVE to filter out cells that are called doublets by either vireo or demuxmix (for relevant samples only) because we have no singlet assignment to a sample for them. We also HAVE to filter out negative cells as they also have no singlet assignment. We will also filter out cells that are called doublets by BOTH scds and scDblFinder (both methods were suggested in the demuxafy paper) that are not already in the vireo/demuxmix doublet pool. The percentages for those doublets calls per capture are shown below.

lapply(sceLst, function(sce){
  scds_scdf <- table(sce$Capture, sce$hybrid_call & sce$scDblFinder.class == "doublet" &
                       !sce$ProcessedHTO == "Doublet")
  vireo_dmm <- table(sce$Capture, sce$ProcessedHTO == "Doublet")
  
  if(length(levels(sce$Capture)) < 4) {
    data.frame(scds_scdf_dbl = scds_scdf[,2]/scds_scdf[,1]*100,
               vireo_dmm_dbl = vireo_dmm[,2]/vireo_dmm[,1]*100) -> tab
  } else {
    data.frame(scds_scdf_dbl = scds_scdf[,2]/scds_scdf[,1]*100) -> tab
  }
  
  tab %>% knitr::kable()
})

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

First, filter out wheeze and CSLD samples as we will focus only on CF.

sceLst <- lapply(sceLst, function(sce){
  sce[, !sce$Disease %in% c("Wheeze", "CSLD")]
})

p <- lapply(1:length(sceLst), function(i){
  sce <- sceLst[[i]]
  colData(sce) %>%
    data.frame() %>%
    ggplot(aes(x = sample.id, fill = Donor)) +
    geom_bar(position = position_dodge()) + 
    geom_text(aes(label = ..count..),
              stat = "count", vjust = 1.5, 
              colour = "black", size = 2) +
    theme(legend.position = "none") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
    ggtitle(glue("batch {i-1}"))
})

wrap_plots(p, ncol = 2)

Past versions of unnamed-chunk-9-1.png

Then, filter out doublets.

sceLst <- lapply(sceLst, function(sce){
  sce[, !(sce$ProcessedHTO %in% c("Doublet") | 
                                  (sce$hybrid_call & sce$scDblFinder.class == "doublet"))]
})

p <- lapply(1:length(sceLst), function(i){
  sce <- sceLst[[i]]
  colData(sce) %>%
    data.frame() %>%
    ggplot(aes(x = sample.id, fill = Donor)) +
    geom_bar(position = position_dodge()) + 
    geom_text(aes(label = ..count..),
              stat = "count", vjust = 1.5, 
              colour = "black", size = 2) +
    theme(legend.position = "none") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
    ggtitle(glue("batch {i-1}"))
})

wrap_plots(p, ncol = 2)

Past versions of unnamed-chunk-10-1.png

Finally, filter out negatives cells.

sceLst <- lapply(sceLst, function(sce){
  sce[, !(sce$sample.id %in% "Negative")]
})

p <- lapply(1:length(sceLst), function(i){
  sce <- sceLst[[i]]
  colData(sce) %>%
    data.frame() %>%
    ggplot(aes(x = sample.id, fill = Donor)) +
    geom_bar(position = position_dodge()) + 
    geom_text(aes(label = ..count..),
              stat = "count", vjust = 1.5, 
              colour = "black", size = 2) +
    theme(legend.position = "none") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
    ggtitle(glue("batch {i-1}"))
})

wrap_plots(p, ncol = 2)

Past versions of unnamed-chunk-11-1.png

Cell counts per sample after filtering

p <- lapply(sceLst, function(sce){
  colData(sce) %>%
    data.frame() %>%
    mutate(donor.fct = factor(Donor, labels = 1:length(unique(Donor)))) %>%
    ggplot(aes(x = sample.id, fill = donor.fct)) +
    geom_bar(position = position_dodge()) + 
    geom_text(aes(label = ..count..),
              stat = "count", vjust = 1.5, 
              colour = "black", size = 2) +
    theme(legend.position = "none") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))
})

wrap_plots(p, ncol = 2)

Past versions of unnamed-chunk-12-1.png

QC per sample after filtering

p <- lapply(sceLst, function(sce){
  colData(sce) %>%
    data.frame() %>% 
    mutate(donor.fct = factor(Donor, labels = 1:length(unique(Donor)))) %>%
    ggplot(aes(x = sample.id, y = sum,
               fill = donor.fct)) +
    geom_violin() +
    scale_y_log10() +
    theme(legend.position = "none") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
    labs(x = "Sample",
         y = "Library size (log10 scale)") -> p1
  
  colData(sce) %>%
    data.frame() %>% 
    mutate(donor.fct = factor(Donor, labels = 1:length(unique(Donor)))) %>%
    ggplot(aes(x = sample.id, y = detected,
               fill = donor.fct)) +
    geom_violin() +
    scale_y_log10() +
    theme(legend.position = "none") +
    theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) +
    labs(x = "Sample",
         y = "No. genes detected (log10 scale)") -> p2
  
  p1 + p2
})

wrap_plots(p, ncol = 1) &
  theme(axis.title.y = element_text(size = 9))

Past versions of unnamed-chunk-13-1.png

Save data

batches <- str_extract(files, "batch[0-6]")

sapply(1:length(sceLst), function(i){
  out <- here("data",
              paste0("C133_Neeland_", batches[i]),
              "data", 
              "SCEs", 
              glue("C133_Neeland_{batches[i]}.doublets_filtered.SCE.rds"))
  if(!file.exists(out)) saveRDS(sceLst[[i]], out)
  fs::file_chmod(out, "664")
  if(any(str_detect(fs::group_ids()$group_name, 
                    "oshlack_lab"))) fs::file_chown(out, 
                                                    group_id = "oshlack_lab")
})

[[1]]
NULL

[[2]]
NULL

[[3]]
NULL

[[4]]
NULL

[[5]]
NULL

[[6]]
NULL

[[7]]
NULL