get_gtex_reference_data.R

# gtex

# custom sketching function
SketchData_Metadata <- function(
  object,
  assay = NULL,
  ncells = 5000L,
  sketched.assay = 'sketch',
  method = c('LeverageScore', 'Uniform'),
  var.name = "leverage.score",
  over.write = FALSE,
  seed = 123L,
  cast = 'dgCMatrix',
  verbose = TRUE,
  meta.col = NULL,
  min.cells = 50,
  target.category = NULL,
  ...
) {
  assay <- assay[1L] %||% DefaultAssay(object = object)
  assay <- match.arg(arg = assay, choices = Assays(object = object))
  method <- match.arg(arg = method)
  if (sketched.assay == assay) {
    abort(message = "Cannot overwrite existing assays")
  }
  if (sketched.assay %in% Assays(object = object)) {
    if (sketched.assay == DefaultAssay(object = object)) {
      DefaultAssay(object = object) <- assay
    }
    object[[sketched.assay]] <- NULL
  }
  if (!over.write) {
    var.name <- Seurat:::CheckMetaVarName(object = object, var.name = var.name) # REMOVE THE SEURAT::: WHEN DONE
  }
  
  if (method == 'LeverageScore') {
    if (verbose) {
      message("Calcuating Leverage Score")
    }
    object <- LeverageScore(
      object = object,
      assay = assay,
      var.name = var.name,
      over.write = over.write,
      seed = seed,
      verbose = FALSE,
      ...
    )
  } else if (method == 'Uniform') {
    if (verbose) {
      message("Uniformly sampling")
    }
    object[[var.name]] <- 1
  }
  
  leverage.score <- object[[var.name]]
  # Modified leverage score to add metadata 
  leverage.score.meta <- data.frame(leverage.score, Label = object@meta.data[[meta.col]])
  layers.data <- Layers(object = object[[assay]], search = 'data')
  cells <- lapply(
    X = seq_along(along.with = layers.data),
    FUN = function(i, seed) {
      set.seed(seed = seed)
      lcells <- Cells(x = object[[assay]], layer = layers.data[i])
      if (length(x = lcells) < ncells) {
        return(lcells)
      }
      return(sample(
        x = lcells,
        size = ncells,
        prob = leverage.score[lcells,]
      ))
    },
    seed = seed
  )
  sampled_cells <- unlist(cells)
  # MODIFYING TO SAMPLE AT LEAST MIN.CELL OF EACH CELL TYPE
  # Check each category for min.cells
  sampled_meta <- leverage.score.meta[sampled_cells, ]
  min_target_category_cells <- ceiling(0.25 * ncells)  # 25% of total cells
  for(label in unique(leverage.score.meta$Label)) {
    count <- sum(sampled_meta$Label == label, na.rm = TRUE)
    if(count < min.cells) {
      to_add <- min.cells - count
    } else {
      to_add <- 0
    }
    # Check if this is the target category and if it needs more cells to reach 25%
    if (!is.null(target.category) && label == target.category && count < min_target_category_cells) {
      to_add <- max(to_add, min_target_category_cells - count)
      print(paste0("adding ", to_add, " cells to the target categorty"))
    }
    
    if (to_add > 0) {
      print(paste(label, "needs additional cells:", to_add))
      supplement_cells <- subset(leverage.score.meta, Label == label)
      # Remove cells already in sampled_cells
      supplement_cells <- supplement_cells[!(rownames(supplement_cells) %in% sampled_cells), ]
      # Order by score and take the top 'to_add' cells
      supplement_cells <- supplement_cells[order(-supplement_cells[[var.name]]), ]
      supplement_cells <- head(supplement_cells, n = to_add)
      # Add the new cells to sampled_cells
      sampled_cells <- c(sampled_cells, rownames(supplement_cells))
    }
  }
  print(length(sampled_cells)) 
  print(length(unique(sampled_cells)))
  sketched <- suppressWarnings(expr = subset(
    x = object[[assay]],
    cells = sampled_cells, # changed this to sampled_cells
    layers = Layers(object = object[[assay]], search = c('counts', 'data'))
  ))
  for (lyr in layers.data) {
    try(
      expr = VariableFeatures(object = sketched, method = "sketch", layer = lyr) <-
        VariableFeatures(object = object[[assay]], layer = lyr),
      silent = TRUE
    )
  }
  if (!is.null(x = cast) && inherits(x = sketched, what = 'Assay5')) {
    sketched <- CastAssay(object = sketched, to = cast, ...)
  }
  Key(object = sketched) <- Key(object = sketched.assay, quiet = TRUE)
  object[[sketched.assay]] <- sketched
  DefaultAssay(object = object) <- sketched.assay
  return(object)
}

