cds_preprocess.Rmd

---
title: "CellCode 10x scRNAseq analysis"
author: "Matthew Wither"
date: "9/1/23"
output: html_document
editor_options: 
  chunk_output_type: console
---
  
Run startup
```{r}

setwd("~/CellCode")

#Run startup script
source("startup.R")

#Load data - generated by Sriram
cds <- readRDS("Robjects/cds_calledconditions_cellcycle_071123.rds")

```

Process cds
```{r}
#Compute number of genes expressed per cell
cds <- detect_genes(cds, min_expr = 0.1)

# add a UMI column into pData
pData(cds)$UMI <- Matrix::colSums(exprs(cds))

#Get expression of mitochondrial genes per cell
mito_genes <- startsWith(rowData(cds)$gene_short_name, 'mt-')
cds_mito <- cds[mito_genes,]

#Add mitochondrial UMIs to pData
pData(cds)$UMI_mito <- Matrix::colSums(exprs(cds_mito))

# Add column for percent mitochondrial reads
pData(cds)$frac_mito <- pData(cds)$UMI_mito / pData(cds)$UMI

pData1 <- as.data.frame(pData(cds))

#Plot UMI vs # genes expressed per cell
ggplot(pData1, aes(x = num_genes_expressed, y = UMI)) + 
  geom_point()

ggplot(pData1, aes(x = UMI, y = UMI_mito)) + 
  #stat_density_2d(aes(fill = after_stat(level)), geom = "polygon") +
  #geom_hex(bins = 50) + scale_fill_continuous(type = "viridis") +
  geom_point(aes(color = frac_mito)) + scale_color_continuous(type = "viridis") +
  theme_bw()

ggplot(subset(pData1, UMI > 500 & UMI_mito < 500), aes(x = UMI, y = UMI_mito)) + 
  #stat_density_2d(aes(fill = after_stat(level)), geom = "polygon") +
  #geom_hex(bins = 50) + scale_fill_continuous(type = "viridis") +
  geom_point(aes(color = frac_mito)) + scale_color_continuous(type = "viridis") +
  theme_bw()

#Visualize fraction of mitochondrial reads by condition
ggplot(pData1, aes(x = condition_norep, y = frac_mito, color = condition_norep)) + 
  geom_boxplot(show.legend = F) + 
  geom_hline(aes(yintercept = 0.075)) +
  labs(title = "Fraction Mitochondrial reads")


#Filter cells based on mitochondrial reads
valid_cells <- row.names(subset(pData1, 
                                frac_mito < 0.075 & #high mitochondrial reads
                                !is.na(condition_norep)))

cds <- cds[,valid_cells]

pData1 <- as.data.frame(pData(cds))


#Filter genes - remove genes expressed in < 20% of cells from all base conditions
genes_by_cond <- pbmclapply(split(pData1, pData1$condition_norep), function(x) {
  
  cds0 <- cds[,row.names(x)]
  cds0 <- detect_genes(cds0, min_expr = 0.1)
  
  #Filter out genes expressed in too few cells
  expressed_genes <- row.names(subset(fData(cds0), num_cells_expressed >= 0.2 * dim(cds0)[2]))
  
  return(expressed_genes)

})

expressed_genes <- unique(unname(unlist(genes_by_cond)))

cds <- cds[expressed_genes,]
cds <- detect_genes(cds, min_expr = 0.1)

pData1 <- as.data.frame(pData(cds))

#Filter cells round 2 - determined empirically
ggplot(pData1, aes(x = num_genes_expressed, y = frac_mito)) +
  #geom_point() +
  geom_hex(bins = 50) + scale_fill_continuous(type = "viridis") +
  geom_abline(aes(intercept = -0.01, slope = 0.000035), color = "red")

#Set threshold frac_mito based on num_genes_expressed for filtering out low quality cells
pData1$quality_score <- pData1$num_genes_expressed*0.000035 - 0.01
pData1$quality <- "High"
pData1$quality[which(pData1$frac_mito > pData1$quality_score)] <- "Low"


ggplot(pData1, aes(x = num_genes_expressed, y = frac_mito)) +
  geom_point(aes(color = quality)) +
  #geom_hex(bins = 50) + scale_fill_continuous(type = "viridis") +
  geom_abline(aes(intercept = -0.01, slope = 0.000035), color = "red")

#Now subset only the high quality cells
colData(cds)$quality <- pData1$quality
high_qual_cells <- row.names(subset(pData1, quality == "High"))

cds <- cds[,high_qual_cells]
pData1 <- as.data.frame(colData(cds))

cds <- preprocess_cds(cds, method = "PCA", num_dim = 100)
cds <- reduce_dimension(cds, preprocess_method = "PCA")

colData(cds)$Cond <- paste(pData1$cytokine, pData1$concentration, sep= "_")
pData1 <- as.data.frame(colData(cds))

#re-name conditions for ease
list_of_base_conditions <- unique(pData1$base)

pData1$base_new <- NA
pData1$base_new[which(pData1$base == list_of_base_conditions[[1]])] <- "CD3"
pData1$base_new[which(pData1$base == list_of_base_conditions[[2]])] <- "Rest"
pData1$base_new[which(pData1$base == list_of_base_conditions[[3]])] <- "CD3/28"
pData1$base_new[which(pData1$base == list_of_base_conditions[[4]])] <- "CD3/28/IL-12"

pData1$unique <- paste(pData1$base_new, pData1$Cond, sep= "_")

colData(cds)$base_new <- pData1$base_new
colData(cds)$unique <- pData1$unique

#Save umap coordinates for plotting with ggplot
umap_coords <- as.data.frame(reducedDims(cds)$UMAP)
colData(cds)$umap_1 <- umap_coords$V1
colData(cds)$umap_2 <- umap_coords$V2
rm(umap_coords)

#Create a unique cell identifier since there might be issues with the barcodes from multiple replicates and can cause issues with downstream functions
rownames(colData(cds)) <- paste("cell", 1:length(pData1$barcode), sep = "_")
colData(cds)$cell_id <- paste("cell", 1:length(pData1$barcode), sep = "_")


#Cluster cells
#Increasing k reduces the number of clusters
set.seed(42)
cds <- cluster_cells(cds, reduction_method = "UMAP", k = 100, cluster_method = "louvain")

#Add clusters to colData
colData(cds)$cluster <- clusters(cds)

plot_cells(cds, group_label_size = 8)

#Save processed cds for downstream analysis
saveRDS(cds, "Robjects/cds_processed.rds")

```

