Set working directory

knitr::opts_knit$set(root.dir = ".")

Libraris, paths, colors

Read in object

## 
##           Astrocytes      Fibroblast-like         Interneurons 
##                 1273                  190                 1567 
## Microglia Macrophage              Neurons     Oligodendrocytes 
##                  424                 7725                  520 
##       Polydendrocyte 
##                  305

Subset

## 
##           Astrocytes      Fibroblast-like Microglia Macrophage 
##                 1273                  190                  424 
##     Oligodendrocytes       Polydendrocyte 
##                  520                  305
## 
## Interneurons      Neurons 
##         1567         7725

Re-process subsets

Glia

## PC_ 1 
## Positive:  NPAS3, PITPNC1, RPE65, PREX2, ENSSSCG00000059499 
## Negative:  CALCR, CXCL10, MX2, ENSSSCG00000031640, ENSSSCG00000000657 
## PC_ 2 
## Positive:  PITPNC1, SPARCL1, TRPM3, PREX2, RPE65 
## Negative:  PLP1, MBP, OPALIN, MAG, UGT8 
## PC_ 3 
## Positive:  APOA1, MBP, MAG, PLP1, CLDN11 
## Negative:  CALCR, ENSSSCG00000031640, CCL8, ENSSSCG00000033909, CSF1R 
## PC_ 4 
## Positive:  CEMIP, ADAM12, CFB, EYA1, PRDM6 
## Negative:  NPAS3, RPE65, SPARCL1, PREX2, ENSSSCG00000063471 
## PC_ 5 
## Positive:  PTGDS, OPALIN, PLEKHH1, ENSSSCG00000008115, TTYH2 
## Negative:  SNX22, BCAN, OPCML, GPR17, TNR

# Significant PCs

## [1] 9
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2712
## Number of edges: 86883
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9603
## Number of communities: 4
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2712
## Number of edges: 86883
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9366
## Number of communities: 7
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2712
## Number of edges: 86883
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9150
## Number of communities: 9
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2712
## Number of edges: 86883
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8937
## Number of communities: 9
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2712
## Number of edges: 86883
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8724
## Number of communities: 9
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2712
## Number of edges: 86883
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8513
## Number of communities: 9
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2712
## Number of edges: 86883
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8346
## Number of communities: 11
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2712
## Number of edges: 86883
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8184
## Number of communities: 11
## Elapsed time: 0 seconds

UMAP plot

# Save

Tree

# Marker exploration ## Find all markers

Top variable genes

# print top variable genes
top20 <- pigs.glia@assays$RNA@var.features[1:20]
top20
##  [1] "CXCL10"             "CLDN5"              "ENSSSCG00000002749"
##  [4] "APOA1"              "ISG15"              "CEMIP"             
##  [7] "B2M"                "GBP1"               "GPR17"             
## [10] "ENSSSCG00000037358" "ENSSSCG00000032591" "IGFBP7"            
## [13] "ENSSSCG00000038912" "ENSSSCG00000010044" "ENSSSCG00000060622"
## [16] "ENSSSCG00000038719" "SERPINE1"           "CD163"             
## [19] "CCL8"               "MGP"
## 
##    0    1    2    3    4    5    6    7    8    9   10 
##  611  508 1080  279  240  178 1933  249  113  361   24
## 
##   0   1   2   3   4   5   6   7   8   9  10 
##  34  62  84  58  43  43 697 100  50 198  24

Cells per cluster

n_cells <- FetchData(pigs.glia, 
                     vars = c("ident", "treatment")) %>%
  dplyr::count(ident,treatment) %>%
  tidyr::spread(ident, n)
n_cells
##   treatment   0   1   2   3   4   5   6   7  8  9 10
## 1     Ecoli 902 395 394 318 217 130 106 101 63 50 36

Cluster Annotation

Cluster 0 - astrocyte

SLC1A2 - astrocyte PITPNC1 - astrocyte GJA1 - astrocyte LRIG1 - bergmann glia / astrocyte

