---
title: "Figure Panels Script"
author: "Dana Schwalbe"
date: "10/16/2023"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# load packages:
```{r}
library(Seurat)
library(dplyr)
library(ggplot2)
library(pheatmap)
library(clustifyr)
```


#Load data:
```{r}

    vlmneurons = readRDS("~/VLM_Neurons_filtered.rds")
    serotonergic_neurons = readRDS("~/n07n08_subclusters.rds")
    cartptcluster = readRDS("~/n09_subclusters.rds")
    vlm_integrated_filtered = readRDS("~/VLM_integrated_filtered.rds")
    astrocytes_integrated = readRDS("~/VLM_astrocytes.rds")

    
```


## Figure 1:
```{r}
# Panel H:
    #UMAP with annotations
    pdf(file = "/path/to/save/file.pdf", width = 6, height = 6)
        
    DimPlot(vlm_integrated_filtered, label=T, raster=F,pt.size=0.001,
            label.size =3)+NoLegend()+
      theme(axis.line=element_blank(),axis.ticks = element_blank(),
            axis.title = element_blank(),axis.text = element_blank())
        
        dev.off()
        
        
# Panel I:
        
#Heatmap of DEGs
  # load saved DEGs .csv file:
    allcell_markers = read.csv("/path/to/file.csv")
  # or run following line:
    #allcell_markers = FindAllMarkers(vlm_integrated_filtered,only.pos=T)

 #check max p-value
    max(allcell_markers$p_val_adj)
    #max pval is below 0.01 --> no need to filter.
    
    heatmap_allcellmarkersgenes = allcell_markers$gene
    
    cluster.averages <- AverageExpression(vlm_integrated_filtered, assays = "RNA",
                          features = heatmap_allcellmarkersgenes,group.by="celltypes")
    
    ##### Save averaged expression to csv
    write.csv(cluster.averages,file="/path/to/save/cluster_averages.csv")
    
    ##### Read in expression values as data matrix
    read.csv(file="/path/to/save/cluster_averages.csv",header=TRUE,row.names=1) -> expr_mat
    expr_mat <- data.matrix(expr_mat)

    
    ## Save pdf
    pdf(file = "/path/to/save/file.pdf", width = 4, height = 10)
    
    pheatmap(expr_mat,border_color="NA", scale="row",color=cm.colors(256), 
         show_rownames = F,
         cluster_rows=F,cluster_cols = F)
    
    dev.off()
    
    
# Panel J:
    
    #Cluster tree
        vlm_integrated_filtered = BuildClusterTree(vlm_integrated_filtered,dims=1:30,reorder = T)
        
        pdf(file = "/path/to/save/file.pdf", width = 4, height = 4)
        
        par(cex = .001)
        PlotClusterTree(vlm_integrated_filtered,cex=1000)
        
        dev.off()
        
        
    #N cells in each cell type
        nCells = as.data.frame(table(Idents(vlm_integrated_filtered)))
    
        pdf(file = "/path/to/save/file.pdf", width = 7, height = 4)
        
        ggplot(nCells) + 
          geom_bar(aes(x=Var1, y=Freq,fill=Var1),stat="identity")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
          theme_bw()+labs(y="Number of Cells",x="CellType")+
          geom_text(aes(x=Var1, y=Freq,label = Freq), vjust = -0.2,size=2.5)+
          theme(panel.grid.major = element_blank(),panel.grid.minor = element_blank())
        
        dev.off()
        
    #nUMIs
        pdf(file = "/path/to/save/file.pdf", width = 4, height = 4)
        
        VlnPlot(vlm_integrated_filtered, "nCount_RNA",pt.size=0)+NoLegend()+
          theme(axis.text=element_text(size=5),title = element_text(size=10),axis.title = element_text(size=7))
        
        dev.off() 
        
    #nGenes
        pdf(file = "/path/to/save/file.pdf", width = 4, height = 4)
        
        VlnPlot(vlm_integrated_filtered, "nFeature_RNA",pt.size=0)+NoLegend()+ theme(axis.text=element_text(size=5),title = element_text(size=10),axis.title = element_text(size=7))
        
        dev.off()
        
        
    #Expression of cell type marker genes
        genes = c("Olig1","Olig2","Cspg4","Agt","Cx3cr1","Ranbp3l","Slc47a1","Slc6a13",
                "Kcnj8","Abcc9","Slco1c1","Foxp2","Rbfox3","Syp","Syn1","Kl","Folr1")
        
        pdf(file = "/path/to/save/file.pdf",  width = 5, height = 7) 
        
        VlnPlot(vlm_integrated_filtered,features=genes,stack=T,flip=F)+NoLegend()
        
        dev.off()
        
```



