library(DESeq2)
library(org.Mm.eg.db)
library(tidyverse)

# Read the data
count_table <- read_csv("https://denvirlab.marshall.edu/PHS675-2026/data/M_hf_20wk_th_vs_b6_counts.csv")

# Sample information is encoded in the sample names, which are
# column names in the count table
sample_data <- count_table %>%
  select(!EnsemblID) %>%
  colnames() %>%
  # Convert the column names to a single-column table:
  as_tibble_col(column_name = "Sample") %>%
  # Split the names into separate columns 
  separate(Sample, into=c("Sex", "Strain", "Diet", "Age", "Id", "Tissue"), remove = FALSE, sep = '-') %>%
  # Make the strain a factor, as it's our experimental variable
  mutate(Strain = factor(Strain, levels=c("B6", "TH")))

# Create the DESeqDataSet, and then run DESeq on it:
dds <- DESeqDataSetFromMatrix(count_table %>% as.data.frame(), colData = sample_data, design=~Strain, tidy=T) %>%
  DESeq

# Make the MA plot:
png("ma_plot.png", width=960, height=960, res=144)
plotMA(dds)
dev.off()

# Perform a variance-stabilizing transform on the data
vsd <- vst(dds)

# Create a PCA plot, using the variance-stabilizing transformed data
png("pca_plot.png", width=960, height=960, res=144)
plotPCA(vsd, intgroup="Strain")
dev.off()

# Extract the results and convert to a tidyverse table
deseq_results <- results(dds) %>% as_tibble(rownames="EnsemblID")

# Significantly upregulated genes
deseq_results %>% filter(padj < 0.1 & log2FoldChange > 0)

# Significantly downregulated genes
deseq_results %>% filter(padj < 0.1 & log2FoldChange < 0)

# Genes with more-than-minimal expression
expressed_results <- deseq_results  %>% filter(!is.na(padj))

# Volcano plot:

expressed_results %>% 
  # Add a column called "Regulation" for coloring points
  # "Upregulated" are padj < 0.1 and fold change > 2 (log2FoldChange > 1)
  # "Downregulated" are padj < 0.1 and fold change < 0.5 (log2FoldChange < -1) 
  mutate(Regulation = case_when(
    padj < 0.1 & log2FoldChange > 1 ~ "Up",
    padj < 0.1 & log2FoldChange < -1 ~ "Down",
    TRUE ~ "NS"
  )) %>%
  # Plot with log2FoldChange on x-axis and -log(padj) on y-axis
  ggplot(aes(log2FoldChange, -log10(padj))) +
  # Small points colored by regulation:
  geom_point(aes(color=Regulation), size = 0.5) + 
  # Manually define color scheme:
  scale_colour_manual(values=c("Up"="red", "Down"="blue", "NS"="black")) + 
  # Refine labels for axes and title:
  labs(
    x="Log fold change",
    y="Negative log adjusted p-value",
    title="Expression profiling of TH and B6 mice fed a high-fat diet")

# Get gene symbols from AnnotationDbi:
gene_map <- mapIds(org.Mm.eg.db, keys=deseq_results$EnsemblID, keytype="ENSEMBL", column = "SYMBOL", multiVals = "first") %>%
  enframe(name = "EnsemblID", value="GeneName")

# Get the top 20 differentially expressed genes:
# Add the gene symbols to the expressed results table:
top20 <- expressed_results %>% inner_join(gene_map) %>%
  # Only consider genes with a name:
  filter(! is.na(GeneName)) %>%
  # Only consider statistically significant genes at padj < 0.1:
  filter(padj < 0.1) %>%
  # Get the top 20 genes by absolute value of log fold change:
  slice_max(abs(log2FoldChange), n=20) %>%
  # Order by log fold change:
  arrange(log2FoldChange)

# Output to CSV file for creating a table:
top20 %>% select(GeneName, log2FoldChange, padj) %>%
  write_csv("Top20SignificantByFoldChange.csv")

# Create heatmap
png("heatmap_top20.png", width=960, height=960, res=144)

# Take the top 20 genes by fold change and add in the variance-stabilized expression:
top20 %>% inner_join(assay(vsd) %>% as_tibble(rownames="EnsemblID")) %>%
  # Extract just the gene name column and the VST expression columns:
  select(GeneName, starts_with("M-")) %>%
  # Make the gene name column a factor and set the ordering (levels) to the current order:
  mutate(GeneName = factor(GeneName, levels = GeneName)) %>%
  # Pivot the expression columns
  # This will replace all the sample columns with two columns, one for sample name
  # and one for expression.
  # Each row will be replaced by multiple rows (one per sample)
  pivot_longer(! GeneName, names_to = "Sample", values_to = "Expression") %>%
  # To calculate the z-score we need to do this by gene, so group the table by gene:
  group_by(GeneName) %>%
  # Compute the z-score for the expression (expression-mean_expression)/sd_expression:
  mutate(Expression_zscore = (Expression - mean(Expression))/sd(Expression)) %>%
  # Manipulate the sample name so it just has strain and id:
  separate(Sample, into=c("Sex", "Strain", "Diet", "Age", "ID", "Tissue"), sep='-') %>%
  mutate(Sample = paste(Strain, ID, sep='-')) %>%
  # Finally, create the plot: the x-axis (columns) are samples,
  # the y-axis (rows) are genes, and the colors are expression z-scores:
  ggplot(aes(Sample, GeneName, fill=Expression_zscore)) + geom_tile() +
  # Red-Blue diverging palette:
  scale_fill_distiller(palette="RdBu") + 
  # Rotate sample names so they don't collide:
  theme(axis.text.x = element_text(angle=90, hjust=1, vjust=0.5)) +
  # Refine labels
  labs(x="Sample", y="Gene", fill="Expression z-score")
dev.off()

sessionInfo()