Create cds for each base condition
```{r}

cds_processed <- readRDS("Robjects/cds_processed.rds")
pData <- as.data.frame(pData(cds_processed))

list_of_bases <- unique(pData$base_new)

cds_by_base <- pbmclapply(list_of_bases, FUN = function(x) {
  if (endsWith(x, "IL-12")) {
    #Add CD3/28_IL-12 1ng/ml as null
    good_cells <- row.names(subset(pData,base_new == x))
    null_cells <- row.names(subset(pData,unique == "CD3/28_IL-12_1ng/mL"))
    good_cells <- c(good_cells, null_cells)
    
  } else if (x == "CD3") {
    good_cells <- row.names(subset(pData,base_new == x))
    CD28_null_cells <- row.names(subset(pData,unique == "CD3/28_Null_"))
    good_cells <- c(good_cells, CD28_null_cells)
    
  } else {
    good_cells <- row.names(subset(pData,base_new == x))
    
  }
 
  cds0 <- cds_processed[,good_cells]
  
  cds0 <- preprocess_cds(cds0, method = "PCA", num_dim = 100)
  
  #Regress out cell cycle phase
  #cds0 <- align_cds(cds0, num_dim = 100, alignment_group = "Phase")
  
  cds0 <- reduce_dimension(cds0, preprocess_method = "PCA")
  cds0 <- cluster_cells(cds0, reduction_method = "UMAP", k = 100, cluster_method = "louvain")
  colData(cds0)$cluster <- clusters(cds0)
  
  umap_coords <- as.data.frame(reducedDims(cds0)$UMAP)
  colData(cds0)$umap_1 <- umap_coords$V1
  colData(cds0)$umap_2 <- umap_coords$V2
  
  if (endsWith(x, "IL-12")) {
     pData0 <- as.data.frame(pData(cds0))
     pData0$unique[which(pData0$base_new == "CD3/28")] <- paste0(x,"_Null_")
     pData0$Cond[which(pData0$base_new == "CD3/28")] <- "Null_"
     pData0$base_new[which(pData0$base_new == "CD3/28")] <- x
     
     colData(cds0)$unique <- pData0$unique
     colData(cds0)$Cond <- pData0$Cond
     colData(cds0)$base_new <- pData0$base_new
  }
  
  if (x == "CD3") {
     pData0 <- as.data.frame(pData(cds0))
     pData0$Cond[which(pData0$base_new == "CD3/28")] <- "CD28"
     pData0$base_new[which(pData0$base_new == "CD3/28")] <- x

     colData(cds0)$Cond <- pData0$Cond
     colData(cds0)$base_new <- pData0$base_new
  }
  
  return(cds0)
}, mc.cores = 4)

names(cds_by_base) <- list_of_bases

saveRDS(cds_by_base, "Robjects/cds_by_base.rds")

```

Create cds of only null conditions to compare base conditions only
```{r}

cds_processed <- readRDS("Robjects/cds_processed.rds")
pData1 <- as.data.frame(pData(cds_processed))

null_cells <- row.names(subset(pData1,unique %in% c("CD3/28_IL-12_1ng/mL", "CD3_Null_", "CD3/28_Null_", "Rest_Null_")))
cds_null <- cds_processed[,null_cells]

```

Create cds of Hi doses only
```{r}

cds_hi_dose <- lapply(cds_by_base, function(cds) {
  pData <- as.data.frame(colData(cds))
  split_by_cyt <- split(pData, pData$cytokine)
  pData2 <- bind_rows(pblapply(split_by_cyt, function(x) {
    x$dose <- as.numeric(gsub("([0-9]+).*$", "\\1", x$concentration))
    if(any(!is.na(x$dose))) {
      max_dose <- max(x$dose)
      x <- subset(x, dose == max_dose)
    }
    return(x)
  }))
  
  cds_out <- cds[,row.names(pData2)]
  return(cds_out)
  
})


```