## Figure 2:
```{r}

# Panel A:
    #UMAP with annotations
        pdf(file = "/path/to/save/file.pdf",
            width = 5, 
            height = 5) 
        
        DimPlot(vlmneurons, label=T, raster=FALSE,pt.size=0.3,label.size=6)+NoLegend()
        
        dev.off()


#Panel B:
    #Cluster tree
        vlmneurons = BuildClusterTree(vlmneurons,dims=1:20)
        
        pdf(file = "/path/to/save/file.pdf", width = 4, height = 4)
        
        par(cex = .001)
        PlotClusterTree(vlmneurons,cex=1000)
        
        dev.off()
        
        
    
    #N cells in each subtype
        nCells = as.data.frame(table(vlmneurons$annotationsnumeric))
        
        pdf(file = "/path/to/save/file.pdf", width = 7, height = 4)
        
        ggplot(nCells) + 
          geom_bar(aes(x=Var1, y=Freq,fill=Var1),stat="identity")+
          theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
        labs(y="Number of Cells",x="CellType")+ geom_text(aes(x=Var1, y=Freq,label = Freq), 
                                                          vjust = -0.2,size=2.5)
        
        dev.off()
    
    #nUMIs
        pdf(file = "/path/to/save/file.pdf", width = 4, height = 4)
        
        VlnPlot(vlmneurons, "nCount_RNA",pt.size=0)+NoLegend()+
          theme(axis.text=element_text(size=5),title = element_text(size=10),
                axis.title = element_text(size=7))
        
        dev.off()
        
        
    #nGenes
         pdf(file = "/path/to/save/file.pdf",
             width = 4, height = 4)
        
        VlnPlot(vlmneurons, "nFeature_RNA",pt.size=0)+NoLegend()+
          theme(axis.text=element_text(size=5),title = element_text(size=10),
                axis.title = element_text(size=7))
        
        dev.off()
        

# Panel C:
    #Dot plot of cluster marker genes
        color=cm.colors(256)
        
        pdf(file = "/path/to/save/file.pdf",
            width = 10, 
            height = 4) 
        neuronmarkers = c("Car8","Calb1","Arhgef33","Pcp2","Selenos","Etv1","Gabra6","Fat2","Il16","Cbln3","Ntn1","Kit","Megf10","Chst9","Cnpy1","Pdzrn3","Satb2","Mlip","Npas2","Ankrd33b","Gpc3","A630012P03Rik","Mmrn1","Popdc3","Trpc4","Hydin","Gm29536","Crhbp","Ptprq","Gm31135","Trh","Cpne7","Tph2","1700042O10Rik","Slc6a4","Ddc","Cpa6","Gm43154","Maoa","Cartpt","Dbh","Th","Atp6v1f","Tuba4a","Tuba1b","Mif","Gli3","Gm12128","Uts2b","Onecut1","Glra3","Kcnh8","Gpr149","Pax8","Pax2","Esr1","Ebf3","Ebf2","Prdm6","Crybg3","Angpt1","Tmem132cos","Gxylt2","Gm29683","Nxph2","Lypd6b","Ttn","Rbm20","Maf","E130114P18Rik","Necab1","Dach2","4930438E09Rik","Mecom","Gm12649","Shox2","Epb41l4a","4930511M06Rik","C1ql2","Prph","Chat","Slc5a7","Gm10754","Creb5","Col27a1","Cdh23","Lama1","Cpne9","Atp1a2","Slc1a3","Plpp3","Aldh1a1","Gpr37l1")
          
          DotPlot(vlmneurons, features=unique(neuronmarkers),dot.scale=4,cols=c(color[1],color[256]))+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1,size=7),axis.text.y = element_text(size=7),axis.title = element_text(size=8),legend.text = element_text(size=6),legend.title = element_text(size=6))
          
          dev.off()
          
          
          
# Panel D:
    multi_gene_list = c("Slc6a5","Gad2","Slc32a1","Gad1","Slc17a7","Slc17a6","Dbh","Slc6a2","Tph2","Slc6a4","Fev","Chat","Slc18a3")
    
    cluster.averages_multi_gene_list <- AverageExpression(vlmneurons, assays = "RNA",
                          features = multi_gene_list,group.by="annotationsnumeric")
    
    ##### Save averaged expression to csv
    write.csv(cluster.averages_multi_gene_list,file="/path/to/save/cluster_averages_multigenelist_neurons.csv")
    
    ##### Read in expression values as data matrix
    read.csv(file="/path/to/save/cluster_averages_multigenelist_neurons.csv",header=TRUE,row.names=1) -> expr_mat
    expr_mat <- data.matrix(expr_mat)
    
    colnames(expr_mat) = c("n01","n02","n03","n04","n05","n06","n07","n08","n09","n10","n11","n12","n13","n14","n15","n16","n17","n18","n19","n20","n21","n22","n23")
  
      pdf(file = "/path/to/save/file.pdf",width = 5, height = 10) 
   
       pheatmap(expr_mat,border_color="NA", scale="row",color=cm.colors(256), 
         show_rownames = T,
         cluster_rows=T,cluster_cols = T)
       
    dev.off()

```