# Number of cells per condition
n_cells[,c(1,2)]
##   treatment   0
## 1     Ecoli 902
# UMAP with only cluster 0
DimPlot(object = subset(pigs.glia, seurat_clusters == "0"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = "#E69F00")

VlnPlot(pigs.glia,
        features = cluster0$gene[1:10],
        cols = colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

## Cluster 1 - oligodendrocyte PLP1 - oligodendrocyte MBP - oligodendrocyte OPALIN - oligodendrocyte ENPP2 - choroid plexus

# Number of cells per condition
n_cells[,c(1,3)]
##   treatment   1
## 1     Ecoli 395
DimPlot(object = subset(pigs.glia, seurat_clusters == "1"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = "#56B4E9")

VlnPlot(pigs.glia,
        features = cluster1$gene[1:10],
        cols = colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

## Cluster 2 - microglia CALCR - neuron CSF1R - microglia CCL8 - microglia C1QC - microglia

# Number of cells per condition
n_cells[,c(1,4)]
##   treatment   2
## 1     Ecoli 394
DimPlot(object = subset(pigs.glia, seurat_clusters == "2"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = "#009E73")

VlnPlot(pigs.glia,
        features = cluster2$gene[1:10],
        cols = colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

## Cluster 3 - astrocyte RPE65 - astrocyte SLC1A2 - astrocyte PITPNC1 - astrocyte / neuron PREX2 - PREX2

# Number of cells per condition
n_cells[,c(1,4)]
##   treatment   2
## 1     Ecoli 394
DimPlot(object = subset(pigs.glia, seurat_clusters == "3"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = "#F0E442")

VlnPlot(pigs.glia,
        features = cluster3$gene[1:10],
        cols = colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

## Cluster 4 - polydendrocyte PDGFRA - polydendrocyte LHFPL3 - polydendrocyte SNX22 - polydendrocyte

# Number of cells per condition
n_cells[,c(1,5)]
##   treatment   3
## 1     Ecoli 318
DimPlot(object = subset(pigs.glia, seurat_clusters == "4"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = "#0072B2")

VlnPlot(pigs.glia,
        features = cluster4$gene[1:10],
        cols = colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

## Cluster 5 - fibroblast-like CEMIP - fibroblast-like ADAM12 - fibroblast-like BMP6 - fibroblast-like

# Number of cells per condition
n_cells[,c(1,6)]
##   treatment   4
## 1     Ecoli 217
DimPlot(object = subset(pigs.glia, seurat_clusters == "5"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = "#D55E00")

VlnPlot(pigs.glia,
        features = cluster5$gene[1:10],
        cols = colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

## Cluster 6 - neuron KCNIP4 - purkinjeneuron LRRC7 - neuron_rora CNKSR2 - neuron

# Number of cells per condition
n_cells[,c(1,7)]
##   treatment   5
## 1     Ecoli 130
DimPlot(object = subset(pigs.glia, seurat_clusters == "6"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = "#CC79A7")

VlnPlot(pigs.glia,
        features = cluster6$gene[1:10],
        cols = colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

## Cluster 7 - poldendrocyte GPR17 - poldendrocyte FYN - poldendrocyte CDCP1 - neuron ENPP6 - poldendrocyte

# Number of cells per condition
n_cells[,c(1,8)]
##   treatment   6
## 1     Ecoli 106
DimPlot(object = subset(pigs.glia, seurat_clusters == "7"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = "#666666")

VlnPlot(pigs.glia,
        features = cluster7$gene[1:10],
        cols = colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

## Cluster 8 - fibroblast-like C4A - microglia MGP - fibroblast-like POSTN - fibroblast-like LCN2 - endothelial

# Number of cells per condition
n_cells[,c(1,9)]
##   treatment   7
## 1     Ecoli 101
DimPlot(object = subset(pigs.glia, seurat_clusters == "8"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = "#AD7700")

VlnPlot(pigs.glia,
        features = cluster8$gene[1:10],
        cols = colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

## Cluster 9 - immune response CLDN5 - endothelial CXCL10 - microglia GBP1 - macrophages APOA1 - ependyma

# Number of cells per condition
n_cells[,c(1,10)]
##   treatment  8
## 1     Ecoli 63
DimPlot(object = subset(pigs.glia, seurat_clusters == "9"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = "#1C91D4")

VlnPlot(pigs.glia,
        features = cluster9$gene[1:10],
        cols = colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

## Cluster 10 - noise? MBP - oligodendrocyte PDGFRA - polydendrocyte BCAN - astrocyte OLIG1 - polydendrocyte SLC35F1 - neuron

# Number of cells per condition
n_cells[,c(1,11)]
##   treatment  9
## 1     Ecoli 50
DimPlot(object = subset(pigs.glia, seurat_clusters == "10"),
        reduction = "umap", 
        label = TRUE,
        label.box = TRUE,
        label.size = 3,
        repel = TRUE,
        cols = "#007756")

VlnPlot(pigs.glia,
        features = cluster10$gene[1:10],
        cols = colors,
        stack = TRUE,
        flip = TRUE,
        split.by = "seurat_clusters")

## Assign identities

# Rename all identities
pigs.glia.annotated <- RenameIdents(object = pigs.glia, 
                               "0" = "Astrocytes",
                               "1" = "Oligodendrocytes",
                               "2" = "Microglia",
                               "3" = "Astrocytes",
                               "4" = "Polydendrocytes",
                               "5" = "Fibroblast",
                               "6" = "Neurons",
                               "7" = "Polydendrocytes",
                               "8" = "Fibroblast",
                               "9" = "Immune response",
                               "10" = "Noise?")
pigs.glia.annotated$annotated_clusters <- factor(Idents(pigs.glia.annotated),
                                             levels = c("Astrocytes",
                                                        "Fibroblast",
                                                        "Microglia",
                                                        "Neurons",
                                                        "Oligodendrocytes",
                                                        "Polydendrocytes", 
                                                        "Immune response", 
                                                        "Noise?"))
Idents(pigs.glia.annotated) <- "annotated_clusters"

UMAP - annotated

u1 <- DimPlot(pigs.glia.annotated,
              group.by = "annotated_clusters",
              cols = colors,
              raster = FALSE,
              label = FALSE) +
  ggtitle("Ecoli pig brain")
u1

## png 
##   2
# umap percent.mt split by sample
FeaturePlot(pigs.glia.annotated, 
            reduction = "umap", 
            features = "CCL8",
            split.by = "sample",
            min.cutoff = 'q10',
            label = TRUE)

n_cells <- FetchData(pigs.glia.annotated, 
                     vars = c("ident", "treatment")) %>%
  dplyr::count(ident,treatment) %>%
  tidyr::spread(ident, n)
n_cells
##   treatment Astrocytes Fibroblast Microglia Neurons Oligodendrocytes
## 1     Ecoli       1220        193       394     106              395
##   Polydendrocytes Immune response Noise?
## 1             318              50     36