#copy seurat object to filter on prior to DE testing Astsaver2.seurat <- Astsaver.seurat #first remove samples that have no counterpart, and younger than 4 month samples, and isogenic pairs with < 30 cells in each sample Astsaver2.seurat <- SetIdent(Astsaver2.seurat, value = "Age") Astsaver2.seurat <- subset(Astsaver2.seurat, idents = c("1m", "2m", "3m"), invert = TRUE) Astsaver2.seurat <- SetIdent(Astsaver2.seurat, value = "SampleName2") Astsaver2.seurat <- subset(Astsaver2.seurat, idents = c("4mFB_337VV3B12_12", "4mFB_406WW1_12", "6mFB_337VV3B12_14", "6mFB_406WW1_14"), invert = TRUE) Astsaver2.seurat <- subset(Astsaver2.seurat, idents = c("4mFB_337VM2_12", "4mFB_337VV2_12", "4mFB_337VM3_12", "4mFB_337VV3F02_12", "4mFB_406RR1_12", "4mFB_406RW1_12", "6mFB_337VM2_14", "6mFB_337VV2_14" ), invert = TRUE) ####subsample samples that need cells removed to match their respective mut/ctrl sample ###4mFB_337VV1_05: remove 24 cells sampled.cells_4mFB_337VV1_05 <- subset(Astsaver2.seurat, idents = "4mFB_337VV1_05") set.seed(111) sampled.cells_4mFB_337VV1_05 <- sample(x = colnames(sampled.cells_4mFB_337VV1_05), size = 24, replace = F) ###4mFB_337VM1_07: remove 27 cells sampled.cells_4mFB_337VM1_07 <- subset(Astsaver2.seurat, idents = "4mFB_337VM1_07") set.seed(111) sampled.cells_4mFB_337VM1_07 <- sample(x = colnames(sampled.cells_4mFB_337VM1_07), size = 27, replace = F) ###4mFB_337VV1_13: remove 132 cells sampled.cells_4mFB_337VV1_13 <- subset(Astsaver2.seurat, idents = "4mFB_337VV1_13") set.seed(111) sampled.cells_4mFB_337VV1_13 <- sample(x = colnames(sampled.cells_4mFB_337VV1_13), size = 132, replace = F) ###4mFB_337VV2_06: remove 18 cells sampled.cells_4mFB_337VV2_06 <- subset(Astsaver2.seurat, idents = "4mFB_337VV2_06") set.seed(111) sampled.cells_4mFB_337VV2_06 <- sample(x = colnames(sampled.cells_4mFB_337VV2_06), size = 18, replace = F) ###6mFB_337VV1_07: remove 302 cells sampled.cells_6mFB_337VV1_07 <- subset(Astsaver2.seurat, idents = "6mFB_337VV1_07") set.seed(111) sampled.cells_6mFB_337VV1_07 <- sample(x = colnames(sampled.cells_6mFB_337VV1_07), size = 302, replace = F) ###6mFB_337VV1_15: remove 96 cells sampled.cells_6mFB_337VV1_15 <- subset(Astsaver2.seurat, idents = "6mFB_337VV1_15") set.seed(111) sampled.cells_6mFB_337VV1_15 <- sample(x = colnames(sampled.cells_6mFB_337VV1_15), size = 96, replace = F) ###6mFB_337VV2_08: remove 126 cells sampled.cells_6mFB_337VV2_08 <- subset(Astsaver2.seurat, idents = "6mFB_337VV2_08") set.seed(111) sampled.cells_6mFB_337VV2_08 <- sample(x = colnames(sampled.cells_6mFB_337VV2_08), size = 126, replace = F) ###6mFB_337VV3F02_14: remove 11 cells sampled.cells_6mFB_337VV3F02_14 <- subset(Astsaver2.seurat, idents = "6mFB_337VV3F02_14") set.seed(111) sampled.cells_6mFB_337VV3F02_14 <- sample(x = colnames(sampled.cells_6mFB_337VV3F02_14), size = 11, replace = F) ###6mFB_406RW1_06: remove 17 cells sampled.cells_6mFB_406RW1_06 <- subset(Astsaver2.seurat, idents = "6mFB_406RW1_06") set.seed(111) sampled.cells_6mFB_406RW1_06 <- sample(x = colnames(sampled.cells_6mFB_406RW1_06), size = 17, replace = F) ###6mFB_406RW1_14: remove 176 cells sampled.cells_6mFB_406RW1_14 <- subset(Astsaver2.seurat, idents = "6mFB_406RW1_14") set.seed(111) sampled.cells_6mFB_406RW1_14 <- sample(x = colnames(sampled.cells_6mFB_406RW1_14), size = 176, replace = F) ###8mFB_406RR1_10: remove 91 cells sampled.cells_8mFB_406RR1_10 <- subset(Astsaver2.seurat, idents = "8mFB_406RR1_10") set.seed(111) sampled.cells_8mFB_406RR1_10 <- sample(x = colnames(sampled.cells_8mFB_406RR1_10), size = 91, replace = F) ####remove sampled cells cells_to_remove <- c(sampled.cells_4mFB_337VV1_05, sampled.cells_4mFB_337VM1_07, sampled.cells_4mFB_337VV1_13, sampled.cells_4mFB_337VV2_06, sampled.cells_6mFB_337VV1_07, sampled.cells_6mFB_337VV1_15, sampled.cells_6mFB_337VV2_08, sampled.cells_6mFB_337VV2_14, sampled.cells_6mFB_337VV3F02_14, sampled.cells_6mFB_406RW1_06, sampled.cells_6mFB_406RW1_14, sampled.cells_8mFB_406RR1_10) Astsaver2.seurat <- subset(Astsaver2.seurat, cells = cells_to_remove, invert = TRUE) #clean up environment rm(list=ls(pattern="^sampled.cells")) rm(cells_to_remove) rm(plot1) rm(plot2) Astsaver2.seurat <- SetIdent(Astsaver2.seurat, value = "Genotype2") Ast_DEGs_MAST <- FindMarkers(object = Astsaver2.seurat, ident.1 = c("337VM", "406RW"), ident.2 = c("337VV", "406RR"), logfc.threshold = 0, min.pct = 0, test.use = "MAST") Ast_DEGs_MAST <- Ast_DEGs_MAST[order(Ast_DEGs_MAST$avg_logFC, decreasing = TRUE),] Ast_DEGs_MAST$p_val_adj_BH <- p.adjust(p = Ast_DEGs_MAST$p_val, method = "BH") Ast_DEGs_MAST$p_val_adj_BY <- p.adjust(p = Ast_DEGs_MAST$p_val, method = "BY") Ast_DEGs_MAST_sig <- Ast_DEGs_MAST[Ast_DEGs_MAST$p_val_adj_BH <= 0.05,] write.csv(Ast_DEGs_MAST_sig, file ='Ast_DEGs_MAST_sig.txt') #Volcano Plot library(EnhancedVolcano) keyvals <- ifelse( Ast_DEGs_MAST$avg_logFC < -0.05 & Ast_DEGs_MAST$p_val_adj_BH < 0.05, 'royalblue', ifelse(Ast_DEGs_MAST$avg_logFC > 0.05 & Ast_DEGs_MAST$p_val_adj_BH < 0.05, 'red2', 'darkgrey')) keyvals[is.na(keyvals)] <- 'darkgrey' names(keyvals)[keyvals == 'red2'] <- 'high' names(keyvals)[keyvals == 'darkgrey'] <- 'mid' names(keyvals)[keyvals == 'royalblue'] <- 'low' #plot without labels EnhancedVolcano(Ast_DEGs_MAST, lab = rownames(Ast_DEGs_MAST), x = 'avg_logFC', y = 'p_val_adj_BH', #xlim = c(-0.4, 0.4), #ylim = c(0,6), xlab = bquote('avg log fold change'), ylab = bquote('-log10(p adj) '), axisLabSize = 28, title = 'astrocytes all samples', pCutoff = 0.05, FCcutoff = 0.05, pointSize = 5, colCustom = keyvals, colAlpha = 1, labSize = 0, gridlines.minor = FALSE) #plot with labels of cholesterol-related DEGs EnhancedVolcano(Ast_DEGs_MAST, lab = rownames(Ast_DEGs_MAST), x = 'avg_logFC', y = 'p_val_adj_BH', #xlim = c(-0.4, 0.4), #ylim = c(0,6), xlab = bquote('avg log fold change'), ylab = bquote('-log10(p adj) '), axisLabSize = 28, title = 'astrocytes all samples', selectLab = c("FDFT1", "FDPS", "HMGCR", "HMGCS1", "IDI1", "MSMO1", "SQLE", "INSIG1", "ACAT2", "LDLR", "APOE"), pCutoff = 0.05, FCcutoff = 0.05, pointSize = 5, colCustom = keyvals, colAlpha = 1, labSize = 8, boxedLabels = TRUE, labFace = 'bold', drawConnectors = TRUE, widthConnectors = 1.0, gridlines.minor = FALSE)