# Figure 3:
```{r}

#Panel A:
    Idents(vlmneurons)=vlmneurons$ID
    cells_SpinalTRAP = subset(vlmneurons, idents="RVLM,SpinalTRAP+")
    Idents(vlmneurons)=vlmneurons$annotationsnumeric
    
    pdf(file = "/path/to/save/file.pdf",
        width = 5, 
        height = 5) 
    
    DimPlot(vlmneurons, label=T,raster=FALSE,pt.size=0.3,label.size=4,
            cells.highlight = colnames(cells_SpinalTRAP),cols.highlight = "#E76BF3")+NoLegend()+
      theme(axis.line=element_blank(),axis.ticks = element_blank(),
            axis.title = element_blank(),axis.text = element_blank())
    
      dev.off()
     
       
# Panel B:
    pdf(file = "/path/to/save/file.pdf",width = 7, height = 5)  
       
    cluster = vlmneurons$annotationsnumeric
    ID = vlmneurons$ID
    data <- data.frame(cluster,ID)
    # Stacked
    ggplot(data) +
      geom_bar(aes(fill=ID, x=cluster),position="stack")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
    
    dev.off()
    
    
# Panel C:
     Idents(vlmneurons)=vlmneurons$ID
     cells_SpinalTRAP = subset(vlmneurons, idents="RVLM,SpinalTRAP+")
     Idents(vlmneurons)=vlmneurons$annotationsnumeric
     
      data_new =(table(cells_SpinalTRAP$annotationsnumeric)/1012)*100
      
      data_new <- data.frame(
  group=names(data_new),
  value=as.numeric(data_new)
)
      
   pdf(file = "/path/to/save/file.pdf",width = 7, height = 5)  

      
      ggplot(data_new, aes(group,value)) + geom_col()+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+ylab("percent")


        dev.off()
        
# Panel E:
    fig3e_markers = c("C1ql2","Lhx4","Lhx3","Vsx2","Shox2")
        
    pdf(file = "/path/to/save/file.pdf", width = 4, height = 10)
    
    DotPlot(vlmneurons, features=unique(fig3e_markers),
            dot.scale=4.5,cols = c(color[1],color[256]))+
          theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1,size=6),
                axis.text.y = element_text(size=7),axis.title = element_text(size=8),
                legend.text = element_text(size=6),legend.title = element_text(size=6))
    
        dev.off()

```


