madzafv/snRNASeqRA

1---2title: "RA snRNA-Seq data processing and clustering"3output:4  html_document:5    theme: united6    df_print: kable7date: 'Compiled: `r format(Sys.Date(), "%B %d, %Y")`'8---9 10```{r}11library(dplyr)12library(Seurat)13library(SeuratObject)14library(ggplot2)15library(gridExtra)16library(grid)17 18getwd();19workingDir = "C:/Users/madza/Experiments/RAsnRNA-Seq/madzafvShinnyApp1"20setwd(workingDir); 21 22# 1) Load Cell Ranger filtered matrix23data1 <- Read10X(data.dir = "C:/Users/madza/Experiments/RAsnRNA-Seq/filtered_feature_bc_matrix_run1")24data2 <- Read10X(data.dir = "C:/Users/madzaExperiments/RAsnRNA-Seq/filtered_feature_bc_matrix_run2")25 26seu1 <- CreateSeuratObject(counts = data1, project = "run1")27seu2 <- CreateSeuratObject(counts = data2, project = "run2")28 29seu1 <- RenameCells(seu1, add.cell.id = "run1")30seu2 <- RenameCells(seu2, add.cell.id = "run2")31 32pbmc_recalc <- merge(seu1, seu2)33 34# 2) Normalize 35pbmc_recalc <- NormalizeData(36  pbmc_recalc,37  normalization.method = "LogNormalize",38  scale.factor         = 10000,39  margin               = 1,40  verbose              = TRUE41)42 43# 3) Find most variable genes 44pbmc_recalc <- FindVariableFeatures(45  pbmc_recalc,46  selection.method    = "vst",47  nfeatures           = 2000,48  loess.span          = 0.3,49  clip.max            = "auto",50  mean.function       = ExpMean,  51  dispersion.function = LogVMR,   52  verbose             = TRUE53)54 55# 4) ScaleData on all genes56all.genes <- rownames(pbmc_recalc)57 58pbmc_recalc <- ScaleData(59  pbmc_recalc,60  features       = all.genes,61  vars.to.regress = NULL,62  model.use       = "linear",63  use.umi         = FALSE,64  do.scale        = TRUE,65  do.center       = TRUE,66  block.size      = 1000,67  verbose         = TRUE68)69 70# 5) RunPCA71pbmc_recalc <- RunPCA(72  pbmc_recalc,73  features        = VariableFeatures(pbmc_recalc),74  npcs            = 50,75  rev.pca         = FALSE,76  weight.by.var   = TRUE,77  verbose         = TRUE,78  ndims.print     = 1:5,79  nfeatures.print = 30,80  seed.use        = 42,81  reduction.name  = "pca",82  reduction.key   = "PC_"83)84 85# 6) FindNeighbors86pbmc_recalc <- FindNeighbors(87  pbmc_recalc,88  reduction = "pca",89  dims      = 1:1090)91 92# 7) FindClusters93pbmc_recalc <- FindClusters(94  pbmc_recalc,95  resolution = 0.996)97 98# 8) RunUMAP99set.seed(42)100 101pbmc_recalc <- RunUMAP(102  pbmc_recalc,103  dims                 = 1:10,104  reduction            = "pca",105  assay                = "RNA",106  slot                 = "data",107  umap.method          = "uwot",108  return.model         = FALSE,109  n.neighbors          = 30,110  n.components         = 2,111  metric               = "cosine",112  learning.rate        = 1,113  min.dist             = 0.3,114  spread               = 1,115  set.op.mix.ratio     = 1,116  local.connectivity   = 1,117  repulsion.strength   = 1,118  negative.sample.rate = 5,119  uwot.sgd             = FALSE,120  angular.rp.forest    = FALSE,121  densmap              = FALSE,122  dens.lambda          = 2,123  dens.frac            = 0.3,124  dens.var.shift       = 0.1,125  verbose              = TRUE,126  reduction.name       = "umap",127  reduction.key        = "UMAP_"128)129 130# 9) Rename cluster identities in the Seurat object131pbmc_recalc <- RenameIdents(132  pbmc_recalc,133  `0`  = "Astrocytes_1",134  `1`  = "Astrocytes_2",135  `2`  = "Excit_5",136  `3`  = "Excit_2.1",137  `4`  = "Olig",138  `5`  = "Excit_4.1",139  `6`  = "Excit_4.2",140  `7`  = "Inhib_1",141  `8`  = "Excit_1",142  `9`  = "Inhib_3",143  `10` = "Inhib_4",144  `11` = "Inhib_2",145  `12` = "Excit_2.2",146  `13` = "Excit_3",147  `14` = "Microg",148  `15` = "Endot",149  `16` = "Inhib_5"150)151 152# 9) Plot the UMAP153DimPlot(154  pbmc_recalc,155  reduction = "umap",156  label     = TRUE,157  repel = TRUE,158  label.box = TRUE159)160 161# 10) Compute markers using renamed identities162pbmc_recalc.markers <- FindAllMarkers(163  pbmc_recalc,164  only.pos = TRUE165)166 167# 11) ATP1A3 Dotplot168 169pbmc_recalc$celltype <- as.character(Idents(pbmc_recalc))170 171cluster_order_top_to_bottom <- c(172  "Endotelial", "Microglia", "Oligodendrocytes",173  "Astrocytes_2", "Astrocytes_1",174  "Inhib_5", "Inhib_4", "Inhib_3", "Inhib_2", "Inhib_1",175  "Excit_5", "Excit_4.2", "Excit_4.1", "Excit_3", "Excit_2.2", "Excit_2.1", "Excit_1"176)177 178pbmc_recalc$celltype <- factor(179  pbmc_recalc$celltype,180  levels = rev(cluster_order_top_to_bottom)181)182 183p <- DotPlot(184  pbmc_recalc,185  features = "ATP1A3",186  group.by = "celltype"187) +188  theme(189    legend.position = "bottom",190    aspect.ratio = 5,191    legend.title = element_text(size = 10),192    legend.text = element_text(size = 8)193  ) +194  guides(195    size = guide_legend(title = "% cells"),196    color = guide_colorbar(title = "Avg. expr")197  ) +198  labs(x = NULL, y = NULL) +199  theme(200    legend.position = "right",201    aspect.ratio = 5,202    legend.title = element_text(size = 10),203    legend.text = element_text(size = 8)204  )205 206dot_tbl <- p$data[, c("id", "features.plot", "pct.exp", "avg.exp", "avg.exp.scaled")]207 208dot_tbl$id <- factor(dot_tbl$id, levels = cluster_order_top_to_bottom)209 210dot_tbl <- dot_tbl[order(dot_tbl$id), ]211 212colnames(dot_tbl) <- c("Cell type", "Gene", "% cells", "Avg. exp", "Scaled avg. exp")213 214dot_tbl[, c("% cells", "Avg. exp", "Scaled avg. exp")] <- round(215  dot_tbl[, c("% cells", "Avg. exp", "Scaled avg. exp")],216  2217)218 219tbl_img <- dot_tbl[, c("Cell type", "% cells", "Avg. exp")]220 221table_grob <- tableGrob(222  tbl_img,223  rows = NULL,224  theme = ttheme_minimal(225    core = list(fg_params = list(cex = 0.8)),226    colhead = list(fg_params = list(fontface = "bold", cex = 0.9))227  )228)229```230 231