# Gtex data from https://gtexportal.org/home/singleCellOverviewPage
#gtex.obj <- sceasy::convertFormat("/brahms/mollag/meta_analysis/get_data/gtex/gtex.h5ad", from = "anndata", to = "seurat")
gtex <- open_matrix_anndata_hdf5("/brahms/mollag/meta_analysis/get_data/gtex/GTEx_8_tissues_snRNAseq_atlas_071421.public_obs.h5ad")
gtex_obs <- LoadH5ADobs("/brahms/mollag/meta_analysis/get_data/gtex/GTEx_8_tissues_snRNAseq_atlas_071421.public_obs.h5ad")
gtex.obj <- CreateSeuratObject(gtex, meta.data = gtex_obs)

#gtex.obj[["RNA"]] <- as(gtex.obj[["RNA"]], "Assay5")
#gtex.obj[["RNA"]]$counts <- write_matrix_dir(gtex.obj[["RNA"]]$counts, "/brahms/mollag/meta_analysis/BPCells_directories/gtex_test_counts_2")
#gtex.obj[["RNA"]]$data <- write_matrix_dir(gtex.obj[["RNA"]]$data, "/brahms/mollag/meta_analysis/BPCells_directories/gtex_test_data_2")

gtex.obj <- NormalizeData(gtex.obj)
gtex.obj <- FindVariableFeatures(gtex.obj)
gtex.obj <- ScaleData(gtex.obj)
gtex.obj <- RunPCA(gtex.obj)
gtex.obj <- RunUMAP(gtex.obj, dims = 1:50)

gtex.obj$granular_cell_type <- gtex.obj$`Granular cell type`
gtex.obj$cell_type_broad <- gtex.obj$`Cell types level 3`
DimPlot(gtex.obj, group.by = "granular_cell_type", label = T, repel = T, alpha = 0.1, raster = F) + NoLegend()
DimPlot(gtex.obj, group.by = "cell_type_broad", label = T, repel = T, alpha = 0.1, raster = F) + NoLegend()
gtex.obj$cell_type <- gtex.obj$granular_cell_type
gtex.obj$sex <- gtex.obj$Sex
gtex.obj$tissue_general <- gtex.obj$tissue
gtex.obj$batch_id <- paste0("gtex_", gtex.obj$batch)
gtex.obj$batch_donor <- paste0(gtex.obj$batch, "_", gtex.obj$individual)
#gtex.obj$donor_id <- gtex.obj$individual
gtex$donor_id <- gtex$`Participant ID`
gtex$dataset_id <- gtex$`Sample ID`

saveRDS(gtex.obj, "/brahms/mollag/meta_analysis/objects/gtex_objs/gtex_obj.rds")
summary(gtex.obj$nCount_RNA)

gtex.meta <- gtex.obj@meta.data
write.csv(gtex.meta, "/brahms/mollag/meta_analysis/metadata_files/full_object_metadata/gtex_full_meta.csv")

# Sketching # WE DONT DO THIS IN THE PREPROCESSING ANYMORE 
#============================================
gtex.obj <- NormalizeData(gtex.obj)
gtex.obj <- FindVariableFeatures(gtex.obj)
gtex.obj <- SketchData(object = gtex.obj, ncells = 100000, method = "LeverageScore", sketched.assay = "sketch")
saveRDS(gtex.obj, "/brahms/mollag/meta_analysis/objects/gtex_objs/gtex_obj_sketch.rds")



# Custom Sketching
gtex.obj <- SketchData_Metadata(object = gtex.obj, ncells = 50000, meta.col = "cell_type", min.cells = 50, 
                                method = "LeverageScore", sketched.assay = "sketch")
saveRDS(gtex.obj, "/brahms/mollag/meta_analysis/objects/gtex_objs/gtex_obj_custom_sketch.rds")