generate_cellpool.Rmd

---
title: "Generate cell pool"
output: html_document
---

```{r,message=FALSE,warning=FALSE}
library(splatter)
library(Seurat)
library(tidyverse)
```

## Splatter

Following codes generate the cell pools we used.

```{r, eval=FALSE}
mnnCT7_seurat <- readRDS("/storage/holab/Nelson/mnnCT7.rds")
count_mnnCT7 = mnnCT7_seurat@assays$RNA@counts
params <- splatEstimate(as.matrix(count_mnnCT7))
params
saveRDS(params, "/storage/holab/linxy/DCATS/mnnCT7params.rds")
```

```{r,eval=FALSE}
set.seed(123)
## simulate cells pool setting
cluster_num = 8
#prob = rep(1/cluster_num, cluster_num)
prob = c(rep(0.8/cluster_num, 7),0.2)
batch_size = cluster_num * 2000  ## generate 3000 cells for each cluster
de_prob = rep(NA, cluster_num)
#de_prob[1:round(cluster_num/3)] = runif(round(cluster_num/3), 0.1, 0.2)
#de_prob[(round(cluster_num/3)+1):cluster_num] = runif((cluster_num - round(cluster_num/3)), 0.5, 1)
de_prob[1:round(cluster_num/2)] = runif(round(cluster_num/2), 0.1, 0.2)
#de_prob[1:round(cluster_num/2)] = sample(c(0.13, 0.15), round(cluster_num/2), replace = TRUE)
de_prob[(round(cluster_num/2)+1):cluster_num] = runif((cluster_num - round(cluster_num/2)), 0.5, 1)
setresolu = 0.5
```

```{r}
par = "mnnCT7"  # can be select from "origin", "mnnCT7"
if (par == "origin"){
  params <- newSplatParams()
} else if (par == "mnnCT7"){
  #params = readRDS("./data/ovarian_params.rds")
  params = readRDS("/storage/holab/linxy/DCATS/mnnCT7params.rds")
}
```

```{r,eval=FALSE}
## simulate cells pool process
set.seed(123)
param.groups <- setParams(params, batchCells = batch_size*prob, nGenes = 1000)
#param.groups <- setParams(params, batchCells = batch_size, nGenes = 1000)
sim <- splatSimulateGroups(param.groups, group.prob = prob, de.prob = de_prob, de.facLoc = 0.01, verbose = FALSE)
sim_mat <- counts(sim)
origLabels = sim@colData@listData$Group
```

```{r}
# Save multiple objects
save(sim_mat, origLabels, file = str_c("/storage/holab/linxy/DCATS/cells_pool", as.character(cluster_num), "_splatter_large.RData"))
#save(sim_mat, origLabels, file = str_c("D:/Data/DCATS/cells_pool", as.character(cluster_num), "_splatter_test.RData"))
```

```{r}
load(str_c("/storage/holab/linxy/DCATS/cells_pool", as.character(cluster_num), "_splatter.RData"))
#load(str_c("./data/cells_pool", as.character(cluster_num), ".RData"))
```

```{r}
seuratObj <- CreateSeuratObject(counts = counts(sim), project="Splatter")
seuratObj <- NormalizeData(seuratObj, verbose = FALSE)
seuratObj <- FindVariableFeatures(seuratObj, selection.method = "vst", nfeatures = 2000, verbose = FALSE)
seuratObj <- ScaleData(seuratObj)
seuratObj <- RunPCA(seuratObj, features = VariableFeatures(object = seuratObj))
seuratObj <- FindNeighbors(seuratObj, dims = 1:30, verbose = FALSE)
seuratObj <- FindClusters(seuratObj, resolution = setresolu, verbose = FALSE)
seuratObj <- RunUMAP(seuratObj, dims = 1:30, verbose = FALSE)
seuratObj[['orig.ident']] <- origLabels
seuratObj[['orig.ident']] <- sim@colData@listData$Group
DimPlot(seuratObj, reduction = "umap")
DimPlot(seuratObj, reduction = "umap", group.by = "orig.ident")
```

## Splatter-discrete

