5_prepareSCVelo.R

library(SeuratWrappers)
library(Seurat)
library(SeuratDisk)

ts <- c("FB_10-A4_T", "FB_12-F1_T", "FB_16-B1_T", "FB_18-F2_T", "FB_4-B3_T", "FB_6-D1_T", "FB_11-C0_T", "FB_17-E3_T", "FB_5-C1_T")
veldata <- list()
velobj <- list()
for (i in ts) {
    veldata[[i]] <- ReadVelocity(file = paste0("Birocchi_SciTraMed2022/velocyto_loom/", i, ".loom"))
    velobj[[i]] <- as.Seurat(x = veldata[[i]])
    velobj[[i]] <- RenameCells(object = velobj[[i]], new.names = gsub(pattern = "merged", replacement = "", x = colnames(velobj[[i]])))
    velobj[[i]] <- RenameCells(object = velobj[[i]], new.names = gsub(pattern = "x$", perl = T, replacement = "", x = colnames(velobj[[i]])))
    velobj[[i]] <- RenameCells(object = velobj[[i]], new.names = gsub(pattern = ":", perl = T, replacement = "_", x = colnames(velobj[[i]])))
    velobj[[i]] <- RenameCells(object = velobj[[i]], new.names = gsub(pattern = "_", perl = T, replacement = "-", x = colnames(velobj[[i]])))
    velobj[[i]] <- RenameCells(object = velobj[[i]], new.names = gsub(pattern = "-T-", perl = T, replacement = "_", x = colnames(velobj[[i]])))
    velobj[[i]]@meta.data$orig.ident <- i
}

velobj_ctrl <- merge(merge(velobj[["FB_10-A4_T"]], velobj[["FB_16-B1_T"]]), velobj[["FB_4-B3_T"]])
velobj_high <- merge(merge(velobj[["FB_11-C0_T"]], velobj[["FB_5-C1_T"]]), velobj[["FB_6-D1_T"]])
velobj_low <- merge(merge(velobj[["FB_12-F1_T"]], velobj[["FB_18-F2_T"]]), velobj[["FB_17-E3_T"]])
velo_full <- merge(merge(velobj_ctrl, velobj_high), velobj_low)

