Monocle_E175.qmd

---
title: Clustering and differential expression analysis using Monocle
---

```{r}
library(monocle3)
library(rjson)
source("./utils.R")
```

Data loading and preparation

```{r}
numTopMarkers = 500
minClusterSize = 0

datasetName = "E175"
datasetPath= paste("./data/", datasetName, sep='')

outDir = paste(datasetPath, '/monocle/', sep='')
outDirDefault = paste(outDir, 'default/', sep='')
outDirCelltypist = paste(outDir, 'celltypist/', sep='')

if (!dir.exists(outDirDefault)) {
  dir.create(outDirDefault, recursive = TRUE)
}
if (!dir.exists(outDirCelltypist)) {
  dir.create(outDirCelltypist, recursive = TRUE)
}

numClusterThreshold = 1/10
numClusterCelltypist = fromJSON(file=paste(datasetPath, '/nclusters.json', sep=''))$nclusters
minNumClusterCelltypist = round(numClusterCelltypist - numClusterThreshold*numClusterCelltypist)
maxNumClusterCelltypist = round(numClusterCelltypist + numClusterThreshold*numClusterCelltypist)
```

```{r}
matrix <- read.csv(paste(datasetPath,'CorticalCells_GSM2861514_E175_cleaned.csv', sep='/'), header=TRUE, row.names=1)
cells <- read.csv(paste(datasetPath, 'barcodes.csv', sep='/'), header=FALSE, row.names=1)
genes <- read.csv(paste(datasetPath, 'genes.csv', sep='/'), header=FALSE, row.names=1)
matrix <- as.matrix(matrix)
cds<- new_cell_data_set(
  matrix,
  cell_metadata = cells,
  gene_metadata = genes
)
```

Studying feature variance

```{r}
cds <- preprocess_cds(cds, num_dim = 10, method = "PCA") # normalization
plot_pc_variance_explained(cds)
```

Clustering with default parameters

```{r}
cds <- cluster_cells(cds, reduction_method = "PCA")
cat(paste("Number of clusters:", length(levels(clusters(cds, reduction_method = "PCA")))))
```

Save the default clustering results

```{r}
defaultLabelsDf = data.frame(clusters(cds, reduction_method = "PCA"))
defaultLabelsDf$cell <- rownames(defaultLabelsDf)
rownames(defaultLabelsDf) <- 1:nrow(defaultLabelsDf)
colnames(defaultLabelsDf)[1] <- "computed_id"
defaultLabelsDf$computed_id <- as.numeric(defaultLabelsDf$computed_id)
write_clustering(outDirDefault, defaultLabelsDf, "cell", "computed_id")
```


Differential expression on default clustering results

```{r}
defaultMarkers <- top_markers(
  cds,
  group_cells_by="cluster",
  genes_to_test_per_group = 3000,
  reduction_method='PCA',
  cores=10
)
```

Save markers

```{r}
write_markers(outDirDefault, defaultMarkers, "gene_id", "cell_group", "marker_test_p_value", FALSE, numTopMarkers)
```

Clustering tuning resolution according to celltypist

```{r}
cat("Clusters to find:")
numClusterCelltypist
```

```{r}
cds <- preprocess_cds(cds, num_dim = 10, method = "PCA") # normalization
plot_pc_variance_explained(cds)
```

```{r}
maxResolution = 3
minResolution = 0
repeat{
  resolution = (maxResolution + minResolution) / 2
  cat(paste("Trying resolution ", resolution, '\n', sep=''))
  cds <- cluster_cells(cds, reduction_method = "PCA", resolution=resolution)
  numClusters = length(levels(clusters(cds, reduction_method = "PCA")))
  cat(paste("Got ", numClusters, ' clusters', '\n', sep=''))
  if (numClusters >= minNumClusterCelltypist & numClusters <= maxNumClusterCelltypist){
    break
  }
  else if (numClusters < minNumClusterCelltypist){
    # increase the resolution
    minResolution = resolution
  }
  else{
    # reduce the resolution
    maxResolution = resolution
  }
}
```


Save the tuned clustering results

```{r}
celltypistLabelsDf = data.frame(clusters(cds, reduction_method = "PCA"))
celltypistLabelsDf$cell <- rownames(celltypistLabelsDf)
rownames(celltypistLabelsDf) <- 1:nrow(celltypistLabelsDf)
colnames(celltypistLabelsDf)[1] <- "computed_id"
celltypistLabelsDf$computed_id <- as.numeric(celltypistLabelsDf$computed_id)
write_clustering(outDirCelltypist, celltypistLabelsDf, "cell", "computed_id")
```

Differential expression on tuned clustering results

```{r}
celltypistMarkers <- top_markers(
  cds,
  group_cells_by="cluster",
  genes_to_test_per_group = 3000,
  reduction_method='PCA',
  cores=10
)
```

Save markers

```{r}
write_markers(outDirCelltypist, celltypistMarkers, "gene_id", "cell_group", "marker_test_p_value", FALSE, numTopMarkers)
```

```{r}
sessionInfo()
```