Replication - Linear and/or logistic regression - Mayo Clinic dataset
Juliana Beker Godoy,
MSc
Graduate Program in Bioinformatics
Professional and Technological Education Sector
Federal University of Parana (UFPR)
Curitiba, Brazil
Graduate Program in Bioinformatics
Professional and Technological Education Sector
Federal University of Parana (UFPR)
Curitiba, Brazil
2025-06-02
Bulk RNA-Seq | Brain region: temporal cortex | 258 unique samples
file_path <- "C:/Users/beker/OneDrive/Documentos/Mestrado/GitHub/Cytokines/replication/T028/RNAseq_Harmonization_Mayo_combined_metadata.csv"
pheno_data <- read.csv(file_path)
pheno_data$diagnosis = factor(pheno_data$diagnosis, levels = c("control", "Alzheimer Disease", "pathological aging", "progressive supranuclear palsy" ))# upload list of cytokines
file_path <- "C:/Users/beker/OneDrive/Documentos/Mestrado/GitHub/Cytokines/list_cytokines/T023/list_cytokines_T023.xlsx"
cytokines <- read_excel(file_path)
cytokines <- subset(cytokines, select = -family) # I removed the 'family' column here because it had NA values
ensembls = cytokines$ensembl## Number of samples: 258Correlation
to_show = colnames(res_test$matrix_pvalue)
transpose = T
show_only_significant = F; signif_cutoff = c("***","**","*")
matrix_rsquared = res_test$matrix_rsquared
matrix_pvalue = res_test$matrix_pvalue # final matrix with the pvalues
matrix_rsquared_to_plot = matrix_rsquared[,to_show]
matrix_pvalue_to_plot = matrix_pvalue[,to_show]
# Saving xlsx
write.xlsx(matrix_pvalue_to_plot, file = "C:/Users/beker/OneDrive/Documentos/Mestrado/GitHub/Cytokines/replication/T029/nom_p_values_MAYO_T029.xlsx", row.names = TRUE)# Heatmap
to_show = colnames(res_test$matrix_pvalue)
transpose = T
show_only_significant = T; signif_cutoff = c("***","**","*")
matrix_rsquared = res_test$matrix_rsquared
matrix_pvalue = res_test$matrix_pvalue # final matrix with the pvalues
matrix_rsquared_to_plot = matrix_rsquared[,to_show]
matrix_pvalue_to_plot = matrix_pvalue[,to_show]
# Adjust P-values by each phenotype separately.
adj_matrix_pvalue_to_plot = matrix_pvalue_to_plot
for(i in 1:ncol(matrix_pvalue_to_plot)){
adj_matrix_pvalue_to_plot[,i] = p.adjust(matrix_pvalue_to_plot[,i], method = "bonferroni")
}
adj_matrix_pvalue_to_plot.signif <- symnum(x = as.matrix(adj_matrix_pvalue_to_plot), corr = FALSE, na = FALSE,
cutpoints = c(0, 0.001, 0.01, 0.05, 0.1, 1),
symbols = c("***", "**", "*", ".", " "))
log_matrix_pvalue_to_plot = -log10(matrix_pvalue_to_plot)
dimnames(log_matrix_pvalue_to_plot) = dimnames(log_matrix_pvalue_to_plot)
if(show_only_significant){
if(is.numeric(signif_cutoff)){
to_keep = colSums(adj_matrix_pvalue_to_plot <= signif_cutoff) > 0
}else{
to_keep = rep(F,ncol(adj_matrix_pvalue_to_plot.signif))
for(cut_i in signif_cutoff){
to_keep = to_keep | colSums(adj_matrix_pvalue_to_plot.signif == cut_i) > 0 # change for the significance you want
}
}
log_matrix_pvalue_to_plot = log_matrix_pvalue_to_plot[,to_keep]
adj_matrix_pvalue_to_plot.signif = adj_matrix_pvalue_to_plot.signif[,to_keep]
}
matrix_pvalue_to_plot_labels = formatC(log_matrix_pvalue_to_plot, format = "f", digits = 2)
log_matrix_pvalue_to_plot_t = t(log_matrix_pvalue_to_plot)
if(transpose){
log_matrix_pvalue_to_plot_t = t(log_matrix_pvalue_to_plot_t)
matrix_pvalue_to_plot_labels = t(matrix_pvalue_to_plot_labels)
adj_matrix_pvalue_to_plot.signif = t(adj_matrix_pvalue_to_plot.signif)
}
# Colored by -log10(pvalue)
# Numbers inside cell = -log10(pvalue): nominal
new_column_names <- c(
"Braak" = "Braak stages",
"thal" = "Thal amyloid stages",
"diagnosis" = "Diagnosis"
)
rownames(log_matrix_pvalue_to_plot_t) <- new_column_names[rownames(log_matrix_pvalue_to_plot_t)]
heatmap_plot <- Heatmap(log_matrix_pvalue_to_plot_t, name = "-log10(P-value)",
cell_fun = function(j, i, x, y, width, height, fill) {
if(as.character(t(adj_matrix_pvalue_to_plot.signif)[i,j]) == " "){
grid.text( t(matrix_pvalue_to_plot_labels)[i,j], x, y,
gp = gpar(fontsize = 8))
}else{
grid.text(paste0( t(matrix_pvalue_to_plot_labels)[i,j],"\n", t(adj_matrix_pvalue_to_plot.signif)[i,j] ), x, y,
gp = gpar(fontsize = 8))
}
},
col = colorRampPalette(rev(brewer.pal(n = 7, name ="RdYlBu")))(100),
column_names_rot = 40,
column_title = "Mayo Clinic bulk RNA-seq data",
row_names_side = "right", show_row_names = T,
cluster_rows = F, cluster_columns = F,
column_names_gp = gpar(fontsize = 9),
row_names_gp = gpar(fontsize = 9),
border = T,
show_row_dend = F, show_column_dend = F, rect_gp = gpar(col = "white", lwd = 1))
ht_draw <- draw(heatmap_plot, padding = unit(c(1, 8, 1, 2), "mm")) # Ajustar as margens (top, right, bottom, left)
# Save to PDF
pdf(file = "C:/Users/beker/OneDrive/Documentos/Mestrado/GitHub/Cytokines/replication/T029/lr_expr_bulk_Mayo_T029.pdf", width = 12, height = 1.8)
print(ht_draw)
dev.off()## png
## 2# Save to PNG
png(file = "C:/Users/beker/OneDrive/Documentos/Mestrado/GitHub/Cytokines/replication/T029/lr_expr_bulk_Mayo_T029.png", width = 12, height = 1.8, units = "in", res = 300)
print(ht_draw)
dev.off()## png
## 2Top results
Top result by covariate.
Session info
## R version 4.3.2 (2023-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 11 x64 (build 26100)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=Portuguese_Brazil.utf8 LC_CTYPE=Portuguese_Brazil.utf8
## [3] LC_MONETARY=Portuguese_Brazil.utf8 LC_NUMERIC=C
## [5] LC_TIME=Portuguese_Brazil.utf8
##
## time zone: America/Sao_Paulo
## tzcode source: internal
##
## attached base packages:
## [1] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] performance_0.12.2 lmerTest_3.1-3 lme4_1.1-35.5
## [4] Matrix_1.6-5 openxlsx_4.2.6.1 readxl_1.4.3
## [7] RColorBrewer_1.1-3 circlize_0.4.16 ComplexHeatmap_2.18.0
## [10] ggsignif_0.6.4 ggeasy_0.1.4 ggpubr_0.6.0
## [13] lubridate_1.9.3 forcats_1.0.0 stringr_1.5.1
## [16] purrr_1.0.2 tidyr_1.3.1 tibble_3.2.1
## [19] tidyverse_2.0.0 ggplot2_3.5.0 readr_2.1.5
## [22] rstatix_0.7.2 dplyr_1.1.4
##
## loaded via a namespace (and not attached):
## [1] rlang_1.1.3 magrittr_2.0.3 clue_0.3-65
## [4] GetoptLong_1.0.5 matrixStats_1.3.0 compiler_4.3.2
## [7] png_0.1-8 vctrs_0.6.5 reshape2_1.4.4
## [10] pkgconfig_2.0.3 shape_1.4.6.1 crayon_1.5.3
## [13] fastmap_1.2.0 backports_1.4.1 magick_2.8.4
## [16] utf8_1.2.4 rmarkdown_2.27 tzdb_0.4.0
## [19] nloptr_2.1.1 xfun_0.46 cachem_1.1.0
## [22] jsonlite_1.8.8 highr_0.11 broom_1.0.6
## [25] parallel_4.3.2 cluster_2.1.4 R6_2.5.1
## [28] bslib_0.8.0 stringi_1.8.3 car_3.1-2
## [31] boot_1.3-28.1 jquerylib_0.1.4 cellranger_1.1.0
## [34] numDeriv_2016.8-1.1 Rcpp_1.0.12 iterators_1.0.14
## [37] knitr_1.48 IRanges_2.36.0 splines_4.3.2
## [40] timechange_0.3.0 tidyselect_1.2.1 rstudioapi_0.16.0
## [43] abind_1.4-5 yaml_2.3.10 doParallel_1.0.17
## [46] codetools_0.2-19 lattice_0.21-9 plyr_1.8.9
## [49] withr_3.0.1 evaluate_0.24.0 zip_2.3.1
## [52] pillar_1.9.0 carData_3.0-5 DT_0.33
## [55] foreach_1.5.2 stats4_4.3.2 insight_0.20.3
## [58] generics_0.1.3 S4Vectors_0.40.2 hms_1.1.3
## [61] munsell_0.5.1 scales_1.3.0 minqa_1.2.7
## [64] glue_1.7.0 tools_4.3.2 crosstalk_1.2.1
## [67] colorspace_2.1-0 nlme_3.1-163 cli_3.6.2
## [70] fansi_1.0.6 gtable_0.3.5 sass_0.4.9
## [73] digest_0.6.36 BiocGenerics_0.48.1 htmlwidgets_1.6.4
## [76] rjson_0.2.21 htmltools_0.5.8.1 lifecycle_1.0.4
## [79] GlobalOptions_0.1.2 MASS_7.3-60