get_czi_reference_data.R

# CZI Sketched Object
# ==================================
## some code from Yuhan ##

library(Seurat)
library(BPCells)
library(tidyverse)
library(rhdf5)
library(dplyr)
library(Azimuth)

# Started with downloading h5ads from czi w python api 
## in "cellxgene_python_download.ipynb"
### did 3 different ways, going with the tissue anndata version 

setwd("/brahms/mollag/meta_analysis/get_data/czi_cellcensus_h5ads_by_tissue/")
h5ad_path <- list.files(pattern = ".h5ad", recursive = TRUE) 

# Turn h5ads into an object

mat.list <- list()
check.list <- list()
failed.list <- list()

for (i in 1:length(h5ad_path)){
  file.i <- h5ad_path[i]
  name <- gsub(".h5ad", "", file.i)
  print(paste('file', i, ' ', file.i))
  worked <- tryCatch({
    mat.i <- open_matrix_anndata_hdf5(path = file.i , group = 'X'  )
    worked <<- TRUE
  }, 
  error=function(cond) {
    print(paste("errored out at", file.i))
    failed.list <<- c(failed.list, file.i)
    worked <<- FALSE
  }
  )
  print(worked)
  if (worked == FALSE) next
  if (  as.sparse(mat.i[, 1])@x[1]%%1 != 0 ) { 
    mat.i <- open_matrix_anndata_hdf5(path = file.i , group = '/raw/X'  )
    print("had to access raw matrix")
  } 
  # because these formats don't have feature or cell names
  var <- h5read(file = file.i, name = "var")
  features <- var$feature_name
  rownames(mat.i) <- features
  colnames(mat.i) <- paste0(colnames(mat.i), "_", name)
  mat.list[[name]] <- mat.i
  check.list[[i]] <- as.sparse(mat.i[, 1])@x[1]%%1==0 
}
all(unlist(check.list))
length(failed.list)


meta.list <- list()
for (i in 1:length(h5ad_path)){
  file.i <- h5ad_path[i]
  name <- gsub(".h5ad", "", file.i)
  print(paste('file', i, ' ',file.i))
  
  meta.i <-  LoadH5ADobs(file.i)
  print("opened metadata")
  
  meta.i$batch_id <- paste0('batch_',i)
  meta.i$batch_code <- name
  
  # make similar cell names as loop above to allow for metadata transfer
  rownames(meta.i) <- paste0(rownames(meta.i), "_", name)
  
  meta.list[[name]] <- meta.i
}


identical(names(mat.list), names(meta.list))


features.list <- lapply(mat.list, function(x) rownames(x))
features.union <- Reduce(union, features.list)

y <- mat.list
for (i in seq_along(along.with = y)) {
  print(i)
  missing_row <- setdiff(x = features.union, y =  rownames(y[[i]]))
  if (length(x = missing_row) > 0) {
    print("setting missing rows")
    zero_i <- SparseEmptyMatrix(nrow = length(x = missing_row), 
                                ncol = ncol(x = y[[i]]), colnames =  colnames( y[[i]]), 
                                rownames = missing_row)
    zero_i <- as(object = zero_i, Class = "IterableMatrix")
    if (y[[i]]@transpose) {
      y[[i]] <- transpose_storage_order(matrix =  y[[i]])
    }
    y[[i]] <- rbind(y[[i]], zero_i)[features.union, 
    ]
  }
}


mat_all <- do.call(cbind, y)

# meta combining 
common.names <- Reduce(intersect,lapply(meta.list, function(x) colnames(x)))
meta.list.filter <- lapply(meta.list, function(x) x[,common.names])

meta_df <- bind_rows(meta.list.filter)


meta_df$donor_id <- gsub('_','-', meta_df$donor_id )
meta_df$cell_type <- gsub('_','-', meta_df$cell_type )
meta_df$batch_donor <- paste0(meta_df$batch_id,'_', meta_df$donor_id)

ncol(mat_all) # should equal 
nrow(meta_df)


write_matrix_dir(mat = mat_all, dir = "/brahms/mollag/meta_analysis/BPCells_directories/all_human_matrix_v3")
write.csv(meta_df, "/brahms/mollag/meta_analysis/metadata_matching_all_human_matrix_v3.csv")


# Seurat object
# ===================================
mat <- open_matrix_dir("/brahms/mollag/meta_analysis/BPCells_directories/all_human_matrix_v3")
meta_df <- read.csv("/brahms/mollag/meta_analysis/metadata_matching_all_human_matrix_v3.csv", row.names = 1)
# bpcells isnt saving the colnames: 
colnames(mat) <- colnames(mat_all)
rownames(mat) <- rownames(mat_all)
czi <- CreateSeuratObject(mat, meta.data = meta_df)
#saveRDS(czi, "/brahms/mollag/meta_analysis/objects/czi_objects/czi_unique_human_cells_v1.rds")

# DEALING WITH NOT UNIQUE CELLS 
# https://chanzuckerberg.github.io/cellxgene-census/cellxgene_census_docsite_data_release_info.html#errata

duplicates <- read.csv('/brahms/mollag/meta_analysis/obs_duplicate_primary.csv')
duplicate_cells <- rownames(czi@meta.data[czi$soma_joinid %in% duplicates$soma_joinid, ])
czi.unique <- subset(czi, cells = duplicate_cells, invert = T) 
saveRDS(czi.unique, "/brahms/mollag/meta_analysis/objects/czi_objects/czi_unique_human_cells.rds")
# this has 33515318 cells

# QC 

czi_subset <- subset(czi.unique, subset = nCount_RNA > 500 & nCount_RNA < 30000 & nFeature_RNA > 100)  
# this reduced to 29963419 cells
czi.meta <- czi_subset@meta.data
write.csv(czi.meta, "/brahms/mollag/meta_analysis/metadata_files/full_object_metadata/czi_full_meta.csv")




# make filtered BPCells matrix 

# Sketching
#============================================
czi_subset <- NormalizeData(czi_subset)
# normalize smart seq data by length - not doing for now 

#saveRDS(czi_subset, "/brahms/mollag/meta_analysis/objects/czi_objects/czi_cells_subset_unique.rds")
czi_subset <- readRDS("/brahms/mollag/meta_analysis/objects/czi_objects/czi_cells_subset_unique.rds")
# REMOVE THESE DATASET IDS because they are repeated in the AzimuthData 
# b07e5164-baf6-43d2-bdba-5a249d0da879 - pancreas muraro # NOT IN CENSUS CURRENTLY, BUT IF IT ENDS UP ON THERE REMOVE IT
# d4e69e01-3ba2-4d6b-a15d-e7048f78f22e - heart Litvinukova
# a43aa46b-bd16-47fe-bc3e-19a052624e79 - liver macparland

czi_subset <- subset(czi_subset, subset = dataset_id %in% c("d4e69e01-3ba2-4d6b-a15d-e7048f78f22e", "a43aa46b-bd16-47fe-bc3e-19a052624e79"), invert = T)
# removed 491261 cells mostly from heart
#filter out gtex
czi_subset <- subset(czi_subset, dataset_id == "4ed927e9-c099-49af-b8ce-a2652d069333", invert = T)
saveRDS(czi_subset, "/brahms/mollag/meta_analysis/objects/czi_objects/czi_cells_subset_unique_filtered.rds")