CorticalCells_GSM2861514_E175-analysis.qmd

---
title: "Cortical cells DGE E17.5 Data-set Anaysis"
---

```{r, warning=FALSE}
library(ggplot2)
library(tibble)
library(zeallot)
library(COTAN)

options(parallelly.fork.enable = TRUE)

outDir <- "Data/Yuzwa_MouseCortex/"

setLoggingLevel(1)
setLoggingFile(file.path(outDir, "CorticalCells_GSM2861514_E175-analysis.log"))
```

Read the dataset from file

```{r}
dataset <- read.csv(file.path(outDir <- "Data/Yuzwa_MouseCortex/"
, "GSM2861514_E175_Only_Cortical_Cells_DGE.txt"),
       header = TRUE, sep = "\t", strip.white = TRUE,row.names = 1)
```

## Cleaning

Crate the COTAN object

```{r}
sampleCondition <- "CorticalCells_GSM2861514_E175"
cc175Obj <- COTAN(raw = dataset)
cc175Obj <- initializeMetaDataset(cc175Obj,
                                  GEO = "GSM2861511_E135",
                                  sequencingMethod = "DropSeq",
                                  sampleCondition = sampleCondition)
```

Inspect cells' sizes

```{r}
cellSizePlot(cc175Obj)
```

Drop cells with too many ritz reads as they are probably duplets

```{r}
cellsSizeThr <- 6500
cc175Obj <- addElementToMetaDataset(cc175Obj, "Cells size threshold", cellsSizeThr)

cells_to_rem <- getCells(cc175Obj)[getCellsSize(cc175Obj) > cellsSizeThr]
cc175Obj <- dropGenesCells(cc175Obj, cells = cells_to_rem)

cellSizePlot(cc175Obj)
```

Inspect the number of expressed genes per cell

```{r}
genesSizePlot(cc175Obj)
```

Drop cells with too high genes expession as they are probably duplets

```{r}
genesSizeThr <- 3000
cc175Obj <- addElementToMetaDataset(cc175Obj, "Num genes threshold", genesSizeThr)

numExprGenes <- getNumExpressedGenes(cc175Obj)
cells_to_rem <- names(numExprGenes)[numExprGenes > genesSizeThr]
cc175Obj <- dropGenesCells(cc175Obj, cells = cells_to_rem)

genesSizePlot(cc175Obj)
```

Check number of mitocondrial genes expressed in each cell

```{r}
mitGenesPattern <- "^mt-"
getGenes(cc175Obj)[grep(mitGenesPattern, getGenes(cc175Obj))]

c(mitPlot, mitSizes) %<-%
  mitochondrialPercentagePlot(cc175Obj, genePrefix = mitGenesPattern)

plot(mitPlot)
```

Cells with a too high percentage of mitocondrial genes are likely dead (or at the last problematic) cells. So we drop them!

```{r}
mitPercThr <- 5.0
cc175Obj <- addElementToMetaDataset(cc175Obj, "Mitoc. perc. threshold", mitPercThr)

cells_to_rem <- rownames(mitSizes)[mitSizes[["mit.percentage"]] > mitPercThr]

cc175Obj <- dropGenesCells(cc175Obj, cells = cells_to_rem)

c(mitPlot, mitSizes) %<-%
  mitochondrialPercentagePlot(cc175Obj, genePrefix = mitGenesPattern)

plot(mitPlot)
```

Check no further outliers after all the culling

```{r}
cellSizePlot(cc175Obj)
genesSizePlot(cc175Obj)
```

### Clean: round 1

```{r}
cc175Obj <- clean(cc175Obj)

c(pcaCellsPlot, pcaCellsData, genesPlot,
  UDEPlot, nuPlot, zoomedNuPlot) %<-% cleanPlots(cc175Obj)

plot(pcaCellsPlot)
plot(genesPlot)

cc175Obj <- addElementToMetaDataset(cc175Obj, "Num drop B group", 0)
```

B group contains high number of hemoglobin genes: so they are not interesting

```{r}
cells_to_rem <- rownames(pcaCellsData)[pcaCellsData[["groups"]] == "B"]

cc175Obj <- dropGenesCells(cc175Obj, cells = cells_to_rem)
```

