PBMC3_filtering.qmd

---
title: Filtering of PBMC3 using COTAN
---

Library import

```{r}
library(dplyr)
library(COTAN)
library(Seurat)
library(tibble)
library(ggplot2)
library(zeallot)
library(DropletUtils)
```

Settings

```{r}
datasetName = 'PBMC3'
datasetFolder = './data/'

inDir  = paste(datasetFolder, datasetName, '/raw/10X/', sep='')
outDir = paste(datasetFolder, datasetName, '/filtered/', sep='')
dir10X = paste(outDir, '10X/', sep='')

if (!dir.exists(outDir)) {
  dir.create(outDir, recursive = TRUE, showWarnings = FALSE)
}

setLoggingLevel(2)
setLoggingFile(paste(outDir, "logfile.log", sep=""))
options(parallelly.fork.enable = TRUE)
```

Data loading

```{r}
dataset = Read10X(data.dir = inDir, strip.suffix = TRUE)
dataset = dataset[[1]]
sampleCond <- datasetName
PBMC3 <- COTAN(raw = dataset)
PBMC3 <- initializeMetaDataset(
  PBMC3,
  GEO = paste("10X ", datasetName, sep=""),
  sequencingMethod = "10X",
  sampleCond = sampleCond
)
```

Inspect cells’ sizes

```{r}
cellSizePlot(PBMC3)
```

Drop cells with too many reads as they are probably doublets

```{r}
cellsSizeThr <- 20000
PBMC3 <- addElementToMetaDataset(PBMC3, "Cells size threshold", cellsSizeThr)

cellsToRem <- getCells(PBMC3)[getCellsSize(PBMC3) > cellsSizeThr]
PBMC3 <- dropGenesCells(PBMC3, cells = cellsToRem)

cellSizePlot(PBMC3, splitPattern = "-", numCol = 2)
```

Inspect the number of expressed genes per cell

```{r}
genesSizePlot(PBMC3, splitPattern = "-", numCol = 2)
```

Drop cells with too high genes expression as they are probably doublets

```{r}
geneSizeThr <- 3500
PBMC3 <- addElementToMetaDataset(PBMC3, "Num genes threshold", geneSizeThr)

numExprGenes <- getNumExpressedGenes(PBMC3)
cellsToRem <- names(numExprGenes)[numExprGenes > geneSizeThr]
PBMC3 <- dropGenesCells(PBMC3, cells = cellsToRem)

genesSizePlot(PBMC3, splitPattern = "-", numCol = 2)
```

Check number of mithocondrial genes expressed in each cell

```{r}
mitGenesPattern <- "^[Mm][Tt]-"
getGenes(PBMC3)[grep(mitGenesPattern, getGenes(PBMC3))]
```

```{r}
c(mitPlot, mitSizes) %<-%
  mitochondrialPercentagePlot(PBMC3, genePrefix = mitGenesPattern,
                              splitPattern = "-", numCol = 2)

plot(mitPlot)
```

We drop cells with a too high percentage of mitocondrial genes (are likely dead)

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

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

PBMC3 <- dropGenesCells(PBMC3, cells = cellsToRem)

c(mitPlot, mitSizes) %<-%
  mitochondrialPercentagePlot(PBMC3, genePrefix = mitGenesPattern,
                              splitPattern = "-", numCol = 2)

plot(mitPlot)
```

Check number of ribosomial genes expressed in each cell

```{r}
ribGenesPattern <- "^RP[SL]\\d+"
getGenes(PBMC3)[grep(ribGenesPattern, getGenes(PBMC3))]
```

```{r}
c(ribPlot, ribSizes) %<-%
  mitochondrialPercentagePlot(PBMC3, genePrefix = ribGenesPattern,
                              splitPattern = "-", numCol = 2)

plot(ribPlot)
```

Check no further outliers after all the culling

```{r}
cellSizePlot(PBMC3, splitPattern = "-", numCol = 2)
```

```{r}
genesSizePlot(PBMC3, splitPattern = "-", numCol = 2)
```

Cleaning, round 1

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

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

plot(pcaCellsPlot)
```

```{r}
plot(genesPlot)
```

```{r}
PBMC3 <- addElementToMetaDataset(PBMC3, "Num drop B group", 0)
```

```{r}
plot(UDEPlot)
```

```{r}
plot(nuPlot)
```

```{r}
plot(zoomedNuPlot)  
```

```{r}
yset=0.16
nuDf <- data.frame("nu" = sort(getNu(PBMC3)), "n" = seq_along(getNu(PBMC3)))
PBMC3 <- addElementToMetaDataset(PBMC3, "Threshold low UDE cells:", yset)
cellsToRem <-rownames(nuDf)[nuDf[["nu"]] < yset]
PBMC3 <- dropGenesCells(PBMC3, cells = cellsToRem)
```

Cleaning, round 2

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

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

plot(pcaCellsPlot)
```

```{r}
plot(pcaCellsData)
```

```{r}
plot(genesPlot)
```

```{r}
plot(UDEPlot)
```

```{r}
plot(nuPlot)
```

```{r}
plot(zoomedNuPlot)
```


```{r}
plot(cellSizePlot(PBMC3, splitPattern = "-", numCol = 2))
```

```{r}
plot(genesSizePlot(PBMC3, splitPattern = "-", numCol = 2))
```

Save the filtered dataset

```{r}
if (!dir.exists(dir10X)) {
  write10xCounts(dir10X, getRawData(PBMC3))
}
```

```{r}
saveRDS(PBMC3, file = paste0(outDir, sampleCond, ".cotan.RDS"))
```

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