simon.R

suppressPackageStartupMessages({
  library(Seurat)
  library(Matrix)
  library(ggplot2)
  library(patchwork)
  library(stringr)
  library(pbapply)
})

{
  base.path <- "/mnt/datadisk//avi/analyses/simon/OTD/outs/filtered_feature_bc_matrix/"
  #base.path <- "/mnt/datadisk//avi/analyses/simon/OTT/outs/filtered_feature_bc_matrix/"
  matrices <- Read10X(data.dir = base.path)

  keep.cells <- intersect(
    colnames(matrices$`Gene Expression`), 
    colnames(matrices$`Antibody Capture`)
  )
  length(keep.cells)
  
  object <- CreateSeuratObject(
    counts = matrices$`Gene Expression`,
    min.cells = 10, min.features = 200
  )
  object[["ADT"]] <- CreateAssayObject(
    counts = matrices$`Antibody Capture`[-c(157:161), Cells(object)]
  )
  object[["HTO"]] <- CreateAssayObject(
    counts = matrices$`Antibody Capture`[157:161, Cells(object)]
  )
  
  VlnPlot(object, features = c("nCount_RNA", "nFeature_RNA",
                               "nCount_ADT", "nFeature_ADT"), 
          log = TRUE, pt.size = 0.001, ncol = 4)
  object
}

#saveRDS(object, "/mnt/datadisk/avi/analyses/simon/rds/OTT.rds")
object <- readRDS("/mnt/datadisk/avi/analyses/simon/rds/OTD.rds")
object

{ #RNA analyses
  VlnPlot(object, features = c("nCount_RNA", "nFeature_RNA"))
  object <- subset(object, subset = nCount_RNA < 20000 & nFeature_RNA < 4000)
  
  VlnPlot(object, features = c("nCount_ADT", "nFeature_ADT"))
  object <- subset(object, subset = nFeature_ADT > 120)
  object
  
  DefaultAssay(object) <- "RNA"
  object <- SCTransform(object)

  object <- RunPCA(object, npcs = 25)
  print(object[["pca"]], dims = 1:5, nfeatures = 5)
  VizDimLoadings(object, dims = 1:2, reduction = "pca")
  DimHeatmap(object, dims=5, cells=500, balanced = T)
  DimPlot(object, reduction = "pca")
  
  object <- RunUMAP(object, dims = 1:25)
  object <- FindNeighbors(object, dims = 1:25)
  object <- FindClusters(object)
  DimPlot(object, label = T)
  
  FeaturePlot(object, features = c("CD19"), order = T, reduction = "umap")
}

{ #ADT analyses
  DefaultAssay(object) <- "ADT"
  object <- NormalizeData(object, normalization.method = "CLR", margin = 2)
  VariableFeatures(object) <- rownames(object[["ADT"]])
  object <- ScaleData(object, assay = "ADT", verbose = F, do.scale = F)
  
  object <- RunPCA(object, reduction.name = "apca", 
                   reduction.key = "apca_", 
                   npcs = 15)
  object <- FindNeighbors(object, reduction = "apca", 
                          dims = 1:15)
  object <- FindClusters(object, verbose = T)
  object <- RunUMAP(object, reduction.name = "aumap", 
                    reduction.key = "aumap_", 
                    dims = 1:15, reduction = "apca")
  DimPlot(object, label = T)
  FeaturePlot(object, features = c("rna_CD19"), order = T, reduction = "aumap")
}

{ #HTO
  DefaultAssay(object) <- "HTO"
  object <- NormalizeData(object,   normalization.method = "CLR", margin = 2, verbose = F)
  VariableFeatures(object) <- rownames(object[["HTO"]]@counts)
  object <- ScaleData(object, assay = "HTO")
  
  object <- HTODemux(object, assay = "HTO", positive.quantile = 0.99, verbose = F)
  Idents(object) <- "HTO_classification.global"
  VlnPlot(object, features = "nCount_HTO", pt.size = 0.1, log = TRUE)
  
  object <- subset(object, idents = "Negative", invert = TRUE)

  object <- RunPCA(object, reduction.name = "hto.pca", reduction.key = "HPC_", verbose = F)
  object <- RunUMAP(object, reduction = "hto.pca", dims = 1:3, 
                    reduction.name = "hto.umap", reduction.key = "HUMAP_", verbose = F)
  DimPlot(object, reduction = "hto.umap", label = T, group.by = "hash.ID")
  table(object$hash.ID)
  
  rownames(object@assays$ADT)
  FeaturePlot(object, features = paste0("adt_", c("CD16", "CD14.1", "CD11b",
                                                  "CD4.1", "CD8", "CD19.1")), 
              reduction = "aumap")
}

{
  objects <- c(
    "otd" = readRDS("/mnt/datadisk/avi/analyses/simon/rds/OTD.rds"),
    "ott" = readRDS("/mnt/datadisk/avi/analyses/simon/rds/OTT.rds")
  )
  p1 <- FeaturePlot(objects$otd, features = paste0("adt_", c("CD16", "CD14.1", "CD11b",
                                                  "CD4.1", "CD8", "CD19.1")), 
              reduction = "aumap")
  p2 <- FeaturePlot(objects$ott, features = paste0("adt_", c("CD16", "CD14.1", "CD11b",
                                                  "CD4.1", "CD8", "CD19.1")), 
              reduction = "aumap")
  p1 | p2
}