wdir <- "Birocchi_SciTraMed2022/results"
macro_obj <- readRDS(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_final.rds", sep = "/"))
colnames(macro_obj)[1:5]

# write.csv(Cells(macro_obj), file = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_cellID_obs.csv", sep = "/"), row.names = FALSE)
# write.csv(Embeddings(macro_obj, reduction = "umap"), file = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_cell_embeddings.csv", sep = "/"))
# colnames(macro_obj@meta.data)
# write.csv(macro_obj@meta.data[,c("orig.ident", "RNA_Group", "RNA_snn_res.0.6", "RNA_snn_res.1", "RNA_snn_res.1.2")],
#           file = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", "MacrophagesFiltered_fastMNN_mdata.csv", sep = "/"))
# 
# #write.csv(Cells(obj), file = "cellID_obs.csv")
# macro_obj_embedding <- Embeddings(macro_obj, reduction = "umap")
# head(macro_obj_embedding)
# 
# #write.csv(Embeddings(obj, reduction = "umap"), file = "cell_embeddings.csv")
# #my_clusters <- obj@meta.data$RNA_snn_res.0.6
# macro_obj_clu <- FetchData(object = macro_obj, vars = "RNA_snn_res.0.6")
# colnames(macro_obj_clu) <- c("Seurat_RNA_snn_res.0.6")
# head(macro_obj_clu)
# #write.csv(obj@meta.data$seurat_clusters, file = "clusters.csv")
# macro_obj_groups <- FetchData(object = macro_obj, "RNA_Group")
# colnames(macro_obj_groups) <- c("RNA_Group")
# head(macro_obj_groups)
# 
# cells <- intersect(Cells(velo_full), Cells(macro_obj))
# length(Cells(velo_full))
# length(Cells(macro_obj))
# length(cells)
# velo_full <- subset(velo_full, cells = cells)
# 
# velo_full[["RNA"]] <- velo_full[["spliced"]]
# #########velo_full[["RNA"]] <- macro_obj@assays$RNA
# # adding clusters and groups
# velo_full <- AddMetaData(object = velo_full, metadata = macro_obj_clu)
# velo_full <- AddMetaData(object = velo_full, metadata = macro_obj_groups)
# 
# colnames(velo_full@meta.data)
# 
# names(velo_full@assays)
# str(velo_full)
# # creating fake umap slots in order to create the reduction slot associated
# velo_full <- SCTransform(velo_full)
# velo_full <- RunPCA(velo_full)
# velo_full <- RunUMAP(velo_full, dims = 1:30, n.components = 2)
# velo_full <- FindNeighbors(velo_full, dims = 1:30)
# velo_full <- FindClusters(velo_full)
# DefaultAssay(velo_full) <- "RNA"
# # restore original UMAPs
# head(velo_full@reductions$umap@cell.embeddings)
# head(macro_obj_embedding[rownames(velo_full@reductions$umap@cell.embeddings),])
# velo_full@reductions$umap@cell.embeddings <- macro_obj_embedding[rownames(velo_full@reductions$umap@cell.embeddings),]
# SaveH5Seurat(velo_full, filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2","MacrophagesFiltered_velo.h5Seurat", sep = "/"), overwrite = T)
# Convert(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2","MacrophagesFiltered_velo.h5Seurat", sep = "/"), dest = "h5ad", overwrite = T)

# Groups
Idents(macro_obj) <- "RNA_Group"
velobj_groups <- list("AB_Ctrl" = velobj_ctrl, "CD_High" = velobj_high, "EF_Low" = velobj_low)
for (group in names(velobj_groups)) {
    print(group)
    velo_sub <- velobj_groups[[group]]
    macro_obj_sub <- subset(macro_obj, idents = group)
    macro_obj_sub_embedding <- Embeddings(macro_obj_sub, reduction = "umap")
    head(macro_obj_sub_embedding)
    macro_obj_sub_clu <- FetchData(object = macro_obj_sub, vars = "RNA_snn_res.0.6")
    head(macro_obj_sub_clu)
    cells_sub <- intersect(Cells(velo_sub), Cells(macro_obj_sub))
    length(Cells(velo_sub))
    length(Cells(macro_obj_sub))
    length(cells_sub)
    velo_sub <- subset(velo_sub, cells = cells_sub)
    velo_sub[["RNA"]] <- velo_sub[["spliced"]]
    # adding clusters
    velo_sub <- AddMetaData(object = velo_sub, metadata = macro_obj_sub_clu)
    # creating fake umap slots in order to create the reduction slot associated
    velo_sub <- SCTransform(velo_sub)
    velo_sub <- RunPCA(velo_sub)
    velo_sub <- RunUMAP(velo_sub, dims = 1:30)
    velo_sub <- FindNeighbors(velo_sub, dims = 1:30)
    velo_sub <- FindClusters(velo_sub)
    DefaultAssay(velo_sub) <- "RNA"
    # restore original UMAPs
    head(velo_sub@reductions$umap@cell.embeddings)
    head(macro_obj_sub_embedding[rownames(velo_sub@reductions$umap@cell.embeddings),])
    velo_sub@reductions$umap@cell.embeddings <- macro_obj_sub_embedding[rownames(velo_sub@reductions$umap@cell.embeddings),]
    SaveH5Seurat(velo_sub, filename = paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", paste0("MacrophagesFiltered-", group, "_velo.h5Seurat"), sep = "/"), overwrite = T)
    Convert(paste(wdir, "MacrophagesFiltered", "01-fastMNN_2", paste0("MacrophagesFiltered-", group, "_velo.h5Seurat"), sep = "/"), dest = "h5ad", overwrite = T)
}