library(readr)
library(org.Hs.eg.db)
library(Hmisc)
library(corrplot)
library(zeitzeiger)
library(tidyr)
library(ggplot2)

Data were retrieved from the GTEx Portal and filtered to include only samples from donors with minimal clinical agony (Hardy Scale 1 or 2) to ensure high transcriptomic integrity (Lonsdale et al. 2013). Further filtering retained only genes with an expression level of >0.5 TPM in at least 50% of samples for downstream analysis

hypothalamus_data <- as.data.frame(read_delim("gene_tpm_v11_brain_hypothalamus.gct.gz",
delim = "\t", escape_double = FALSE,
trim_ws = TRUE, skip = 2))
sample_metadata <- read.delim("GTEx_Analysis_v11_Annotations_SampleAttributesDS.txt")
subject_pheno <- read.delim("GTEx_Analysis_v11_Annotations_SubjectPhenotypesDS.txt")
fast_death_subjects <- subject_pheno[subject_pheno$DTHHRDY %in% c(1, 2), "SUBJID"]
hypothalamus_data <- hypothalamus_data[!grepl("^ENSG", hypothalamus_data$Description), ]
hypothalamus_data <- hypothalamus_data[,-1]
hypothalamus_data <- aggregate(. ~ Description, data = hypothalamus_data, FUN = sum)
rownames(hypothalamus_data) <- hypothalamus_data[,1]
hypothalamus_fd <- names(hypothalamus_data)[grepl(paste(fast_death_subjects, collapse="|"), names(hypothalamus_data))]
hypothalamus_data <- hypothalamus_data[, hypothalamus_fd]
is_expressed <- rowMeans(hypothalamus_data > 0.5) > 0.5
hypothalamus_test <- hypothalamus_data[is_expressed, ]

Initially, we estimated the correlations between components of the dopaminergic system and genes involved in circadian rhythms. The analysis was conducted using the corrplot R package (v. 0.95) (Wei et al. 2021).

circadian_genes <- select(org.Hs.eg.db,
keys = "GO:0007623",
keytype = "GOALL",
columns = "SYMBOL")
# Genes demonstrating the most robust association with circadian rhythmicity were selected for further analysis, namely:
circadian_genes <- circadian_genes[circadian_genes$EVIDENCEALL %in% 
                                     c("EXP", #experimental
                                    "IDA", #direct assay
                                    "IPI", #physical interactions
                                    "IMP", #mutant phenotype
                                    "IGI", #genetic interaction
                                    "IEP"),] #or expression pattern
circadian_list <- circadian_genes$SYMBOL
DA_list <- c("DRD1", "DRD2", "DRD4", "DRD5", "TH", "DDC", "SLC18A2", "COMT", "MAOA", "MAOB") # DRD3 and SLC6A3 (DAT) were excluded from this analysis due to low expression levels. These genes were omitted during the initial filtration step.
combined_data <- t(rbind(hypothalamus_test[DA_list, ], hypothalamus_test[circadian_list, ]))
res <- rcorr(as.matrix(combined_data), type = "spearman")
cor_sub <- res$r[rownames(res$r) %in% DA_list, colnames(res$r) %in% circadian_list]
p_sub <- res$P[rownames(res$r) %in% DA_list, colnames(res$r) %in% circadian_list]
corrplot(cor_sub,
method = "color",
type = "full",
p.mat = p_sub,
sig.level = 0.05,     
insig = "blank",      
addCoef.col = "black",
number.digits = 1,
number.cex = 0.7,
font=1,
tl.col = "black",
cl.lim = c(-1, 1),
col = colorRampPalette(c( "magenta4", "lightyellow", "orange2"))(200))
title(main = "Correlations between expression of genes involved in dopamine signaling and circadian genes", 
      font.main = 1,
      line = 1.2,   
      cex.main = 1.3)

Since the time of death is not provided in the GTEx metadata, we estimated biological time directly from the transcriptomic profiles using the zeitzeiger R package (v. 2.1.3) (Hughey, Hastie, and Butte 2016). We applied an unsupervised circadian reconstruction approach to predict the phase for each sample. Data visualization was performed using the ggplot2 R package (v. 3.5.2) (Wickham 2016).

x <- t(hypothalamus_test[rownames(hypothalamus_test) %in% circadian_list, ]) 
pca_res <- prcomp(x)
pc1 <- pca_res$x[,1]
pc2 <- pca_res$x[,2]
pseudotime <- atan2(pc2, pc1)
pseudotime <- (pseudotime + pi) / (2 * pi) * 24
plot_df <- data.frame(
Pseudotime = pseudotime,
log2(t(hypothalamus_test[DA_list, , drop=FALSE])+1)
) %>%
pivot_longer(cols = -Pseudotime, names_to = "Gene", values_to = "Expression")
ggplot(plot_df, aes(x = Pseudotime, y = Expression, color = Gene)) +
geom_smooth(method = "gam", formula = y ~ s(x, bs = "cc"), size = 1.5) +
scale_x_continuous(breaks = seq(0, 24, by = 4), limits = c(0, 24)) +
theme_minimal() +
theme(axis.title = element_text(size = 14), axis.text = element_text(size = 12)) +
labs(title = "Gene Expression along Reconstructed Circadian Phase",
x = "Pseudotime (Hours 0-24)",
y = "Relative Expression (Z-score)") +
theme(legend.position = "bottom")+
scale_color_manual(values = c("gold2", "coral", "forestgreen", "skyblue2", "cyan3", "red2", "black", "yellow2", "darkmagenta", "brown"))

Hughey, Jacob J., Trevor Hastie, and Atul J. Butte. 2016. “ZeitZeiger: Supervised Learning for High-Dimensional Data from an Oscillatory System.” Nucleic Acids Research 44 (8): e80. https://doi.org/10.1093/nar/gkw030.
Lonsdale, John, Jeffrey Thomas, Mike Salvatore, Rebecca Phillips, Edmund Lo, Saboor Shad, Richard Hasz, et al. 2013. “The Genotype-Tissue Expression (GTEx) Project.” Nature Genetics 45 (6): 580–85. https://doi.org/10.1038/ng.2653.
Wei, T, V Simko, M Levy, Y Xie, Y Jin, and J Zemla. 2021. “R Package ‘Corrplot’: Visualization of a Correlation Matrix [Internet]. 2021 Jul.” 2021.
Wickham, Hadley. 2016. Ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org.