```{r,eval=FALSE}
set.seed(123)
## simulate cells pool setting
cluster_num = 8
prob = rep(1/cluster_num, cluster_num)
batch_size = cluster_num * 1500
de_prob = rep(NA, cluster_num)
#de_prob[1:round(cluster_num/3)] = runif(round(cluster_num/3), 0.1, 0.2)
#de_prob[(round(cluster_num/3)+1):cluster_num] = runif((cluster_num - round(cluster_num/3)), 0.5, 1)
#de_prob[1:round(cluster_num/2)] = runif(round(cluster_num/2), 0.1, 0.2)
de_prob[1:round(cluster_num/2)] = runif(round(cluster_num/2), 0.2, 0.4)
de_prob[(round(cluster_num/2)+1):cluster_num] = runif((cluster_num - round(cluster_num/2)), 0.5, 1)
setresolu = 0.5
```

```{r}
par = "mnnCT7"  # can be select from "origin", "mnnCT7"
if (par == "origin"){
  params <- newSplatParams()
} else if (par == "mnnCT7"){
  #params = readRDS("./data/ovarian_params.rds")
  params = readRDS("D:/Data/DCATS/simulation/mnnCT7params.rds")
}
```

```{r,eval=FALSE}
## simulate cells pool process
set.seed(123)
#param.groups <- setParams(params, batchCells = batch_size*prob, nGenes = 1000)
param.groups <- setParams(params, batchCells = batch_size, nGenes = 1000)
sim <- splatSimulateGroups(param.groups, group.prob = prob, de.prob = de_prob, de.facLoc = 0.01, verbose = FALSE)
sim_mat <- counts(sim)
origLabels = sim@colData@listData$Group
```

```{r}
# Save multiple objects
#save(sim_mat, origLabels, file = str_c("D:/Data/DCATS/cells_pool", as.character(cluster_num), "_splatter_discrete.RData"))
#save(sim_mat, origLabels, file = str_c("D:/Data/DCATS/cells_pool", as.character(cluster_num), "_splatter_test.RData"))
```

```{r}
load(str_c("D:/Data/DCATS/cells_pool", as.character(cluster_num), "_splatter_discrete.RData"))
#load(str_c("./data/cells_pool", as.character(cluster_num), ".RData"))
```

```{r}
seuratObj <- CreateSeuratObject(counts = sim_mat, project="Splatter")
seuratObj <- NormalizeData(seuratObj, verbose = FALSE)
seuratObj <- FindVariableFeatures(seuratObj, selection.method = "vst", nfeatures = 2000, verbose = FALSE)
seuratObj <- ScaleData(seuratObj)
seuratObj <- RunPCA(seuratObj, features = VariableFeatures(object = seuratObj))
seuratObj <- FindNeighbors(seuratObj, dims = 1:30, verbose = FALSE)
seuratObj <- FindClusters(seuratObj, resolution = setresolu, verbose = FALSE)
seuratObj <- RunUMAP(seuratObj, dims = 1:30, verbose = FALSE)
seuratObj[['orig.ident']] <- origLabels
DimPlot(seuratObj, reduction = "umap")
DimPlot(seuratObj, reduction = "umap", group.by = "orig.ident")
table(seuratObj$orig.ident)
```

## real world data

```{r}
cluster_num = 7
```

```{r}
data_dir1 = "./data/real_world/Exp6_Kang/GSM2560248"
list.files(data_dir1) # Should show barcodes.tsv, genes.tsv, and matrix.mtx
expression_matrix1 <- Read10X(data.dir = data_dir1)
seuratObj1 = CreateSeuratObject(counts = expression_matrix1)
seuratObj1 = AddMetaData(object = seuratObj1, metadata = rep("ctrl1", dim(seuratObj1)[2]), col.name = 'sample')

data_dir2 = "./data/real_world/Exp6_Kang/GSM2560249"
list.files(data_dir2) # Should show barcodes.tsv, genes.tsv, and matrix.mtx
expression_matrix2 <- Read10X(data.dir = data_dir2)
seuratObj2 = CreateSeuratObject(counts = expression_matrix2)
seuratObj2 = AddMetaData(object = seuratObj2, metadata = rep("ctrl2", dim(seuratObj2)[2]), col.name = 'sample')
```

```{r}
cluster_info = read.delim("./data/real_world/Exp6_Kang/GSE96583_batch2.total.tsne.df.tsv") %>% 
  rownames_to_column(var = "barcode")
cluster_info %>% group_by(stim, cell) %>% 
  summarise(n = n())
summary(as.factor(cluster_info$cell))
```

```{r}
ggplot(cluster_info, aes(x = tsne1, y = tsne2, color = cell)) +
  geom_point() +
  facet_grid(. ~ stim)
```

