library("dplyr")
library("openxlsx")
library("readxl")
library("stringr")
library("ggplot2")

#####Directories#####
input_dir <- "/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_Bresesti/Analysis_MM/GitLab_scripts/Input"
output_dir <- "/beegfs/scratch/ric.squadrito/ric.squadrito/90-935462466_scRNAseq_Bresesti/Analysis_MM/GitLab_scripts/Output"

################################################################################
# R script to generate data and figures for manuscript figures 1B, 1C, and 1D.
#
# This script reads RNA-seq and miRNA-seq data from Excel files, performs
# normalization and data transformation, and generates output tables and plots.
#
# Figures generated/data produced:
#   - Figure 1B: Heatmap data for selected marker genes (Excel table)
#   - Figure 1C.sx: Pairwise volcano plots for RNAseq data
#   - Figure 1C.dx: Pairwise MA plots for miRNA data
#   - Figure 1D: Heatmap data for selected miRNA families (Excel table)
#
# Input data files:
#   - Figure 1B: "miRNA_QIAseq_1509_QIAseqUltraplexRNA_181342.xlsx" (Sheet 3: umis.genes.polyA-mouse)
#   - Figure 1C.sx: "miRNA_QIAseq_1509_181342_edgeR_results.xlsx"
#   - Figure 1C.dx & 1D: "miRNA_QIAseq_1510_173308.all_samples.summary.xlsx" (Sheet 2: miRNA_piRNA)
#   - Figure 1D: "miRNA_Family.xlsx" (Sheet 1)
#
################################################################################


#### Figure 1B ####
df1 <- as.data.frame(read_excel(paste0(input_dir, "/miRNA_QIAseq_1509_QIAseqUltraplexRNA_181342.xlsx"), sheet=3, col_names = T, skip=1)) # sheet: umis.genes.polyA-mouse
df1 <-df1[,-c(1,3:6)] # keep only gene name and UMI counts
df1[,-1] <-apply(df1[,-1],2,function(x){x/sum(x)*1000000}) # UMI normalized by CPM
gene_vector <- c("Cd19", "Ms4a1", "Fcer2a", "Ighm", "Cd8a", "Xcr1", "Itgae", "Itgax", 
                 "Ccr2", "Itgam", "Mgl2", "Cd68", "Vcam1", "Csf1r", "Adgre1", 
                 "Siglec1", "Hmox1", "Timd4", "Vsig4", "Clec4f", "Marco", "Pecam1", 
                 "Tek", "Lyve1", "Stab2")
df2 <- df1[df1$gene%in%gene_vector,] %>% arrange(factor(gene, levels=gene_vector))
rownames(df2) <- df2[,1]
df2.scaled <- as.data.frame(t(scale(t(df2[-1])))) # Zscore normalization. Scaled only works on columns, so need to transform
openxlsx::write.xlsx(df2.scaled, paste0(output_dir, "/Figure_1B_table.xlsx"), rowNames=T)
# the df was exported and used to create a heatmap using Graphpad


#### Figure 1C.sx ####

dge_res <- list()
ff <- paste0(input_dir, "/miRNA_QIAseq_1509_181342_edgeR_results.xlsx")

for (contr in getSheetNames(ff)) {
  dge_res[[contr]] <- read.xlsx(ff, rowNames = T, sheet = contr)
}

sample_order <- c("RPM", "cDC1", "cDC2", "B220+", "LSEC", "KC")

# Plot DGE volcano
diff_genes_union <- c()
diff_genes_inter <- c()
pva_thr <- 0.001
logfc_thr <- 2
for (contr in names(dge_res)) {
  gg <- row.names(subset(dge_res[[contr]], PValue < pva_thr & abs(logFC) > logfc_thr))
  diff_genes_union <- union(diff_genes_union, gg)
  diff_genes_inter <- intersect(diff_genes_inter, gg)
}
length(diff_genes_union)
length(diff_genes_inter)

full.t <- data.frame()
for (contr in names(dge_res)) {
  df <- dge_res[[contr]]
  df$Sample1 <- str_split_fixed(contr, "-", 2)[,1]
  df$Sample2 <- str_split_fixed(contr, "-", 2)[,2]
  full.t <- rbind(full.t, df)
}
head(full.t)

# Volcano plot
data <- full.t
data$pval <- -log10(data$PValue)
data <- mutate(data, color = case_when(data$logFC > logfc_thr & data$pval > -log10(pva_thr) ~ "Overexpressed",
                                       data$logFC < -logfc_thr & data$pval > -log10(pva_thr) ~ "Underexpressed",
                                       abs(data$logFC) < logfc_thr | data$pval < -log10(pva_thr) ~ "NonSignificant"))
data$Sample1 <- factor(data$Sample1, levels = sample_order)
data$Sample2 <- factor(data$Sample2, levels = sample_order)
num_genes <- data.frame()
for (contr in names(dge_res)) {
  s1 <- str_split_fixed(contr, "-", 2)[,1]
  s2 <- str_split_fixed(contr, "-", 2)[,2]
  dd <- subset(data, Sample1 == s1 & Sample2 == s2 & PValue < pva_thr)
  n_up <- nrow(subset(dd, logFC > logfc_thr))
  n_down <- nrow(subset(dd, logFC < -logfc_thr))
  print(nrow(dd))
  df <- data.frame("Sample1" = s1,
                   "Sample2" = s2,
                   "label" = paste(paste0("Down: ", n_down), paste0("Up: ", n_up), sep = "               "))
  num_genes <- rbind(num_genes, df)
}
num_genes$Sample1 <- factor(num_genes$Sample1, levels = sample_order)
num_genes$Sample2 <- factor(num_genes$Sample2, levels = sample_order)

