#Load the required library

library('Seurat')
## Loading required package: SeuratObject
## Loading required package: sp
## 
## Attaching package: 'SeuratObject'
## The following object is masked from 'package:base':
## 
##     intersect
library('ggplot2')
library('plyr')
library('dplyr')
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:plyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library('viridis')
## Loading required package: viridisLite
library('scCustomize')
## scCustomize v2.0.1
## If you find the scCustomize useful please cite.
## See 'samuel-marsh.github.io/scCustomize/articles/FAQ.html' for citation info.
library('ggthemes')
library('plotly')
## 
## Attaching package: 'plotly'
## The following objects are masked from 'package:plyr':
## 
##     arrange, mutate, rename, summarise
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout

#load the rds files for Figure Generation

LMN_Final = readRDS("LMN_Final_Remake_Jan8.rds")
LMN_CST_InjVUninj <- FindMarkers(LMN_Final,
                                 subset.ident = "CST",
                                 group.by = "injury.status",
                                 ident.1 = "injured",
                                 ident.2 = "uninjured",
                                 logfc.threshold = 0.2,
                                 min.pct = 0.01)
## For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
## (default method for FindMarkers) please install the presto package
## --------------------------------------------
## install.packages('devtools')
## devtools::install_github('immunogenomics/presto')
## --------------------------------------------
## After installation of presto, Seurat will automatically use the more 
## efficient implementation (no further action necessary).
## This message will be shown once per session
write.csv(x = LMN_CST_InjVUninj, 
          file = "LMN_CST_InjvUninj.csv",
          quote = FALSE)


LMN_HYP_InjVUninj <- FindMarkers(LMN_Final,
                                 subset.ident = "HYP",
                                 group.by = "injury.status",
                                 ident.1 = "injured",
                                 ident.2 = "uninjured",
                                 logfc.threshold = 0.2,
                                 min.pct = 0.01)


write.csv(x = LMN_HYP_InjVUninj, 
          file = "LMN_HYP_InjvUninj.csv",
          quote = FALSE)

LMN_RN_InjVUninj <- FindMarkers(LMN_Final,
                                 subset.ident = "RN",
                                 group.by = "injury.status",
                                 ident.1 = "injured",
                                 ident.2 = "uninjured",
                                 logfc.threshold = 0.2,
                                 min.pct = 0.01)


write.csv(x = LMN_RN_InjVUninj, 
          file = "LMN_RN_InjvUninj.csv",
          quote = FALSE)


LMN_DP_InjVUninj <- FindMarkers(LMN_Final,
                                subset.ident = "DP",
                                group.by = "injury.status",
                                ident.1 = "injured",
                                ident.2 = "uninjured",
                                logfc.threshold = 0.2,
                                min.pct = 0.01)


write.csv(x = LMN_DP_InjVUninj, 
          file = "LMN_DP_InjvUninj.csv",
          quote = FALSE)

LMN_MED_InjVUninj <- FindMarkers(LMN_Final,
                                subset.ident = "MED",
                                group.by = "injury.status",
                                ident.1 = "injured",
                                ident.2 = "uninjured",
                                logfc.threshold = 0.2,
                                min.pct = 0.01)


write.csv(x = LMN_MED_InjVUninj, 
          file = "LMN_MED_InjvUninj.csv",
          quote = FALSE)

LMN_ALL_InjVUninj <- FindMarkers(LMN_Final,
                                 group.by = "injury.status",
                                 ident.1 = "injured",
                                 ident.2 = "uninjured",
                                 logfc.threshold = 0.2,
                                 min.pct = 0.01)


write.csv(x = LMN_ALL_InjVUninj, 
          file = "LMN_ALL_InjvUninj.csv",
          quote = FALSE)

Idents(LMN_Final) <- "seurat_clusters"
levels(LMN_Final)
##  [1] "0"  "1"  "4"  "11" "16" "5"  "20" "18" "15" "2"  "17" "19" "12" "3"  "8" 
## [16] "6"  "14" "7"  "9"  "10" "13"
LMN_Final_1vs0 <- FindMarkers(LMN_Final,
                                 ident.1 = "1",
                                 ident.2 = "0",
                                 logfc.threshold = 0.2,
                                 min.pct = 0.01)


write.csv(x = LMN_Final_1vs0, 
          file = "LMN_Final_1vs0.csv",
          quote = FALSE)
