###homozygous mutant: make heatmaps for genes that are deregulated in isogenic pairs from the following categories: #cholesterol-related genes and fatty acid related genes. ###add new metadata column with ctrl/mut info "ctrl" "het" "hom" #add ctrl/mut CellsMeta <- Astsaver.seurat[[]] CellsMeta["Genotype3"] <- CellsMeta["Genotype"] CellsMeta[CellsMeta$Genotype3 == "337VM",]$Genotype3 <- "het" CellsMeta[CellsMeta$Genotype3 == "406RW",]$Genotype3 <- "het" CellsMeta[CellsMeta$Genotype3 == "406WW",]$Genotype3 <- "hom" CellsMeta[CellsMeta$Genotype3 == "wildtype",]$Genotype3 <- "ctrl" #randomly sample new matadata to verify samplenames are correct and repeat until you are convinced library(dplyr) sample_n(CellsMeta, 1) #now extract SampleName column and add to seurat object CellsMetaTrim <- subset(CellsMeta, select = c("Genotype3")) Astsaver.seurat <- AddMetaData(Astsaver.seurat, CellsMetaTrim) ###make seurat object with good only samples, including the homozygous mutant Astsaver.goodonly.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 Astsaver.goodonly.seurat <- SetIdent(Astsaver.goodonly.seurat, value = "Age") Astsaver.goodonly.seurat <- subset(Astsaver.goodonly.seurat, idents = c("1m", "2m", "3m"), invert = TRUE) Astsaver.goodonly.seurat <- SetIdent(Astsaver.goodonly.seurat, value = "SampleName2") Astsaver.goodonly.seurat <- subset(Astsaver.goodonly.seurat, idents = c("4mFB_337VV3B12_12", "6mFB_337VV3B12_14"), invert = TRUE) Astsaver.goodonly.seurat <- subset(Astsaver.goodonly.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) #then keep only concordant pairs Astsaver.goodonly.seurat <- subset(Astsaver.goodonly.seurat, idents = c("4mFB_337VM1_05", "4mFB_337VV1_05", "4mFB_337VM2_06", "4mFB_337VV2_06", "6mFB_337VM1_07", "6mFB_337VV1_07", "6mFB_337VM3_14", "6mFB_337VV3F02_14", "6mFB_406RW1_06", "6mFB_406RR1_06", "8mFB_406RW1_10", "8mFB_406RR1_10", "6mFB_406WW1_14", "4mFB_406WW1_12")) ####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(Astsaver.goodonly.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_337VV2_06: remove 18 cells sampled.cells_4mFB_337VV2_06 <- subset(Astsaver.goodonly.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(Astsaver.goodonly.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_337VV3F02_14: remove 11 cells sampled.cells_6mFB_337VV3F02_14 <- subset(Astsaver.goodonly.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(Astsaver.goodonly.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) ###8mFB_406RR1_10: remove 91 cells sampled.cells_8mFB_406RR1_10 <- subset(Astsaver.goodonly.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_337VV2_06, sampled.cells_6mFB_337VV1_07, sampled.cells_6mFB_337VV3F02_14, sampled.cells_6mFB_406RW1_06, sampled.cells_8mFB_406RR1_10) Astsaver.goodonly.seurat <- subset(Astsaver.goodonly.seurat, cells = cells_to_remove, invert = TRUE) Astsaver.goodonly.seurat <- SetIdent(Astsaver.goodonly.seurat, value = "Genotype2") #clean up environment rm(list=ls(pattern="^sampled.cells")) rm(cells_to_remove) ##prepare and plot heatmapls library(gplots) library(pheatmap) ############################## ### cholesterol-related DEGs #create a vector of genes of interest chol_genes <- c("CYP51A1", "DHCR7", "DHCR24", "FDFT1", "FDPS", "HMGCR", "HMGCS1", "HSD17B7", "IDI1", "LSS", "MVD", "MSMO1", "SC5D", "SQLE", "TM7SF2", "INSIG1", "SREBF2", "ACAT2", "STARD4", "LDLR", "APOE") #create an object only with the genes from your list chol_genes <- subset(Astsaver.goodonly.seurat, features = chol_genes) #calculate average expression (average will be calculated according to the identity of your object) chol_genes <- SetIdent(chol_genes, value = "Genotype3") expCholGenes <- AverageExpression(chol_genes)$RNA #scaling the rows cal_z_score <- function(x){(x - mean(x)) / sd(x)} expCholGenes <- t(apply( expCholGenes , 1, cal_z_score)) #re-order rows and columns expCholGenes <- expCholGenes[c("DHCR7", "DHCR24", "FDFT1", "FDPS", "HMGCR", "HMGCS1", "HSD17B7", "IDI1", "LSS", "MVD", "MSMO1", "SC5D", "SQLE", "INSIG1", "SREBF2", "ACAT2", "STARD4", "LDLR"),c(2,1,3)] paletteLength <- 50 myColor <- colorRampPalette(c("blue", "white", "red"))(paletteLength) myBreaks <- c(seq(min(expCholGenes), 0, length.out=ceiling(paletteLength/2) + 1), seq(max(expCholGenes)/paletteLength, max(expCholGenes), length.out=floor(paletteLength/2))) pheatmap(expCholGenes, color = myColor, breaks = myBreaks, cluster_rows = FALSE, cluster_cols = FALSE, show_rownames = TRUE, show_colnames = TRUE, fontsize_col=10) ############################## ### fatty acid biosynthesis genes fa_genes <- c("FASN","INSIG1","FADS1","MGST1","PLP1","SCD","FADS2","ELOVL6") ## Create an object only with the genes from your list fa_genes <- subset(Astsaver.goodonly.seurat, features = fa_genes) fa_genes <- SetIdent(fa_genes, value = "Genotype3") expFaGenes <- AverageExpression(fa_genes)$RNA ### Scaling the rows expFaGenes <- t(apply( expFaGenes , 1, cal_z_score)) #re-order the columns expFaGenes <- expFaGenes[,c(2,1,3)] paletteLength <- 50 myColor <- colorRampPalette(c("blue", "white", "red"))(paletteLength) myBreaks <- c(seq(min(expFaGenes), 0, length.out=ceiling(paletteLength/2) + 1), seq(max(expFaGenes)/paletteLength, max(expFaGenes), length.out=floor(paletteLength/2))) pheatmap(expFaGenes, color = myColor, breaks = myBreaks, cluster_rows = FALSE, cluster_cols = FALSE, show_rownames = TRUE, show_colnames = TRUE, fontsize_col=10)