### Clean: round 2

```{r}
cc175Obj <- clean(cc175Obj)

c(pcaCellsPlot, pcaCellsData, genesPlot,
  UDEPlot, nuPlot, zoomedNuPlot) %<-% cleanPlots(cc175Obj)

plot(pcaCellsPlot)
plot(genesPlot)

cc175Obj <- addElementToMetaDataset(cc175Obj, "Num drop B group", 1)
```

Visualize if all is ok:

```{r}
plot(pcaCellsData)

plot(UDEPlot)

plot(nuPlot)
plot(zoomedNuPlot)
```

Drop very low UDE cells as they are likely outliers

```{r}
lowUDEThr <- 0.2
cc175Obj <- addElementToMetaDataset(cc175Obj, "Low UDE threshold", lowUDEThr)

nuDf <- data.frame("nu" = sort(getNu(cc175Obj)), "n" = seq_along(getNu(cc175Obj)))

cells_to_rem <- rownames(nuDf)[nuDf[["nu"]] < lowUDEThr]
cc175Obj <- dropGenesCells(cc175Obj, cells = cells_to_rem)
```

Final cleaning to check all is OK

```{r}
cc175Obj <- clean(cc175Obj)

c(pcaCellsPlot, pcaCellsData, genesPlot,
  UDEPlot, nuPlot, zoomedNuPlot) %<-% cleanPlots(cc175Obj)

plot(pcaCellsPlot)
plot(genesPlot)
plot(UDEPlot)
plot(nuPlot)

plot(cellSizePlot(cc175Obj))
plot(genesSizePlot(cc175Obj))
```

Calculate genes' COEX

```{r, eval=FALSE, echo=TRUE}
cc175Obj <- proceedToCoex(cc175Obj, calcCoex = TRUE, cores = 12,
                          saveObj = TRUE, outDir = outDir)

```

## GDI

```{r, echo=FALSE,eval=TRUE}
cc175Obj <- readRDS(file.path(outDir, paste0(sampleCondition, ".cotan.RDS")))
```

```{r}
gdiData <- calculateGDI(cc175Obj)

genesToLabel <- head(rownames(gdiData[order(gdiData[["GDI"]],
                                            decreasing = TRUE), ]), n = 10L)

genesToLabel

gdiPlot <- GDIPlot(cc175Obj, GDIIn = gdiData, GDIThreshold = 1.4,
                   genes = list("Top 10 GDI genes" = genesToLabel))

plot(gdiPlot)
```

Save the COTAN object

```{r, eval=FALSE, echo=TRUE}
saveRDS(cc175Obj, file = file.path(outDir, paste0(sampleCondition, ".cotan.RDS")))
```

## Consistent Transcript Cohorts (clustering)

```{r, eval=FALSE, echo=TRUE}
c(splitClusters, splitCoexDF) %<-%
  cellsUniformClustering(cc175Obj, GDIThreshold = 1.4, cores = 13,
                         saveObj = TRUE, outDir = outDir)


cc175Obj <- addClusterization(cc175Obj, clName = "split",
                              clusters = splitClusters,
                              coexDF = splitCoexDF, override = TRUE)

table(splitClusters)
```

```{r, eval=FALSE,echo=TRUE}
saveRDS(cc175Obj, file = file.path(outDir, paste0(sampleCondition, ".cotan.RDS")))
```

```{r, eval=FALSE, echo=TRUE}
c(mergedClusters, mergedCoexDF) %<-%
  mergeUniformCellsClusters(cc175Obj, clusters = splitClusters,
                            GDIThreshold = 1.4, cores = 13,
                            saveObj = TRUE, outDir = outDir)

cc175Obj <- addClusterization(cc175Obj, clName = "merge",
                              clusters = mergedClusters,
                              coexDF = mergedCoexDF,
                              override = TRUE)

table(mergedClusters)
```

```{r, eval=FALSE}
saveRDS(cc175Obj, file = file.path(outDir, paste0(sampleCondition, ".cotan.RDS")))
```

------------------------------------------------------------------------

```{r}
Sys.time()
```

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