recluster microglia
Kennedi Todd & Kimberly Olney
3/24/2022
Notes about reclustering
https://github.com/satijalab/seurat/issues/1621 https://github.com/satijalab/seurat/issues/2340 https://github.com/satijalab/seurat/issues/1883
Setup
knitr::opts_knit$set(
root.dir = ".")library(dplyr)
library(dittoSeq)
library(ggrepel)
library(ggtree)
library(parallel)
library(plotly) # 3D plot
library(Seurat) # Idents()
library(SeuratDisk) # SaveH5Seurat()
library(tibble) # rownnames_to_column
library(harmony) # RunHarmony()
#options(mc.cores = detectCores() - 1)Read in object
pigs.merged <- readRDS("../../rObjects/brain_annotated.rds")
Idents(pigs.merged) <- pigs.merged$merged_clusters
DefaultAssay(pigs.merged) <- "RNA"
# natural-log
pigs.merged <- NormalizeData(pigs.merged, verbose = FALSE)
pigs.merged## An object of class Seurat
## 42773 features across 8934 samples within 2 assays
## Active assay: RNA (21805 features, 0 variable features)
## 1 other assay present: SCT
## 3 dimensional reductions calculated: pca, harmony, umapAnnotated UMAP
# UMAP
u1 <- DimPlot(object = pigs.merged,
reduction = "umap",
label = TRUE,
label.box = TRUE,
label.size = 4,
repel = TRUE)
u1
microglias before re-cluster
# subset
microglia <- subset(pigs.merged, merged_clusters == "microglia")
# UMAP of microglia only
microglia_colors <- c("gold")
u1 <- DimPlot(object = microglia,
reduction = "umap",
label = TRUE,
label.box = TRUE,
label.size = 4,
repel = TRUE,
cols = microglia_colors)
u1
# show microglia cells by treatment
u2 <- DimPlot(object = microglia,
reduction = "umap",
label = TRUE,
label.box = TRUE,
label.size = 4,
repel = TRUE,
group.by = "treatment",
cols = treatment_colors)
u2
## png
## 2## null device
## 1## pdf
## 2## null device
## 1Quality checks
Number of cells
# Visualize the number of cell counts per sample
data <- as.data.frame(table(microglia$sample))
colnames(data) <- c("sample","frequency")
ncells <- ggplot(data, aes(x = sample, y = frequency, fill = sample)) +
geom_col() +
theme_classic() +
geom_text(aes(label = frequency),
position=position_dodge(width=0.9),
vjust=-0.25) +
scale_fill_manual(values = sample_colors) +
# scale_y_continuous(breaks = seq(0,30000, by = 5000), limits = c(0,30000)) +
ggtitle("Cells per sample") +
theme(legend.position = "none") +
theme(axis.text.x = element_text(angle = 45, hjust=1))
ncells
## Top transcripts
df <- data.frame(row.names = rownames(microglia))
df$rsum <- rowSums(x = microglia, slot = "counts")
df$gene_name <- rownames(df)
df <- df[order(df$rsum,decreasing = TRUE),]
head(df, 10)## rsum gene_name
## KK-MALAT1 8629 KK-MALAT1
## ARHGAP24 3459 ARHGAP24
## UBE2E2 3375 UBE2E2
## LRMDA 3112 LRMDA
## CALCR 3093 CALCR
## PRKCB 2841 PRKCB
## GAB2 2563 GAB2
## SLC8A1 2544 SLC8A1
## RN18S 2318 RN18S
## AGAP1 2301 AGAP1Top variable genes
# Identify the most variable genes
microglia <- FindVariableFeatures(microglia,
selection.method = "vst", # default vst
nfeatures = 2000, # default 2000
verbose = FALSE)## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -2.472## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.30103## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 1.8489e-14# view top variable genes
top40 <- head(VariableFeatures(microglia), 40)
top40## [1] "SULF1" "CEMIP" "NRG1"
## [4] "RORA" "TMSB10" "REV1"
## [7] "RBMS3" "NKAIN2" "CTNND2"
## [10] "DLG2" "LSAMP" "BMP6"
## [13] "NRXN3" "ADAM12" "LRP1B"
## [16] "TCF4" "RBFOX1" "AUTS2"
## [19] "ENSSSCG00000038719" "PBX1" "DAPK1"
## [22] "ENSSSCG00000026043" "ROBO2" "OPCML"
## [25] "CADM2" "PARD3" "NTM"
## [28] "BNC2" "SVIL" "FBXL7"
## [31] "NLGN1" "NRXN1" "ENSSSCG00000037775"
## [34] "RPS8" "SLC1A2" "GREB1L"
## [37] "NRG3" "NAV2" "GPC5"
## [40] "RN18S"# plot variable features with labels
VarFeatPlot <- VariableFeaturePlot(microglia, cols = c("gray47","red"))
VarFeatPlotLabel <- LabelPoints(plot = VarFeatPlot,
points = top40, repel = TRUE, fontface="italic",
xnudge = 0, ynudge = 0, max.overlaps = 12)
VarFeatPlotLabel
## Median absolute deviation
Remove outliers with log-library size greater than 3 median absolute deviations (MADs) or below the median log-library size.
# MAD example
log.lib <- as.numeric(log10(microglia$nCount_RNA))
med <- median(log.lib)
abs.dev <- abs(log.lib - med)
mad <- median(abs.dev)
mad <- mad * 1.4826 # multiply by consistency cutoff
mad## [1] 0.2874665# stats function
mad <- mad(log.lib)
mad## [1] 0.2874665# remove outliers greater than 3 MADs
remove <- abs(log.lib - median(log.lib)) / mad(log.lib) > 3
table(remove)## remove
## FALSE TRUE
## 591 2Re-cluster
Split, transform
# split object by sample
Idents(microglia) <- microglia$sample
microglia.split <- SplitObject(microglia, split.by = "sample")
# SCTransform and regress percent.mt
microglia.split <- lapply(microglia.split, SCTransform, vars.to.regress = "percent.mt")## Calculating cell attributes from input UMI matrix: log_umi## Variance stabilizing transformation of count matrix of size 6812 by 166## Model formula is y ~ log_umi## Get Negative Binomial regression parameters per gene## Using 2000 genes, 166 cells##
|
| | 0%
|
|================== | 25%
|
|=================================== | 50%
|
|==================================================== | 75%
|
|======================================================================| 100%## Found 29 outliers - those will be ignored in fitting/regularization step## Second step: Get residuals using fitted parameters for 6812 genes##
|
| | 0%
|
|===== | 7%
|
|========== | 14%
|
|=============== | 21%
|
|==================== | 29%
|
|========================= | 36%
|
|============================== | 43%
|
|=================================== | 50%
|
|======================================== | 57%
|
|============================================= | 64%
|
|================================================== | 71%
|
|======================================================= | 79%
|
|============================================================ | 86%
|
|================================================================= | 93%
|
|======================================================================| 100%## Computing corrected count matrix for 6812 genes##
|
| | 0%
|
|===== | 7%
|
|========== | 14%
|
|=============== | 21%
|
|==================== | 29%
|
|========================= | 36%
|
|============================== | 43%
|
|=================================== | 50%
|
|======================================== | 57%
|
|============================================= | 64%
|
|================================================== | 71%
|
|======================================================= | 79%
|
|============================================================ | 86%
|
|================================================================= | 93%
|
|======================================================================| 100%## Calculating gene attributes## Wall clock passed: Time difference of 7.081976 secs## Determine variable features## Place corrected count matrix in counts slot## Regressing out percent.mt## Centering data matrix## Set default assay to SCT## Calculating cell attributes from input UMI matrix: log_umi## Variance stabilizing transformation of count matrix of size 6961 by 86## Model formula is y ~ log_umi## Get Negative Binomial regression parameters per gene## Using 2000 genes, 86 cells##
|
| | 0%
|
|================== | 25%
|
|=================================== | 50%
|
|==================================================== | 75%
|
|======================================================================| 100%## Found 42 outliers - those will be ignored in fitting/regularization step## Second step: Get residuals using fitted parameters for 6961 genes##
|
| | 0%
|
|===== | 7%
|
|========== | 14%
|
|=============== | 21%
|
|==================== | 29%
|
|========================= | 36%
|
|============================== | 43%
|
|=================================== | 50%
|
|======================================== | 57%
|
|============================================= | 64%
|
|================================================== | 71%
|
|======================================================= | 79%
|
|============================================================ | 86%
|
|================================================================= | 93%
|
|======================================================================| 100%## Computing corrected count matrix for 6961 genes##
|
| | 0%
|
|===== | 7%
|
|========== | 14%
|
|=============== | 21%
|
|==================== | 29%
|
|========================= | 36%
|
|============================== | 43%
|
|=================================== | 50%
|
|======================================== | 57%
|
|============================================= | 64%
|
|================================================== | 71%
|
|======================================================= | 79%
|
|============================================================ | 86%
|
|================================================================= | 93%
|
|======================================================================| 100%## Calculating gene attributes## Wall clock passed: Time difference of 5.297404 secs## Determine variable features## Place corrected count matrix in counts slot## Regressing out percent.mt## Centering data matrix## Set default assay to SCT## Calculating cell attributes from input UMI matrix: log_umi## Variance stabilizing transformation of count matrix of size 1939 by 45## Model formula is y ~ log_umi## Get Negative Binomial regression parameters per gene## Using 1939 genes, 45 cells##
|
| | 0%
|
|================== | 25%
|
|=================================== | 50%
|
|==================================================== | 75%
|
|======================================================================| 100%## Found 23 outliers - those will be ignored in fitting/regularization step## Second step: Get residuals using fitted parameters for 1939 genes##
|
| | 0%
|
|================== | 25%
|
|=================================== | 50%
|
|==================================================== | 75%
|
|======================================================================| 100%## Computing corrected count matrix for 1939 genes##
|
| | 0%
|
|================== | 25%
|
|=================================== | 50%
|
|==================================================== | 75%
|
|======================================================================| 100%## Calculating gene attributes## Wall clock passed: Time difference of 3.296428 secs## Determine variable features## Place corrected count matrix in counts slot## Regressing out percent.mt## Centering data matrix## Set default assay to SCT## Calculating cell attributes from input UMI matrix: log_umi## Variance stabilizing transformation of count matrix of size 7569 by 296## Model formula is y ~ log_umi## Get Negative Binomial regression parameters per gene## Using 2000 genes, 296 cells##
|
| | 0%
|
|================== | 25%
|
|=================================== | 50%
|
|==================================================== | 75%
|
|======================================================================| 100%## Found 93 outliers - those will be ignored in fitting/regularization step## Second step: Get residuals using fitted parameters for 7569 genes##
|
| | 0%
|
|==== | 6%
|
|========= | 12%
|
|============= | 19%
|
|================== | 25%
|
|====================== | 31%
|
|========================== | 38%
|
|=============================== | 44%
|
|=================================== | 50%
|
|======================================= | 56%
|
|============================================ | 62%
|
|================================================ | 69%
|
|==================================================== | 75%
|
|========================================================= | 81%
|
|============================================================= | 88%
|
|================================================================== | 94%
|
|======================================================================| 100%## Computing corrected count matrix for 7569 genes##
|
| | 0%
|
|==== | 6%
|
|========= | 12%
|
|============= | 19%
|
|================== | 25%
|
|====================== | 31%
|
|========================== | 38%
|
|=============================== | 44%
|
|=================================== | 50%
|
|======================================= | 56%
|
|============================================ | 62%
|
|================================================ | 69%
|
|==================================================== | 75%
|
|========================================================= | 81%
|
|============================================================= | 88%
|
|================================================================== | 94%
|
|======================================================================| 100%## Calculating gene attributes## Wall clock passed: Time difference of 10.12541 secs## Determine variable features## Place corrected count matrix in counts slot## Regressing out percent.mt## Centering data matrix## Set default assay to SCT# Choose the features to use when integrating multiple datasets.
# will use nfeatures as 3000 as defined by running SCTransform above
var.features <- SelectIntegrationFeatures(object.list = microglia.split,
nfeatures = 3000)
# Merge the split object
microglia.sct.merged <- merge(x = microglia.split[[1]],
y = c(microglia.split[[2]], microglia.split[[3]], microglia.split[[4]]),
project = "LPS Pigs microglias")
# Define the variable features
VariableFeatures(microglia.sct.merged) <- var.features
# Run PCA on the merged object
microglia.sct.merged <- RunPCA(object = microglia.sct.merged, assay = "SCT")## PC_ 1
## Positive: CALCR, UBE2E2, ENSSSCG00000023479, GPM6B, ARHGAP24, PDE3B, GAB2, DOCK4, ANK1, SLC8A1
## NAV3, FGF1, IL1RAPL2, HK2, ENSSSCG00000033041, ENSSSCG00000012238, TCF7L2, PDE7A, ABR, EPB41L2
## ENSSSCG00000015645, PIK3AP1, P2RY12, SMAP2, TMCC3, CLCN5, APP, ST18, INPP5D, SLCO2B1
## Negative: KK-MALAT1, NOL11, FOXP1, PRKCE, HIVEP2, ZBTB16, FNBP1, HSP90AA1, DOCK1, ITPR1
## MED13L, ENSSSCG00000033262, CYTH1, ENSSSCG00000050772, CADM1, PARD3, NEDD4L, SORBS1, MED13, SFMBT2
## SNX24, OSBPL8, SSBP3, BCAS3, ANTXR2, ZBTB20, CHST11, ENSSSCG00000009327, ENSSSCG00000031462, RALGAPA1
## PC_ 2
## Positive: SNX24, RBPJ, ASAP1, PTPN1, B4GALT5, ZEB2, CMIP, EXT1, CIITA, SNX9
## MS4A7, DENND5A, ANKHD1, PFKFB3, PLEK, LCP1, KK-MALAT1, ENSSSCG00000033909, TNS3, IRAK2
## KDM6A, JARID2, STAB1, GAB2, DIAPH2, KMT2C, PICALM, CD53, CDK12, RNF2
## Negative: ZBTB20, NEGR1, FRMD4A, TLN2, ENSSSCG00000011121, ENSSSCG00000050772, CAMTA1, DST, GPC5, NFIA
## PARD3, SORBS1, CALCR, DPYSL2, ENSSSCG00000048207, SLC1A3, JMJD1C, ENSSSCG00000009327, ANK1, ENSSSCG00000049691
## BCAS3, FOXP1, NEDD4L, SRPK2, OPHN1, FOXN3, DYNC1H1, ENSSSCG00000051610, SRGAP1, PDE4D
## PC_ 3
## Positive: ENSSSCG00000017146, MX2, FOXP1, NLRC5, VAV3, ZEB2, HERC6, OAS2, EPSTI1, PARP14
## LRMDA, GNAQ, ARHGAP22, ANTXR2, ITSN1, ENSSSCG00000009240, ENTPD1, QKI, ENSSSCG00000033089, ITPR1
## SNX24, ARHGAP15, ARHGAP25, PARP12, AP3B1, FYB1, HERC5, RBM47, DNAJC13, MX1
## Negative: ASAP1, ST18, RNF2, ENSSSCG00000048263, NAV3, C1orf21, FGF1, MAML3, PAFAH1B1, ENSSSCG00000009327
## ST3GAL3, SLC7A1, EPS8, ZEB1, SBNO2, HK2, KLF12, ENSSSCG00000033041, PPARG, PIK3CD
## CAMKMT, TOR1AIP1, ITGA9, TPM3, STAT5B, PTK2B, CMSS1, ENSSSCG00000047147, DNAJC7, MICU2
## PC_ 4
## Positive: TCF7L2, ASAP1, CIITA, ENSSSCG00000033089, HERC6, PLD1, ENSSSCG00000017146, RNF2, NIBAN1, ENSSSCG00000048263
## EPSTI1, MAML3, NEBL, PTK2B, PARP14, ZNF146, ENSSSCG00000030801, MX2, ETV6, EPS8
## KLF12, PRKCE, STK3, PIK3C2A, MX1, PIK3CD, ZBTB20, C3, SNTB1, SLC11A2
## Negative: DOCK2, IFNAR2, XYLT1, ENSSSCG00000000530, ARHGAP25, HPCAL1, BLNK, FYB1, RPS6KA3, PRKCB
## ENSSSCG00000028461, TGFBR2, ARID3A, ARHGAP22, RBM47, ENSSSCG00000033041, ZMYND8, SIN3A, CST3, ARHGAP18
## EHBP1L1, ZMIZ1, DGKD, DOCK8, RCSD1, ZEB2, ST6GAL1, SNX29, ENSSSCG00000028035, RRBP1
## PC_ 5
## Positive: RBPJ, AGAP1, MERTK, LRMDA, SNX24, STAB1, SLC7A7, TTC7B, UBASH3B, SLC9A9
## WDFY3, NCOR2, VAV3, PDGFC, NCOA2, CADM1, GMDS, TLN2, FRMD4B, DST
## HPCAL1, CTBP2, EXT1, NEK6, PPP1R37, FRMD4A, CSF1R, DENND1A, ARHGAP18, ARHGEF10L
## Negative: ARHGAP15, RABGAP1L, PTPRC, FRYL, ZEB1, FOXN3, NLRC5, RUNX1, MBNL1, KLF12
## XYLT1, STAT1, P2RY6, CNOT4, TGFBR2, MSI2, FNBP1, LRP8, ITPR2, ZHX2
## RHOH, ENSSSCG00000033089, PHF20L1, NFKB1, STX11, PARP14, SNX29, SNTB1, PDE7A, ZC3HAV1No integration
# get significant PCs
stdv <- microglia.sct.merged[["pca"]]@stdev
sum.stdv <- sum(microglia.sct.merged[["pca"]]@stdev)
percent.stdv <- (stdv / sum.stdv) * 100
cumulative <- cumsum(percent.stdv)
co1 <- which(cumulative > 90 & percent.stdv < 5)[1]
co2 <- sort(which((percent.stdv[1:length(percent.stdv) - 1] -
percent.stdv[2:length(percent.stdv)]) > 0.1),
decreasing = T)[1] + 1
min.pc <- min(co1, co2)
min.pc## [1] 6# run umap
microglia.nointergration <- RunUMAP(microglia.sct.merged, dims = 1:min.pc, reduction = "pca")## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session## 19:18:04 UMAP embedding parameters a = 0.9922 b = 1.112## 19:18:04 Read 593 rows and found 6 numeric columns## 19:18:04 Using Annoy for neighbor search, n_neighbors = 30## 19:18:04 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 19:18:04 Writing NN index file to temp file /tmp/RtmpGS27MT/file104a52102ed67
## 19:18:04 Searching Annoy index using 1 thread, search_k = 3000
## 19:18:04 Annoy recall = 100%
## 19:18:05 Commencing smooth kNN distance calibration using 1 thread
## 19:18:07 Initializing from normalized Laplacian + noise
## 19:18:07 Commencing optimization for 500 epochs, with 20562 positive edges
## 19:18:09 Optimization finished# cluster
microglia.nointergration <- FindNeighbors(object = microglia.nointergration, dims = 1:min.pc) ## Computing nearest neighbor graph
## Computing SNN# Determine the clusters for various resolutions
microglia.nointergration <- FindClusters(object = microglia.nointergration,resolution = seq(0.1,0.8,by=0.1))## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 593
## Number of edges: 19402
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9000
## Number of communities: 1
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 593
## Number of edges: 19402
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8302
## Number of communities: 2
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 593
## Number of edges: 19402
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7669
## Number of communities: 2
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 593
## Number of edges: 19402
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7242
## Number of communities: 3
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 593
## Number of edges: 19402
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.6942
## Number of communities: 5
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 593
## Number of edges: 19402
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.6722
## Number of communities: 5
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 593
## Number of edges: 19402
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.6503
## Number of communities: 5
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 593
## Number of edges: 19402
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.6283
## Number of communities: 5
## Elapsed time: 0 secondsIdents(microglia.nointergration) <- "SCT_snn_res.0.1"
microglia.nointergration$seurat_clusters <- microglia.nointergration$SCT_snn_res.0.6UMAP no integration
u1 <- dittoDimPlot(object = microglia.nointergration,
var = "seurat_clusters",
reduction.use = "umap",
do.label = TRUE,
labels.highlight = FALSE)
u1
u2 <- dittoDimPlot(object = microglia.nointergration,
var = "seurat_clusters",
reduction.use = "umap",
do.label = TRUE,
split.by = "treatment",
labels.highlight = FALSE)
u2
# show microglia cells by treatment
u3 <- DimPlot(object = microglia.nointergration,
reduction = "umap",
label = FALSE,
repel = TRUE,
group.by = "treatment",
cols = treatment_colors)
u3
# show microglia cells by treatment
#remove(u1,u2,u3)Cells per cluster
# Cells per treatment per cluster
treatment_ncells <- FetchData(microglia.nointergration,
vars = c("ident", "treatment")) %>%
dplyr::count(ident,treatment) %>%
tidyr::spread(ident, n)
treatment_ncells## treatment 0
## 1 LPS 341
## 2 Saline 252# Cells per sample per cluster
sample_ncells <- FetchData(microglia.nointergration,
vars = c("ident", "sample")) %>%
dplyr::count(ident,sample) %>%
tidyr::spread(ident, n)
sample_ncells## sample 0
## 1 10.Saline 86
## 2 12.LPS 296
## 3 4.R.Saline 166
## 4 8.R.LPS 45# treatment
b1 <- microglia.nointergration@meta.data %>%
group_by(seurat_clusters, treatment) %>%
dplyr::count() %>%
group_by(seurat_clusters) %>%
dplyr::mutate(percent = 100*n/sum(n)) %>%
ungroup() %>%
ggplot(aes(x=seurat_clusters,y=percent, fill=treatment)) +
geom_col() +
scale_fill_manual(values = treatment_colors) +
ggtitle("Percentage of treatment per cluster")
b1
# sample
b2 <- microglia.nointergration@meta.data %>%
group_by(seurat_clusters, sample) %>%
dplyr::count() %>%
group_by(seurat_clusters) %>%
dplyr::mutate(percent = 100*n/sum(n)) %>%
ungroup() %>%
ggplot(aes(x=seurat_clusters,y=percent, fill=sample)) +
geom_col() +
#scale_fill_manual(values = sample_colors) +
ggtitle("Percentage of sample per cluster")
b2
#remove(b1, b2)Integrate with harmony
# Run harmony to harmonize over samples
microglia.integrated <- RunHarmony(object = microglia.sct.merged,
group.by.vars = "sample",
assay.use = "SCT",
reduction = "pca",
plot_convergence = TRUE)## Harmony 1/10## Harmony 2/10## Harmony 3/10## Harmony 4/10## Harmony 5/10## Harmony 6/10## Harmony converged after 6 iterations## Warning: Invalid name supplied, making object name syntactically valid. New
## object name is Seurat..ProjectDim.SCT.harmony; see ?make.names for more details
## on syntax validity
Find significant PCs
First metric
# Determine percent of variation associated with each PC
stdv <- microglia.integrated[["pca"]]@stdev
sum.stdv <- sum(microglia.integrated[["pca"]]@stdev)
percent.stdv <- (stdv / sum.stdv) * 100
# Calculate cumulative percents for each PC
cumulative <- cumsum(percent.stdv)
# Determine which PC exhibits cumulative percent greater than 90% and
# and % variation associated with the PC as less than 5
co1 <- which(cumulative > 90 & percent.stdv < 5)[1]
co1## [1] 45Second metric
# Determine the difference between variation of PC and subsequent PC
co2 <- sort(which(
(percent.stdv[1:length(percent.stdv) - 1] -
percent.stdv[2:length(percent.stdv)]) > 0.1),
decreasing = T)[1] + 1
# last point where change of % of variation is more than 0.1%.
co2## [1] 6Usually, we would choose the minimum of these two metrics as the PCs covering the majority of the variation in the data.
# Minimum of the two calculation
min.pc <- min(co1, co2)
min.pc## [1] 6Use min.pc we just calculated to generate the clusters. We can plot the elbow plot again and overlay the information determined using our metrics:
# Create a dataframe with values
plot_df <- data.frame(pct = percent.stdv,
cumu = cumulative,
rank = 1:length(percent.stdv))
# Elbow plot to visualize
ggplot(plot_df, aes(cumulative, percent.stdv, label = rank, color = rank > min.pc)) +
geom_text() +
geom_vline(xintercept = 90, color = "grey") +
geom_hline(yintercept = min(percent.stdv[percent.stdv > 5]), color = "grey") +
theme_bw()## Warning in min(percent.stdv[percent.stdv > 5]): no non-missing arguments to min;
## returning Inf
Run UMAP and Cluster
# Run UMAP
microglia.integrated <- RunUMAP(microglia.integrated,
dims = 1:min.pc,
reduction = "pca",
n.components = 3) # set to 3 to use with VR## 19:18:15 UMAP embedding parameters a = 0.9922 b = 1.112## 19:18:15 Read 593 rows and found 6 numeric columns## 19:18:15 Using Annoy for neighbor search, n_neighbors = 30## 19:18:15 Building Annoy index with metric = cosine, n_trees = 50## 0% 10 20 30 40 50 60 70 80 90 100%## [----|----|----|----|----|----|----|----|----|----|## **************************************************|
## 19:18:15 Writing NN index file to temp file /tmp/RtmpGS27MT/file104a57049adff
## 19:18:15 Searching Annoy index using 1 thread, search_k = 3000
## 19:18:15 Annoy recall = 100%
## 19:18:16 Commencing smooth kNN distance calibration using 1 thread
## 19:18:17 Initializing from normalized Laplacian + noise
## 19:18:17 Commencing optimization for 500 epochs, with 20562 positive edges
## 19:18:20 Optimization finished# Determine the K-nearest neighbor graph
microglia.integrated <- FindNeighbors(object = microglia.integrated,dims = 1:min.pc)## Computing nearest neighbor graph
## Computing SNN# Determine the clusters for various resolutions
microglia.unannotated <- FindClusters(object = microglia.integrated,resolution = seq(0.1,0.8,by=0.1))## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 593
## Number of edges: 19402
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9000
## Number of communities: 1
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 593
## Number of edges: 19402
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8302
## Number of communities: 2
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 593
## Number of edges: 19402
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7669
## Number of communities: 2
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 593
## Number of edges: 19402
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7242
## Number of communities: 3
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 593
## Number of edges: 19402
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.6942
## Number of communities: 5
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 593
## Number of edges: 19402
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.6722
## Number of communities: 5
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 593
## Number of edges: 19402
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.6503
## Number of communities: 5
## Elapsed time: 0 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 593
## Number of edges: 19402
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.6283
## Number of communities: 5
## Elapsed time: 0 secondsIdents(microglia.unannotated) <- "SCT_snn_res.0.6"
# save
saveRDS(microglia.unannotated,"../../rObjects/microglia_unannotated.rds")microglias after re-cluster
UMAP
DefaultAssay(microglia.unannotated) <- "RNA"
#microglia.unannotated <- NormalizeData(microglia.unannotated, verbose = FALSE)
microglia.unannotated$seurat_clusters <- microglia.unannotated$SCT_snn_res.0.6
u1 <- dittoDimPlot(object = microglia.unannotated,
var = "seurat_clusters",
reduction.use = "umap",
do.label = TRUE,
labels.highlight = FALSE)
u1
u2 <- dittoDimPlot(object = microglia.unannotated,
var = "seurat_clusters",
reduction.use = "umap",
do.label = TRUE,
split.by = "treatment",
labels.highlight = FALSE)
u2
# show microglia cells by treatment
u3 <- DimPlot(object = microglia.unannotated,
reduction = "umap",
label = FALSE,
repel = TRUE,
group.by = "treatment",
cols = treatment_colors)
u3
remove(u1,u2,u3)Cluster tree
cluster_colors <- c("gold","firebrick1","dodgerblue","green",
"cyan","chocolate4","gray40","purple", "blue")
microglia.unannotated <- BuildClusterTree(object = microglia.unannotated,
dims = 1:min.pc,
reorder = FALSE,
reorder.numeric = FALSE)
tree <- microglia.unannotated@tools$BuildClusterTree
tree$tip.label <- paste0("Cluster ", tree$tip.label)
p <- ggtree::ggtree(tree, aes(x, y)) +
scale_y_reverse() +
ggtree::geom_tree() +
ggtree::theme_tree() +
ggtree::geom_tiplab(offset = 1) +
ggtree::geom_tippoint(color = cluster_colors[1:length(tree$tip.label)], shape = 16, size = 5) +
coord_cartesian(clip = 'off') +
theme(plot.margin = unit(c(0,2.5,0,0), 'cm'))## Scale for 'y' is already present. Adding another scale for 'y', which will
## replace the existing scale.p
Cells per cluster
# Cells per treatment per cluster
treatment_ncells <- FetchData(microglia.unannotated,
vars = c("ident", "treatment")) %>%
dplyr::count(ident,treatment) %>%
tidyr::spread(ident, n)
treatment_ncells## treatment 0 1 2 3 4
## 1 LPS 88 111 57 45 40
## 2 Saline 95 17 60 46 34# Cells per sample per cluster
sample_ncells <- FetchData(microglia.unannotated,
vars = c("ident", "sample")) %>%
dplyr::count(ident,sample) %>%
tidyr::spread(ident, n)
sample_ncells## sample 0 1 2 3 4
## 1 10.Saline 34 6 17 27 2
## 2 12.LPS 77 96 42 41 40
## 3 4.R.Saline 61 11 43 19 32
## 4 8.R.LPS 11 15 15 4 NA# treatment
b1 <- microglia.unannotated@meta.data %>%
group_by(seurat_clusters, treatment) %>%
dplyr::count() %>%
group_by(seurat_clusters) %>%
dplyr::mutate(percent = 100*n/sum(n)) %>%
ungroup() %>%
ggplot(aes(x=seurat_clusters,y=percent, fill=treatment)) +
geom_col() +
scale_fill_manual(values = treatment_colors) +
ggtitle("Percentage of treatment per cluster")
b1
# sample
b2 <- microglia.unannotated@meta.data %>%
group_by(seurat_clusters, sample) %>%
dplyr::count() %>%
group_by(seurat_clusters) %>%
dplyr::mutate(percent = 100*n/sum(n)) %>%
ungroup() %>%
ggplot(aes(x=seurat_clusters,y=percent, fill=sample)) +
geom_col() +
#scale_fill_manual(values = sample_colors) +
ggtitle("Percentage of sample per cluster")
b2
remove(b1, b2)Conserved markers
# set levels
microglia.unannotated$treatment <- factor(microglia.unannotated$treatment,
levels = c("LPS","Saline"))
microglia.unannotated$SCT_snn_res.0.6 <- factor(microglia.unannotated$SCT_snn_res.0.6)
# initialize df
conserved.markers <- data.frame()
all.clusters <- unique(microglia.unannotated$SCT_snn_res.0.6)
# loop through each cluster
for (i in all.clusters) {
# print the cluster you're on
print(i)
# find conserved marker in cluster vs all other clusters
markers <- FindConservedMarkers(microglia.unannotated,
ident.1 = i, # subset to single cluster
grouping.var = "treatment", # compare by treatment
only.pos = TRUE, # only positive markers
min.pct = 0.1, # default
logfc.threshold = 0.25, # default
test.use = "MAST"
)
# skip if none
if(nrow(markers) == 0) {
next
}
# make rownames a column
markers <- rownames_to_column(markers, var = "gene")
# make cluster number a column
markers$cluster <- i
# add to final table
conserved.markers <- rbind(conserved.markers, markers)
}## [1] "3"## Testing group Saline: (3) vs (4, 0, 2, 1)##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## Testing group LPS: (3) vs (1, 0, 2, 4)##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## [1] "4"## Testing group Saline: (4) vs (3, 0, 2, 1)
##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## Testing group LPS: (4) vs (1, 0, 2, 3)##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## [1] "0"## Testing group Saline: (0) vs (3, 4, 2, 1)
##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## Testing group LPS: (0) vs (1, 2, 3, 4)##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## [1] "2"## Testing group Saline: (2) vs (3, 4, 0, 1)
##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## Testing group LPS: (2) vs (1, 0, 3, 4)##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## [1] "1"## Testing group Saline: (1) vs (3, 4, 0, 2)
##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## Testing group LPS: (1) vs (0, 2, 3, 4)##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!# create delta_pct
conserved.markers$Saline_delta_pct <- abs(conserved.markers$Saline_pct.1 -
conserved.markers$Saline_pct.2)
conserved.markers$LPS_delta_pct <- abs(conserved.markers$LPS_pct.1 -
conserved.markers$LPS_pct.2)
# more stringent filtering
markers.strict <- conserved.markers[
conserved.markers$Saline_delta_pct > summary(conserved.markers$Saline_delta_pct)[5],]
markers.strict <- conserved.markers[
conserved.markers$LPS_delta_pct > summary(conserved.markers$LPS_delta_pct)[5],]
markers.strict$gene_name <- markers.strict$gene
markers.strict$row.num <- 1:nrow(markers.strict)
# compare
table(conserved.markers$cluster)##
## 0 1 2 3 4
## 225 207 231 474 382table(markers.strict$cluster)##
## 0 1 2 3 4
## 19 9 11 266 75saveRDS(conserved.markers, paste0("../../rObjects/brain_conserved_markers_", cell, ".rds"))
# subset
cluster0 <- markers.strict[markers.strict$cluster == 0,]
cluster1 <- markers.strict[markers.strict$cluster == 1,]
cluster2 <- markers.strict[markers.strict$cluster == 2,]
cluster3 <- markers.strict[markers.strict$cluster == 3,]
cluster4 <- markers.strict[markers.strict$cluster == 4,]Cluster Annotations
Cluster 0
VlnPlot(microglia.unannotated,
features = cluster0$gene[1:19],
split.by = "seurat_clusters",
cols = cluster_colors,
flip = TRUE,
stack = TRUE)## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##
## This message will be shown once per session.
### Cluster 1
VlnPlot(microglia.unannotated,
features = cluster1$gene[1:9],
split.by = "seurat_clusters",
cols = cluster_colors,
flip = TRUE,
stack = TRUE)
### Cluster 2
VlnPlot(microglia.unannotated,
features = cluster2$gene[1:11],
split.by = "seurat_clusters",
cols = cluster_colors,
flip = TRUE,
stack = TRUE)
Cluster 3
VlnPlot(microglia.unannotated,
features = cluster3$gene[1:21],
split.by = "seurat_clusters",
cols = cluster_colors,
flip = TRUE,
stack = TRUE)
VlnPlot(microglia.unannotated,
features = cluster3$gene[21:40],
split.by = "seurat_clusters",
cols = cluster_colors,
flip = TRUE,
stack = TRUE)
### Cluster 4
VlnPlot(microglia.unannotated,
features = cluster4$gene[1:20],
split.by = "seurat_clusters",
cols = cluster_colors,
flip = TRUE,
stack = TRUE)
VlnPlot(microglia.unannotated,
features = cluster4$gene[21:40],
split.by = "seurat_clusters",
cols = cluster_colors,
flip = TRUE,
stack = TRUE)
# microglia markers
# Annotation
# umap with annotations Differential expression
LPS vs Saline within each cluster
cell_types <- unique(microglia.unannotated$SCT_snn_res.0.6)
microglia.unannotated$cell_type <- microglia.unannotated$SCT_snn_res.0.6
DE.df <- data.frame()
for (i in cell_types) {
print(i)
cluster <- subset(microglia.unannotated, cell_type == i)
Idents(cluster) <- cluster$treatment
markers <- FindMarkers(object = cluster,
ident.1 = "LPS",
ident.2 = "Saline",
only.pos = FALSE, # default
min.pct = 0.10, # default
test.use = "MAST",
verbose = TRUE,
assay = "RNA")
if(nrow(markers) == 0) {
next
}
markers$cluster <- i
DE.df <- rbind(DE.df, markers)
}## [1] "3"##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## [1] "4"##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## [1] "0"##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## [1] "2"##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!## [1] "1"##
## Done!## Combining coefficients and standard errors## Calculating log-fold changes## Calculating likelihood ratio tests## Refitting on reduced model...##
## Done!# change column names
colnames(DE.df)[c(3,4)] <- c("percent_LPS","percent_Saline")
# add gene name column
DE.df$gene <- rownames(DE.df)
# change row names
rownames(DE.df) <- 1:nrow(DE.df)
# calculate percent difference
DE.df$percent_difference <- abs(DE.df$percent_LPS - DE.df$percent_Saline)
# reorder columns
DE.df <- DE.df[,c(6,7,1,5,2,3,4,8)]
# write table
write.table(DE.df, "../../results/recluster/microglia/LPS_vs_Saline_DEGs.tsv", sep = "\t",
quote = FALSE, row.names = FALSE)Shiny App
Cleanup object
metadata <- microglia.unannotated@meta.data
microglia.unannotated@meta.data <- metadata
microglia.unannotated@assays$SCT@meta.features <- metadata
microglia.unannotated@assays$RNA@meta.features <- metadataOutput folder
# load libraries
library(ShinyCell)
# make shiny folder
DefaultAssay(microglia.unannotated) <- "RNA"
Idents(microglia.unannotated) <- "merged_clusters"
sc.config <- createConfig(microglia.unannotated)
makeShinyApp(microglia.unannotated, sc.config, gene.mapping = TRUE,
shiny.title = "recluster microglia")