# Figure 4:
```{r}
# Panel A
    # violin plot:

    genesf4 = c("Tacr1","Sst","Oprm1","Penk","Cdh9","Slc17a6","Slc6a5", "Slc32a1", "Gad2", "Gpc3", "Glra3","Trpc4")
    
    pdf(file = "/path/to/save/file.pdf",   
        width = 5,
        height = 7)
    
    VlnPlot(vlmneurons,features=genesf4,stack=T,flip=F)+NoLegend()
    
  
    dev.off()
```


# Figure 5:
```{r}
# Panel A:

    genesf5a = rev(c("Slc6a4","Tph2","Trh","Cpa6","Cpne7","Gm43154"))
  
    pdf(file = "/path/to/save/file.pdf",
        width = 4, 
        height = 9) 
       
    VlnPlot(vlmneurons, features=genesf5a,stack=T)+NoLegend()

    dev.off()
    
    
# Panel D:    
    #UMAP, subclusters with annotations
        pdf(file = "/path/to/save/file.pdf",
            width = 4, 
            height = 4)
    
        DimPlot(serotonergic_neurons, label=T, raster=FALSE,pt.size=0.3,label.size=4)+NoLegend()+theme(axis.line=element_blank(),axis.ticks = element_blank(),axis.title = element_blank(),axis.text = element_blank())
        
        dev.off()
        
#Panel E:
        
    #UMAP, SpinalTRAP+ samples highlighted
       
      Idents(serotonergic_neurons)=serotonergic_neurons$ID
        cellssero_SpinalTRAP = subset(serotonergic_neurons, idents="RVLM,SpinalTRAP+")
        Idents(serotonergic_neurons)=serotonergic_neurons$subclusters_numeric
        
    
     pdf(file = "/path/to/save/file.pdf",
            width = 4, 
            height = 4)
    
        DimPlot(serotonergic_neurons, label=F, raster=FALSE,pt.size=0.3,
                cells.highlight = colnames(cellssero_SpinalTRAP), cols.highlight = "#E76BF3")+
          NoLegend()+theme(axis.line=element_blank(),axis.ticks = element_blank(),
                           axis.title = element_blank(),axis.text = element_blank())
        
        dev.off()
        
        # pie chart values:
        (table(cellssero_SpinalTRAP$subclusters_numeric)/364)*100     
        
# Panel F:
              
      sero_markers = c("Rps5","Rpl32","Rps16","Rps29","Sncg","Thsd7b","Ldb2",
                       "Csgalnact1","Zfp536","Il1rapl2","Megf11","Ryr3","Kcnt2",
                       "Zeb2","Hs3st4","Tmem132c","Tmem132cos","Pou6f2","Tafa2",
                       "Adgrl2","Gm30382","Pcdh11x","Gm2516")
          
      pdf(file = "/path/to/save/file.pdf", width = 10, height = 4)
      
      DotPlot(serotonergic_neurons, features=unique(sero_markers),
              dot.scale=4.5,cols = c(color[1],color[256]))+
            theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1,size=6),
                  axis.text.y = element_text(size=7),axis.title = element_text(size=8),
                  legend.text = element_text(size=6),legend.title = element_text(size=6))
      
          dev.off()
        
#Panel G:
        sero_genes_of_interst = c("Trh","Cpne7","Cpa6","Gm43154","Tac1","Slc17a8","Hcrtr1"
                 ,"Calcr","Adra1a","Gad1","Gad2")
          
      pdf(file = "/path/to/save/file.pdf", width = 6, height = 4)
      
          VlnPlot(serotonergic_neurons, features=sero_genes_of_interst,stack=T)+NoLegend()
      
      
          dev.off()
        

```