vol <- ggplot(data, aes(x = logFC, y = pval, colour = color)) +
  theme_bw(base_size = 12) +
  theme(legend.position = "none") +
  xlab("log2 Fold Change") +
  ggtitle(paste0("DGE results: pval < ", pva_thr, " & abs(logFC) > ", logfc_thr)) +
  ylab(expression(-log[10]("adjusted p-value"))) +
  geom_hline(yintercept = -log10(pva_thr), colour = "darkgray") +
  geom_vline(xintercept = logfc_thr, colour = "darkgray") +
  geom_vline(xintercept = -logfc_thr, colour = "darkgray") +
  geom_point(size = 2, alpha = 0.8, na.rm = T) +
  scale_color_manual(name = "Expression",
                     values = c(Overexpressed = "red", # #CD4F39
                                Underexpressed = "blue", # #008B00
                                NonSignificant = "darkgray")) +
  geom_text(data = num_genes, mapping = aes(x = -1, y = -3, label = label, color = "black")) +
  facet_grid(Sample1 ~ Sample2, scales = "free")
ggsave(filename = paste(output_dir, "edgeR_DGE_res_volcano_Fig.3C.sx.pdf", sep = "/"), plot = vol, width = 15, height = 16)
# "edgeR_DGE_res_volcano_Fig.3C.sx.pdf" was imported in illustrator and merged with the plot generated by the code below

#### Figure 1C.dx ####
df1 <- as.data.frame(read_excel(paste0(input_dir, "/miRNA_QIAseq_1510_173308.all_samples.summary.xlsx"), sheet=2, col_names=T)) # sheet: miRNA_piRNA
df1 <- df1[,1:7] # UMI 
rownames(df1)<-df1$miRNA
df2 <- apply(df1[,-1], 2, function(x) log(x)) # Log counts
df2 <- as.data.frame(limma::normalizeCyclicLoess(df2, weights = NULL, span=0.7, iterations = 5, method = "pairs")) # Cyclic loess normalization
colnames(df2)<-c("B cell","RPM","cDC1","cDC2","LSEC","KC") # Rename columns
upper.panel<-function(x, y, ...){
  points(((x+y)/2),(x-y), cex=0.6, pch=19, col="grey")
  above <- (x-y-1)*((x+y)/2-5) > 5 & ((x+y)/2)>5
  points(((x+y)/2)[above],(x-y)[above], col="red", cex=0.6, pch=19)
  below <- (x-y+1)*((x+y)/2-5) < -5 & ((x+y)/2)>5
  points(((x+y)/2)[below],(x-y)[below], col="blue", cex=0.6, pch=19)
} # function for MA plot
pdf(file = paste(output_dir, "MA_plots_miRNA_Fig_1C.dx.pdf", sep = "/"), width=9, height=9)
pairs(df2[,1:6], lower.panel = NULL, upper.panel = upper.panel, 
      ylim=c(-8.5,8.5), xlim=c(2,14), cex.labels = 2) # pairwise plot
dev.off()

#### Figure 1D ####
df1 <- as.data.frame(read_excel(paste0(input_dir, "/miRNA_QIAseq_1510_173308.all_samples.summary.xlsx"), sheet=2, col_names=T)) # sheet: miRNA_piRNA
df1 <- df1[-which(grepl("piR",df1[,1])),1:7] # UMI and piRNA removing
df1[,1] <- gsub("/.*","",df1[,1]) # leave only first miRNA for ambiguous entries
df1[,-1] <- apply(df1[,-1],2,function(x){x/sum(x, na.rm=T)*1000000}) # UMI normalized by CPM
df.families <- as.data.frame(read_excel(paste0(input_dir, "/miRNA_Family.xlsx"), sheet=1 ,col_names=T)) # miRNA families from miRBase
df.families <- df.families[df.families[,3]==10090,c(4,1)] # select mouse entries
df.families[,1] <- sub(pattern = "p.*", replacement="p", x = df.families[,1])
df2 <- merge(df.families,df1,by=1)
df2 <- aggregate(df2[,-c(1:2)], by = df2["miR family"], FUN=sum, na.rm=T) # sum counts by family
gene_vector <- c("miR-150-5p","miR-25-3p/32-5p/92-3p/363-3p/367-3p","miR-142-3p.1",
                 "miR-17-5p/20-5p/93-5p/106-5p","miR-191-5p",
                 "miR-15-5p/16-5p/195-5p/322-5p/497-5p","miR-26-5p","miR-138-5p",
                 "miR-223-3p","miR-342-3p","miR-22-3p","miR-192-5p/215-5p",
                 "miR-125-5p/351-5p","miR-126-3p.1","miR-199-3p")
df3 <- subset(df2, df2$`miR family`%in% gene_vector) %>% arrange(factor(`miR family`, levels=gene_vector))
rownames(df3) <- df3[,1]
df3.scaled <- as.data.frame(t(scale(t(df3[-1])))) # Zscore normalization. Scaled only works on columns, so need to transform
openxlsx::write.xlsx(df2.scaled, paste0(output_dir, "/Figure_1D_table.xlsx"), rowNames=T)
# The df was exported and used to create a heatmap using Graphpad