Idents(LMN_Final) <- "manual.clusters.1"
levels(LMN_Final)
##  [1] "CST"           "PVH"           "LH"            "RN"           
##  [5] "EW"            "PAG"           "LC"            "DP"           
##  [9] "Raphe / Inhib" "Inhib"         "HB-1-Daam2"    "HB-2-Pard3b"  
## [13] "HB-3-Col19a1"  "HB-4-Nox4"     "HB-5-Pappa"    "HB-6-Other"
DimPlot_scCustom(seurat_object = LMN_Final, figure_plot = TRUE, label = F)
## Warning: Removed 1 rows containing missing values (`geom_point()`).

ggsave("Figure_1_Cluster.png",height = 4,width = 6,dpi = 300,units = "in")
## Warning: Removed 1 rows containing missing values (`geom_point()`).

#Figure_1_B

Idents(LMN_Final) <- "seurat_clusters"

DimPlot(LMN_Final, reduction = "umap", label = FALSE, pt.size = 1, 
        label.size = 0, group.by =  "injury.status", cols = c("red","green"),alpha = 0.3, split.by = "injury.status",
        order = FALSE
) + NoLegend() + NoAxes() 

ggsave("Figure_1_B.png",height = 4,width = 8,dpi = 300,units = "in")
DimPlot(LMN_Final, reduction = "umap", label = FALSE, pt.size = 1, 
        label.size = 0, group.by =  "injury.status", cols = c("red","green"),alpha = 0.3,
        order = FALSE
) + NoLegend() + NoAxes() 

ggsave("Figure_1_B_2.png",height = 4,width = 4,dpi = 300,units = "in")

#For Figure 1C - Compare Creb5 between injured and uninjured

Idents(LMN_Final) <- "manual.clusters.1"
# Subsetting for each cluster
CST <- subset(LMN_Final, idents = "CST")
DP <- subset(LMN_Final, idents = "DP")
RN <- subset(LMN_Final, idents = "RN")
HB <- subset(LMN_Final, idents = c("HB-1-Daam2","HB-2-Pard3b","HB-3-Col19a1","HB-4-Nox4","HB-5-Pappa","HB-6-Other"))
PVH <- subset(LMN_Final, idents = "PVH")
LH <- subset(LMN_Final, idents = "LH")
PAG <- subset(LMN_Final, idents = "PAG")
LC <- subset(LMN_Final, idents = "LC")
Raphe_Inhib <- subset(LMN_Final, idents = "Raphe / Inhib")
pal <- viridis(n = 10, option = "D")
FeaturePlot_scCustom(seurat_object = LMN_Final,pt.size = 1, features = "Creb5", colors_use = pal,split.by = "injury.status", order = T) & NoAxes() + theme(axis.text = element_text(size = 20),text = element_text (size = 20))
## 
## NOTE: FeaturePlot_scCustom uses a specified `na_cutoff` when plotting to
## color cells with no expression as background color separate from color scale.
## Please ensure `na_cutoff` value is appropriate for feature being plotted.
## Default setting is appropriate for use when plotting from 'RNA' assay.
## When `na_cutoff` not appropriate (e.g., module scores) set to NULL to
## plot all cells in gradient color palette.
## 
## -----This message will be shown once per session.-----

ggsave("Figure_1_C_full.png",height = 4,width = 8,dpi = 300,units = "in")
#Making plot for clusters seprately 
FeaturePlot_scCustom(seurat_object = CST,pt.size = 1, features = "Creb5", colors_use = pal,split.by = "injury.status", order = T) & NoAxes() & NoLegend()+ theme(axis.text = element_text(size = 20),text = element_text (size = 20))

ggsave("Figure_1_C_CST.png",height = 4,width = 8,dpi = 300,units = "in")
FeaturePlot_scCustom(seurat_object = DP,pt.size = 1, features = "Creb5", colors_use = pal,split.by = "injury.status", order = T) & NoAxes() & NoLegend()+ theme(axis.text = element_text(size = 20),text = element_text (size = 20))

ggsave("Figure_1_C_DP.png",height = 4,width = 8,dpi = 300,units = "in")
FeaturePlot_scCustom(seurat_object = HB,pt.size = 1, features = "Creb5", colors_use = pal,split.by = "injury.status", order = T) & NoAxes() & NoLegend()+ theme(axis.text = element_text(size = 20),text = element_text (size = 20))

ggsave("Figure_1_C_HB.png",height = 4,width = 8,dpi = 300,units = "in")
FeaturePlot_scCustom(seurat_object = RN,pt.size = 1, features = "Creb5", colors_use = pal,split.by = "injury.status", order = T) & NoAxes() & NoLegend()+ theme(axis.text = element_text(size = 20),text = element_text (size = 20))

