Load Packages

library(tximport) library(DESeq2) library(biomaRt) library(vsn) library(dplyr) library(pheatmap) library(gplots) library(RColorBrewer) library(ComplexHeatmap) library(readr) library(circlize) library(EnhancedVolcano) library(ggrepel) library(edgeR)

Set Directory

dir <- “~/Desktop/Feature-Counts/PV1SN3/All/Longitudinal/”

Read in Meta Data

file<- paste0(dir, “All-Genes-All-AD-Male-Longitudinal-Meta.csv”) samples <- read.csv(file, stringsAsFactors = FALSE, header = TRUE, sep = “,”)

Assign Sample Variables as Factor

samples$Age<-factor(samples$Age, levels = c(“6”,“9”,“12”,“15”,“18”))

Set Up Variable Data Frame

sampleTable <- data.frame(sampleName = samples$Name, Age = samples$Age)

Load Data

data<-read.csv(paste0(dir, “All-Genes-All-AD-Male-Longitudinal-Counts.csv”), sep = “,”, stringsAsFactors = F, header = T)

Read data into DESeq2

data.1<-as.matrix(data[,2:18]) rownames(data.1)<-data[,1] dds <- DESeqDataSetFromMatrix(countData = data.1, colData = sampleTable, design= ~ Age)

Drop Low Read Counts

smallestGroupSize <- 3 keep <- rowSums(counts(dds) >=10) >= smallestGroupSize dds <- dds[keep,]

Run DESeq

dds_age <- DESeq(dds)

DESeq Results

Repeat for each age comparison to baseline [coef values correspond to each comparison from resultsNames: Age_9_vs_6, Age_12_vs_6, Age_15_vs_6, etc.]

resultsNames(dds_age) resLFC <- lfcShrink(dds_age, coef = 5, type = “apeglm”)

resLFC_ord <- resLFC[order(resLFC$padj),]

Write csv file for complete results

write.csv(as.data.frame(resLFC_ord), file = paste0(dir,“Males-AD-Longitudinal-Age-18v6.csv”), quote = FALSE, row.names = T, col.names = T, sep = “,”)

res_file<- paste0(dir, “Males-AD-Longitudinal-Age-18v6.csv”) all_data <- read.csv(res_file, stringsAsFactors = FALSE, header = TRUE, sep = “,”)

Generate Volcano Plot

ppi=300 png(paste0(dir, “18v6_M_ADvWT_volcanoPlot_DEG_Males_AD”, “.png”), width=7.5ppi, height=6ppi, res=ppi) EnhancedVolcano(all_data, lab = NA, x=“log2FoldChange”,y=“padj”, #selectLab = labels, xlim=c(-4,4),ylim=c(0.5,5), axisLabSize=10, title= “AD 18 vs. 6 month old mice”, titleLabSize=10, subtitle = NULL, captionLabSize = 6, pCutoff=0.05, FCcutoff = 2, labSize = 2.5, legendPosition = “right”, legendLabSize = 6, #shape=c(6,6,19,16), caption = “FC cutoff, 2; p-value cutoff, 0.05”) dev.off()

Generate Heatmap

Subset @ 5% adjusted pval sig. level

res_order_FDR_05 <- subset(res_ord, padj <0.05) nrow(res_order_FDR_05)

Take the regularized log

rld <- rlog(dds_age, blind = FALSE)

mat2 <- assay(rld)[rownames(res_order_FDR_05),] head(mat2)

Scale matrix by each col. values

mat_scaled = t(apply(mat2, 1, scale)) print(mat_scaled)

Set up colors for heatmap

col = colorRamp2(c(-3, 0, 3), c(“blue”, “white”, “red”)) cols1 <- brewer.pal(11, “Paired”)

Subset coldata for samples in untx and ex groups

colData_sub <- colData(dds_age)

Set up annotation bar for samples

ha1 = HeatmapAnnotation(Age = colData_sub$Age, col = list(Age = c(“6” = “#66BF39”, “9” = “#106B03”, “12” = “#FEFF45”, “15” = “#FFD505”, “18” = “#FF8000”), show_legend = TRUE))

Set up column annotation labels (samples)

ha = columnAnnotation(x = anno_text(colData_sub$sampleName, which=“column”, rot = 45, gp = gpar(fontsize = 10)))

Generate heatmap object

ht1 = Heatmap(mat_scaled, name = “Expression”, col = col, top_annotation = c(ha1), bottom_annotation = c(ha), show_row_names = FALSE, show_column_names = FALSE, show_row_dend = FALSE)

Plot the heatmap

png(paste0(dir, “18v6M-AD-heatmap_DEG_males”, “.png”), width=7.5ppi, height=6ppi, res=ppi) draw(ht1, row_title = “Genes”, column_title = “Hierarchical clustering of DEGs (padj<0.05)”) dev.off()