Bioinformatics Final Project
Mia Kovacs
2025-04-25
Introduction
As people age, it is common to see signs of decline across multiple dimensions, such as cognitive and motor function and an overall vulnerability to stress. This multifaceted age related decline is classified as frailty syndrome, which is characterized by weakness, frequent falls, cognitive decline, increased hospitalization and death, and increased vulnerability to stressors due to a loss of biological reserves (Dacomo et al., 2024). Although it is not a disease itself, it is linked to a variety of serious diseases including hypertension, congestive heart failure, diabetes mellitus, stroke, chronic lung disease, and cancer, as well as neurological disorders such dementia, Alzheimer’s disease, and brain atrophy (Dacomo et al., 2024). Fortunately, frailty is reversible, if treated early (Dacomo et al., 2024). Given both the potential for remediation and the severe diseases linked to frailty, identifying biomarkers associated with frailty may help improve lifespan and quality of life for aging populations through early identification and precision medicine. The study Transcriptomic and metabolomic changes might predict frailty in SAMP8 mice aims to do this by investigating potential frailty biomarkers using mice models.
Given the close relationship between frailty and age related conditions, I hypothesized that there would be increased expression of Alzheimer’s related genes in the differentially expressed genes in SAMP8 mice. Within the original study, mice at the 6 month time point had the largest quantity of differentially expressed genes, so for this bioinformatic analysis, only mice from this time point were chosen. To further investigate this potential relationship, first analysis of differentially expressed genes was conducted, including DESeq2 analysis, the creation of PCA plots, volcano plots, and heatmaps. Following this, Alzheimer’s Disease related genes were analyzed also using heatmaps and volcano plots, as well as gene count plots. Additionally, figures elucidating potential shortcomings of the original study’s bioinformatic pipeline were created to further explain certain results.
Load Libraries
# loading required packages
pacman::p_load(DESeq2, dplyr, tidyr, ggplot2, genekitr, pheatmap,
RColorBrewer, EnhancedVolcano, clusterProfiler, DOSE,
ggnewscale, enrichplot, knitr, cowplot, msigdbr, pathview,
patchwork, ggplotify)Data Wrangling
Loading countdata
Given roadblocks encountered within the initial bioinformatic pipeline, countdata for this project was obtained directly from the original paper. This data was adjusted to only include data from 6 month old mice. Additionally, a seperate meta data table was created for annotation purposes.
# Read in countdata
countdata <- read.csv("SAMP8_SAMR1_six_months.csv",
header=TRUE,
row.names=1)
# Transform countdata to a data matrix
countdata <- as.matrix(countdata)# Read in meta data file that contains information on treatment
coldata <- read.csv(file = "meta_SAMP8_SAMR1_six_months.csv",
header = T,
row.names = 1)
rownames(coldata)## [1] "SAMP8_cortex_six_months1" "SAMP8_cortex_six_months2"
## [3] "SAMP8_cortex_six_months3" "SAMP8_cortex_six_months4"
## [5] "SAMP8_hippocampus_six_months1" "SAMP8_hippocampus_six_months2"
## [7] "SAMP8_hippocampus_six_months3" "SAMR1_cortex_six_months1"
## [9] "SAMR1_cortex_six_months2" "SAMR1_cortex_six_months3"
## [11] "SAMR1_hippocampus_six_months1" "SAMR1_hippocampus_six_months2"
## [13] "SAMR1_hippocampus_six_months3" "SAMR1_hippocampus_six_months4"# Ensure that column names and row names match up
all(colnames(countdata) == rownames(coldata))## [1] TRUEDifferential Expression Analysis
Create DESeq Object
To continue in this analysis, first a DESeq object of differentially expressed genes needs to be created using DESeq2. Treatment within this study includes the mouse model, age, and tissue sample (either cortex or hippocampus).
# Create dds object from countdata and coldata
dds <- DESeqDataSetFromMatrix(countData = countdata,
colData = coldata,
design = ~ replicates + treatment)Pre-filtering
After creating the initial DESeq object, the data is filtered to remove genes with low numbers of reads.
# Filter out rows with less than 20 total reads
keep <- rowSums(counts(dds)) >= 20
dds <- dds[keep,]Run Differential expression analysis
Following the creation of the DESeq object and filtering of the data, the differential expression analysis can be run.
dds <- DESeq(dds)
#saveRDS(dds, file = "dds.RDS")
#dds <- DESeq(dds, test="Wald", sfType="poscount")Pairwise comparison
After running the differential expression analysis, pairwise comparison can be conducted to learn more regarding the data.
# Want to compare SAMP8 to SAMR1 cortexes, as the model used is like the treatment. Will do this again for hippocampus
SAMP8_SAMR1_CTX <- results(dds, contrast = c("treatment", "SAMP8_cortex", "SAMR1_cortex"))
SAMP8_SAMR1_HIP <- results(dds, contrast = c("treatment", "SAMP8_hippocampus", "SAMR1_hippocampus"))Data visualization
PCA plot of top differentially expressed genes
PCA plots can be a very useful tool in analyzing data, taking large data sets and present findings in a clear visual manner. With PCA plots, clustering can give insight into overall patterns in data, as well as potential batch effects. Within the original paper, one of the transcriptomic figures created was a PCA plot of the top differentially expressed genes. This figure is a recreation of this, using the top 60 differentially expressed genes.
# Create PCA plot of normalized counts
# Creat pcaData object, include strain and tissue
vstcounts <- vst(dds, blind=TRUE)
pcaData <- plotPCA(vstcounts,
intgroup = c("strain", "tissue"),
ntop = 60,
returnData = TRUE)
# returns a numeric vector that is used in axis labeling
percentVar <- round(100 * attr(pcaData, "percentVar"))
# Creat blank graph
ggplot(pcaData, aes(x = PC1, y = PC2, color = strain, shape = tissue)) +
# Add geoms
geom_point(size = 4, alpha = 0.8) +
# Add labels
xlab(paste0("PC1: ", percentVar[1], "% variance")) +
ylab(paste0("PC2: ", percentVar[2], "% variance")) +
labs(title = "PCA of Top 60 Differentially Expressed Genes") +
# Clean up
theme_grey() +
theme(
legend.position = "right",
plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5)
) +
scale_color_brewer(palette = "Set1")
Figure 1. PCA Plot of Top 60 DE Genes PCA plot of the top 60 differentially expressed genes within SAMP8 and SAMR1 mouse cortexes and hippocampi. There is a 49% PC1 variance and 21% PC2 variance.
PCA plot of all differentially expressed genes
Selecting only the top differentially expressed genes in a data set is likely to force clustering to appear anyways, so a more comprehensive view including all genes should be created to investigate if there are more broad patterns within gene expression.
# Create PCA plot of normalized counts
# Creat pcaData object, include strain and tissue
pcaData <- plotPCA(vstcounts,
intgroup = c("strain", "tissue"),
returnData = TRUE)
# returns a numeric vector that is used in axis labeling
percentVar <- round(100 * attr(pcaData, "percentVar"))
# Creat blank graph
ggplot(pcaData, aes(x = PC1, y = PC2, color = strain, shape = tissue)) +
# Add geoms
geom_point(size = 4, alpha = 0.8) +
# Add labels
xlab(paste0("PC1: ", percentVar[1], "% variance")) +
ylab(paste0("PC2: ", percentVar[2], "% variance")) +
labs(title = "PCA of VST-normalized Counts") +
# Clean up
theme_grey() +
theme(
legend.position = "right",
plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5)
) +
scale_color_brewer(palette = "Set1")
Figure 2. PCA Plot of All DE Genes PCA plot of all differentially expressed genes within SAMP8 and SAMR1 mouse cortexes and hippocampi. There is a 25% PC1 variance and 23% PC2 variance.
Plot of number of reads per sample
Within the PCA plots there was somewhat weak clustering, as well as multiple outliers. To investigate potential sources of this, a plot of reads per sample can provide insight into potential errors in earlier pipelines.
# Shorten names
colnames(countdata) <- gsub("SAMP8_cortex_six_months", "SAMP8_cortex", colnames(countdata))
colnames(countdata) <- gsub("SAMP8_hippocampus_six_months", "SAMP8_hippocampus", colnames(countdata))
colnames(countdata) <- gsub("SAMR1_cortex_six_months", "SAMR1_cortex", colnames(countdata))
colnames(countdata) <- gsub("SAMR1_hippocampus_six_months", "SAMR1_hippocampus", colnames(countdata))
# Create data frame to make ggplot graph
bar_dataframe <- data.frame(Sample = colnames(countdata),
Total_Counts = colSums(countdata),
Group = factor(c(rep("SAMP_cortex", 4), rep("SAMP8_hippocampus", 3),
rep("SAMR1_cortex", 3), rep("SAMR1_hippocampus", 4))))
# Create blank graph
ggplot(data = bar_dataframe,
mapping = aes(x = Sample,
y = Total_Counts,
fill = Group)) +
# Add geoms
geom_bar(stat = "identity",
color = "black",
size = 0.2) +
# Clean up and add labels
theme_grey() +
theme(axis.text.x = element_text(angle = 45,
hjust = 1,
size = 8),
plot.title = element_text(hjust = 0.5, face = "bold")) +
labs(x = "Replicates",
y = "Total Counts",
title = "Number of Reads per Sample") +
# Select colors
scale_fill_manual(values = c("tomato", "lightblue", "forestgreen", "orchid"))
Figure 3. Bar Plot of Number of Reads per Sample There appears to be a high level of variation between the number of reads each sample had, with SAMP8_cortex1 having the largest number of reads at 30 million, and SAMP8_hippocampus3 having the smallest.
Merged Results
Based on the results of the PCA plots, further analysis will be conducted specifically using hippocampus data, as there was more clear clustering among these two groups. Given this, a new data frame was created that included the differentially expressed genes for the hippocampus and new columns including if the genes were unchanged, up, or down regulated were created.
SAMP8_SAMR1_HIP <- merge(as.data.frame(SAMP8_SAMR1_HIP),
as.data.frame(counts(dds, normalized=TRUE)),
by="row.names",
sort=FALSE)
names(SAMP8_SAMR1_HIP)[1] <- "Gene_id"
#Add a column to dataframe to specify up or down regulated
SAMP8_SAMR1_HIP$diffexpressed <- "Unchanged"
#if log2FC > 1 and padj < 0.05 set as upregulated
SAMP8_SAMR1_HIP$diffexpressed[SAMP8_SAMR1_HIP$padj < 0.05 & SAMP8_SAMR1_HIP$log2FoldChange > 1] <- "Upregulated"
#if log2FC > -1 and padj < 0.05 set as downregulated
SAMP8_SAMR1_HIP$diffexpressed[SAMP8_SAMR1_HIP$padj < 0.05 & SAMP8_SAMR1_HIP$log2FoldChange < -1] <- "Downregulated"GeneID conversion
The gene information in its current form is difficult to understand and will not allow for further analysis, so the gene IDs were converted to symbol and entrez IDs.
# Remove the version from Ensembl gene IDs
SAMP8_SAMR1_HIP$Ensembl_id <- sub("\\.\\d+$", "", SAMP8_SAMR1_HIP$Gene_id)
# Extract the cleaned Ensembl IDs
ids <- SAMP8_SAMR1_HIP$Ensembl_id
# Convert to gene symbols and Entrez IDs
gene_symbol_map <- transId(ids,
transTo = c("symbol", "entrez"),
org = "mouse", # switch to mouse
unique = TRUE) # ensures one-to-one mapping
SAMP8_SAMR1_HIP <- merge(SAMP8_SAMR1_HIP,
gene_symbol_map,
by.x = "Ensembl_id",
by.y = "input_id")
str(SAMP8_SAMR1_HIP)## 'data.frame': 18492 obs. of 25 variables:
## $ Ensembl_id : 'AsIs' chr "ENSMUSG00000000001" "ENSMUSG00000000028" "ENSMUSG00000000031" "ENSMUSG00000000037" ...
## $ Gene_id : 'AsIs' chr "ENSMUSG00000000001.5" "ENSMUSG00000000028.16" "ENSMUSG00000000031.17" "ENSMUSG00000000037.18" ...
## $ baseMean : num 56.65 3.41 1.47 3.66 85.43 ...
## $ log2FoldChange : num 0.395 -0.572 3.668 4.304 0.127 ...
## $ lfcSE : num 0.424 1.526 2.392 1.463 0.353 ...
## $ stat : num 0.934 -0.375 1.533 2.941 0.359 ...
## $ pvalue : num 0.35051 0.70791 0.12521 0.00327 0.71925 ...
## $ padj : num 0.908 1 0.692 0.207 1 ...
## $ SAMP8_cortex_six_months1 : num 26.16 6.83 0.61 1.34 123.49 ...
## $ SAMP8_cortex_six_months2 : num 58.132 3.185 0.796 0.796 88.393 ...
## $ SAMP8_cortex_six_months3 : num 38.51 0 1.75 1.75 96.28 ...
## $ SAMP8_cortex_six_months4 : num 30.78 4.29 0 2.86 83.74 ...
## $ SAMP8_hippocampus_six_months1: num 42.84 8.57 17.14 0 94.25 ...
## $ SAMP8_hippocampus_six_months2: num 83.5 0 0 37.1 92.8 ...
## $ SAMP8_hippocampus_six_months3: num 123.6 0 0 0 61.8 ...
## $ SAMR1_cortex_six_months1 : num 42.01 5.09 0 1.91 84.02 ...
## $ SAMR1_cortex_six_months2 : num 53.18 1.22 0 2.44 89.85 ...
## $ SAMR1_cortex_six_months3 : num 54.11 1.59 0 1.59 87.53 ...
## $ SAMR1_hippocampus_six_months1: num 57.4 4.5 0 0 78.8 ...
## $ SAMR1_hippocampus_six_months2: num 64.1 0 0 0 37.4 ...
## $ SAMR1_hippocampus_six_months3: num 77.625 4.168 0.347 1.389 95.512 ...
## $ SAMR1_hippocampus_six_months4: num 41.13 8.23 0 0 82.25 ...
## $ diffexpressed : chr "Unchanged" "Unchanged" "Unchanged" "Unchanged" ...
## $ symbol : chr "Gnai3" "Cdc45" "H19" "Scml2" ...
## $ entrezid : chr "14679" "12544" "14955" "107815" ...#write.csv(SAMP8_SAMR1_HIP,
#"SAMP8_SAMR1_HIP_full_matrix_042525.csv",
#row.names = FALSE)Volcano plot
Volcano plots are useful tools in visualizing up and down regulated genes in samples, as well as identifying specific genes that are more differentially expressed.
# Select genes to label
library(ggrepel)
# Chose 3 related to Alzheimers, 3 upregulated, and 3 downregulated
genes_to_label <- c("Cox5a", "Sdhd", "Ndufa9", "Vip", "Rnase4", "P2rx3", "Il11", "Pr29", "Nopchap1")
# If genes_to_lable in dataframe, add to volcano plot
SAMP8_SAMR1_HIP$label <- ifelse(SAMP8_SAMR1_HIP$symbol %in% genes_to_label,
SAMP8_SAMR1_HIP$symbol,
NA)# Create volcano plot
# Make blank graph
ggplot(SAMP8_SAMR1_HIP) +
# Add geoms
geom_point(aes(x = log2FoldChange, y = -log10(padj),
color = diffexpressed, alpha = 0.8)) +
# Add gene labels
geom_text_repel(aes(x = log2FoldChange,
y = -log10(padj),
label = label),
size = 3.5,
max.overlaps = Inf) +
# Add labels
theme_grey() +
ggtitle("Volcano Plot of Differentially Expressed Genes") +
xlab("Log2 fold change") +
ylab("-Log10(padj)") +
# Clean up
theme(legend.position = "right",
plot.title = element_text(size = rel(1.25), hjust = 0.5),
axis.title = element_text(size = rel(1.25))) +
scale_color_manual(name = "Differential Expression",
values = c("red", "grey", "blue"),
labels = c("Downregulated", "Not Significant", "Upregulated")) +
scale_x_continuous(limits = c(-20, 20)) +
coord_cartesian(ylim = c(0, 15), clip = "off") +
# Remove alpha legend
guides(alpha = "none")
Figure 4. Volcano Plot of DE Genes There were a total of 29 upregulated genes and 55 downregulated genes. Labeled are several genes of interest including two upregulated genes, two downregulated genes, and two Alzheimer’s disease related genes.
Heatmaps
Another useful tool in visualizing broad patterns of expression is heat maps. For this portion of the analysis, three different heatmaps were created. The first replicated the Heatmap present within the original paper, which includes only the top 60 differentially expressed genes. Following this, a second heatmap including all differentially expressed genes was created to see broader patterns present in the data.
# Create matrix to make heatmap of top 60 DE genes
# Convert to data frame
deg_matrix <- data.frame(SAMP8_SAMR1_HIP)
# Select top DE genes
de_top <- deg_matrix[(deg_matrix$baseMean > 50) & (abs(deg_matrix$log2FoldChange) > 1.495),]
# Number of rows in de_top should equal 60, for top 60 DE genes
summary(de_top)## Ensembl_id Gene_id baseMean log2FoldChange
## Length:60 Length:60 Min. : 50.12 Min. :-4.916
## Class :AsIs Class :AsIs 1st Qu.: 62.54 1st Qu.:-2.203
## Mode :character Mode :character Median : 82.38 Median :-1.684
## Mean : 611.24 Mean :-1.337
## 3rd Qu.: 120.63 3rd Qu.:-1.509
## Max. :23537.42 Max. : 2.691
## lfcSE stat pvalue padj
## Min. :0.2672 Min. :-6.096 Min. :0.0000000 Min. :0.0000014
## 1st Qu.:0.5142 1st Qu.:-3.114 1st Qu.:0.0006647 1st Qu.:0.0976711
## Median :0.6554 Median :-2.477 Median :0.0046327 Median :0.2318238
## Mean :0.9205 Mean :-1.399 Mean :0.0404551 Mean :0.3073441
## 3rd Qu.:1.2403 3rd Qu.:-1.308 3rd Qu.:0.0537154 3rd Qu.:0.5477974
## Max. :2.9817 Max. : 4.585 Max. :0.3155471 Max. :0.8844648
## SAMP8_cortex_six_months1 SAMP8_cortex_six_months2 SAMP8_cortex_six_months3
## Min. : 0.73 Min. : 11.15 Min. : 0.00
## 1st Qu.: 34.07 1st Qu.: 56.84 1st Qu.: 38.51
## Median : 57.23 Median : 82.22 Median : 60.39
## Mean : 2870.87 Mean : 307.38 Mean : 214.37
## 3rd Qu.: 93.56 3rd Qu.: 152.20 3rd Qu.: 115.10
## Max. :155721.60 Max. :7509.04 Max. :6595.08
## SAMP8_cortex_six_months4 SAMP8_hippocampus_six_months1
## Min. : 0.00 Min. : 0.00
## 1st Qu.: 44.91 1st Qu.: 20.35
## Median : 75.15 Median : 42.84
## Mean : 240.68 Mean : 101.39
## 3rd Qu.: 113.26 3rd Qu.: 77.11
## Max. :7037.46 Max. :1619.32
## SAMP8_hippocampus_six_months2 SAMP8_hippocampus_six_months3
## Min. : 0.00 Min. : 0.00
## 1st Qu.: 18.56 1st Qu.: 15.45
## Median : 46.39 Median : 30.90
## Mean : 173.13 Mean : 129.14
## 3rd Qu.: 96.27 3rd Qu.: 92.70
## Max. :5349.21 Max. :3198.29
## SAMR1_cortex_six_months1 SAMR1_cortex_six_months2 SAMR1_cortex_six_months3
## Min. : 3.82 Min. : 0.00 Min. : 12.73
## 1st Qu.: 50.29 1st Qu.: 46.61 1st Qu.: 46.55
## Median : 87.52 Median : 70.60 Median : 84.35
## Mean : 1373.13 Mean : 121.40 Mean : 444.79
## 3rd Qu.: 220.40 3rd Qu.: 102.69 3rd Qu.: 206.10
## Max. :72151.13 Max. :1139.95 Max. :12653.04
## SAMR1_hippocampus_six_months1 SAMR1_hippocampus_six_months2
## Min. : 9.00 Min. : 0.00
## 1st Qu.: 37.69 1st Qu.: 36.05
## Median : 59.64 Median : 96.13
## Mean : 903.32 Mean : 188.75
## 3rd Qu.: 87.20 3rd Qu.: 166.90
## Max. :44754.33 Max. :1404.61
## SAMR1_hippocampus_six_months3 SAMR1_hippocampus_six_months4 diffexpressed
## Min. : 5.21 Min. : 0.00 Length:60
## 1st Qu.: 52.10 1st Qu.: 41.12 Class :character
## Median : 73.98 Median : 82.25 Mode :character
## Mean : 873.70 Mean : 615.30
## 3rd Qu.: 252.76 3rd Qu.: 123.38
## Max. :29095.91 Max. :27644.34
## symbol entrezid label
## Length:60 Length:60 Length:60
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
## # Order differentially expressed genes by log fold change
de_top <- de_top[order(de_top$log2FoldChange,
decreasing = TRUE),]
# Select relevant columns
heatmap_matrix <- select(de_top,
symbol,
SAMP8_cortex_six_months1:SAMR1_hippocampus_six_months4)
heatmap_matrix <- data.frame(heatmap_matrix, row.names = 1)# Create annotation file for heatmap
# Read in meta data file
heatmap_meta <- read.csv("meta_SAMP8_SAMR1_six_months.csv",
header = TRUE)
# Assign sample names as row names
rownames(heatmap_meta) <- heatmap_meta$Sample_ID
# Select replicate and treatment
heatmap_meta <- select(heatmap_meta, treatment)
# Re-level factors
heatmap_meta$treatment <- factor(heatmap_meta$treatment,
levels = c("SAMR1_cortex", "SAMR1_hippocampus", "SAMP8_cortex", "SAMP8_hippocampus"))
# Reorder columns to match figure in paper
ordered_samples <- rownames(heatmap_meta)[order(heatmap_meta$treatment)]
heatmap_matrix <- heatmap_matrix[, ordered_samples]
heatmap_meta <- heatmap_meta[ordered_samples, , drop = FALSE]# Pick colors
heatmap_colors <- rev(brewer.pal(10, "RdBu"))
# Plot with pheatmap
pheatmap(heatmap_matrix,
color = heatmap_colors,
annotation_col = heatmap_meta,
show_colnames = FALSE,
show_rownames = FALSE,
cluster_cols = FALSE,
clustering_distance_cols = "euclidean",
clustering_method = "complete",
scale = "row",
main = "Heatmap of Top 60 DE Genes") 
Figure 5. Heatmap of Top 60 DE Genes Up and down regulation patterns vary across each of the groups, with hippocampus samples exhibiting stronger patterns in regulation. SAMP8 hippocampus samples tended to have more upregulated regions when those same regions were down regulated in SAMR1, with the inverse also being true. However, these patterns are still weak.
# Create heatmap of all genes to compare DE expression broadly across mouse models
# Convert to data frame
deg_norm_matrix <- data.frame(SAMP8_SAMR1_HIP)
# Order differentially expressed genes by log fold change
deg_norm_matrix <- deg_norm_matrix[order(deg_norm_matrix$log2FoldChange,
decreasing = TRUE),]
# Select relevant columns
heatmap_nor_matrix <- select(deg_norm_matrix,
symbol,
SAMP8_cortex_six_months1:SAMR1_hippocampus_six_months4)
# Make rownames unique to fix error
heatmap_nor_matrix$symbol <- make.unique(as.character(heatmap_nor_matrix$symbol))
rownames(heatmap_nor_matrix) <- heatmap_nor_matrix$symbol
heatmap_nor_matrix <- data.frame(heatmap_nor_matrix[, -1])# Pick colors
heatmap_colors <- rev(brewer.pal(10, "RdBu"))
# Plot with pheatmap
pheatmap(heatmap_nor_matrix,
color = heatmap_colors,
annotation_col = heatmap_meta,
show_colnames = FALSE,
show_rownames = FALSE,
cluster_cols = FALSE,
clustering_distance_cols = "euclidean",
clustering_method = "complete",
scale = "row",
main = "Heatmap of All genes")
Figure 6. Heatmap of All DE Genes Similarly to Figure
5., more patterns can be seen emerging within the hippocampus samples;
however, broadly there are minimal patterns observed in the heat
map.
#write.csv(deg_matrix, "Deg_Matrix.csv", row.names = FALSE)Alzheimer Gene Analysis
This portion of the analysis primarily focuses on Alzheimer’s disease related genes within the data set. A list of these genes was obtained by cross referencing the differentially expressed genes with a list of 166 Alzheimer related genes from KEGG. Of these 166 genes, 148 were present within this data set. From here, a heat map and volcano plot these specific genes was created.
# Create data frame containing only genes in KEGG_ALZHEIMERS_DISEASE
# Read in alzheimer's reference
alz_ref <- read.csv("Alzheimers_Ref.csv")
# Conver deg_matrix to upper case to match alz_ref
deg_upper <- deg_matrix |>
mutate(symbol = toupper(trimws(as.character(symbol))))
# Use semi join to create new data table
alz_genes <- semi_join(deg_upper, alz_ref, by = "symbol")# Create heatmap of alzheimer genes to see if there are different expression patterns
# Convert to data frame
alz_matrix <- data.frame(alz_genes)
# Order differentially expressed genes by log fold change
alz_matrix <- alz_matrix[order(alz_matrix$log2FoldChange,
decreasing = TRUE),]
# Select relevant columns
heatmap_alz_matrix <- select(alz_matrix,
symbol,
SAMP8_cortex_six_months1:SAMR1_hippocampus_six_months4)
heatmap_alz_matrix <- data.frame(heatmap_alz_matrix, row.names = 1)# Pick colors
heatmap_colors <- rev(brewer.pal(10, "RdBu"))
# Plot with pheatmap
pheatmap(heatmap_alz_matrix,
color = heatmap_colors,
annotation_col = heatmap_meta,
show_colnames = FALSE,
show_rownames = FALSE,
cluster_cols = FALSE,
clustering_distance_cols = "euclidean",
clustering_method = "complete",
scale = "row",
main = "Heatmap of Alzheimer's Related Genes")
Figure 7. Heatmap of Alzheimer’s Disease Related Genes When looking specifically at Alheimer’s disease related genes, there is even more variation in the heatmap, suggesting minimal relation between these genes and different mouse models.
# Create volcano plot of alz genes
ggplot(alz_genes) +
geom_point(aes(x = log2FoldChange, y = -log10(padj), color = diffexpressed, alpha = 0.8)) +
geom_text_repel(aes(x = log2FoldChange,
y = -log10(padj),
label = label),
size = 3.5,
max.overlaps = Inf) +
theme_grey() +
ggtitle("Volcano Plot of Alzheimer's Related Genes") +
xlab("Log2 fold change") +
ylab("-Log10(padj)") +
theme(legend.position = "right",
plot.title = element_text(size = rel(1.25), hjust = 0.5),
axis.title = element_text(size = rel(1.25))) +
scale_color_manual(name = "Differential Expression",
values = c("red", "grey", "blue"),
labels = c("Downregulated", "Not Significant", "Upregulated")) +
scale_x_continuous(limits = c(-20, 20)) +
coord_cartesian(ylim = c(0, 15), clip = "off") +
# Remove alpha legend
guides(alpha = "none")
Figure 8. Volcano Plot of Alzheimer’s Disease Related Genes There was virtually no differential expression of Alzheimer’s disease related genes, with all of the genes tending to be more downregulated.
Plots of gene counts by treatment
Lastly, multiple specific genes that were upregulated, downregulated, or Alzheimer’s disease related genes were selected for further analysis using plot counts. These plot counts can provide information on how specific gene expression compares across treatments, and can provide further illustration of why certain genes are up or down regulated.
# Create plot count of upregulated gene
# Creat plotCounts object from dds
Vip <- plotCounts(dds,
gene="ENSMUSG00000019772.6",
intgroup = "treatment",
returnData = TRUE) |>
# Create blank graph
ggplot(aes(x = treatment, y = count, fill = treatment)) +
# Add geoms
geom_boxplot(outlier.shape = NA, color = "black") +
geom_jitter(shape = 16, position = position_jitter(0.2)) +
# Add labels
labs(title = "Vip Expression",
y = "Counts",
x = NULL) +
# Clean up
theme_grey() +
theme(plot.title = element_text(hjust = 0.5),
axis.text = element_text(angle = 45,
hjust = 1,
size = 8),
axis.title = element_text(size = 10)) +
scale_fill_manual(name = "Treatment",
values = c("tomato", "lightblue", "forestgreen", "orchid"))
Vip
Figure 9. Plot Count of Vip As reflected in the volcano plot, Vip appears to be upregulated in SAMP8 mice compared to SAMR1 mice. Vip counts within SAMP8 mice were around 90 and 95 for the cortex and hippocampus respectively, and around 55 and 20 within SAMR1 mice.
# Create plot count of downregulated gene
# Creat plotCounts object from dds
Il11 <- plotCounts(dds,
gene="ENSMUSG00000004371.16",
intgroup = "treatment",
returnData = TRUE) |>
# Create blank graph
ggplot(aes(x = treatment, y = count, fill = treatment)) +
# Add geoms
geom_boxplot(outlier.shape = NA, color = "black") +
geom_jitter(shape = 16, position = position_jitter(0.2)) +
# Add labels
labs(title = "Il11 Expression",
y = "Counts",
x = NULL) +
# Clean up
theme_grey() +
theme(plot.title = element_text(hjust = 0.5),
axis.text = element_text(angle = 45,
hjust = 1,
size = 8),
axis.title = element_text(size = 10)) +
scale_fill_manual(name = "Treatment",
values = c("tomato", "lightblue", "forestgreen", "orchid"))
Il11
Figure 10. Plot Count of Il11 Il11 appears to be downregulated in SAMP8 mice compared to SAMR1 mice. Il11 counts within SAMP8 mice were around 1.25 for the cortex and 1 for the hippocampus, and around 1 to 3 for SAMR1 mice. However, it appears that this gene is largely downregulated in both SAMP8 and SAMR1, with these being remarkably low counts.
# Create plot count of alzheimer gene
# Creat plotCounts object from dds
Cox5a <- plotCounts(dds,
gene="ENSMUSG00000000088.8",
intgroup = "treatment",
returnData = TRUE) |>
# Create blank graph
ggplot(aes(x = treatment, y = count, fill = treatment)) +
# Add geoms
geom_boxplot(outlier.shape = NA, color = "black") +
geom_jitter(shape = 16, position = position_jitter(0.2)) +
# Add labels
labs(title = "Cox5a Expression",
y = "Counts",
x = NULL) +
# Clean up
theme_grey() +
theme(plot.title = element_text(hjust = 0.5),
axis.text = element_text(angle = 45,
hjust = 1,
size = 8),
axis.title = element_text(size = 10)) +
scale_fill_manual(name = "Treatment",
values = c("tomato", "lightblue", "forestgreen", "orchid"))
Cox5a
Figure 11. Plot Count of Cox5a Cox5a is an Alzheimer’s disease related gene that appears to be upregulated in SAMP8 mice compared to SAMR1 mice. Vip counts within SAMP8 mice were around 155 and 220 for the cortex and hippocampus respectively, and around 115 and 140 within SAMR1 mice.
# Create plot count of alzheimer gene
# Creat plotCounts object from dds
Ndufa9 <- plotCounts(dds,
gene="ENSMUSG00000000399.11",
intgroup = "treatment",
returnData = TRUE) |>
# Create blank graph
ggplot(aes(x = treatment, y = count, fill = treatment)) +
# Add geoms
geom_boxplot(outlier.shape = NA, color = "black") +
geom_jitter(shape = 16, position = position_jitter(0.2)) +
# Add labels
labs(title = "Ndufa9 Expression",
y = "Counts",
x = NULL) +
# Clean up
theme_grey() +
theme(plot.title = element_text(hjust = 0.5),
axis.text = element_text(angle = 45,
hjust = 1,
size = 8),
axis.title = element_text(size = 10)) +
scale_fill_manual(name = "Treatment",
values = c("tomato", "lightblue", "forestgreen", "orchid"))
Ndufa9
Figure 12. Plot Count of Ndufa9 Ndufa9 is another Alzheimer’s disease related gene that appears to be have similar levels of expression in SAMP8 mice compared to SAMR1 mice, with SAMR1 mice having slightly higher expression, particularly in the hippocampus. Ndufa9 counts within SAMP8 mice were around 155 for the cortex and hippocampus, and around 160 in the cortex and 185 in the hippocampus within SAMR1 mice.
# Remove legend from graphs
Vip_noleg <- Vip + theme(legend.position = "none")
Il11_noleg <- Il11 + theme(legend.position = "none")
Cox5a_noleg <- Cox5a + theme(legend.position = "none")
Ndufa9_noleg <- Ndufa9 + theme(legend.position = "none")
# Add legend in again to the side
legend <- get_legend(Vip
+ theme(legend.position = "right"))
# Create multi-paneled figure of box plots
boxplot_grid <- plot_grid(Vip_noleg, Il11_noleg, Cox5a_noleg, Ndufa9_noleg,
labels = c("A", "B", "C", "D"),
align = "v",
ncol = 2)
boxplot_grid <- plot_grid(boxplot_grid, legend, rel_widths = c(4, 2))
boxplot_grid
Figure 13. Multipanneled Figure of Plot Counts of All Selected Genes When looking at the expression of each of these genes in a multipanneled figure, it is easier to see that the expression of Alzheimer’s related genes varies less across treatments compared to known up and down regulated genes.
Conclusion
Frailty, although not a disease itself, is often linked to a number of age related conditions, including neurodegenerative diseases. To further study frailty, mouse models can be utilized including the SAMP8 accelerated age model and the SAMR1 control model. Given this, and that the original paper utilized tissue samples from the hippocampus and cortex, I hypothesized that there would be increased expression of Alzheimer’s disease related genes within the SAMP8 mouse models when compared to the SAMR1 mice.
However, the findings of this bioinformatic analysis demonstrated that there likely is no relation between Alzheimer’s disease related genes and SAMP8 aged model mice. My analysis utilized a KEGG gene set of 166 Alzheimer’s disease related genes, and of these 166 genes, 148 were present in the SAMP8 and SAMR1 mouse data set. When looking at Figure 7.—a heatmap of the expression of these genes—there is no clear pattern in up or down regulation of these genes, with low levels of deregulation being observed across SAMP8 and SAMR1 cortex and hippocampus samples. This lack of differential expression can further be visualized in Figure 8. This volcano plot shows low levels of differential expression across all genes, with all of the genes being slightly downregulated. Furthermore, when comparing gene plot counts of two specific Alzheimer’s disease related genes (Cox5a and Ndufa9) to upregulated and downregulated genes, it continues to be evident that expression of Alzheimer genes remains primarily consistent across SAMP8 and SAMR1 mouse models.
It is worth noting though, that there were several potential limitations and shortcomings of this data set. A primary one being an uneven distribution of reads for each sample. When looking at Figure 3, the amount of reads per sample varied widely, with SAMP8_cortex1 having the largest number of reads at 30 million, and SAMP8_hippocampus3 having the smallest. This large difference in read was most likely due to over amplificaition during PCR or overloading of samples. This increased number of reads can also explain outliers present within the PCA plots in Figures 1. and 2.
Furthermore, based on this bioinformatic analysis, some of the findings of the original paper appear to be overstated. Within the original paper, heatmaps and PCA plots were created only using top differentially expressed genes, meaning that it is likely that patterns of deregulation would be observed anyways. When looking at heatmaps and PCA plots of all differentially expressed genes, as shown in Figures 2. and 6., there is little clustering or patterns in differential expression. While the heatmaps specify that the top 60 differentially expressed genes were chosen, it is not specified how many genes were chosen for the PCA plots, so I assumed it was this same 60 genes in recreating these figures.
When comparing my recreation of the heatmap and PCA plot of differentially expressed genes at 6 months of age to the originals within the paper, there are notable differences between the figures. Both the PCA plot and heatmap of the original paper contained much clearer patterns of clustering and differential expression across SAMP8 and SAMR1. However, this may be due less to failings of the original authors and more as a result of my own limited experience with generating figures in R. To improve upon these results, future studies should limit potential errors in sequencing and more analysis of all differentially expressed genes—not just top genes—should be included. For my own analysis, I would be interested in exploring gene pathways, as well as other age related diseases to see if there may still be a link between frailty mouse models and disease specific biomarkers.
Session info
sessionInfo()## R version 4.4.2 (2024-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sequoia 15.5
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] fstcore_0.10.0 ggplotify_0.1.2
## [3] patchwork_1.3.0 pathview_1.46.0
## [5] msigdbr_10.0.1 cowplot_1.1.3
## [7] knitr_1.50 enrichplot_1.26.6
## [9] ggnewscale_0.5.1 DOSE_4.0.1
## [11] clusterProfiler_4.14.6 EnhancedVolcano_1.24.0
## [13] ggrepel_0.9.6 RColorBrewer_1.1-3
## [15] pheatmap_1.0.12 genekitr_1.2.8
## [17] ggplot2_3.5.2 tidyr_1.3.1
## [19] dplyr_1.1.4 DESeq2_1.46.0
## [21] SummarizedExperiment_1.36.0 Biobase_2.66.0
## [23] MatrixGenerics_1.18.1 matrixStats_1.5.0
## [25] GenomicRanges_1.58.0 GenomeInfoDb_1.42.3
## [27] IRanges_2.40.1 S4Vectors_0.44.0
## [29] BiocGenerics_0.52.0
##
## loaded via a namespace (and not attached):
## [1] rstudioapi_0.17.1 jsonlite_2.0.0 magrittr_2.0.3
## [4] ggtangle_0.0.6 farver_2.1.2 rmarkdown_2.29
## [7] fs_1.6.6 zlibbioc_1.52.0 vctrs_0.6.5
## [10] memoise_2.0.1 europepmc_0.4.3 RCurl_1.98-1.17
## [13] ggtree_3.14.0 htmltools_0.5.8.1 S4Arrays_1.6.0
## [16] progress_1.2.3 SparseArray_1.6.2 gridGraphics_0.5-1
## [19] sass_0.4.10 bslib_0.9.0 plyr_1.8.9
## [22] cachem_1.1.0 igraph_2.1.4 lifecycle_1.0.4
## [25] pkgconfig_2.0.3 Matrix_1.7-3 R6_2.6.1
## [28] fastmap_1.2.0 gson_0.1.0 GenomeInfoDbData_1.2.13
## [31] digest_0.6.37 geneset_0.2.7 aplot_0.2.5
## [34] colorspace_2.1-1 AnnotationDbi_1.68.0 RSQLite_2.3.9
## [37] org.Hs.eg.db_3.20.0 labeling_0.4.3 urltools_1.7.3
## [40] httr_1.4.7 polyclip_1.10-7 abind_1.4-8
## [43] compiler_4.4.2 bit64_4.6.0-1 withr_3.0.2
## [46] BiocParallel_1.40.2 viridis_0.6.5 DBI_1.2.3
## [49] ggforce_0.4.2 R.utils_2.13.0 MASS_7.3-65
## [52] DelayedArray_0.32.0 tools_4.4.2 ape_5.8-1
## [55] zip_2.3.2 fst_0.9.8 R.oo_1.27.0
## [58] glue_1.8.0 nlme_3.1-168 GOSemSim_2.32.0
## [61] grid_4.4.2 ggvenn_0.1.10 reshape2_1.4.4
## [64] fgsea_1.32.4 generics_0.1.3 gtable_0.3.6
## [67] R.methodsS3_1.8.2 data.table_1.17.0 hms_1.1.3
## [70] tidygraph_1.3.1 xml2_1.3.8 XVector_0.46.0
## [73] pillar_1.10.2 stringr_1.5.1 babelgene_22.9
## [76] yulab.utils_0.2.0 splines_4.4.2 tweenr_2.0.3
## [79] treeio_1.30.0 lattice_0.22-7 bit_4.6.0
## [82] tidyselect_1.2.1 GO.db_3.20.0 locfit_1.5-9.12
## [85] Biostrings_2.74.1 gridExtra_2.3 xfun_0.52
## [88] graphlayouts_1.2.2 KEGGgraph_1.66.0 stringi_1.8.7
## [91] UCSC.utils_1.2.0 lazyeval_0.2.2 ggfun_0.1.8
## [94] yaml_2.3.10 pacman_0.5.1 evaluate_1.0.3
## [97] codetools_0.2-20 ggraph_2.2.1 tibble_3.2.1
## [100] qvalue_2.38.0 Rgraphviz_2.50.0 graph_1.84.1
## [103] cli_3.6.4 munsell_0.5.1 jquerylib_0.1.4
## [106] Rcpp_1.0.14 triebeard_0.4.1 png_0.1-8
## [109] XML_3.99-0.18 parallel_4.4.2 assertthat_0.2.1
## [112] blob_1.2.4 prettyunits_1.2.0 bitops_1.0-9
## [115] viridisLite_0.4.2 tidytree_0.4.6 scales_1.3.0
## [118] openxlsx_4.2.8 purrr_1.0.4 crayon_1.5.3
## [121] rlang_1.1.6 fastmatch_1.1-6 KEGGREST_1.46.0