ggsave("Figure_1_C_RN.png",height = 4,width = 8,dpi = 300,units = "in")
LMN_Final_DotPlot2 <- read.csv("LMN_Final_Dotplot2.csv", sep = ",")
DotPlot_scCustom (LMN_Final, features = LMN_Final_DotPlot2$Marker, dot.scale = 6, dot.min = .01,colors_use = pal) +
  coord_flip() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1),
        axis.text = element_text(size = 15))
## Warning: Removed 444 rows containing missing values (`geom_point()`).

ggsave("Figure_1_Dotplor.png",height = 14,width = 9,dpi = 300,units = "in",bg = "white")
## Warning: Removed 444 rows containing missing values (`geom_point()`).

#Volcano plot

deg_data <- read.csv("LMN_ALL_InjvUninj.csv")
Volcano_filtered_data <- deg_data %>%
  filter(!grepl("Rik", X), !grepl("^Gm", X),!grepl("Grik", X))
p = ggplot(Volcano_filtered_data, aes(x=avg_log2FC, y=-log10(p_val_adj), color=ifelse(avg_log2FC > 2, "green", ifelse(avg_log2FC < -2, "red", "grey"))))+
  geom_point(aes(text=paste("Gene_name:", X, "<br>FC:", avg_log2FC, "<br>P-value:", -log10(p_val_adj))), alpha=0.5) +
  scale_color_identity() +
  geom_hline(yintercept=-log10(0.05), color="red", linetype="dashed") +
  geom_vline(xintercept=0, color="grey", linetype="dashed") +
  geom_vline(xintercept=0.5, color="blue", linetype="dashed") +
  geom_vline(xintercept=-0.5, color="blue", linetype="dashed") +
  labs(title="",
       x="Average log2 Fold Change", y="Adjusted P-value (log scale)") +
  theme_few() +
  scale_x_continuous(limits=c(-6, 6), breaks=seq(-6, 6, by=1))
## Warning in geom_point(aes(text = paste("Gene_name:", X, "<br>FC:", avg_log2FC,
## : Ignoring unknown aesthetics: text
p

p_interactive=ggplotly(p)
htmlwidgets::saveWidget(p_interactive, file="VolcanoPlot_figure_1.html")
#MA Plot
# Add a new column for color based on avg_log2FC
average_expression <- rowMeans(Volcano_filtered_data[,c("pct.1", "pct.2")])
Volcano_filtered_data$color <- ifelse(Volcano_filtered_data$avg_log2FC > 0.5, "green",
                         ifelse(Volcano_filtered_data$avg_log2FC < -0.5, "red", "black"))
# Create an MA plot
ma_plot <- ggplot(Volcano_filtered_data, aes(x=average_expression, y=avg_log2FC, color=color)) +
  geom_point(alpha=0.3) +
  scale_color_identity() +  # This will use the colors specified in the 'color' column
  geom_hline( yintercept =  0, color="gray", linetype="dashed") +
  geom_hline(yintercept = 0.5, color="blue", linetype="dashed") +
  geom_hline(yintercept = -0.5, color="blue", linetype="dashed") +
  xlab("Average Expression") + ylab("Log2 Fold Change") +
  #ggtitle("MA Plot") +
  theme_few()  + scale_y_continuous(limits=c(-6, 6), breaks=seq(-6, 6, by=0.5))
# Print the plot
print(ma_plot)

#Histogram

library(ggplot2)
library(cowplot)
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggthemes':
## 
##     theme_map
deg_data <- Volcano_filtered_data
# Define common aesthetics with continuous fill based on avg_log2FC
common_aes <- aes(x=avg_log2FC, fill=avg_log2FC)

# Define the color scale
color_scale <- scale_fill_gradient2(low = "red", mid = "black", high = "green", 
                                    midpoint = 0, limit = c(-5, 5),
                                    breaks = c(-0.5, 0.5))

# Define common theme
common_theme <- theme_few() + theme(text = element_text(size = 16)) +
                theme(axis.text = element_text(size=12),
                      axis.title = element_text(size=14, face="bold"),
                      plot.title = element_text(size=16, hjust = 0.5))

# Adjust the size of the labels here
label_size <- 4.5  # Increase this value to make the labels larger

