Replication - Linear and/or logistic regression - Mount Sinai Brain Bank (MSBB) 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 data | Brain region: frontal lobe | 304 unique samples
file_path <- "C:/Users/beker/OneDrive/Documentos/Mestrado/GitHub/Cytokines/replication/T026/RNAseq_Harmonization_MSBB_combined_metadata.csv"
pheno_data <- read.csv(file_path)# 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: 304Correlation
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]# 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",
"CERAD" = "CERAD score",
"CDR" = "CDR scale",
"plaqueMean" = "Plaque mean"
)
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 = "Mount Sinai bulk RNA-seq data",
row_names_side = "right", show_row_names = T,
cluster_rows = F, cluster_columns = F,
column_names_gp = gpar(fontsize = 9, fontface = "italic"),
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, 4, 1, 1), "mm")) # Ajustar as margens (top, right, bottom, left)
# Save to PDF
pdf(file = "C:/Users/beker/OneDrive/Documentos/Mestrado/GitHub/Cytokines/replication/T027/lr_expr_bulk_MSBB_T027.pdf", width = 6.4, height = 1.95)
print(ht_draw)
dev.off()## png
## 2# Save to PNG
png(file = "C:/Users/beker/OneDrive/Documentos/Mestrado/GitHub/Cytokines/replication/T027/lr_expr_bulk_MSBB_T027.png", width = 6.4, height = 1.95, 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 magick_2.8.4 backports_1.4.1
## [16] utf8_1.2.4 rmarkdown_2.27 tzdb_0.4.0
## [19] nloptr_2.1.1 bit_4.0.5 xfun_0.46
## [22] cachem_1.1.0 jsonlite_1.8.8 highr_0.11
## [25] broom_1.0.6 parallel_4.3.2 cluster_2.1.4
## [28] R6_2.5.1 bslib_0.8.0 stringi_1.8.3
## [31] car_3.1-2 boot_1.3-28.1 jquerylib_0.1.4
## [34] cellranger_1.1.0 numDeriv_2016.8-1.1 Rcpp_1.0.12
## [37] iterators_1.0.14 knitr_1.48 IRanges_2.36.0
## [40] splines_4.3.2 timechange_0.3.0 tidyselect_1.2.1
## [43] rstudioapi_0.16.0 abind_1.4-5 yaml_2.3.10
## [46] doParallel_1.0.17 codetools_0.2-19 plyr_1.8.9
## [49] lattice_0.21-9 withr_3.0.1 evaluate_0.24.0
## [52] zip_2.3.1 pillar_1.9.0 carData_3.0-5
## [55] DT_0.33 foreach_1.5.2 stats4_4.3.2
## [58] insight_0.20.3 generics_0.1.3 vroom_1.6.5
## [61] S4Vectors_0.40.2 hms_1.1.3 munsell_0.5.1
## [64] scales_1.3.0 minqa_1.2.7 glue_1.7.0
## [67] tools_4.3.2 crosstalk_1.2.1 colorspace_2.1-0
## [70] nlme_3.1-163 cli_3.6.2 fansi_1.0.6
## [73] gtable_0.3.5 sass_0.4.9 digest_0.6.36
## [76] BiocGenerics_0.48.1 htmlwidgets_1.6.4 rjson_0.2.21
## [79] htmltools_0.5.8.1 lifecycle_1.0.4 GlobalOptions_0.1.2
## [82] bit64_4.0.5 MASS_7.3-60