# Figure 6:
```{r}

# Panel A:
        pdf(file = "/path/to/save/file.pdf",
            width = 4, 
            height = 9) 
           
        VlnPlot(vlmneurons, features=c("Cartpt","Th","Dbh"),stack=T)+NoLegend()
    
        dev.off()
        
        
# Panel B
        
      #Original UMAP, zoomed in on cluster 9 --> feature plots of Pnmt, Slc6a2, SpinalTRAP+
      
      cluster9 = subset(vlmneurons, idents = "n09")
      
      pdf(file = "/path/to/save/file.pdf",
              width = 4,
              height = 4) 
      
      FeaturePlot(cluster9,"Pnmt", order=T)
      
          dev.off()
          
        pdf(file = "/path/to/save/file.pdf",
              width = 4,
              height = 4) 
      
      FeaturePlot(cluster9,"Slc6a2", order=T)
      
          dev.off()
          
          
          
      Idents(cluster9) = cluster9$ID
      spinal_trap_cluster9 = subset(cluster9, idents="RVLM,SpinalTRAP+")
      Idents(cluster9) = cluster9$annotationsnumeric
      
      pdf(file = "/path/to/save/file.pdf",
              width = 4,
              height = 4) 
      
      
         DimPlot(cluster9, label=F,raster=FALSE,pt.size=0.3,
                  cells.highlight = colnames(spinal_trap_cluster9),cols.highlight = "#E76BF3")+NoLegend()+
            theme(axis.line=element_blank(),axis.ticks = element_blank(),
                  axis.title = element_blank(),axis.text = element_blank())
          
          dev.off()
          
          
# Panel C:
      #UMAP with annotations
    pdf(file = "/path/to/save/file.pdf",   
        width = 4, 
        height = 4)
    
    DimPlot(cartptcluster, label=T, raster=FALSE,pt.size=0.3,label.size=4)+NoLegend()+theme(axis.line=element_blank(),axis.ticks = element_blank(),axis.title = element_blank(),axis.text = element_blank())
    
    dev.off()
    

# Panel D:
    
      Idents(cartptcluster) = cluster9$ID
      cartpt_C1pos = subset(cartptcluster, idents="RVLM,C1pos")
      Idents(cartptcluster) = cartptcluster$subclusters_numeric
      
      pdf(file = "/path/to/save/file.pdf",
              width = 4,
              height = 4) 
      
      
         DimPlot(cartptcluster, label=F,raster=FALSE,pt.size=0.3,
                  cells.highlight = colnames(cartpt_C1pos),cols.highlight = "#00BF7D")+NoLegend()+
            theme(axis.line=element_blank(),axis.ticks = element_blank(),
                  axis.title = element_blank(),axis.text = element_blank())
          
          dev.off()
          
          
      # pie chart values:
      (table(cartpt_C1pos$subclusters_numeric)/259)*100     


# Panel E:
          
      Idents(cartptcluster) = cartptcluster$ID
      spinaltrappos_cartpt_subclusters = subset(cartptcluster, idents="RVLM,SpinalTRAP+")
      Idents(cartptcluster) = cartptcluster$subclusters_numeric
      
      pdf(file = "/path/to/save/file.pdf",
              width = 4,
              height = 4) 
      
      
         DimPlot(cartptcluster, label=F,raster=FALSE,pt.size=0.3,
                  cells.highlight = colnames(spinaltrappos_cartpt_subclusters),cols.highlight = "#E76BF3")+NoLegend()+
            theme(axis.line=element_blank(),axis.ticks = element_blank(),
                  axis.title = element_blank(),axis.text = element_blank())
          
          dev.off()
          
      #pie chart values:
      (table(spinaltrappos_cartpt_subclusters$subclusters_numeric)/105)*100     
          

#Panel F:
          
      pdf(file = "/path/to/save/file.pdf", width = 6, height = 4)

      VlnPlot(cartptcluster, features=c("Cartpt","Th","Dbh","Pnmt",
                                        "Phox2b","Slc17a6","Slc6a2",
                                        "Slc18a2","Penk","Npy","Adcyap1",
                                        "Tacr1"),stack=T)+
        NoLegend()+ theme(axis.text=element_text(size=7),
                    title = element_text(size=10),axis.title = element_text(size=9))

      dev.off()
      
          
        cartpt_markers = c("9630002D21Rik","Chrm3","A330102I10Rik","Gm13832","Chat","Galr1","Gm44618","Lama1","Fbln5","Htr2c","Gabrb2","Gm20713","4930419G24Rik","A230004M16Rik","Htr2a","Dscaml1","Aoah","Gm13944","Pstpip1","Arhgap15","Ryr3","Gm12128","Myl1","Lypd1","St8sia6","Eya1","Calcr","Tshz1","Sema3d","1700016P03Rik","Gm12027","Htr4","Slc5a7","Pde7b","Prkg2","Pdk4","9530026P05Rik","Hs3st2","Col8a1","Nell1os","Lncbate1","Nell1","Npy")
    
      pdf(file = "/path/to/save/file.pdf", width = 6, height = 4)
    
      DotPlot(cartptcluster, features=unique(cartpt_markers),
            dot.scale=4.5,cols=c(color[1],color[256]))+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1,size=6)
                               ,axis.text.y = element_text(size=7),axis.title = element_text(size=8),
                               legend.text = element_text(size=6),legend.title = element_text(size=6))
    
        dev.off()
                            
          
```