# Lower plot (zoomed-in version)
p1 <- ggplot(deg_data, common_aes) + 
  geom_histogram(binwidth=0.5, color="white", alpha=0.2, show.legend=FALSE) +
  color_scale +
  geom_vline(xintercept = 0, color="black", linetype="dashed", linewidth=0.7) +
  geom_vline(xintercept = c(-0.5, 0.5), color="grey", linetype="dashed", linewidth=0.7, alpha=0.5) +
  geom_text(stat='bin', aes(label=ifelse(..count.. >= 1, ..count.., "")), vjust=-0.8, size=label_size, binwidth=0.5, color="black", 
            label.padding = unit(0.2, "lines"), label.size = 0.25, label.r = unit(0.15, "lines"), label.color = NA, fill="white") +
  xlab("Log2 Fold Change") + ylab("Number of Genes") +
  scale_x_continuous(breaks=seq(-5, 5, by=0.5), limits=c(-5, 5)) + 
  coord_cartesian(ylim = c(0, 1000)) + 
  common_theme
## Warning in geom_text(stat = "bin", aes(label = ifelse(..count.. >= 1, ..count.., : Ignoring unknown parameters: `label.padding`, `label.size`, `label.r`,
## `label.colour`, and `fill`
# Upper plot (for the range with fewer genes)
p2 <- ggplot(deg_data, common_aes) + 
  geom_histogram(binwidth=0.5, color="white", alpha=0.2, show.legend=FALSE) +
  color_scale +
  geom_vline(xintercept = 0, color="black", linetype="dashed", linewidth=0.7) +
  geom_vline(xintercept = c(-0.5, 0.5), color="grey", linetype="dashed", linewidth=0.7, alpha=0.5) +
  geom_text(stat='bin', aes(label=ifelse(..count.. >= 1000, ..count.., "")), vjust=-1.5, size=label_size, binwidth=0.5, color="black", 
            label.padding = unit(0.2, "lines"), label.size = 0.25, label.r = unit(0.15, "lines"), label.color = NA, fill="white") +
  coord_cartesian(ylim = c(2000, 20000)) + 
  common_theme +
  xlab("Log2 Fold Change") + ylab("Number of Genes") + 
  scale_x_continuous(breaks=seq(-5, 5, by=0.5), limits=c(-5, 5)) + 
  scale_y_continuous(breaks = seq(0,20000,by=2000)) +
  annotate("text", x=-1.5, y=14000, label="Downregulated", size=4, color="red", angle= 0) +
  annotate("text", x=1.5, y=14000, label="Upregulated", size=4, color="green", angle = 0) +
  annotate("text", x=0, y=18000, label="No Change", size=4, color="blue", angle = 0)
## Warning in geom_text(stat = "bin", aes(label = ifelse(..count.. >= 1000, : Ignoring unknown parameters: `label.padding`, `label.size`, `label.r`,
## `label.colour`, and `fill`
p3 <- p2 + 
  theme(axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank()) 
p4 = p1 +
  theme(axis.title.y = element_blank())

# Combine the two plots with p1 at the bottom
combined_plot <- plot_grid(p3, p4, ncol=1, align='v', rel_heights = c(2/3, 1/3))
## Warning: The dot-dot notation (`..count..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(count)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The following aesthetics were dropped during statistical transformation: fill
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?
## The following aesthetics were dropped during statistical transformation: fill
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?
## Warning: Removed 2 rows containing missing values (`geom_bar()`).
## Warning: The following aesthetics were dropped during statistical transformation: fill
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?
## The following aesthetics were dropped during statistical transformation: fill
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?
## Warning: Removed 2 rows containing missing values (`geom_bar()`).
print(combined_plot)

#heatmap

top_50_genes <- head(deg_data[order(deg_data$p_val_adj),], 50)

heatmap_plot <- ggplot(top_50_genes, aes(x = reorder(top_50_genes[[1]], avg_log2FC), y = 1, fill = avg_log2FC)) + 
  geom_tile(color = "white") +  # Adding white border to the tiles
  geom_text(aes(label = sprintf("%.2f", avg_log2FC)), vjust = 0.5, hjust = 0.5, color = "black", angle = 90) +
  scale_fill_gradient2(low = "#b35806", mid = "white", high = "#542788", midpoint = 0, limits = c(-4, 4)) +
  labs(x = "", y = "", fill = "Log2FC") +
  theme_few() + theme(text = element_text(size = 14)) +
  theme(axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_text(angle = 90, vjust = 0.5),
        axis.title.x = element_text(angle = 90, vjust = 0.5))
heatmap_plot
## Warning: Use of `top_50_genes[[1]]` is discouraged.
## ℹ Use `.data[[1]]` instead.
## Use of `top_50_genes[[1]]` is discouraged.
## ℹ Use `.data[[1]]` instead.