```{r}
ctrl1 = cluster_info %>% 
  filter(cell %in% c("CD4 T cells", "CD14+ Monocytes", "FCGR3A+ Monocytes", "B cells", "CD8 T cells", "NK cells", "Dendritic cells", "Megakaryocytes")[1:cluster_num]) %>% 
  filter(stim == "ctrl")
ctrl2 = cluster_info %>% 
  filter(cell %in% c("CD4 T cells", "CD14+ Monocytes", "FCGR3A+ Monocytes", "B cells", "CD8 T cells", "NK cells", "Dendritic cells", "Megakaryocytes")[1:cluster_num]) %>% 
  filter(stim == "stim")

mat1 = subset(seuratObj1, cells = ctrl1$barcode)@assays$RNA@counts %>% 
  as.matrix()
mat2 = subset(seuratObj2, cells = ctrl2$barcode)@assays$RNA@counts %>% 
  as.matrix()
```

```{r}
sim_mat = cbind(mat1, mat2)
barcodesL = colnames(sim_mat)
selected_cells = cluster_info %>% 
  filter(barcode %in% barcodesL)
selected_cells[match(barcodesL, selected_cells$barcode),]
origLabels = selected_cells$cell
save(sim_mat, origLabels, file = str_c("./data/cells_pool", as.character(cluster_num), "_realWorld.RData"))
```

## GCTB

```{r}
GCTB5.integrated = readRDS("./data/real_world/GCTB5.integrated.1217.RDS")
```

```{r}
DimPlot(GCTB5.integrated, group.by = "seurat_clusters",
        reduction = "umap", split.by ="orig.ident", label = T, ncol = 2 )
DimPlot(GCTB5.integrated, group.by = "seurat_clusters",
        reduction = "umap", label = T, ncol = 2 )
DimPlot(GCTB5.integrated, group.by = "celltype",
        reduction = "umap",  label = T)
```

## SPARSim with batch

```{r}
library(SPARSim)
library(tidyverse)
library(Seurat)
```

```{r}
cluster_num = 8
batch_size = cluster_num * 1500
setresolu = 0.2
batch_num = 4
cell.prop = c(1/3, 1/3, 1/3)
de.prop = c(runif(round(cluster_num/3), 0.01, 0.2), runif((cluster_num - round(cluster_num/3)), 0., 1))
options(future.globals.maxSize = 15000 * 1024^2)
```

```{r}
set.seed(123)
# Load preset
data(Chu_param_preset)
cluster_name = c("A", "B", "C", "D", "E", "F", "G")
sim_resultL = list()
for (i in 1:cluster_num){
  ## cluster Info
  cluster_params <- SPARSim_create_simulation_parameter(
    intensity = Chu_param_preset[[i]]$intensity, 
    variability = Chu_param_preset[[i]]$variability, 
    library_size = sample(Chu_param_preset[[i]]$lib_size, size = round(batch_size*cell.prop[i])*batch_num, replace = TRUE),
    condition_name = str_c("cluster", cluster_name[i]))
  
  # Batch effect
  batch_set = list(batch1 = SPARSim_create_batch(name = "Lane1", distribution = "normal", param_A = 0, param_B = 1))
  batch_sample_association = rep("Lane1", round(batch_size*cell.prop[i]))
  for (j in 2:batch_num){
    batch_set[[str_c("batch", as.character(j))]] = SPARSim_create_batch(name = str_c("Lane", as.character(j)), distribution = "normal", param_A = j-2, param_B = 1)
    batch_sample_association = c(batch_sample_association, rep(str_c("Lane", as.character(j)), round(batch_size*cell.prop[1])))
}
  batch_params <- SPARSim_create_batch_parameter(batch_set = batch_set,  batch_sample = batch_sample_association)
  
  # Run SPARSim simulation using the just created simulation parameter, the 'parameter' object should be a list even there is only one condi
  sim_result <- SPARSim_simulation(dataset_parameter = list(cluster_params), batch_parameter = batch_params)
  sim_resultL[[str_c("cluster", cluster_name[i])]] = sim_result
}

```

```{r}
sim_result = matrix(NA, nrow = dim(sim_resultL[[1]]$count_matrix)[1], ncol = 0)
for (i in 1:length(sim_resultL)){
  sim_result = cbind(sim_result, sim_resultL[[i]]$count_matrix)
}
```