# Figure 7:
```{r}
#Panel A:
    
    pdf(file = "/path/to/save/file.pdf",
        width = 9, 
        height = 4) 
       
    VlnPlot(vlmneurons, features=c("Eya1"))+NoLegend()

    dev.off()
    
# Panel B:
    
        pdf(file = "/path/to/save/file.pdf",
            width = 4,
            height = 4) 
    
    FeaturePlot(cartptcluster,"Eya1", order=T, label=T)
    
        
        dev.off()
  
# Panel c:
      pdf(file = "/path/to/save/file.pdf",
        width = 16,
        height = 4) 

FeaturePlot(cartptcluster,c("Eya1","Slc6a2"),blend=T, order = T)
    
    dev.off()

```


# Figure 8:
```{r}

#Panel A:
        pdf(file = "/path/to/save/file.pdf",
            width = 9, 
            height = 4) 
           
        VlnPlot(vlmneurons, features=c("Hk2"))+NoLegend()
    
        dev.off()
        
        
#Panel B:
    pdf(file = "/path/to/save/file.pdf",
        width = 4,
        height = 4) 

    FeaturePlot(cartptcluster,"Hk2", order=T, label=T)

    
    dev.off()
```


# Suppplemental Figure 2:
```{r}

#Panel A:
    #UMAP by batch
    pdf(file = "/path/to/save/file.pdf",   
        width = 4, 
        height = 4) 
    
    DimPlot(vlm_integrated_filtered, label=F,shuffle=T, raster=FALSE,pt.size=0.01,
          group.by = "batch")+theme(axis.line=element_blank(),axis.ticks = element_blank(),axis.title = element_blank(),axis.text = element_blank())
    
    dev.off()
    
#Panel B:
  #Batch counts in data
  cluster = vlm_integrated_filtered$celltypes
  batch = vlm_integrated_filtered$batch
  data <- data.frame(cluster,batch)
    # Stacked
    pdf(file = "/path/to/save/file.pdf", width = 4, height = 4)
    
   ggplot(data) + 
    geom_bar(aes(fill=batch, x=cluster),position="stack")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+theme(legend.text = element_text(size=6),legend.title = element_text(size=7))
    
    dev.off()
    
#Panel C:
    #Correlation matrix:
    matrix_rvlm = vlm_integrated_filtered[["RNA"]]@data
    
    meta_rvlm = as.data.frame(vlm_integrated_filtered$celltypes)
    colnames(meta_rvlm)=c("celltypes")
    
    vlm_integrated_filtered = FindVariableFeatures(vlm_integrated_filtered)
    
    vargenes = VariableFeatures(vlm_integrated_filtered)

     new_ref_matrix <- average_clusters(
      mat = matrix_rvlm,
      metadata = meta_rvlm, 
      cluster_col = "celltypes",
      if_log = TRUE                    
    )
    
    res <- clustify(
      input = matrix_rvlm,
      metadata = meta_rvlm, 
      cluster_col = "celltypes", 
      ref_mat = new_ref_matrix, 
      query_genes = vargenes 
    )
    
    pdf(file = "/path/to/save/file.pdf",width = 4, height = 4)

    color=cm.colors(256)

    plot_cor_heatmap(cor_mat = res,col=c("white",color[1],color[256]))

    dev.off()
    

# Panel E:    
    pdf(file = "/path/to/save/file.pdf", width = 4, height = 4)
    
        DimPlot(astrocytes_integrated, label=T, raster=FALSE,pt.size=0.3,label.size=4)+NoLegend()+theme(axis.line=element_blank(),axis.ticks = element_blank(),axis.title = element_blank(),axis.text = element_blank())
    
    dev.off()

    color=cm.colors(256)

    DefaultAssay(astrocytes_integrated)="RNA"

    pdf(file = "/path/to/save/file.pdf", width = 4, height = 4)

  DotPlot(astrocytes_integrated, features=c("Agt","Gfap","Slc1a3","Rbfox1","Prkg1","Mybpc1","Ust","Cdh4","Nkain2","Slc24a2","Tmeff2","Pex5l","St18","Ablim2","6330411D24Rik","Myoc","Ccdc148","Gm11099","Gm4951","Gm50237","Oasl2","Iigp1","Gbp6"),dot.scale=4,cols = c(color[1],color[256]))+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1,size=7),axis.text.y = element_text(size=7),axis.title = element_text(size=8),legend.text = element_text(size=6),legend.title = element_text(size=6))
  
    dev.off()   
    

# Panel F:
    cluster = vlmneurons$annotationsnumeric
    batch = vlmneurons$batch
    data <- data.frame(cluster,batch)
    
    pdf("/path/to/save/file.pdf",
        width = 4,
        height = 4) 
    
    ggplot(data) + 
        geom_bar(aes(fill=batch, x=cluster),position="stack")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+theme(legend.text = element_text(size=6),legend.title = element_text(size=7))
    
    dev.off()
    
# Panel G:
    pdf(file = "/path/to/save/file.pdf", width = 6, height = 4)
    
    cluster = serotonergic_neurons$subclusters_numeric
    batch = serotonergic_neurons$batch
    data <- data.frame(cluster,batch)
    
    ggplot(data) + 
        geom_bar(aes(fill=batch, x=cluster),position="stack")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+theme(legend.text = element_text(size=6),legend.title = element_text(size=7))
    
        dev.off()
        
# Panel H:
    pdf(file = "/path/to/save/file.pdf", width = 6, height = 4)
    
    cluster = cartptcluster$subclusters_numeric
    batch = cartptcluster$batch
    data <- data.frame(cluster,batch)
    ggplot(data) + 
        geom_bar(aes(fill=batch, x=cluster),position="stack")+theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+theme(legend.text = element_text(size=6),legend.title = element_text(size=7))

    dev.off() 
    
    
    
```