```{r}
seuratObj <- CreateSeuratObject(counts = sim_mat, project="SPARSim")
seuratObj <- NormalizeData(seuratObj, verbose = FALSE)
seuratObj <- FindVariableFeatures(seuratObj, selection.method = "vst", nfeatures = 2000, verbose = FALSE)
seuratObj <- ScaleData(seuratObj)
seuratObj <- RunPCA(seuratObj, features = VariableFeatures(object = seuratObj))
seuratObj <- FindNeighbors(seuratObj, dims = 1:30, verbose = FALSE)
seuratObj <- FindClusters(seuratObj, resolution = setresolu, verbose = FALSE)
seuratObj <- RunUMAP(seuratObj, dims = 1:30, verbose = FALSE)
seuratObj[['orig.ident']] = origLabels
DimPlot(seuratObj, reduction = "umap", group.by = "orig.ident")
```

```{r}
sim_mat = sim_result
origLabels = colnames(sim_result) %>% str_extract("cluster[ABC]") %>% tail()
batchInfo = rep(batch_sample_association, 3) %>% str_replace("Lane", "Batch")
save(sim_mat, origLabels, batchInfo, file = str_c("./data/cells_pool", as.character(cluster_num), "_SPARSim.RData"))
```

## SPARSim without batch

```{r}
library(SPARSim)
library(tidyverse)
library(Seurat)
```

```{r}
cluster_num = 4
batch_size = cluster_num * 1500
setresolu = 0.2
cell.prop = rep(1/cluster_num, cluster_num)
de.prop = c(runif(round(cluster_num/3), 0.1, 0.2), runif((cluster_num - round(cluster_num/3)), 0.2, 0.4))
options(future.globals.maxSize = 15000 * 1024^2)
```

```{r}
source("functionsV2.r")
```

```{r}
set.seed(123)
# Load preset
data(Chu_param_preset)
cluster_name = c("A", "B", "C", "D", "E", "F", "G")
cluster_params = list()
cluster_params[["clusterA"]] = SPARSim_create_simulation_parameter(
    intensity = Chu_param_preset$Chu_C1$intensity, 
    variability = Chu_param_preset$Chu_C1$variability, 
    library_size = sample(Chu_param_preset$Chu_C1$lib_size, size = round(batch_size*cell.prop[1]), replace = TRUE),
    condition_name = "clusterA")

for (clu in 2:cluster_num) {
  DE_multiplier = c(runif(n = 1000*de.prop[clu]/2, min = 0.05, max = 1), runif(n = 1000*de.prop[clu]/2, min = 1, max = 20))
  genes_num = length(cluster_params$clusterA$intensity)
  add_foldchange = rep(1, genes_num)
  add_foldchange[sample(1:genes_num, 1000*de.prop[clu])] = DE_multiplier
  cluster_params[[str_c("cluster", cluster_name[clu])]] = SPARSim_create_DE_genes_parameter(
                                  sim_param = cluster_params[[1]], 
                                  fc_multiplier = add_foldchange, 
                                  N_cells = 1000,
                                  condition_name = str_c("cluster", cluster_name[clu]))
}

SPARSim_result <- SPARSim_simulation(cluster_params)
sendEmail("Cell pool generation is done!!")
```

```{r}
#seuratObj <- CreateSeuratObject(counts = sim_mat, project="SPARSim")
seuratObj <- CreateSeuratObject(counts = SPARSim_result$count_matrix, project="SPARSim")
seuratObj <- NormalizeData(seuratObj, verbose = FALSE)
seuratObj <- FindVariableFeatures(seuratObj, selection.method = "vst", nfeatures = 2000, verbose = FALSE)
seuratObj <- ScaleData(seuratObj)
seuratObj <- RunPCA(seuratObj, features = VariableFeatures(object = seuratObj))
seuratObj <- FindNeighbors(seuratObj, dims = 1:30, verbose = FALSE)
seuratObj <- FindClusters(seuratObj, resolution = setresolu, verbose = FALSE)
seuratObj <- RunUMAP(seuratObj, dims = 1:30, verbose = FALSE)
DimPlot(seuratObj, reduction = "umap")
```

```{r}
sim_mat = SPARSim_result$count_matrix
origLabels = colnames(sim_mat) %>% str_extract("cluster.")
save(sim_mat, origLabels, file = str_c("./data/cells_pool", as.character(cluster_num), "_SPARSim.RData"))
```