OSCR

Identification of suicide brain transcriptomic signatures using meta-analysis of multiple cohorts.

Code ↔ Paper

23 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 23 matches
  1. [1] § Materials and methods › Differential expression analysis in the primary cohorts ↔ Cell_expression_imputation_analysis.R, lines 1713 to 1777 · score 0.87 · filterByExpr, DGEList, lmFit, topTable, eBayes, design matrix · Report a mismatch
  2. [2] § Results › Gene set enrichment analysis and classification ↔ AI_gene_classification.ipynb, lines 72 to 96 · score 0.83 · secreted molecules, immune receptors, structural proteins, transcription factors, enzymes, kinase · Report a mismatch
  3. [3] § Results › Generation of datasets ↔ Moderator_calculation.R, lines 79 to 160 · score 0.79 · GSE92538_U133A, GSE92538_U133_PLUS2, dorsolateral prefrontal cortex, brain regions, sex, GSE248260 · Report a mismatch
  4. [4] § Materials and methods › Probe and transcript mapping to gene symbols ↔ Data_preprocessing_analysis/TRANSCRIPT_SUICIDE_PREPR_ANALYSIS_SCRIPT.R, lines 667 to 709 · score 0.77 · BioMart, probe IDs, map probes, gene symbols, Ensembl, Python · Report a mismatch
  5. [5] § Results › Generation of datasets ↔ OLINK_PSY_SUICIDE_SCRIPT.R, lines 2338 to 2373 · score 0.74 · GSE92538_U133A, GSE92538_U133_PLUS2, Single cell, prefrontal cortex, GSE248260, GSE199536 · Report a mismatch
  6. [6] § Materials and methods › Cohort identification ↔ OLINK_PSY_SUICIDE_SCRIPT.R, lines 1656 to 1699 · score 0.74 · genome tiling, MeSH, Expression profiling, Fields, Organism, human · Report a mismatch
  7. [7] § Results › Generation of datasets ↔ Data_preprocessing_analysis/DEMOGRAPHICS_TRANSCRIPT_SCRIPT.R, lines 1193 to 1223 · score 0.73 · aINS, insula, BA11, Cg25, OFC, gyrus · Report a mismatch
  8. [8] § Results › Generation of datasets ↔ Data_preprocessing_analysis/TRANSCRIPT_SUICIDE_PREPR_ANALYSIS_SCRIPT.R, lines 3165 to 3210 · score 0.73 · aINS, insula, BA11, Cg25, OFC, gyrus · Report a mismatch
  9. [9] § Results › Random effects meta-analyses ↔ OLINK_PSY_SUICIDE_SCRIPT.R, lines 2436 to 2503 · score 0.72 · AC159540.14, CX3CR1, P2RY12, FOS, DEPP1, KLF4 · Report a mismatch
  10. [10] § Results › Sensitivity and moderator analyses ↔ OLINK_PSY_SUICIDE_SCRIPT.R, lines 2436 to 2503 · score 0.71 · AC159540.14, CX3CR1, P2RY12, DEPP1, KLF4, SOX9 · Report a mismatch
  11. [11] § Results › Cell deconvolution analysis and cell-specific meta-analysis ↔ OLINK_PSY_SUICIDE_SCRIPT.R, lines 2145 to 2192 · score 0.69 · precursor cells, excitatory neurons, mixed expression, oligodendrocytes, suicide, meta · Report a mismatch
  12. [12] § Results › Generation of datasets ↔ Data_preprocessing_analysis/DEMOGRAPHICS_TRANSCRIPT_SCRIPT.R, lines 995 to 1047 · score 0.65 · major depressive disorder, bipolar disorder, GSE92538 U133A, schizophrenia, diagnoses, DLPFC · Report a mismatch
  13. [13] § Results › Generation of datasets ↔ Moderator_calculation.R, lines 270 to 312 · score 0.64 · major depressive disorder, bipolar disorder, GSE92538 U133A, schizophrenia, GSE5388, GSE208338 · Report a mismatch
  14. [14] § Results › Comparison of brain meta-analysis with blood transcriptome and PSY proteome ↔ OLINK_PSY_SUICIDE_SCRIPT.R, lines 1510 to 1562 · score 0.63 · PLA2G10, AKT1S1, CCL27, RBKS, antidepressant, PSY · Report a mismatch
  15. [15] § Materials and methods › Cohort identification ↔ OLINK_PSY_SUICIDE_SCRIPT.R, lines 1745 to 1789 · score 0.62 · MeSH, Expression profiling, Fields, sequencing, Organism, human · Report a mismatch
  16. [16] § Results › Differential expression analysis in the individual cohorts ↔ Moderator_calculation.R, lines 79 to 160 · score 0.62 · nucleus accumbens, GSE92538_U133_PLUS2, dorsolateral prefrontal cortex, GSE202537, GSE144136, GSE213982 · Report a mismatch
  17. [17] § Materials and methods › Cohort identification ↔ OLINK_PSY_SUICIDE_SCRIPT.R, lines 2338 to 2373 · score 0.61 · GSE92538 U133A, GSE248260, GSE199536, PSY, GSE247998, GSE101521 · Report a mismatch
  18. [18] § Materials and methods › Meta-analysis of differential expression changes ↔ Cell_expression_imputation_analysis.R, lines 2290 to 2369 · score 0.57 · p.adjust, rma.uni, Weights, SE, FDR, model · Report a mismatch
  19. [19] § Materials and methods › Cell deconvolution meta-analysis ↔ Cell_expression_imputation_analysis.R, lines 775 to 818 · score 0.54 · sampled signature matrix, single cell, bulk, CIBERSORTx, GSE144136, GSE213982 · Report a mismatch
  20. [20] § Materials and methods › Bulk RNAseq data preprocessing ↔ Data_preprocessing_analysis/TRANSCRIPT_SUICIDE_PREPR_ANALYSIS_SCRIPT.R, lines 3873 to 3935 · score 0.53 · FastQC, SRA, conda, tools, preprocessed, STAR · Report a mismatch
  21. [21] § Materials and methods › Cohort identification ↔ Moderator_calculation.R, lines 314 to 361 · score 0.53 · major depressive disorder, psychiatric diagnoses, suicide · Report a mismatch
  22. [22] § Materials and methods › Meta-analysis of differential expression changes ↔ OLINK_PSY_SUICIDE_SCRIPT.R, lines 675 to 718 · score 0.52 · rma.uni, Log2FC, Weights, SE, model, meta · Report a mismatch
  23. [23] § Materials and methods › Meta-analysis of differential expression changes ↔ Data_preprocessing_analysis/TRANSCRIPT_SUICIDE_PREPR_ANALYSIS_SCRIPT.R, lines 4204 to 4256 · score 0.50 · causal pathway, status, depression, diagnosis, phenotypes, variables · Report a mismatch

Paper

Loaded from Europe PMC by your browser; OSCR does not store it: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 4,359 lines · 195 KB · GPL-3.0 · 9 matches

  1. # The script is intended viewed in R studio
  2. # Blocks are titled
  3. # "=" symbol was used as an assignment operator
  4. # Data processing and DE in initial cohorts is located in another file "Data_preprocessing_analysis/TRANSCRIPT_SUICIDE_PREPR_ANALYSIS_SCRIPT.R"
  5. # Calculation of cohort-level moderators is located in "Moderator_calculation.R"
  6. # Cell expression deconvolution tests are located in "Cell_expression_imputation_tests.R"
  7. # Cell expression deconvolution final runs and analysis are located in "Cell_expression_imputation_analysis.R"
  10. setwd("/home/aleksandr/Desktop/WORK/OLINK_suicide_PSY_project")
  12. # Setting options
  13. getOption("scipen") # Default number notation is 0
  14. options(scipen=999)
  15. options(stringsAsFactors = FALSE)
  16. options(show.error.messages = TRUE)
  18. ################### Package import ###################
  19. library(R.utils)
  20. library(fun)
  21. library(stringr)
  22. library(dplyr)
  23. library(ggplot2)
  24. library(openxlsx)
  25. library(grid)
  26. library(gdata)
  27. library(RColorBrewer)
  28. library(networkD3)
  29. library(webshot)
  30. library(htmlwidgets)
  31. library(magrittr)
  32. library(igraph)
  33. library(visNetwork)
  34. library(data.table)
  35. library(XML)
  36. library(rvest)
  37. library(RCurl)
  38. library(HGNChelper)
  39. library(stringi)
  40. library(httr)
  41. library(lubridate)
  42. library(rjson)
  43. library(rtracklayer)
  44. library(rstudioapi)
  45. library(tidyr)
  46. library(Gviz)
  47. library(limma)
  48. library(FactoMineR)
  49. library(ggthemes)
  50. library(igraph)
  51. library(RSelenium)
  52. library(lumi)
  53. library(outliers)
  54. library(svglite)
  55. library(scatterplot3d)
  56. library(sva)
  57. library(jsonlite)
  58. library(ggrepel)
  59. library(parallel)
  60. library(bacon)
  61. library(gridExtra)
  62. library(ggplotify)
  63. library(HGNChelper)
  64. library(jetset)
  65. library(GEOquery)
  66. library(clusterProfiler)
  67. library(DOSE)
  68. library(enrichplot)
  69. library(chromoMap)
  70. library(RIdeogram)
  71. library(ggVennDiagram)
  72. library(seqinr)
  73. library(Biostrings)
  74. library(dbparser)
  75. library(LDlinkR)
  76. library(parallel)
  77. library(HardyWeinberg)
  78. library(sva)
  79. library(MatrixEQTL)
  80. library(MASS)
  81. library(metafor)
  82. library(AnnotationDbi)
  83. library(org.Hs.eg.db)
  84. library(RobustRankAggreg)
  86. ################### Defining basic functions ###################
  88. # NOT IN operator
  89. '%!in%' = function(x,y){!('%in%'(x,y))}
  91. # A function to convert list to data frame
  92. list_to_df = function(data_list){
  93. if (length(data_list) > 1){
  94. data_list = do.call(rbind, data_list)
  95. } else {
  96. data_list = data_list[[1]]
  97. }
  98. return(data_list)
  99. }
  101. # A function to replace multiple patterns by multiple replacements in a string
  102. multiple_stri_replacer = function(string, pattern_vector, replacement_vector){
  104. # Pattern_vector and replacement_vector should have the same length
  105. for (i in 1:length(pattern_vector)){
  106. string = stri_replace_all_fixed(str = string, pattern = pattern_vector[i], replacement = replacement_vector[i])
  107. }
  108. return(string)
  109. }
  111. # A function to read text files fast; uses data.table::fread
  112. smart_fread = function(x, ...){
  113. x = as.data.frame(fread(x, nThread = 10, ...))
  114. if ("V1" %in% colnames(x)){
  115. rownames(x) = x$V1
  116. x$V1 = NULL
  117. }
  118. return(x)
  119. }
  121. # A function to detect at least one pattern in a string
  122. multiple_stri_detector = function(string, pattern_vector){
  123. output_list = list()
  124. for (i in 1:length(pattern_vector)){
  125. output_list[[i]] = stri_detect_fixed(str = string, pattern = pattern_vector[i])
  126. }
  127. output_list = do.call(rbind, output_list)
  128. apply(output_list, 2, any)
  129. }
  131. # A function to expand a data frame where several columns contain condensed cells with a specified separator
  132. multiple_expander = function(df, cols_to_expand, pattern){
  133. #
  134. orig_colnames = colnames(df)
  135. df_modif = df[, cols_to_expand, drop = FALSE]
  136. df_const = df[,-cols_to_expand, drop = FALSE]
  137. orig_colnames_modif = colnames(df_modif)
  139. # Running expansion
  140. df_list = list()
  141. for (i in 1:nrow(df_const)){
  142. print(i)
  143. curr_df_const = df_const[i,, drop = FALSE]
  144. curr_df_modif = df_modif[i,, drop = FALSE]
  145. curr_df_modif = apply(curr_df_modif, 2, function(x) unlist(stri_split_fixed(x, pattern = pattern)), simplify = FALSE)
  147. if (length(cols_to_expand) > 1){
  148. curr_df_modif = do.call(cbind, curr_df_modif)
  149. } else {
  150. curr_df_modif = unlist(curr_df_modif)
  151. }
  153. if (is.matrix(curr_df_modif)){
  154. for (b in 1:nrow(curr_df_modif)){
  155. curr_df_const[b, ] = curr_df_const[1, ]
  156. }
  157. } else {
  158. for (b in 1:length(curr_df_modif)){
  159. curr_df_const[b, ] = curr_df_const[1, ]
  160. }
  161. }
  163. curr_df_combined = cbind(curr_df_const, curr_df_modif)
  165. if (length(cols_to_expand) == 1){
  166. colnames(curr_df_combined)[ncol(curr_df_combined)] = orig_colnames[cols_to_expand]
  167. }
  168. df_list[[i]] = curr_df_combined
  169. }
  171. if (length(df_list) > 1){
  172. df_list = do.call(rbind, df_list)
  173. } else {
  174. df_list = df_list[[1]]
  175. }
  177. df_list = df_list[, orig_colnames]
  178. return(df_list)
  179. }
  181. # A function to check gene symbols
  182. # Importing NIH dataset
  183. Homo_Sapiens_Gene_info_NIH = smart_fread("/home/aleksandr/Desktop/WORK/UCSC_ID_MAP/Homo_sapiens.gene_info")
  184. # https://ftp.ncbi.nih.gov/gene/DATA/GENE_INFO/Mammalia/ (replace with an appropriate path)
  185. Homo_Sapiens_Gene_info_NIH_expanded = multiple_expander(df = Homo_Sapiens_Gene_info_NIH, cols_to_expand = 5, pattern = "|")
  186. check_gene_symbol_NIH = function(PRF_gene_symbols, PRF_ref_NIH_expanded, PRF_replace_NA_with_old = FALSE){
  187. PRF_gene_symbols_check = lapply(PRF_gene_symbols, function(x){
  188. if (x %in% PRF_ref_NIH_expanded$Symbol_from_nomenclature_authority){
  189. Curr_gene = x
  190. Approved = TRUE
  191. Suggested.Symbol = x
  192. } else if (x %in% PRF_ref_NIH_expanded$Symbol){
  193. RRF_df = PRF_ref_NIH_expanded[PRF_ref_NIH_expanded$Symbol == x,]
  194. Curr_gene = x
  195. Approved = FALSE
  196. Suggested.Symbol = RRF_df[,"Symbol_from_nomenclature_authority"]
  197. Suggested.Symbol = unique(Suggested.Symbol)
  198. Suggested.Symbol = Suggested.Symbol[Suggested.Symbol != "-"]
  199. if (length(Suggested.Symbol) >= 1){
  200. Suggested.Symbol = Suggested.Symbol[1]
  201. if (Suggested.Symbol == x){
  202. Approved = TRUE
  203. }
  204. } else {
  205. Suggested.Symbol = RRF_df[,"Symbol"]
  206. Suggested.Symbol = unique(Suggested.Symbol)
  207. Suggested.Symbol = Suggested.Symbol[1]
  208. }
  209. } else if (x %in% PRF_ref_NIH_expanded$Synonyms){
  210. RRF_df = PRF_ref_NIH_expanded[PRF_ref_NIH_expanded$Synonyms == x,]
  211. Curr_gene = x
  212. Approved = FALSE
  213. Suggested.Symbol = RRF_df[,"Symbol_from_nomenclature_authority"]
  214. Suggested.Symbol = unique(Suggested.Symbol)
  215. Suggested.Symbol = Suggested.Symbol[Suggested.Symbol != "-"]
  216. if (length(Suggested.Symbol) >= 1){
  217. Suggested.Symbol = Suggested.Symbol[1]
  218. } else {
  219. Suggested.Symbol = RRF_df[,"Symbol"]
  220. Suggested.Symbol = unique(Suggested.Symbol)
  221. Suggested.Symbol = Suggested.Symbol[1]
  222. }
  223. } else {
  224. Curr_gene = x
  225. Approved = FALSE
  226. Suggested.Symbol = NA
  227. if (PRF_replace_NA_with_old){
  228. Suggested.Symbol = x
  229. }
  230. }
  231. Dataset = data.frame(x = Curr_gene, Approved = Approved, Suggested.Symbol = Suggested.Symbol)
  232. })
  233. PRF_gene_symbols_check = list_to_df(PRF_gene_symbols_check)
  234. return(PRF_gene_symbols_check)
  235. }
  237. test_genes_OLINK_linear_models_generalized = function(PREFIX_Pheno_df,
  238. PREFIX_Expression_matrix,
  239. PREFIX_contrast_col_number,
  240. PREFIX_participants_col_number,
  241. PREFIX_contrast_vector,
  242. PREFIX_model.covariates,
  243. PREFIX_Remove_NA_predictors = TRUE,
  244. PREFIX_show_all_coef = FALSE,
  245. PREFIX_round_output = TRUE){
  247. writeLines("Testing all genes")
  249. #Preparing data
  250. colnames(PREFIX_Pheno_df)[PREFIX_contrast_col_number] = "Case_Control"
  251. PREFIX_Pheno_df = PREFIX_Pheno_df[PREFIX_Pheno_df$Case_Control %in% PREFIX_contrast_vector, ]
  252. PREFIX_Pheno_df$Case_Control = factor(PREFIX_Pheno_df$Case_Control, levels = c(PREFIX_contrast_vector[1], PREFIX_contrast_vector[2]))
  254. colnames(PREFIX_Pheno_df)[PREFIX_participants_col_number] = "Participant_ID_FUN"
  256. # getting expression values
  257. PREFIX_Tested_expression = PREFIX_Expression_matrix
  258. PREFIX_Tested_expression = PREFIX_Tested_expression[,colnames(PREFIX_Tested_expression) %in% PREFIX_Pheno_df$Participant_ID_FUN]
  260. # preparing pheno df before model
  261. if (!is.null(PREFIX_model.covariates)){
  262. PREFIX_Pheno_df = PREFIX_Pheno_df[,c("Participant_ID_FUN", "Case_Control", PREFIX_model.covariates)]
  263. } else {
  264. PREFIX_Pheno_df = PREFIX_Pheno_df[,c("Participant_ID_FUN", "Case_Control")]
  265. }
  266. PREFIX_MODEL_LIST = list()
  268. for (PREFIX_i in 1:nrow(PREFIX_Tested_expression)){
  269. PREFIX_Curr_Test = rownames(PREFIX_Tested_expression)[PREFIX_i]
  270. PREFIX_Curr_values = PREFIX_Tested_expression[PREFIX_Curr_Test,]
  271. PREFIX_tmp_df = cbind(PREFIX_Pheno_df, t(PREFIX_Curr_values))
  273. if (PREFIX_Remove_NA_predictors){
  274. writeLines("Removing NA predictors")
  275. PREFIX_Participants_missing_data = apply(PREFIX_tmp_df, 1, function(x){
  276. x = as.character(x)
  277. if (any(is.na(x))){
  278. return(TRUE)
  279. } else{
  280. return(FALSE)
  281. }
  282. })
  283. PREFIX_Participants_missing_data = PREFIX_tmp_df$Participant_ID_FUN[PREFIX_Participants_missing_data]
  284. }
  286. PREFIX_tmp_df = PREFIX_tmp_df[PREFIX_tmp_df$Participant_ID_FUN %!in% PREFIX_Participants_missing_data,]
  288. if (!is.null(PREFIX_model.covariates)){
  289. PREFIX_Model.formula.string = paste0(PREFIX_model.covariates, collapse = " + ")
  290. PREFIX_Model.formula.string = paste0(PREFIX_Curr_Test, " ~ ", "Case_Control", " + ", PREFIX_Model.formula.string)
  291. } else {
  292. PREFIX_Model.formula.string = paste0(PREFIX_Curr_Test, " ~ ", "Case_Control")
  293. }
  295. PREFIX_model_fit = lm(formula = as.formula(PREFIX_Model.formula.string), data = PREFIX_tmp_df)
  296. if (PREFIX_show_all_coef){
  297. PREFIX_model_output = modify_coef_log_model(PREFIX_model_fit, model = "lin.model", rounding = FALSE)
  298. } else {
  299. PREFIX_model_output = modify_coef_log_model(PREFIX_model_fit, model = "lin.model", rounding = FALSE)
  300. PREFIX_model_output = PREFIX_model_output[2,]
  301. }
  302. PREFIX_model_output$Model.formula.string = PREFIX_Model.formula.string
  303. PREFIX_model_output$Gene = PREFIX_Curr_Test
  304. PREFIX_model_output$Contrast = paste0(PREFIX_contrast_vector, collapse = "/")
  305. PREFIX_model_output$Excluded_particip = paste0(PREFIX_Participants_missing_data, collapse = ";")
  306. PREFIX_model_output$Excluded_particip_number = length(PREFIX_Participants_missing_data)
  307. PREFIX_model_output$Df = PREFIX_model_fit$df.residual
  309. # Calculate logFC
  310. PREFIX_expression_values_controls = PREFIX_tmp_df[PREFIX_tmp_df$Case_Control == PREFIX_contrast_vector[1], ncol(PREFIX_tmp_df)]
  311. PREFIX_expression_values_cases = PREFIX_tmp_df[PREFIX_tmp_df$Case_Control == PREFIX_contrast_vector[2], ncol(PREFIX_tmp_df)]
  312. PREFIX_expression_values_controls_exp = 2^PREFIX_expression_values_controls
  313. PREFIX_expression_values_cases_exp = 2^PREFIX_expression_values_cases
  314. Mean_1 = mean(PREFIX_expression_values_controls_exp)
  315. Mean_2 = mean(PREFIX_expression_values_cases_exp)
  316. FC = Mean_2/Mean_1
  317. Log_FC = log2(FC)
  318. PREFIX_model_output$logFC = Log_FC
  320. PREFIX_model_output = PREFIX_model_output[,c("Model.formula.string", "Contrast", "Gene", "Excluded_particip", "Excluded_particip_number",
  321. "Coef.", "logFC", "β","95% CI β","Std. Error β","T.value","P.value", "Df")]
  322. PREFIX_MODEL_LIST[[PREFIX_i]] = PREFIX_model_output
  323. }
  325. if (length(PREFIX_MODEL_LIST) > 1){
  326. PREFIX_MODEL_LIST = do.call(rbind, PREFIX_MODEL_LIST)
  327. } else {
  328. PREFIX_MODEL_LIST = PREFIX_MODEL_LIST[[1]]
  329. }
  331. if (PREFIX_show_all_coef == FALSE){
  332. PREFIX_MODEL_LIST$Adj.P.val = p.adjust(p = PREFIX_MODEL_LIST$P.value, method = "bonferroni")
  333. PREFIX_MODEL_LIST = dplyr::arrange(PREFIX_MODEL_LIST, P.value)
  334. } else {
  335. PREFIX_MODEL_LIST$Adj.P.val = NA
  336. }
  337. colnames(PREFIX_MODEL_LIST) = stri_replace_all_fixed(colnames(PREFIX_MODEL_LIST), pattern = "PREFIX_", replacement = "")
  339. if (PREFIX_round_output){
  340. writeLines("Rounding Output")
  341. PREFIX_MODEL_LIST$β = as.numeric(PREFIX_MODEL_LIST$β)
  342. PREFIX_MODEL_LIST$`Std. Error β` = as.numeric(PREFIX_MODEL_LIST$`Std. Error β`)
  343. PREFIX_MODEL_LIST$P.value = as.numeric(PREFIX_MODEL_LIST$P.value)
  344. PREFIX_MODEL_LIST$`95% CI β` = sapply(PREFIX_MODEL_LIST$`95% CI β`, function(x){
  345. x = stri_split_fixed(x, pattern = " \u2013 ")
  346. x = unlist(x)
  347. x = x[x!=""]
  348. x = as.numeric(x)
  349. x = round(x, digits = 3)
  350. x = paste0(x, collapse = " \u2013 ")
  351. })
  352. PREFIX_MODEL_LIST = PREFIX_MODEL_LIST %>%
  353. mutate_if(is.numeric, function(x) (round(x, digits = 3)))
  354. }
  355. return(PREFIX_MODEL_LIST)
  356. }
  359. modify_coef_log_model = function(x, digits = 3, model = "log_regr", rounding = TRUE){
  360. if (model == "log_regr"){
  361. Summary_model = summary(x)[["coefficients"]]
  362. if (rounding){
  363. Summary_model[,-ncol(Summary_model)] = apply(Summary_model[,-ncol(Summary_model)], 2, function(x) round(x, digits = digits))
  364. }
  365. #Summary_model[,ncol(Summary_model)] = round(Summary_model[,ncol(Summary_model)], digits = 9)
  366. Conf_inter = confint(x, level = 0.95)
  367. if (rounding){
  368. Conf_inter = apply(Conf_inter, 2, function(x) round(x, digits = digits))
  369. }
  370. Conf_inter = mapply(Conf_inter[,1], Conf_inter[,2], FUN = function(x,y){paste0(x, " \u2013 " ,y)})
  371. Exponent_beta_coef = round(exp(Summary_model[,1]), digits = digits)
  372. Table_output = cbind(rownames(Summary_model), Summary_model[,1], Conf_inter, Exponent_beta_coef, Summary_model[,-1])
  373. colnames(Table_output) = c("Coef.", "\U03B2", "95% CI \U03B2", "exp(\U03B2)", "Std. Error \U03B2","Z.value", "P.value")
  374. rownames(Table_output) = NULL
  375. Table_output = as.data.frame(Table_output)
  376. } else {
  377. Summary_model = summary(x)[["coefficients"]]
  378. if (rounding){
  379. Summary_model[,-ncol(Summary_model)] = apply(Summary_model[,-ncol(Summary_model)], 2, function(x) round(x, digits = digits))
  380. }
  381. #Summary_model[,ncol(Summary_model)] = round(Summary_model[,ncol(Summary_model)], digits = 9)
  382. Conf_inter = confint(x, level = 0.95)
  383. if (rounding){
  384. Conf_inter = apply(Conf_inter, 2, function(x) round(x, digits = digits))
  385. }
  386. Conf_inter = mapply(Conf_inter[,1], Conf_inter[,2], FUN = function(x,y){paste0(x, " \u2013 " ,y)})
  387. Table_output = cbind(rownames(Summary_model), Summary_model[,1], Conf_inter, Summary_model[,-1])
  388. colnames(Table_output) = c("Coef.", "\U03B2", "95% CI \U03B2", "Std. Error \U03B2", "T.value", "P.value")
  389. rownames(Table_output) = NULL
  390. Table_output = as.data.frame(Table_output)
  391. }
  392. return(Table_output)
  393. }
  395. characterize_categorical_variable = function(df, Category_1, Category_2, Variable_name, P.val.valid = TRUE, keep_missing = TRUE){
  396. Total_variable = df[, Variable_name]
  398. if (!is.ordered(Variable_name)){
  399. df[, Variable_name] = ordered(Total_variable, levels = sort(unique(Total_variable)))
  400. Total_variable = df[, Variable_name]
  401. }
  403. Contingency_table = table(df[,Variable_name], df[,"variable_to_split"])
  405. if (keep_missing){
  406. Contingency_table = table(df[,Variable_name], df[,"variable_to_split"], exclude = NULL)
  407. }
  409. Contingency_table_percents_1 = (Contingency_table[,Category_1]/(sum(Contingency_table[,Category_1])) * 100) %>% round(., digits = 1)
  410. Contingency_table_percents_2 = (Contingency_table[,Category_2]/(sum(Contingency_table[,Category_2])) * 100) %>% round(., digits = 1)
  411. Report_1 = paste0(row.names(Contingency_table), ": ")
  412. Report_1 = paste0(Report_1, Contingency_table[,Category_1], " (",Contingency_table_percents_1, "%)", collapse = "\n")
  413. Report_2 = paste0(row.names(Contingency_table), ": ")
  414. Report_2 = paste0(Report_2, Contingency_table[,Category_2], " (",Contingency_table_percents_2, "%)", collapse = "\n")
  416. Contingency_table_pval = table(df[,Variable_name], df[,"variable_to_split"])
  417. if (P.val.valid){
  418. tryCatch({
  419. Var.pval <<- chisq.test(Contingency_table_pval)
  420. Var.pval <<- Var.pval$p.value}, warning = function(w){
  421. message = paste0(Variable_name, " has small counts in some groups, using P.val from Fisher's Exact Test")
  422. writeLines(message)
  423. Var.pval <<- fisher.test(Contingency_table_pval)
  424. Var.pval <<- Var.pval$p.value
  425. })
  426. } else {
  427. Var.pval = "Not valid"
  428. }
  430. Output_list = list()
  431. Output_list[[1]] = Report_1
  432. Output_list[[2]] = Report_2
  433. Output_list[[3]] = Var.pval
  434. return(Output_list)
  435. }
  437. characterize_numeric_variable = function(df, Category_1, Category_2, Variable_name, Mention_NAs = TRUE){
  438. Total_variable = df[, Variable_name]
  439. Variable_1 = df[df$variable_to_split == Category_1, Variable_name]
  440. Variable_1_report = describe_vector_numeric(Variable_1, Mention_NAs = Mention_NAs)
  441. Variable_2 = df[df$variable_to_split == Category_2, Variable_name]
  442. Variable_2_report = describe_vector_numeric(Variable_2, Mention_NAs = Mention_NAs)
  443. Shapiro_normality_check = shapiro.test(Total_variable)
  444. Shapiro_normality_check = Shapiro_normality_check$p.value
  445. if (Shapiro_normality_check < 0.05){
  446. message = paste0(Variable_name, " is not normally distributed. Using P-val from Mann Whitney U Test")
  447. writeLines(message)
  448. formula_test = paste0(Variable_name, "~ variable_to_split")
  450. tryCatch(
  451. {
  452. Total_pval = wilcox.test(formula = as.formula(formula_test), data = df)
  453. Total_pval = Total_pval$p.value
  454. },
  455. error = function(cond) {
  456. message(conditionMessage(cond))
  457. message("Returning NA for p-value")
  458. Total_pval <<- NA
  459. }
  460. )
  462. } else {
  463. message = paste0(Variable_name, " is normally distributed. Using P-val from T Test")
  464. writeLines(message)
  465. formula_test = paste0(Variable_name, "~ variable_to_split")
  466. Var_test_check = var.test(formula = as.formula(formula_test), data = df)
  467. Var_test_check = Var_test_check$p.value
  468. if (Var_test_check >= 0.05){
  469. Total_pval = t.test(formula = as.formula(formula_test), data = df, var.equal = TRUE)
  470. } else {
  471. Total_pval = t.test(formula = as.formula(formula_test), data = df, var.equal = FALSE)
  472. }
  473. Total_pval = Total_pval$p.value
  474. }
  475. Output_list = list()
  476. Output_list[[1]] = Variable_1_report
  477. Output_list[[2]] = Variable_2_report
  478. Output_list[[3]] = Total_pval
  479. return(Output_list)
  480. }
  482. describe_vector_numeric = function(x, Mention_NAs, show_as_percent = FALSE, Describe_min_max = TRUE){
  484. if (all(is.na(x))){
  485. Missing_count = length(x[is.na(x)])
  486. Missing_val_percent = (Missing_count/length(x) * 100) %>% round(., digits = 2)
  487. Var_report = paste0("Missing val: ", Missing_count, " (", Missing_val_percent, "%)")
  488. return(Var_report)
  489. }
  491. Mean_var = mean(x, na.rm = TRUE) %>% round(., digits = 2)
  492. SD_var = sd(x, na.rm = TRUE) %>% round(., digits = 2)
  493. Var_report = paste0(Mean_var, " ± ", SD_var)
  495. if (is.na(SD_var)){
  496. Var_report = paste0("Single non-NA value")
  497. }
  499. if (show_as_percent & Var_report != "Single non-NA value"){
  500. Mean_var = Mean_var * 100/Mean_var
  501. SD_var = SD_var * 100/Mean_var
  502. Var_report = paste0(Mean_var, " ± ", SD_var, "%")
  503. }
  505. if (Mention_NAs){
  506. Missing_count = length(x[is.na(x)])
  507. Missing_val_percent = (Missing_count/length(x) * 100) %>% round(., digits = 2)
  508. if (Missing_count > 0){
  509. Var_report = paste0(Var_report, "\nMissing val: ", Missing_count, " (", Missing_val_percent, "%)")
  510. }
  511. }
  514. if (Describe_min_max){
  515. Min_var = min(x, na.rm = TRUE) %>% round(., digits = 2)
  516. Max_var = max(x, na.rm = TRUE) %>% round(., digits = 2)
  517. if (show_as_percent & !is.na(SD_var)){
  518. Min_var = Min_var * 100/Mean_var
  519. Max_var = Max_var * 100/Mean_var
  520. Var_report = paste0(Var_report, "\n", "Min: ", Min_var, "%, Max: ",Max_var, "%")
  521. } else {
  522. Var_report = paste0(Var_report, "\n", "Min: ", Min_var, ", Max: ",Max_var)
  523. }
  524. }
  525. return(Var_report)
  526. }
  528. characterize_dataset_generelized_two_subgroups = function(dataset,
  529. study_char,
  530. contrast_col_number,
  531. participants_col_number,
  532. model_covariates_col_vector,
  533. columns_to_characterise_vector,
  534. Remove_NA_predictors = FALSE,
  535. drop_P = TRUE,
  536. simplif_P = 3){
  538. # Running Raw Statistics
  539. total_number = nrow(dataset)
  540. total_number_message = paste0("Initial dataset includes ", total_number, " participants")
  541. writeLines(total_number_message)
  542. first_row_output = data.frame(Name = total_number_message,Category.1 = "", Category.2 = "", P.value = NA)
  544. # Removing missing data
  545. colnames(dataset)[participants_col_number] = "Participant_ID_FUN"
  546. indeces = c(participants_col_number, contrast_col_number, model_covariates_col_vector)
  547. indeces = unique(indeces)
  548. Model_df = dataset[,indeces]
  549. if (Remove_NA_predictors){
  551. writeLines("Removing participants with missing data")
  552. Participants_missing_data = apply(Model_df, 1, function(x){
  553. x = as.character(x)
  555. if (any(is.na(x))){
  556. return(TRUE)
  557. } else{
  558. return(FALSE)
  559. }
  561. })
  562. Participants_missing_data = Model_df$Participant_ID_FUN[Participants_missing_data]
  564. if (length(Participants_missing_data) < 1){
  565. excluded_number = 0
  567. } else {
  569. excluded_number = length(Participants_missing_data)
  570. }
  572. Participants_excluded_df = dataset[dataset$Participant_ID_FUN %in% Participants_missing_data, ]
  573. dataset = dataset[dataset$Participant_ID_FUN %!in% Participants_missing_data, ]
  574. excluded_message = paste0("Participants with missing data excluded: ", excluded_number,"\n",
  575. "Resulting number of participants: ", nrow(dataset))
  576. second_row_output = data.frame(Name = excluded_message,Category.1 = "", Category.2 = "", P.value = NA)
  577. writeLines(excluded_message)
  579. } else {
  581. writeLines("Exclusion of participants was not performed")
  582. second_row_output = data.frame(Name = "Exclusion of participants was not performed", Category.1 = "", Category.2 = "", P.value = NA)
  583. Participants_excluded_df = NA
  584. }
  586. main_indeces = c(contrast_col_number, columns_to_characterise_vector)
  587. main_indeces = unique(main_indeces)
  588. Main_dataset = list()
  591. # Making categories
  592. Initial_names = colnames(dataset)[main_indeces]
  593. colnames(dataset)[contrast_col_number] = "variable_to_split"
  594. levels_var = dataset[,"variable_to_split"]
  595. levels_var = levels(levels_var)
  596. names(levels_var) = NULL
  597. Category_1 = levels_var[1]
  598. Category_2 = levels_var[2]
  600. for (i in 1:length(main_indeces)){
  601. curr_variable = dataset[, main_indeces[i]]
  602. if (is.factor(curr_variable)){
  603. Characterised_var = characterize_categorical_variable(df = dataset,
  604. Category_1 = Category_1,
  605. Category_2 = Category_2,
  606. Variable_name = colnames(dataset)[main_indeces[i]],
  607. keep_missing = TRUE)
  608. } else if (is.character(curr_variable)){
  609. message = paste0("Bad index: ", main_indeces[i])
  610. print(message)
  611. stop("All columns should be either a Factor or Numeric")
  612. } else if (is.numeric(curr_variable)){
  613. Characterised_var = characterize_numeric_variable(df = dataset,
  614. Category_1 = Category_1,
  615. Category_2 = Category_2,
  616. Variable_name = colnames(dataset)[main_indeces[i]],
  617. Mention_NAs = TRUE)
  618. }
  619. Current_DF = data.frame(Name = Initial_names[i],
  620. Category.1 = Characterised_var[[1]],
  621. Category.2 = Characterised_var[[2]],
  622. P.value = Characterised_var[[3]])
  623. Main_dataset[[i]] = Current_DF
  624. }
  625. Main_dataset = list_to_df(data_list = Main_dataset)
  626. Main_dataset$P.value = round(Main_dataset$P.value, digits = simplif_P)
  628. # Compiling full dataframe
  629. header = data.frame(Name = study_char, Category.1 = "", Category.2 = "", P.value = NA)
  630. Full_table = rbind(header, first_row_output, second_row_output, Main_dataset)
  632. if(drop_P){
  633. Full_table$P.value = NULL
  634. }
  636. output = list()
  637. output$Table = Full_table
  638. output$Excluded = Participants_excluded_df
  639. return(output)
  640. }
  643. # SCAS score calculator
  644. SCAS_calculator = function(df_SCAS_ordered, NA_param = FALSE){
  646. sum_function = function(x){
  647. x = sum(x, na.rm = NA_param)
  648. return(x)
  649. }
  651. Total_score = apply(df_SCAS_ordered[,-c(11, 17, 26, 31, 38, 43)], MARGIN = 1, sum_function)
  653. Social_phobia = apply(df_SCAS_ordered[,c(6, 7, 9, 10, 29, 35)], MARGIN = 1, sum_function)
  655. GAD = apply(df_SCAS_ordered[,c(1, 3, 4, 20, 22, 24)], MARGIN = 1, sum_function)
  657. FOPH = apply(df_SCAS_ordered[,c(2, 18, 23, 25, 33)], MARGIN = 1, sum_function)
  659. Panic_agoraphobia = apply(df_SCAS_ordered[,c(13, 21, 28, 30, 32, 34, 36, 37, 39)], MARGIN = 1, sum_function)
  661. OCD = apply(df_SCAS_ordered[,c(14, 19, 27, 40, 41, 42)], MARGIN = 1, sum_function)
  663. SAD = apply(df_SCAS_ordered[,c(5, 8, 12, 15, 16, 44)], MARGIN = 1, sum_function)
  665. Score_df = cbind(Total_score, Social_phobia, GAD, FOPH, Panic_agoraphobia, OCD, SAD)
  667. colnames(Score_df) = paste0("SCAS_", colnames(Score_df))
  669. Score_df = as.data.frame(Score_df)
  671. return(Score_df)
  673. }
  675. # plot forest for genes
  676. plot_forest_meta_gene = function(gene_name, meta_df){
  678. tmp_df_genes = meta_df[meta_df$Corrected_symbol == gene_name,]
  680. if(nrow(tmp_df_genes)<1){
  681. writeLines("Gene is missing in the data")
  682. return(NULL)
  683. }
  685. tmp_meta_model = rma.uni(yi = avgLog2FC, vi = maxSE^2, data = tmp_df_genes, method = "SJ", weighted = TRUE)
  687. # Example Visualization
  688. weights = fmtx(weights(tmp_meta_model), digits=1)
  689. sav = forest(tmp_meta_model, slab = Study, ilab = weights)
  690. k = nrow(tmp_df_genes)
  691. colp <- "red"
  692. segments(coef(tmp_meta_model), 0, coef(tmp_meta_model), k, col=colp, lty="33", lwd=0.8)
  694. # Add text
  695. par(xpd=NA)
  696. par(cex=sav$cex, font=2)
  697. difference_coords = sav$xlim[2] - sav$xlim[1]
  698. # Headers
  699. text(sav$xlim[1], k+2.5, pos=4, "Cohort")
  700. text(sav$xlim[1]+difference_coords*0.25, k+2.5, pos=4, "Weight %")
  701. text(0, k+2.7, "Log2FC,\n(95% CI)")
  702. segments(sav$ilab.xpos[1]-0.22, k+2.5, sav$ilab.xpos[2]+0.13, k+2.8)
  703. text(sav$xlim[2]-0.10, k+2.7, "Log2FC\n(95% CI)")
  705. # Use a non-bold font for the rest of the text
  706. # Y coord for heterogen
  707. y_low = sav$ylim[1]
  708. y_high = sav$ylim[2]
  709. par(cex=sav$cex, font=1)
  710. text(sav$ilab.xpos[3], 0, "100.0")
  711. text(sav$xlim[1], y_low, pos=4, bquote(paste("Test for heterogeneity: ",
  712. I^2, "=", .(fmtx(tmp_meta_model$I2, digits=2)), ", ",
  713. tau^2, "=", .(fmtx(tmp_meta_model$tau2, digits=2)), ", ",
  714. chi^2, "=", .(fmtx(tmp_meta_model$QE, digits=2)),
  715. ", df=", .(tmp_meta_model$k - tmp_meta_model$p), ", ",
  716. .(fmtp(tmp_meta_model$QEp, digits=2, pname="P", add0=TRUE, equal=TRUE)))))
  717. title(paste0("Gene: ", gene_name))
  718. }
  720. plot_forest_simple_gene = function(gene_name, meta_df, study_column_char = "Study"){
  721. tmp_df_genes = meta_df[meta_df$Corrected_symbol == gene_name,]
  722. if(nrow(tmp_df_genes)<1){
  723. writeLines("Gene is missing in the data")
  724. return(NULL)
  725. }
  726. forest(tmp_df_genes$avgLog2FC, tmp_df_genes$maxSE**2, slab = tmp_df_genes[,study_column_char], xlab = paste0("Gene: ", gene_name))
  727. }
  729. # A function to make Venn diagrams from named lists (uses ggplot2 and ggVennDiagram)
  730. make_Venn_digram_list = function(named_list, plot_full_path = NULL, label_text_size=4, ...){
  732. # Customizable Venn diagram
  733. venn = Venn(named_list)
  734. data = process_data(venn)
  735. data@region$full_lable = sapply(data@region$count, function(x){
  736. Number = x
  737. Percent = x/sum(data@region$count)
  738. Percent = Percent*100
  739. Percent = round(Percent, digits = 1)
  740. Label_full = paste0(Number, "\n","(", Percent, "%)")
  741. return(Label_full)
  742. })
  744. # to see available shapes plot_shapes()
  745. plot = ggplot() +
  747. # 1. region count layer
  748. geom_sf(aes(fill = id), data = venn_region(data), ...) +
  750. # 2. set edge layer
  751. geom_sf(color="black", size = 0.5, data = venn_setedge(data), show.legend = FALSE, ...) +
  753. # 3. set label layer
  754. geom_sf_text(aes(label = name), data = venn_setlabel(data), size = label_text_size, ...) +
  756. # 4. region label layer
  757. geom_sf_label(aes(label = full_lable), data = venn_region(data), size = label_text_size, ...) +
  758. scale_fill_brewer(...) +
  759. scale_x_continuous(expand = expansion(mult = .2)) +
  760. theme(
  761. axis.title = element_blank(),
  762. axis.text = element_blank(),
  763. axis.ticks = element_blank(),
  764. axis.line = element_blank(),
  765. legend.background = element_blank(),
  766. legend.key = element_blank(),
  767. legend.position = "none",
  768. panel.grid.major = element_blank(),
  769. panel.grid.minor = element_blank(),
  770. plot.background = element_blank(),
  771. panel.background = element_blank(),
  772. plot.margin = margin(1,1,1,1, "cm")
  773. )
  774. print(plot)
  776. if (is.character(plot_full_path)){
  778. if (stri_detect_fixed(plot_full_path, ".pdf")){
  779. pdf(plot_full_path, width = 10, height = 10)
  780. print(plot)
  781. dev.off()
  782. }
  784. if (stri_detect_fixed(plot_full_path, ".png")){
  785. png(plot_full_path, width = 1024, height = 1024)
  786. print(plot)
  787. dev.off()
  788. }
  789. }
  790. }
  792. run_enrichment_GO_KEGG_gene_set = function(genes, universe, categories_to_show = 30, folder, plot_name_pref){
  794. # genes should be submitted as Entrez IDs
  795. if (length(genes) < 10) {
  796. Minsize = length(genes)
  797. } else {
  798. Minsize = 10
  799. }
  801. # IDs should be characters
  802. genes = as.character(genes)
  803. universe = as.character(universe)
  805. # Running GO analysis
  806. GO_String = c("BP", "MF", "CC")
  807. GO_Enrichment = lapply(GO_String, function(x){
  808. result = clusterProfiler::enrichGO(gene = genes, OrgDb = 'org.Hs.eg.db', ont = x, universe = universe, minGSSize = Minsize)
  809. return(result)
  810. })
  811. KEGG_Enrichment = clusterProfiler::enrichKEGG(gene = genes, organism = 'hsa' , universe = universe, minGSSize = Minsize)
  813. # Saving GO results
  814. dir.create(folder)
  815. combined_enrich = list()
  816. for (i in 1:length(GO_String)){
  817. filename = paste0(plot_name_pref, "_GO_", GO_String[i], ".xlsx")
  818. filename = paste0(folder, "/", filename)
  819. Curr_GO = GO_Enrichment[[i]]
  820. Curr_GO_readable = setReadable(Curr_GO, 'org.Hs.eg.db', 'ENTREZID')
  821. Curr_GO_result = Curr_GO@result
  822. Curr_GO_result$geneSymbol = Curr_GO_readable@result$geneID
  823. combined_enrich[[i]] = Curr_GO_result
  824. openxlsx::write.xlsx(x = Curr_GO_result, file = filename, overwrite = TRUE)
  825. }
  827. # Saving KEGG results
  828. if (!is.null(KEGG_Enrichment)){
  829. KEGG_Enrichment_readable = setReadable(KEGG_Enrichment, 'org.Hs.eg.db', 'ENTREZID')
  830. KEGG_Enrichment_result = KEGG_Enrichment@result
  831. KEGG_Enrichment_result$geneSymbol = KEGG_Enrichment_readable@result$geneID
  832. combined_enrich[[4]] = KEGG_Enrichment_result
  833. openxlsx::write.xlsx(x = KEGG_Enrichment_result, file = paste0(folder, "/", plot_name_pref, "_KEGG.xlsx"), overwrite = TRUE)
  834. } else {
  835. combined_enrich[[4]] = NA
  836. }
  837. # Saving full results
  838. openxlsx::write.xlsx(x = combined_enrich, file = paste0(folder, "/", plot_name_pref, "_enrichment_combined_sheet.xlsx"), overwrite = TRUE)
  840. # Making plots
  841. GO_Enrichment_readable = lapply(GO_Enrichment, function(x) setReadable(x, 'org.Hs.eg.db', 'ENTREZID'))
  843. if (!is.null(KEGG_Enrichment)){
  844. GO_Enrichment_readable[[4]] = KEGG_Enrichment_readable
  845. } else {
  846. GO_Enrichment_readable[[4]] = NA
  847. }
  849. GO_String = c("BP", "MF", "CC", "KEGG")
  850. gc()
  852. for (i in 1:length(GO_String)){
  853. Curr_GO = GO_Enrichment_readable[[i]]
  855. if (is.na(Curr_GO)){
  856. writeLines(paste0(filename, " has NO results and analysis was not conducted"))
  857. next
  858. }
  860. filename = paste0(plot_name_pref, "_GO_", GO_String[i])
  861. filename = paste0(folder, "/", filename)
  863. if (GO_String[i] == "KEGG"){
  864. # Figures for KEGG
  865. dimensions = c(19*nrow(summary(Curr_GO))/categories_to_show, 11*nrow(summary(Curr_GO))/categories_to_show)
  867. if (any(dimensions < 1)){
  868. dimensions = dimensions*10
  869. }
  871. if (nrow(summary(Curr_GO)) < 1){
  872. writeLines(paste0(filename, " has NO significant results"))
  873. next
  874. }
  876. # Barplot
  877. pdf(paste0(filename, "_bar.pdf")
  878. , width = 10, height = 10)
  879. print(barplot(Curr_GO, showCategory = categories_to_show, font.size = 10))
  880. dev.off()
  882. pdf(paste0(filename, "_bar_full.pdf")
  883. , width = 10, height = 10*nrow(summary(Curr_GO))/categories_to_show)
  884. print(barplot(Curr_GO, showCategory = nrow(summary(Curr_GO)), font.size = 10))
  885. dev.off()
  887. # Dotplot
  888. pdf(paste0(filename, "_dot.pdf")
  889. , width = 10, height = 10)
  890. print(dotplot(Curr_GO, showCategory= categories_to_show, font.size = 10))
  891. dev.off()
  893. pdf(paste0(filename, "_dot_full.pdf")
  894. , width = 10, height = dimensions[2])
  895. print(dotplot(Curr_GO, showCategory = nrow(summary(Curr_GO)), font.size = 10))
  896. dev.off()
  898. # Gene-Concept Network
  899. pdf(paste0(filename, "_gene_conc.pdf")
  900. , width = 19, height = 11)
  901. print(cnetplot(Curr_GO, circular = FALSE, colorEdge = TRUE, color_category='firebrick', color_gene='steelblue', layout = "kk",
  902. showCategory = categories_to_show, cex_gene = 0.7, cex_label_gene = 0.7, shadowtext = "gene", max.overlaps = 1000, force = 3, force_pull = 0.5, max.time = 2))
  903. dev.off()
  905. tryCatch({
  906. pdf(paste0(filename, "_gene_conc_full.pdf")
  907. , width = dimensions[1], height = dimensions[2])
  908. print(cnetplot(Curr_GO, circular = FALSE, colorEdge = TRUE, color_category='firebrick', color_gene='steelblue', layout = "kk",
  909. showCategory = nrow(summary(Curr_GO)), cex_gene = 0.7, cex_label_gene = 0.7, shadowtext = "gene", max.overlaps = 1000, force = 3, force_pull = 0.5, max.time = 2))
  910. dev.off()
  911. }, error = function(e) {writeLines(paste0("Full gene conc plot is not available for ", filename))})
  913. } else {
  914. # Figures for GO
  915. dimensions = c(19*nrow(summary(Curr_GO))/categories_to_show, 11*nrow(summary(Curr_GO))/categories_to_show)
  917. if (any(dimensions < 1)){
  918. dimensions = dimensions*10
  919. }
  921. if (nrow(summary(Curr_GO)) < 1){
  922. writeLines(paste0(filename, " has NO significant results"))
  923. next
  924. }
  926. # Barplot
  927. pdf(paste0(filename, "_bar.pdf")
  928. , width = 10, height = 10)
  929. print(barplot(Curr_GO, showCategory = categories_to_show, font.size = 10))
  930. dev.off()
  932. pdf(paste0(filename, "_bar_full.pdf")
  933. , width = 10, height = dimensions[2])
  934. print(barplot(Curr_GO, showCategory = nrow(summary(Curr_GO)), font.size = 10))
  935. dev.off()
  937. # Dotplot
  938. pdf(paste0(filename, "_dot.pdf")
  939. , width = 10, height = 10)
  940. print(dotplot(Curr_GO, showCategory= categories_to_show, font.size = 10))
  941. dev.off()
  943. pdf(paste0(filename, "_dot_full.pdf")
  944. , width = 10, height = dimensions[2])
  945. print(dotplot(Curr_GO, showCategory = nrow(summary(Curr_GO)), font.size = 10))
  946. dev.off()
  948. # Gene-Concept Network
  949. pdf(paste0(filename, "_gene_conc.pdf")
  950. , width = 19, height = 11)
  951. print(cnetplot(Curr_GO, circular = FALSE, colorEdge = TRUE, color_category='firebrick', color_gene='steelblue', layout = "kk",
  952. showCategory = categories_to_show, cex_gene = 0.7, cex_label_gene = 0.7, shadowtext = "gene", max.overlaps = 1000, force = 3, force_pull = 0.5, max.time = 2))
  953. dev.off()
  955. tryCatch({
  956. pdf(paste0(filename, "_gene_conc_full.pdf"),
  957. width = dimensions[1], height = dimensions[2])
  958. print(cnetplot(Curr_GO, circular = FALSE, colorEdge = TRUE, color_category='firebrick', color_gene='steelblue', layout = "kk",
  959. showCategory = nrow(summary(Curr_GO)), cex_gene = 0.7, cex_label_gene = 0.7, shadowtext = "gene", max.overlaps = 1000, force = 3, force_pull = 0.5, max.time = 2))
  960. dev.off()
  961. }, error = function(e) {writeLines(paste0("Full gene conc plot is not available for ", filename))})
  963. # GO induced graph
  964. tryCatch({
  965. pdf(paste0(filename, "_go_graph.pdf"),
  966. width = 19, height = 11)
  967. print(goplot(Curr_GO, showCategory = categories_to_show))
  968. dev.off()
  969. }, error = function(e) {writeLines(paste0("Go graph is not available for ", filename))})
  971. tryCatch({
  972. pdf(paste0(filename, "_go_graph_full.pdf"),
  973. width = dimensions[1], height = dimensions[2])
  974. print(goplot(Curr_GO, showCategory = nrow(summary(Curr_GO))))
  975. dev.off()
  976. }, error = function(e) {writeLines(paste0("Go graph (big) is not available for ", filename))})
  978. # Tree plot
  979. Curr_GO_pairwise = enrichplot::pairwise_termsim(Curr_GO)
  980. dimensions = c(19*nrow(Curr_GO_pairwise@termsim)/categories_to_show, 11*nrow(Curr_GO_pairwise@termsim)/categories_to_show)
  981. if (any(dimensions < 1)){
  982. dimensions = dimensions*10
  983. }
  985. tryCatch({
  986. pdf(paste0(filename, "_tree.pdf"),
  987. width = 19, height = 11)
  988. print(treeplot(Curr_GO_pairwise, showCategory = categories_to_show, nCluster = 10))
  989. dev.off()
  990. }, error = function(e) {writeLines(paste0("Tree plot is not available for ", filename))})
  992. tryCatch({
  993. pdf(paste0(filename, "_tree_full.pdf"),
  994. width = dimensions[1], height = dimensions[2])
  995. print(treeplot(Curr_GO_pairwise, showCategory = nrow(Curr_GO_pairwise@termsim), nCluster = 10))
  996. dev.off()
  997. }, error = function(e) {writeLines(paste0("Tree plot (big) is not available for ", filename))})
  1000. # Enrichment map
  1001. tryCatch({
  1002. pdf(paste0(filename, "_enr_map.pdf"),
  1003. width = 19, height = 11)
  1004. print(emapplot(Curr_GO_pairwise, showCategory = categories_to_show, layout = "fr", cex_label_category = 0.9,
  1005. cex_line = 0.8, repel = TRUE))
  1006. dev.off()
  1007. }, error = function(e) {writeLines(paste0("Enrich. map is not available for ", filename))})
  1009. tryCatch({
  1010. pdf(paste0(filename, "_enr_map_full.pdf"),
  1011. width = dimensions[1], height = dimensions[2])
  1012. print(emapplot(Curr_GO_pairwise, showCategory = nrow(Curr_GO_pairwise@termsim), layout = "fr", cex_label_category = 0.9,
  1013. cex_line = 0.8, repel = TRUE))
  1014. dev.off()
  1015. }, error = function(e) {writeLines(paste0("Enrich. map is not available for ", filename))})
  1017. }
  1018. }
  1020. # Preparing outputs to load into the environment
  1021. GO_Enrichment[[4]] = KEGG_Enrichment
  1023. if (is.null(KEGG_Enrichment)){
  1024. GO_Enrichment[[4]] = NA
  1025. }
  1027. Output = list()
  1028. Output[[1]] = GO_Enrichment
  1029. Output[[2]] = GO_Enrichment_readable
  1030. Output[[3]] = combined_enrich
  1031. names(Output) = c("Enrichment", "Mapped Enrichment", "Combined Results")
  1033. for (i in 1:length(Output)){
  1034. names(Output[[i]]) = GO_String
  1035. }
  1037. return(Output)
  1038. }
  1040. run_enrichment_GO_KEGG_gene_set_only_tables = function(genes, universe, folder, plot_name_pref){
  1042. # genes should be submitted as Entrez IDs
  1043. if (length(genes) < 10) {
  1044. Minsize = length(genes)
  1045. } else {
  1046. Minsize = 10
  1047. }
  1049. # IDs should be characters
  1050. genes = as.character(genes)
  1051. universe = as.character(universe)
  1053. # Running GO analysis
  1054. GO_String = c("BP", "MF", "CC")
  1055. GO_Enrichment = lapply(GO_String, function(x){
  1056. result = clusterProfiler::enrichGO(gene = genes, OrgDb = 'org.Hs.eg.db', ont = x, universe = universe, minGSSize = Minsize)
  1057. return(result)
  1058. })
  1059. KEGG_Enrichment = clusterProfiler::enrichKEGG(gene = genes, organism = 'hsa' , universe = universe, minGSSize = Minsize)
  1061. # Saving GO results
  1062. dir.create(folder)
  1063. combined_enrich = list()
  1064. for (i in 1:length(GO_String)){
  1065. filename = paste0(plot_name_pref, "_GO_", GO_String[i], ".xlsx")
  1066. filename = paste0(folder, "/", filename)
  1067. Curr_GO = GO_Enrichment[[i]]
  1068. Curr_GO_readable = setReadable(Curr_GO, 'org.Hs.eg.db', 'ENTREZID')
  1069. Curr_GO_result = Curr_GO@result
  1070. Curr_GO_result$geneSymbol = Curr_GO_readable@result$geneID
  1071. combined_enrich[[i]] = Curr_GO_result
  1072. openxlsx::write.xlsx(x = Curr_GO_result, file = filename, overwrite = TRUE)
  1073. }
  1075. # Saving KEGG results
  1076. if (!is.null(KEGG_Enrichment)){
  1077. KEGG_Enrichment_readable = setReadable(KEGG_Enrichment, 'org.Hs.eg.db', 'ENTREZID')
  1078. KEGG_Enrichment_result = KEGG_Enrichment@result
  1079. KEGG_Enrichment_result$geneSymbol = KEGG_Enrichment_readable@result$geneID
  1080. combined_enrich[[4]] = KEGG_Enrichment_result
  1081. openxlsx::write.xlsx(x = KEGG_Enrichment_result, file = paste0(folder, "/", plot_name_pref, "_KEGG.xlsx"), overwrite = TRUE)
  1082. } else {
  1083. combined_enrich[[4]] = NA
  1084. }
  1085. # Saving full results
  1086. openxlsx::write.xlsx(x = combined_enrich, file = paste0(folder, "/", plot_name_pref, "_enrichment_combined_sheet.xlsx"), overwrite = TRUE)
  1088. GO_Enrichment_readable = lapply(GO_Enrichment, function(x) setReadable(x, 'org.Hs.eg.db', 'ENTREZID'))
  1091. # Preparing outputs to load into the environment
  1092. GO_Enrichment[[4]] = KEGG_Enrichment
  1094. if (is.null(KEGG_Enrichment)){
  1095. GO_Enrichment[[4]] = NA
  1096. }
  1098. Output = list()
  1099. Output[[1]] = GO_Enrichment
  1100. Output[[2]] = GO_Enrichment_readable
  1101. Output[[3]] = combined_enrich
  1102. names(Output) = c("Enrichment", "Mapped Enrichment", "Combined Results")
  1104. for (i in 1:length(Output)){
  1105. names(Output[[i]]) = GO_String
  1106. }
  1108. return(Output)
  1109. }
  1111. ################### PSY Recall phenotype preprocessing ###################
  1112. DAWBA_DS = read.xlsx("DAWBA PSY June 2021.xlsx", sheet = 1)
  1115. Recall_2_Pheno = list.files()
  1116. Recall_2_Pheno = Recall_2_Pheno[stri_detect_fixed(Recall_2_Pheno, "RECALL_2")]
  1117. Recall_2_Pheno = lapply(Recall_2_Pheno, function(x) read.xlsx(x, sheet = 1, startRow = 1, detectDates = TRUE))
  1118. which(colnames(Recall_2_Pheno[[1]]) == "Medication")
  1119. which(colnames(Recall_2_Pheno[[2]]) == "Medication")
  1120. names(Recall_2_Pheno) = c("Recall_2A", "Recall_2B")
  1121. Columns_Recall_2A = c(1, 4, 5, 6, 7, 8:16, 17, 18,20, 21:24, 25:33, 111:132, 168:178, 179:226, 152:167) # HADS - 152:167
  1122. Columns_Recall_2B = c(1:30, 31:52, 53:63, 64:111)
  1123. Recall_2A_Pheno = Recall_2_Pheno[[1]]
  1124. Recall_2A_Pheno = Recall_2A_Pheno[,Columns_Recall_2A]
  1125. Recall_2B_Pheno = Recall_2_Pheno[[2]]
  1126. Recall_2B_Pheno = Recall_2B_Pheno[,Columns_Recall_2B]
  1127. Recall_2B_Pheno = Recall_2B_Pheno[!is.na(Recall_2B_Pheno$Code),]
  1129. ############ Recall 2A
  1130. Recall_2A_Pheno_small = Recall_2A_Pheno
  1131. Columns_to_select_RC2A = c(1,2, 22:26, 30, 31:51, 52, 53:63, 64:108, 111, 112, 113)
  1132. Recall_2A_Pheno_small = Recall_2A_Pheno_small[,Columns_to_select_RC2A]
  1133. Recall_2A_Pheno_small$DAWBA_DEPBAND = sapply(Recall_2A_Pheno_small$DAWBA.ID, function(x){
  1134. index = which(DAWBA_DS$dawbaID == x)
  1135. output = DAWBA_DS$depband[index]
  1136. if (length(output)<1){return(NA)}
  1137. return(output)
  1138. })
  1139. SCAS_DF_RC2A = Recall_2A_Pheno_small[,42:86]
  1140. SCAS_DF_RC2A = apply(SCAS_DF_RC2A, 2, as.numeric)
  1142. # corrupt total score check
  1143. corrupt_total_score = which(apply(SCAS_DF_RC2A[,-c(11, 17, 26, 31, 38, 43)], 1, function(x) any(is.na(x))))
  1144. corrupt_total_score
  1145. # Rows 2 199 200 336 should be technically excluded
  1147. corrupt_total_score_item_count = apply(SCAS_DF_RC2A[,-c(11, 17, 26, 31, 38, 43)], 1, function(x) sum(is.na(x)))
  1148. corrupt_total_score_item_count[corrupt_total_score_item_count>0]
  1149. # All rows just miss 1 item
  1150. SCAS_scores_RC2A = SCAS_calculator(df_SCAS_ordered = SCAS_DF_RC2A, NA_param = TRUE)
  1151. SCAS_scores_RC2A[corrupt_total_score, ]
  1152. Recall_2A_Pheno_small = cbind(Recall_2A_Pheno_small, SCAS_scores_RC2A)
  1153. Recall_2A_Pheno_small$Code[2] # "PSY0021"
  1155. # item check
  1156. item_numbers = sapply(colnames(SCAS_DF_RC2A), function(x) unlist(stri_split_fixed(x, pattern = "."))[2])
  1157. names(item_numbers) = NULL
  1158. is.unsorted(as.numeric(item_numbers)) # FALSE
  1161. ############ Recall 2B
  1162. Recall_2B_Pheno_small = Recall_2B_Pheno
  1163. Columns_to_select_RC2B = c(1, 2, 22:26, 30, 31:51, 52, 53:63, 64:108, 111)
  1164. Recall_2B_Pheno_small = Recall_2B_Pheno_small[,Columns_to_select_RC2B]
  1165. Recall_2B_Pheno_small$DAWBA_DEPBAND = sapply(Recall_2B_Pheno_small$DAWBA.ID, function(x){
  1166. index = which(DAWBA_DS$dawbaID == x)
  1167. output = DAWBA_DS$depband[index]
  1168. if (length(output)<1){return(NA)}
  1169. return(output)
  1170. })
  1171. SCAS_DF_RC2B = Recall_2B_Pheno_small[,42:86]
  1172. SCAS_DF_RC2B = apply(SCAS_DF_RC2B, 2, as.numeric)
  1174. # corrupt total score check
  1175. corrupt_total_score = which(apply(SCAS_DF_RC2B[,-c(11, 17, 26, 31, 38, 43)], 1, function(x) any(is.na(x))))
  1176. corrupt_total_score
  1177. # ALL rows are OK!
  1179. SCAS_scores_RC2B = SCAS_calculator(df_SCAS_ordered = Recall_2B_Pheno_small[,42:86])
  1180. Recall_2B_Pheno_small = cbind(Recall_2B_Pheno_small, SCAS_scores_RC2B)
  1182. # item check
  1183. item_numbers = sapply(colnames(SCAS_DF_RC2B), function(x) unlist(stri_split_fixed(x, pattern = "."))[2])
  1184. names(item_numbers) = NULL
  1185. is.unsorted(as.numeric(item_numbers)) # FALSE
  1187. ############ Saving phenotypes
  1188. write.csv(Recall_2A_Pheno_small, "Recall_2A_Pheno_small.csv", row.names = FALSE)
  1189. write.csv(Recall_2B_Pheno_small, "Recall_2B_Pheno_small.csv", row.names = FALSE)
  1191. ################### PSY Recall phenotype curation ###################
  1192. PSY_PHENO_RC2A = smart_fread("Recall_2A_Pheno_small.csv")
  1193. PSY_PHENO_RC2B = smart_fread("Recall_2B_Pheno_small.csv")
  1194. intersecting_cols_RC2 = intersect(colnames(PSY_PHENO_RC2A), colnames(PSY_PHENO_RC2B))
  1195. PSY_PHENO_RC_TOTAL = rbind(PSY_PHENO_RC2A[,intersecting_cols_RC2], PSY_PHENO_RC2B[,intersecting_cols_RC2])
  1197. # Curating recall
  1198. PSY_PHENO_RC_TOTAL$Code = toupper(PSY_PHENO_RC_TOTAL$Code)
  1199. PSY_PHENO_RC_TOTAL$`Gender.(0)male.(1)female` = sapply(PSY_PHENO_RC_TOTAL$`Gender.(0)male.(1)female`, function(x){
  1200. if (x == 0){
  1201. return("M")
  1202. } else {
  1203. return("W")
  1204. }
  1205. })
  1206. colnames(PSY_PHENO_RC_TOTAL) = stri_replace_all_fixed(colnames(PSY_PHENO_RC_TOTAL), pattern = "Gender.(0)male.(1)female", replacement = "Gender")
  1207. ############ Saving phenotypes
  1208. write.csv(PSY_PHENO_RC_TOTAL, "PSY_PHENO_RC_TOTAL.csv", row.names = FALSE)
  1210. ################### Curating phenotype medication intake and removing duplicates ###################
  1211. ANXIETY_CURATED_PSY_PHENO_RC_TOTAL = openxlsx::read.xlsx("ANXIETY_PSY_PHENO_RC_TOTAL_CURATED.xlsx", sheet = 1)
  1212. CURATED_PSY_PHENO_RC_TOTAL = PSY_PHENO_RC_TOTAL
  1213. CURATED_PSY_PHENO_RC_TOTAL$Code = toupper(CURATED_PSY_PHENO_RC_TOTAL$Code)
  1215. any(duplicated(CURATED_PSY_PHENO_RC_TOTAL$Code)) # TRUE There are duplicates in phenotypes
  1216. CURATED_PSY_PHENO_RC_TOTAL$Code[duplicated(CURATED_PSY_PHENO_RC_TOTAL$Code)] # "PSY0662" "PSY0867"
  1218. duplicated_recall = CURATED_PSY_PHENO_RC_TOTAL[CURATED_PSY_PHENO_RC_TOTAL$Code %in% c("PSY0662", "PSY0867"),]
  1219. duplicated_recall_orig = Recall_2A_Pheno[Recall_2A_Pheno$Code %in% c("PSY0662", "PSY0867"),]
  1221. # "PSY0662" should be female as indicated in screening "Collected data_Screening.xlsx", height=159
  1222. as.Date(43158, origin = "1899-12-30") # "2018-02-27", timestamp from PSY_recall_1_2_Biobank_masterfile - Updated_Aleksandr_Sokolov_Oct2021.xlsx
  1223. # "PSY0867" # 02/27/2018 is the correct date -> Age 16
  1224. duplicated_recall_selected = duplicated_recall[c(1,3),]
  1226. # curating recall dataset
  1227. CURATED_PSY_PHENO_RC_TOTAL = CURATED_PSY_PHENO_RC_TOTAL[CURATED_PSY_PHENO_RC_TOTAL$Code %!in% duplicated_recall$Code,]
  1228. CURATED_PSY_PHENO_RC_TOTAL = rbind(CURATED_PSY_PHENO_RC_TOTAL, duplicated_recall_selected)
  1230. # curating recall anxiety dataset
  1231. duplicated_recall = ANXIETY_CURATED_PSY_PHENO_RC_TOTAL[ANXIETY_CURATED_PSY_PHENO_RC_TOTAL$Code %in% c("PSY0662", "PSY0867"),]
  1232. duplicated_recall_selected = duplicated_recall[c(1,3),]
  1233. ANXIETY_CURATED_PSY_PHENO_RC_TOTAL = ANXIETY_CURATED_PSY_PHENO_RC_TOTAL[ANXIETY_CURATED_PSY_PHENO_RC_TOTAL$Code %!in% duplicated_recall$Code,]
  1234. ANXIETY_CURATED_PSY_PHENO_RC_TOTAL = rbind(ANXIETY_CURATED_PSY_PHENO_RC_TOTAL, duplicated_recall_selected)
  1237. # check
  1238. any(duplicated(CURATED_PSY_PHENO_RC_TOTAL$Code)) # FALSE
  1239. any(duplicated(ANXIETY_CURATED_PSY_PHENO_RC_TOTAL$Code)) # FALSE
  1242. CURATED_PSY_PHENO_RC_TOTAL$Antidepr_anxiety = sapply(CURATED_PSY_PHENO_RC_TOTAL$Code, function(x){
  1243. if (x %in% ANXIETY_CURATED_PSY_PHENO_RC_TOTAL$Code){
  1244. treatment_status = ANXIETY_CURATED_PSY_PHENO_RC_TOTAL[ANXIETY_CURATED_PSY_PHENO_RC_TOTAL$Code == x, "Antidepr_anxiety"]
  1245. treatment_status = as.numeric(treatment_status)
  1246. } else {
  1247. treatment_status = NA
  1248. }
  1249. return(treatment_status)
  1250. })
  1252. # PSY0021
  1253. CURATED_PSY_PHENO_RC_TOTAL[CURATED_PSY_PHENO_RC_TOTAL$Code == "PSY0021", "Antidepr_anxiety"] = 0
  1256. ################### Initial Variable prep ###################
  1257. any(is.na(CURATED_PSY_PHENO_RC_TOTAL$SUAS.Total.score)) # FALSE
  1260. CURATED_PSY_PHENO_RC_TOTAL$Gender = factor(CURATED_PSY_PHENO_RC_TOTAL$Gender, levels = c("W", "M"))
  1261. CURATED_PSY_PHENO_RC_TOTAL$SUAS_binary = ifelse(CURATED_PSY_PHENO_RC_TOTAL$SUAS.Total.score >= 31,
  1262. "Risk group", "Normal")
  1263. CURATED_PSY_PHENO_RC_TOTAL$SUAS_binary = factor(CURATED_PSY_PHENO_RC_TOTAL$SUAS_binary, levels = c("Normal", "Risk group"))
  1265. # Variables
  1266. CURATED_PSY_PHENO_RC_TOTAL$Gender
  1267. CURATED_PSY_PHENO_RC_TOTAL$Age
  1268. CURATED_PSY_PHENO_RC_TOTAL$BMI
  1269. CURATED_PSY_PHENO_RC_TOTAL$SUAS.Total.score
  1270. CURATED_PSY_PHENO_RC_TOTAL$SUAS_binary
  1273. ################### Importing proteomics data ###################
  1274. # Preprocess OLINK function
  1275. preprocess_OLINK_main = function(PRFX_OLINK_DS){
  1276. # modifying original DS
  1277. PRFX_idx = which(PRFX_OLINK_DS[,1] == "Panel")
  1278. PRFX_OLINK_DS = PRFX_OLINK_DS[PRFX_idx:nrow(PRFX_OLINK_DS),]
  1279. PRFX_OLINK_DS = as.data.frame(t(PRFX_OLINK_DS))
  1280. colnames(PRFX_OLINK_DS) = PRFX_OLINK_DS[1,]
  1281. PRFX_OLINK_DS = PRFX_OLINK_DS[-1,]
  1282. rownames(PRFX_OLINK_DS) = PRFX_OLINK_DS$Assay
  1284. #preparing data
  1285. PRFX_probes = PRFX_OLINK_DS[,c(1:4, (ncol(PRFX_OLINK_DS)-2):ncol(PRFX_OLINK_DS))]
  1286. PRFX_probes = PRFX_probes[-((nrow(PRFX_OLINK_DS)-1):nrow(PRFX_OLINK_DS)),]
  1287. PRFX_probes$`Missing Data freq.` = as.numeric(PRFX_probes$`Missing Data freq.`)
  1288. PRFX_probes$LOD = as.numeric(PRFX_probes$LOD)
  1290. PRFX_sample_status = PRFX_OLINK_DS[(nrow(PRFX_OLINK_DS)-1):nrow(PRFX_OLINK_DS), ]
  1291. PRFX_sample_status = PRFX_sample_status[,colnames(PRFX_sample_status) %!in% colnames(PRFX_probes)]
  1292. PRFX_sample_status = as.data.frame(t(PRFX_sample_status))
  1293. PRFX_sample_status$Sample = rownames(PRFX_sample_status)
  1295. PRFX_expression = PRFX_OLINK_DS[, colnames(PRFX_OLINK_DS) %in% PRFX_sample_status$Sample]
  1296. PRFX_expression = PRFX_expression[-((nrow(PRFX_OLINK_DS)-1):nrow(PRFX_OLINK_DS)),]
  1297. PRFX_expression = as.matrix(PRFX_expression)
  1298. PRFX_matrix_names = rownames(PRFX_expression)
  1299. PRFX_expression = apply(PRFX_expression, 2, as.numeric)
  1300. rownames(PRFX_expression) = PRFX_matrix_names
  1302. #filtering probes and samples
  1303. PRFX_sample_status$QC_result = ifelse(PRFX_sample_status$`QC Warning` == "Pass", "Included", "Excluded")
  1304. PRFX_expression_pass = PRFX_expression[, colnames(PRFX_expression) %in% PRFX_sample_status[PRFX_sample_status$`QC Warning` == "Pass", "Sample"]]
  1306. PRFX_probes$Probe_QC = ifelse(PRFX_probes$`Missing Data freq.` > 0.2, "Excluded", "Pass")
  1307. for (PRFX_I in 1 : nrow(PRFX_probes)){
  1308. if (PRFX_probes$LOD[PRFX_I] > 2.5 & PRFX_probes$Probe_QC[PRFX_I] == "Excluded"){
  1309. PRFX_probes$Probe_QC[PRFX_I] = "High LOD"
  1310. }
  1311. }
  1312. PRFX_passed_probes = PRFX_probes[PRFX_probes$Probe_QC == "Pass", "Assay"]
  1313. PRFX_expression_pass = PRFX_expression_pass[rownames(PRFX_expression_pass) %in% PRFX_passed_probes,]
  1314. PRFX_expression_pass_strict = PRFX_expression_pass
  1316. for (PRFX_I in 1 : nrow(PRFX_expression_pass_strict)){
  1317. PRFX_Curr_Probe = rownames(PRFX_expression_pass_strict)[PRFX_I]
  1318. PRFX_Curr_LOD = PRFX_probes[PRFX_probes$Assay == PRFX_Curr_Probe, "LOD"]
  1319. PRFX_Curr_row = PRFX_expression_pass_strict[PRFX_I,]
  1320. PRFX_Curr_row = ifelse(PRFX_Curr_row >= PRFX_Curr_LOD, PRFX_Curr_row, NA)
  1321. PRFX_expression_pass_strict[PRFX_I,] = PRFX_Curr_row
  1322. }
  1324. PRFX_expression_pass = as.data.frame(PRFX_expression_pass)
  1325. PRFX_expression_pass_strict = as.data.frame(PRFX_expression_pass_strict)
  1326. PRFX_Output = list(Expression_matrix = PRFX_expression_pass,
  1327. Expression_matrix_strict = PRFX_expression_pass_strict,
  1328. Probes = PRFX_probes,
  1329. Samples = PRFX_sample_status)
  1330. return(PRFX_Output)
  1331. }
  1333. # Batch 2021
  1334. OLINK_DS_2021 = openxlsx::read.xlsx("OLINK_DIANA/Proteomics in healthy adolescents_PSY_NPX_no duplicates.xlsx")
  1335. OLINK_DS_2021_prep = preprocess_OLINK_main(PRFX_OLINK_DS = OLINK_DS_2021)
  1336. View(OLINK_DS_2021_prep$Samples)
  1337. View(OLINK_DS_2021_prep$Probes)
  1338. table(OLINK_DS_2021_prep$Samples$QC_result) # 4 - excluded, 313 - included
  1339. OLINK_DS_2021_prep$Samples[OLINK_DS_2021_prep$Samples$QC_result == "Excluded", ]
  1340. table(OLINK_DS_2021_prep$Probes$Probe_QC) # 43 - excluded, 49 - passed
  1341. OLINK_DS_2021_EXPR = OLINK_DS_2021_prep$Expression_matrix
  1343. # Batch 2022
  1344. OLINK_DS_2022 = openxlsx::read.xlsx("OLINK_AS/Proteomic analyses_Plasma_NPX.xlsx")
  1345. OLINK_DS_2022_prep = preprocess_OLINK_main(PRFX_OLINK_DS = OLINK_DS_2022)
  1346. View(OLINK_DS_2022_prep$Samples)
  1347. View(OLINK_DS_2022_prep$Probes)
  1348. table(OLINK_DS_2022_prep$Samples$QC_result) # 5 - excluded, 233 - included
  1349. OLINK_DS_2022_prep$Samples[OLINK_DS_2022_prep$Samples$QC_result == "Excluded", ]
  1350. table(OLINK_DS_2022_prep$Probes$Probe_QC) # 39 - excluded, 10 - High LOD, 43 - passed
  1351. OLINK_DS_2022_EXPR = OLINK_DS_2022_prep$Expression_matrix
  1353. # Check for repeated measurement
  1354. table(colnames(OLINK_DS_2022_EXPR) %in% colnames(OLINK_DS_2021_EXPR)) # 8 samples were analyzed the 2nd time
  1355. over_analyzed_samples = colnames(OLINK_DS_2022_EXPR)[colnames(OLINK_DS_2022_EXPR) %in% colnames(OLINK_DS_2021_EXPR)]
  1356. OLINK_DS_2022_EXPR = OLINK_DS_2022_EXPR[,colnames(OLINK_DS_2022_EXPR) %!in% over_analyzed_samples]
  1358. # Updating phenotype information
  1359. ALL_PSY_OLINK = c(OLINK_DS_2021_prep$Samples$Sample, OLINK_DS_2022_prep$Samples$Sample)
  1360. ALL_PSY_OLINK = unique(ALL_PSY_OLINK)
  1362. CURATED_PSY_PHENO_RC_TOTAL$OLINK_STATUS = NA
  1363. CURATED_PSY_PHENO_RC_TOTAL$OLINK_BATCH = NA
  1364. CURATED_PSY_PHENO_RC_TOTAL$OLINK_PLATE = NA
  1366. for (i in 1:length(ALL_PSY_OLINK)){
  1368. PRF_sample = ALL_PSY_OLINK[i]
  1369. PRF_name_frag = unlist(stri_split_fixed(PRF_sample, pattern = "_"))
  1370. PRF_name_frag = unlist(stri_split_fixed(PRF_name_frag, pattern = "-"))
  1372. # Meta data
  1373. PRF_praticip = PRF_name_frag[1]
  1374. PRF_visit = PRF_name_frag[2]
  1376. # Get Batch info
  1377. if (PRF_sample %in% OLINK_DS_2021_prep$Samples$Sample){
  1378. PRF_OLINK_BATCH = "2021"
  1379. PRF_SAMPLE_DF = OLINK_DS_2021_prep$Samples
  1380. } else if (PRF_sample %in% OLINK_DS_2022_prep$Samples$Sample){
  1381. PRF_OLINK_BATCH = "2022"
  1382. PRF_SAMPLE_DF = OLINK_DS_2022_prep$Samples
  1383. }
  1385. # Get Plate Info and Status
  1386. PRF_OLINK_STATUS = PRF_SAMPLE_DF[PRF_SAMPLE_DF$Sample == PRF_sample, "QC_result"]
  1387. PRF_OLINK_PLATE = PRF_SAMPLE_DF[PRF_SAMPLE_DF$Sample == PRF_sample, "Plate ID"]
  1389. if (stri_detect_fixed(PRF_visit, pattern = "2")){
  1390. # Data from recall
  1391. CURATED_PSY_PHENO_RC_TOTAL[CURATED_PSY_PHENO_RC_TOTAL$Code == PRF_praticip, "OLINK_STATUS"] = PRF_OLINK_STATUS
  1392. CURATED_PSY_PHENO_RC_TOTAL[CURATED_PSY_PHENO_RC_TOTAL$Code == PRF_praticip, "OLINK_BATCH"] = PRF_OLINK_BATCH
  1393. CURATED_PSY_PHENO_RC_TOTAL[CURATED_PSY_PHENO_RC_TOTAL$Code == PRF_praticip, "OLINK_PLATE"] = PRF_OLINK_PLATE
  1394. }
  1396. # Cleaning
  1397. rm(list = c("PRF_sample",
  1398. "PRF_name_frag",
  1399. "PRF_praticip",
  1400. "PRF_visit",
  1401. "PRF_SAMPLE_DF",
  1402. "PRF_OLINK_STATUS",
  1403. "PRF_OLINK_BATCH",
  1404. "PRF_OLINK_PLATE"))
  1406. }
  1408. table(CURATED_PSY_PHENO_RC_TOTAL$OLINK_STATUS)
  1409. table(CURATED_PSY_PHENO_RC_TOTAL$OLINK_BATCH)
  1410. table(CURATED_PSY_PHENO_RC_TOTAL$OLINK_PLATE)
  1412. ################### Participant Selection ###################
  1413. colnames(OLINK_DS_2022_EXPR) = toupper(colnames(OLINK_DS_2022_EXPR))
  1414. colnames(OLINK_DS_2021_EXPR) = toupper(colnames(OLINK_DS_2021_EXPR))
  1415. colnames(OLINK_DS_2022_EXPR) = stri_replace_all_fixed(colnames(OLINK_DS_2022_EXPR), pattern = "-", replacement = "_")
  1416. colnames(OLINK_DS_2021_EXPR) = stri_replace_all_fixed(colnames(OLINK_DS_2021_EXPR), pattern = "-", replacement = "_")
  1418. combined_cols = c(colnames(OLINK_DS_2022_EXPR), colnames(OLINK_DS_2021_EXPR))
  1419. combined_cols = toupper(combined_cols)
  1420. any(duplicated(combined_cols)) # FALSE -> no duplicated tubes
  1421. combined_cols = stri_replace_all_fixed(combined_cols, pattern = "-", replacement = "_")
  1423. recall_samples = combined_cols[multiple_stri_detector(combined_cols, pattern_vector = c("_2_", "_2B_", "_2A_"))]
  1424. recall_particip = sapply(recall_samples, function(x) unlist(stri_split_fixed(x, pattern = "_"))[1])
  1425. PSY_PHENO_RECALL_OLINK_VALID = CURATED_PSY_PHENO_RC_TOTAL[CURATED_PSY_PHENO_RC_TOTAL$Code %in% recall_particip, ]
  1427. # Filtering Samples
  1428. OLINK_DS_2021_EXPR_KEEP = OLINK_DS_2021_EXPR[, colnames(OLINK_DS_2021_EXPR) %in% recall_samples]
  1429. OLINK_DS_2022_EXPR_KEEP = OLINK_DS_2022_EXPR[, colnames(OLINK_DS_2022_EXPR) %in% recall_samples]
  1430. all(recall_samples %in% c(colnames(OLINK_DS_2021_EXPR), colnames(OLINK_DS_2022_EXPR))) # TRUE
  1432. # Check for duplicated individuals
  1433. total_samples_keep_individ = sapply(recall_samples, function(x) unlist(stri_split_fixed(x, pattern = "_"))[1])
  1434. any(duplicated(total_samples_keep_individ)) # -> FALSE, no duplicated persons
  1436. # Filtering phenotypes
  1437. any(duplicated(PSY_PHENO_RECALL_OLINK_VALID$Code)) # FALSE There are NO duplicates in phenotypes
  1438. PSY_PHENO_RECALL_OLINK_FINAL_curated = PSY_PHENO_RECALL_OLINK_VALID
  1440. # QC
  1441. any(duplicated(PSY_PHENO_RECALL_OLINK_FINAL_curated$Code)) # FALSE
  1443. # Finding overlapping proteins
  1444. proteins_passed_2021 = rownames(OLINK_DS_2021_EXPR_KEEP)
  1445. proteins_passed_2022 = rownames(OLINK_DS_2022_EXPR_KEEP)
  1446. proteins_passed_both_batches = intersect(proteins_passed_2021, proteins_passed_2022)
  1448. # Obtaining combined proteomic dataset
  1449. OLINK_DS_2021_EXPR_KEEP_both = OLINK_DS_2021_EXPR_KEEP[proteins_passed_both_batches,]
  1450. OLINK_DS_2022_EXPR_KEEP_both = OLINK_DS_2022_EXPR_KEEP[proteins_passed_both_batches,]
  1451. all(rownames(OLINK_DS_2021_EXPR_KEEP_both) == rownames(OLINK_DS_2022_EXPR_KEEP_both)) # -> TRUE
  1452. OLINK_DS_BOTH_BATCHES_FINAL = cbind(OLINK_DS_2021_EXPR_KEEP_both, OLINK_DS_2022_EXPR_KEEP_both)
  1454. # Adding Sample column
  1455. PSY_PHENO_RECALL_OLINK_FINAL_curated$Matching_sample = sapply(PSY_PHENO_RECALL_OLINK_FINAL_curated$Code, function(x){
  1456. sample = colnames(OLINK_DS_BOTH_BATCHES_FINAL)[stri_detect_fixed(str = colnames(OLINK_DS_BOTH_BATCHES_FINAL), pattern = x)]
  1457. return(sample)
  1458. })
  1460. # QC
  1461. any(duplicated(colnames(OLINK_DS_BOTH_BATCHES_FINAL))) # FALSE
  1462. any(duplicated(colnames(PSY_PHENO_RECALL_OLINK_FINAL_curated$Matching_sample))) # FALSE
  1463. # uncovered participant in phenotype
  1464. colnames(OLINK_DS_BOTH_BATCHES_FINAL)[colnames(OLINK_DS_BOTH_BATCHES_FINAL) %!in% PSY_PHENO_RECALL_OLINK_FINAL_curated$Matching_sample]
  1465. # all participants are covered!
  1467. PSY_PHENO_RECALL_OLINK_FINAL_curated$OLINK_PLATE = as.factor(PSY_PHENO_RECALL_OLINK_FINAL_curated$OLINK_PLATE)
  1468. table(PSY_PHENO_RECALL_OLINK_FINAL_curated$OLINK_PLATE)
  1470. # Reordering proteomics
  1471. OLINK_DS_BOTH_BATCHES_FINAL = OLINK_DS_BOTH_BATCHES_FINAL[,PSY_PHENO_RECALL_OLINK_FINAL_curated$Matching_sample]
  1472. all(PSY_PHENO_RECALL_OLINK_FINAL_curated$Matching_sample == colnames(OLINK_DS_BOTH_BATCHES_FINAL)) # -> TRUE, order is matching
  1474. # Checking missing values in the data
  1475. corrupt_samples_proteomics = apply(OLINK_DS_BOTH_BATCHES_FINAL, 2, function(x) any(is.na(x)))
  1476. OLINK_DS_BOTH_BATCHES_FINAL = OLINK_DS_BOTH_BATCHES_FINAL[,!corrupt_samples_proteomics]
  1477. PSY_PHENO_RECALL_OLINK_FINAL_curated = PSY_PHENO_RECALL_OLINK_FINAL_curated[!corrupt_samples_proteomics, ]
  1479. all(colnames(OLINK_DS_BOTH_BATCHES_FINAL) == PSY_PHENO_RECALL_OLINK_FINAL_curated$Matching_sample) # TRUE
  1481. ######## Saving processed data
  1482. write.csv(OLINK_DS_BOTH_BATCHES_FINAL, "OLINK_DS_BOTH_BATCHES_FINAL.csv", row.names = TRUE) # 404 samples, 43 proteins
  1483. write.csv(PSY_PHENO_RECALL_OLINK_FINAL_curated, "PSY_PHENO_RECALL_OLINK_FINAL_curated.csv", row.names = FALSE)
  1486. ################### SUAS cut-off cpecs ###################
  1488. CURATED_PSY_PHENO_RC_TOTAL$SUAS_binary = ifelse(CURATED_PSY_PHENO_RC_TOTAL$SUAS.Total.score >= 29,
  1489. "Risk group", "Normal")
  1490. CURATED_PSY_PHENO_RC_TOTAL$SUAS_binary = factor(CURATED_PSY_PHENO_RC_TOTAL$SUAS_binary, levels = c("Normal", "Risk group"))
  1493. ################### Adjusting for batch covariates ###################
  1495. ########### Raw Adjustment (no Covariates)
  1496. OLINK_DS_BOTH_BATCHES_FINAL_ADJUSTED = ComBat(OLINK_DS_BOTH_BATCHES_FINAL, batch = PSY_PHENO_RECALL_OLINK_FINAL_curated$OLINK_PLATE)
  1497. OLINK_DS_BOTH_BATCHES_FINAL_ADJUSTED = as.data.frame(OLINK_DS_BOTH_BATCHES_FINAL_ADJUSTED)
  1500. ########### Adjustment with covar (no antidepr)
  1502. # Batch effect
  1503. # Participants must be columns
  1504. model_batch = model.matrix(~ SUAS_binary + Gender + Age, data = PSY_PHENO_RECALL_OLINK_FINAL_curated)
  1505. OLINK_DS_BOTH_BATCHES_FINAL_SUAS_ADJ = ComBat(dat=OLINK_DS_BOTH_BATCHES_FINAL, batch = PSY_PHENO_RECALL_OLINK_FINAL_curated$OLINK_PLATE, mod=model_batch)
  1506. OLINK_DS_BOTH_BATCHES_FINAL_SUAS_ADJ = as.data.frame(OLINK_DS_BOTH_BATCHES_FINAL_SUAS_ADJ)
  1508. ########### Adjustment with covar (including antidepr)
  1510. # Selecting participants that all have SUAS score and information on antidepressants
  1511. PSY_PHENO_RECALL_OLINK_FINAL_ANTIDEPR = PSY_PHENO_RECALL_OLINK_FINAL_curated[!is.na(PSY_PHENO_RECALL_OLINK_FINAL_curated$Antidepr_anxiety),]
  1512. OLINK_DS_BOTH_BATCHES_FINAL_ANTIDEPR = OLINK_DS_BOTH_BATCHES_FINAL[,PSY_PHENO_RECALL_OLINK_FINAL_ANTIDEPR$Matching_sample]
  1513. all(PSY_PHENO_RECALL_OLINK_FINAL_ANTIDEPR$Matching_sample == colnames(OLINK_DS_BOTH_BATCHES_FINAL_ANTIDEPR)) # -> TRUE, order is matching
  1514. PSY_PHENO_RECALL_OLINK_FINAL_ANTIDEPR$Antidepr_anxiety = factor(PSY_PHENO_RECALL_OLINK_FINAL_ANTIDEPR$Antidepr_anxiety, levels = c(0,1))
  1516. # Batch effect
  1517. # Participants must be columns
  1518. PSY_PHENO_RECALL_OLINK_FINAL_ANTIDEPR$OLINK_PLATE = factor(PSY_PHENO_RECALL_OLINK_FINAL_ANTIDEPR$OLINK_PLATE)
  1519. model_batch = model.matrix(~ SUAS_binary + Antidepr_anxiety + Gender + Age, data = PSY_PHENO_RECALL_OLINK_FINAL_ANTIDEPR)
  1521. test_table = model.matrix(~ SUAS_binary + Antidepr_anxiety + Gender + Age + OLINK_PLATE, data = PSY_PHENO_RECALL_OLINK_FINAL_ANTIDEPR)
  1522. test_table = cor(test_table, method = "pearson", use = "complete.obs")
  1524. OLINK_DS_BOTH_BATCHES_FINAL_SUAS_ANTIDEPR_ADJ = ComBat(dat=OLINK_DS_BOTH_BATCHES_FINAL_ANTIDEPR, batch = PSY_PHENO_RECALL_OLINK_FINAL_ANTIDEPR$OLINK_PLATE, mod=model_batch)
  1525. OLINK_DS_BOTH_BATCHES_FINAL_SUAS_ANTIDEPR_ADJ = as.data.frame(OLINK_DS_BOTH_BATCHES_FINAL_SUAS_ANTIDEPR_ADJ)
  1527. ################### Limma differential expression ###################
  1528. PSY_PHENO_RECALL_OLINK_FINAL_curated$OLINK_PLATE = factor(PSY_PHENO_RECALL_OLINK_FINAL_curated$OLINK_PLATE)
  1530. # Analysis with limma
  1531. design.matrix = model.matrix(~ SUAS_binary + Gender + Age, data = PSY_PHENO_RECALL_OLINK_FINAL_curated)
  1532. fit = lmFit(OLINK_DS_BOTH_BATCHES_FINAL_SUAS_ADJ, design.matrix)
  1533. fitE = eBayes(fit, robust = TRUE)
  1534. Top_table = limma::topTable(fit = fitE, coef = 2, adjust.method = "fdr", number = Inf)
  1535. Top_table$Gene = rownames(Top_table) # RBKS CCL27 AKT1S1
  1537. ################### Limma differential expression (antidepress) ###################
  1539. # Analysis with limma
  1540. design.matrix = model.matrix(~ SUAS_binary + Antidepr_anxiety + Gender + Age, data = PSY_PHENO_RECALL_OLINK_FINAL_ANTIDEPR)
  1541. fit = lmFit(OLINK_DS_BOTH_BATCHES_FINAL_SUAS_ANTIDEPR_ADJ, design.matrix)
  1542. fitE = eBayes(fit, robust = TRUE)
  1543. Top_table_antidepr = limma::topTable(fit = fitE, coef = 2, adjust.method = "fdr", number = Inf)
  1544. Top_table_antidepr$Gene = rownames(Top_table_antidepr) # PLA2G10
  1546. # Boxplot for CCL27 in PSY
  1547. probe = "CCL27"
  1548. plot_df = PSY_PHENO_RECALL_OLINK_FINAL_curated
  1550. probe_vals = OLINK_DS_BOTH_BATCHES_FINAL_SUAS_ADJ[probe,plot_df$Matching_sample]
  1551. all(colnames(probe_vals) == plot_df$Matching_sample)
  1552. probe_vals = as.numeric(probe_vals)
  1553. plot_df = cbind(plot_df, probe_vals)
  1555. plot = ggplot(data = plot_df, aes(x = SUAS_binary, y = probe_vals, fill = SUAS_binary)) +
  1556. geom_boxplot(notchwidth = 0.25) +
  1557. scale_fill_manual(values = c("#C44652","#EAC151")) +
  1558. geom_violin(alpha = 0.25) +
  1559. labs(y = "NPX: CCL27", x = "SUAS risk group", fill = "Disease status") +
  1560. ggthemes::theme_clean(base_size = 14)
  1561. plot
  1562. ggsave(file = "CCL27_PSY_box.pdf", plot = plot, width = 15, height = 20, units = "cm")
  1564. # Boxplot for AKT1S1 in PSY
  1565. probe = "AKT1S1"
  1566. plot_df = PSY_PHENO_RECALL_OLINK_FINAL_curated
  1568. probe_vals = OLINK_DS_BOTH_BATCHES_FINAL_SUAS_ADJ[probe,plot_df$Matching_sample]
  1569. all(colnames(probe_vals) == plot_df$Matching_sample)
  1570. probe_vals = as.numeric(probe_vals)
  1571. plot_df = cbind(plot_df, probe_vals)
  1573. plot = ggplot(data = plot_df, aes(x = SUAS_binary, y = probe_vals, fill = SUAS_binary)) +
  1574. geom_boxplot(notchwidth = 0.25) +
  1575. scale_fill_manual(values = c("#C44652","#EAC151")) +
  1576. geom_violin(alpha = 0.25) +
  1577. labs(y = "NPX: AKT1S1", x = "SUAS risk group", fill = "Disease status") +
  1578. ggthemes::theme_clean(base_size = 14)
  1579. plot
  1580. ggsave(file = "AKT1S1_PSY_box.pdf", plot = plot, width = 15, height = 20, units = "cm")
  1582. ################### Saving PSY data ###################
  1584. PSY_results_limma = list(
  1585. "PSY_no_antidepr" = Top_table,
  1586. "PSY_antidep_adj" = Top_table_antidepr
  1588. )
  1590. names_for_list_PSY = c(
  1591. "PSY_no_antidepr",
  1592. "PSY_antidep_adj"
  1593. )
  1595. wb = createWorkbook()
  1597. for (i in 1:2) {
  1598. # Add a new worksheet with the sheet name
  1599. addWorksheet(wb, names_for_list_PSY[i])
  1601. # Write the data frame to the worksheet
  1602. writeData(wb, sheet = names_for_list_PSY[i], PSY_results_limma[[i]])
  1604. # Adjust column widths to fit the text
  1605. setColWidths(wb, sheet = names_for_list_PSY[i], cols = 1:ncol(PSY_results_limma[[i]]), widths = "auto")
  1606. }
  1608. saveWorkbook(wb, "PSY_proteomics_SUAS.xlsx", overwrite = TRUE)
  1610. # Dataset caracterization
  1612. PSY_charact_df = PSY_PHENO_RECALL_OLINK_FINAL_curated
  1614. # variables
  1615. # SUAS_binary + Antidepr_anxiety + Gender + Age
  1616. which(colnames(PSY_charact_df) == "Code") # 1
  1617. which(colnames(PSY_charact_df) == "Gender") # 4
  1618. which(colnames(PSY_charact_df) == "Age") # 3
  1619. which(colnames(PSY_charact_df) == "SUAS_binary") # 93
  1620. which(colnames(PSY_charact_df) == "Antidepr_anxiety") # 92
  1621. PSY_charact_df$Antidepr_anxiety = ifelse(PSY_charact_df$Antidepr_anxiety == 0, "No treatment", "Treatment")
  1622. PSY_charact_df$Antidepr_anxiety = factor(PSY_charact_df$Antidepr_anxiety, levels = c("No treatment", "Treatment"))
  1624. PSY_summary = characterize_dataset_generelized_two_subgroups(dataset = PSY_charact_df,
  1625. study_char = "PSY",
  1626. contrast_col_number = 93,
  1627. participants_col_number = 1,
  1628. model_covariates_col_vector = c(4,3,92),
  1629. columns_to_characterise_vector = c(93,4,3,92),
  1630. Remove_NA_predictors = FALSE,
  1631. drop_P = TRUE)
  1633. PSY_summary = PSY_summary$Table
  1634. PSY_summary = PSY_summary[-3,]
  1635. PSY_summary
  1636. PSY_summary = apply(X = PSY_summary, MARGIN = 2, function(x){
  1637. x = multiple_stri_replacer(string = x,
  1638. pattern_vector = c("Initial dataset includes",
  1639. "SUAS_binary",
  1640. "Antidepr_anxiety",
  1641. "W: ",
  1642. "M: "),
  1643. replacement_vector = c("Dataset includes",
  1644. "Suicide risk group (SUAS)",
  1645. "Treatment for depression or anxiety",
  1646. "Female: ",
  1647. "Male: "))
  1649. })
  1651. PSY_summary = as.data.frame(PSY_summary)
  1652. PSY_summary
  1654. wb = createWorkbook()
  1656. # Add a new worksheet with the sheet name
  1657. addWorksheet(wb, sheetName = "PSY_demographics")
  1659. # Write the data frame to the worksheet
  1660. writeData(wb, sheet = "PSY_demographics", PSY_summary)
  1662. # Adjust column widths to fit the text
  1663. setColWidths(wb, sheet = "PSY_demographics", cols = 1:ncol(PSY_summary), widths = "auto")
  1665. # Adjust row heights
  1666. setRowHeights(wb, sheet = "PSY_demographics", rows = 1:nrow(PSY_summary), heights = "auto")
  1668. saveWorkbook(wb, "PSY_demographics.xlsx", overwrite = TRUE)
  1671. ################### Array Suicide datasets GEO ###################
  1672. # Jul 17 2024
  1673. # Organism - human
  1674. # Study type: ("suicide"[MeSH Terms] OR suicide[All Fields]) AND ("human"[Organism] AND ("Expression profiling by array"[Filter] OR "Expression profiling by genome tiling array"[Filter]))
  1676. # https://www.ncbi.nlm.nih.gov/gds/?term=suicide
  1678. ## GSE146446
  1679. geo_test = getGEO("GSE146446")
  1680. geo_test = geo_test[[1]]
  1681. geo_test_pheno = pData(geo_test)
  1682. geo_test_pheno$`tissue:ch1` # only blood -> no data on suicide
  1683. # GSE146446 is only for antidepressant response
  1685. ## GSE92538
  1686. # seems good
  1687. geo_test = getGEO("GSE92538")
  1688. geo_test = geo_test[[1]]
  1689. geo_test_pheno = pData(geo_test)
  1690. geo_test_pheno$`suicide (1=yes):ch1`
  1691. table(geo_test_pheno$`suicide (1=yes):ch1`) # 95 33
  1692. # GSE92538 seems good
  1694. ## GSE208338
  1695. geo_test = getGEO("GSE208338")
  1696. geo_test = geo_test[[1]]
  1697. geo_test_pheno = pData(geo_test)
  1698. table(geo_test_pheno$`suicide:ch1`, geo_test_pheno$`diagnosis:ch1`)
  1699. # GSE208338 this cohort is good!
  1701. ## GSE5388 and GSE5389
  1702. geo_test = getGEO("GSE5388")
  1703. geo_test = geo_test[[1]]
  1704. geo_test_pheno = pData(geo_test)
  1705. table(geo_test_pheno$`Suicide:ch1`, geo_test_pheno$`Disease_status:ch1`)
  1706. # GSE5388 seems good
  1708. geo_test = getGEO("GSE5389")
  1709. geo_test = geo_test[[1]]
  1710. geo_test_pheno = pData(geo_test)
  1711. table(geo_test_pheno$`Suicide:ch1`, geo_test_pheno$`Disease_status:ch1`)
  1712. # GSE5389 seems good
  1714. ## GSE66937
  1715. geo_test = getGEO("GSE66937")
  1716. geo_test = geo_test[[1]]
  1717. geo_test_pheno = pData(geo_test)
  1718. table(geo_test_pheno$`group:ch1`, geo_test_pheno$`brain region:ch1`)
  1719. # GSE66937 seems good
  1721. ## GSE125681
  1722. geo_test = getGEO("GSE125681")
  1723. geo_test = geo_test[[1]]
  1724. geo_test_pheno = pData(geo_test)
  1725. table(geo_test_pheno$`suicide:ch1`)
  1726. # GSE125681 only suicide and there are no controls...
  1729. ## GSE24095
  1730. geo_test = getGEO("GSE24095")
  1731. geo_test = geo_test[[1]]
  1732. geo_test_pheno = pData(geo_test)
  1733. table(geo_test_pheno$`Cause of death:ch1`)
  1734. table(geo_test_pheno$`Cause of death:ch2`)
  1735. # GSE24095 weird design -> better to avoid?
  1736. # Not a comparison of suicide VS control -> exclude
  1738. ## GSE199536
  1739. geo_test = getGEO("GSE199536")
  1740. geo_test = geo_test[[1]]
  1741. geo_test_pheno = pData(geo_test)
  1742. table(geo_test_pheno$`group:ch1`)
  1743. # GSE199536 seems good
  1745. ## GSE20568
  1746. geo_test = getGEO("GSE20568")
  1747. geo_test = geo_test[[1]]
  1748. geo_test_pheno = pData(geo_test)
  1749. table(geo_test_pheno$`cause of death:ch1`) # 1 suicide -> useless
  1750. # GSE20568 is useless
  1752. ## GSE93112
  1753. geo_test = getGEO("GSE93112")
  1754. geo_test = geo_test[[1]]
  1755. geo_test_pheno = pData(geo_test)
  1756. table(geo_test_pheno$`disease state:ch1`)
  1757. table(geo_test_pheno$`response phenotype, alda scale:ch1`)
  1758. # response to treatment Yes/No
  1759. # GSE93112 is useless
  1761. ## GSE9058
  1762. geo_test = getGEO("GSE9058")
  1763. geo_test = geo_test[[1]]
  1764. geo_test_pheno = pData(geo_test)
  1765. # bad phenotypes and 1 suicide -> useless
  1766. # GSE9058 is useless
  1768. #### Array Cohorts to include: GSE208338, GSE5388, GSE5389, GSE66937, GSE199536, GSE92538
  1771. ################### RNA Seq Suicide datasets GEO ###################
  1773. # ("suicide"[MeSH Terms] OR suicide[All Fields]) AND ("human"[Organism] AND "Expression profiling by high throughput sequencing"[Filter])
  1775. ## GSE102556
  1776. geo_test = getGEO("GSE102556")
  1777. geo_test = geo_test[[1]]
  1778. geo_test_pheno = pData(geo_test)
  1779. table(geo_test_pheno$organism_ch1) # 263 homo sapiens
  1780. table(geo_test_pheno$`tissue:ch1`, geo_test_pheno$`Cause of death:ch1`)
  1781. # GSE102556 is good
  1783. ## GSE202537
  1784. geo_test = getGEO("GSE202537")
  1785. geo_test = geo_test[[1]]
  1786. geo_test_pheno = pData(geo_test)
  1787. table(geo_test_pheno$organism_ch1) # 215 homo sapiens
  1788. table(geo_test_pheno$`tissue:ch1`, geo_test_pheno$`manner of death:ch1`)
  1789. # GSE202537 is good
  1791. ## GSE243356
  1792. geo_test = getGEO("GSE243356")
  1793. geo_test = geo_test[[1]]
  1794. geo_test_pheno = pData(geo_test)
  1795. table(geo_test_pheno$organism_ch1) # 61 homo sapiens
  1796. table(geo_test_pheno$`group:ch1`, geo_test_pheno$`tissue:ch1`)
  1797. # GSE243356 is good
  1799. ## GSE247998
  1800. geo_test = getGEO("GSE247998")
  1801. geo_test = geo_test[[1]]
  1802. geo_test_pheno = pData(geo_test)
  1803. table(geo_test_pheno$organism_ch1) # 100 homo sapiens
  1804. table(geo_test_pheno$`suicide attempt:ch1`)
  1805. # GSE247998 is good
  1807. ## GSE101521
  1808. geo_test = getGEO("GSE101521")
  1809. geo_test = geo_test[[1]]
  1810. geo_test_pheno = pData(geo_test)
  1811. table(geo_test_pheno$organism_ch1) # 27 homo sapiens
  1812. table(geo_test_pheno$`diagnosis:ch1`)
  1813. # GSE101521 is good
  1815. ## GSE157197
  1816. geo_test = getGEO("GSE157197")
  1817. geo_test = geo_test[[1]]
  1818. geo_test_pheno = pData(geo_test)
  1819. table(geo_test_pheno$organism_ch1) # 45 homo sapiens
  1820. table(geo_test_pheno$`diagnosis:ch1`)
  1821. # GSE157197 is useless
  1823. ## GSE213982
  1824. geo_test = getGEO("GSE213982")
  1825. geo_test = geo_test[[1]]
  1826. geo_test_pheno = pData(geo_test)
  1827. table(geo_test_pheno$organism_ch1) # 36 homo sapiens
  1828. table(geo_test_pheno$`group:ch1`)
  1829. # GSE213982 is good
  1831. ## GSE144136
  1832. # GSE144136 is single cell -> skipped for now as difficult to compare to bulk
  1834. ## GSE248260
  1835. geo_test = getGEO("GSE248260")
  1836. geo_test = geo_test[[1]]
  1837. geo_test_pheno = pData(geo_test)
  1838. table(geo_test_pheno$organism_ch1) # 24 homo sapiens
  1839. table(geo_test_pheno$`suicide:ch1`)
  1840. # GSE248260 is good
  1842. ## GSE12297
  1843. geo_test = getGEO("GSE12297")
  1844. geo_test = geo_test[[1]]
  1845. geo_test_pheno = pData(geo_test)
  1846. table(geo_test_pheno$organism_ch1) # 20 homo sapiens
  1847. # only 3 suicide
  1848. # GSE12297 is useless
  1851. ## GSE193417
  1852. geo_test = getGEO("GSE193417")
  1853. geo_test = geo_test[[1]]
  1854. geo_test_pheno = pData(geo_test)
  1855. table(geo_test_pheno$organism_ch1) # 12 homo sapiens
  1856. # only 3 suicide
  1857. # GSE193417 is useless
  1859. ## GSE208438
  1860. # GSE208438 is on cell lines -> useless
  1862. # Single cell
  1863. # GSE144136 (SC) MDD vs control (suicide)
  1864. # GSE213982 (SC) MDD vs control (suicide)
  1866. #### RNA seq cohorts to include: GSE102556, GSE202537, GSE243356, GSE247998, GSE101521, GSE213982, GSE248260
  1868. ################### Selected cohorts ###################
  1870. #### Array Cohorts to include: GSE208338, GSE5388, GSE5389, GSE66937, GSE199536, GSE92538 (limma-trend (standard pipeline))
  1871. #### RNA seq cohorts to include: GSE102556, GSE202537, GSE243356, GSE247998, GSE101521, GSE213982, GSE144136, GSE248260 (limma-voom (RNAseq specific))
  1873. # genome version in records:
  1874. # GSE102556 hg19 n>240
  1875. # GSE243356 most likely hg19 n=61
  1876. # GSE248260 hg19 n=24 # Macedonian/New York State Psychiatric Institute Brain Collection
  1877. # GSE247998 hg19 n=100 # Study participants were adults 20–65 years old and were recruited in the New York metropolitan area by advertising, via the Columbia University Medical Center Portal
  1878. # GSE202537 hg38 n=215
  1879. # GSE101521 hg19 n=86
  1881. # Single cell GSE144136
  1882. # Single cell GSE213982
  1884. # Criteria: More than 5 suicide cases
  1885. # Criteria: Individual signals/counts are available
  1886. # Criteria: Study/Methods allow comparison between suicide and control
  1887. # Criteria: Human samples only
  1889. # Single cell (separate)
  1890. # GSE144136 (SC) MDD vs control (suicide)
  1891. # GSE213982 (SC) MDD vs control (suicide)
  1893. #### Competitor: https://www.sciencedirect.com/science/article/pii/S0924977X21016631?via%3Dihub#sec0002
  1894. # competitor cohorts:
  1895. # GSE5389
  1896. # GSE5388
  1897. # GSE66937
  1898. # GSE102556
  1899. # GSE92538
  1902. #### Inspection of datasets from https://www.sciencedirect.com/science/article/pii/S0022395621000480?via%3Dihub (post-analysis)
  1903. # GSE21138 - no clear evidence for suicide in samples on GEO, paper has potentially matchable variables in the supplement https://www.sciencedirect.com/science/article/pii/S0006899308019410?via%3Dihub#app1
  1904. # GSE92538 - used by us
  1905. # GSE53987 - no clear evidence for suicide in samples on GEO, paper has NO variables (only grouped): https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0121744#pone-0121744-t001
  1906. # GSE54567 - no clear evidence for suicide in samples on GEO, paper does not show variables (shared https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0090980)
  1907. # GSE54568 - no clear evidence for suicide in samples on GEO, paper does not show variables (shared https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0090980)
  1908. # GSE5388 - used by us
  1909. # GSE5389 - used by us
  1910. # GSE54570 - no clear evidence for suicide in samples on GEO, paper does not show variables (shared https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0090980)
  1911. # SMRI 2 and SMRI 7 were not available as www.stanleyresearch.org was not responsive
  1913. # Conclusion -> GSE21138 could be theoretically useful due to 13 suicide cases however does not fit inclusion protocol based on GEO search
  1914. # Other datasets have no readily available infromation on suicide or already used by us
  1916. ################### Importing DE summary stats ###################
  1917. summary_paths_no_covar = list.files(path = "Data_preprocessing_analysis", recursive = TRUE, full.names = TRUE)
  1918. summary_paths_no_covar = summary_paths_no_covar[stri_detect_fixed(summary_paths_no_covar, pattern = "no_covar")]
  1919. summary_paths_no_covar = summary_paths_no_covar[!stri_detect_fixed(summary_paths_no_covar, pattern = "_SV")]
  1920. summary_paths_no_covar = summary_paths_no_covar[!stri_detect_fixed(summary_paths_no_covar, pattern = "_cell_types")]
  1922. summary_paths_with_covar = list.files(path = "Data_preprocessing_analysis", recursive = TRUE, full.names = TRUE)
  1923. summary_paths_with_covar = summary_paths_with_covar[stri_detect_fixed(summary_paths_with_covar, pattern = "Top_table")]
  1924. summary_paths_with_covar = summary_paths_with_covar[!stri_detect_fixed(summary_paths_with_covar, pattern = "_SV")]
  1925. summary_paths_with_covar = summary_paths_with_covar[!stri_detect_fixed(summary_paths_with_covar, pattern = "no_covar")]
  1927. summary_paths_cell_types = list.files(path = "Data_preprocessing_analysis", recursive = TRUE, full.names = TRUE)
  1928. summary_paths_cell_types = summary_paths_cell_types[stri_detect_fixed(summary_paths_cell_types, pattern = "_cell_types")]
  1930. data_list_no_covar = lapply(summary_paths_no_covar, function(x){
  1931. study_name = unlist(stri_split_fixed(x, pattern = "/"))
  1932. study_name = study_name[length(study_name)]
  1933. study_name = stri_replace_all_fixed(study_name, pattern = "_Top_table_no_covar.csv", replacement = "")
  1935. x = smart_fread(x)
  1936. x$Study = study_name
  1937. return(x)
  1938. })
  1939. reference_colnames = colnames(data_list_no_covar[[1]])
  1940. data_list_no_covar = lapply(data_list_no_covar, function(x){
  1942. if ("Gene_symbol_non_hgnc" %!in% colnames(x)){
  1943. x$Gene_symbol_non_hgnc = NA
  1944. }
  1946. x = x[,reference_colnames]
  1947. return(x)
  1948. })
  1951. data_list_with_covar = lapply(summary_paths_with_covar, function(x){
  1952. study_name = unlist(stri_split_fixed(x, pattern = "/"))
  1953. study_name = study_name[length(study_name)]
  1954. study_name = stri_replace_all_fixed(study_name, pattern = "_Top_table.csv", replacement = "")
  1956. x = smart_fread(x)
  1957. x$Study = study_name
  1958. return(x)
  1959. })
  1960. reference_colnames = colnames(data_list_with_covar[[1]])
  1961. data_list_with_covar = lapply(data_list_with_covar, function(x){
  1963. if ("Gene_symbol_non_hgnc" %!in% colnames(x)){
  1964. x$Gene_symbol_non_hgnc = NA
  1965. }
  1967. x = x[,reference_colnames]
  1968. return(x)
  1969. })
  1971. data_list_cell_types = lapply(summary_paths_cell_types, function(x){
  1972. study_name = unlist(stri_split_fixed(x, pattern = "/"))
  1973. study_name = study_name[length(study_name)]
  1974. study_name = stri_replace_all_fixed(study_name, pattern = "_Top_table_no_covar_cell_types.csv", replacement = "")
  1976. x = smart_fread(x)
  1977. x$Study = study_name
  1978. return(x)
  1979. })
  1980. data_list_cell_types = lapply(data_list_cell_types, function(x){
  1982. if ("Gene_symbol_non_hgnc" %!in% colnames(x)){
  1983. x$Gene_symbol_non_hgnc = NA
  1984. }
  1986. x = x[,reference_colnames]
  1987. return(x)
  1988. })
  1990. # SV datasets
  1991. summary_paths_SV = list.files(path = "Data_preprocessing_analysis", recursive = TRUE, full.names = TRUE)
  1992. summary_paths_SV = summary_paths_SV[stri_detect_fixed(summary_paths_SV, pattern = "Top_table")]
  1993. summary_paths_SV = summary_paths_SV[stri_detect_fixed(summary_paths_SV, pattern = "_SV")]
  1995. data_list_SV = lapply(summary_paths_SV, function(x){
  1996. study_name = unlist(stri_split_fixed(x, pattern = "/"))
  1997. study_name = study_name[length(study_name)]
  1998. study_name = stri_replace_all_fixed(study_name, pattern = "_Top_table_SV.csv", replacement = "")
  2000. x = smart_fread(x)
  2001. x$Study = study_name
  2002. return(x)
  2003. })
  2004. reference_colnames = colnames(data_list_SV[[1]])
  2005. data_list_SV = lapply(data_list_SV, function(x){
  2007. if ("Gene_symbol_non_hgnc" %!in% colnames(x)){
  2008. x$Gene_symbol_non_hgnc = NA
  2009. }
  2011. x = x[,reference_colnames]
  2012. return(x)
  2013. })
  2015. ################### Harmonizing gene names ###################
  2017. gene_names_no_covar = lapply(data_list_no_covar, function(x){
  2019. gene_symbols = x$Gene_symbol
  2020. gene_symbols_non_HGNC = x$Gene_symbol_non_hgnc
  2022. gene_symbols_all = c(gene_symbols, gene_symbols_non_HGNC)
  2023. gene_symbols_all = unlist(gene_symbols_all)
  2024. gene_symbols_all = unique(gene_symbols_all)
  2025. gene_symbols_all = unlist(stri_split_fixed(gene_symbols_all, pattern = ";"))
  2026. gene_symbols_all = gene_symbols_all[!is.na(gene_symbols_all)]
  2027. gene_symbols_all = gene_symbols_all[gene_symbols_all != ""]
  2028. gene_symbols_all = gene_symbols_all[gene_symbols_all != " "]
  2029. gene_symbols_all = gene_symbols_all[gene_symbols_all != "NA"]
  2031. return(gene_symbols_all)
  2032. })
  2033. gene_names_no_covar = unlist(gene_names_no_covar)
  2034. gene_names_no_covar = unique(gene_names_no_covar)
  2036. gene_names_with_covar = lapply(data_list_with_covar, function(x){
  2038. gene_symbols = x$Gene_symbol
  2039. gene_symbols_non_HGNC = x$Gene_symbol_non_hgnc
  2041. gene_symbols_all = c(gene_symbols, gene_symbols_non_HGNC)
  2042. gene_symbols_all = unlist(gene_symbols_all)
  2043. gene_symbols_all = unique(gene_symbols_all)
  2044. gene_symbols_all = unlist(stri_split_fixed(gene_symbols_all, pattern = ";"))
  2045. gene_symbols_all = gene_symbols_all[!is.na(gene_symbols_all)]
  2046. gene_symbols_all = gene_symbols_all[gene_symbols_all != ""]
  2047. gene_symbols_all = gene_symbols_all[gene_symbols_all != " "]
  2048. gene_symbols_all = gene_symbols_all[gene_symbols_all != "NA"]
  2050. return(gene_symbols_all)
  2051. })
  2052. gene_names_with_covar = unlist(gene_names_with_covar)
  2053. gene_names_with_covar = unique(gene_names_with_covar)
  2055. gene_names_from_sv = lapply(data_list_SV, function(x){
  2057. gene_symbols = x$Gene_symbol
  2058. gene_symbols_non_HGNC = x$Gene_symbol_non_hgnc
  2060. gene_symbols_all = c(gene_symbols, gene_symbols_non_HGNC)
  2061. gene_symbols_all = unlist(gene_symbols_all)
  2062. gene_symbols_all = unique(gene_symbols_all)
  2063. gene_symbols_all = unlist(stri_split_fixed(gene_symbols_all, pattern = ";"))
  2064. gene_symbols_all = gene_symbols_all[!is.na(gene_symbols_all)]
  2065. gene_symbols_all = gene_symbols_all[gene_symbols_all != ""]
  2066. gene_symbols_all = gene_symbols_all[gene_symbols_all != " "]
  2067. gene_symbols_all = gene_symbols_all[gene_symbols_all != "NA"]
  2069. return(gene_symbols_all)
  2070. })
  2071. gene_names_from_sv = unlist(gene_names_from_sv)
  2072. gene_names_from_sv = unique(gene_names_from_sv)
  2074. gene_names_full = c(gene_names_no_covar, gene_names_with_covar, gene_names_from_sv)
  2075. gene_names_full = unique(gene_names_full)
  2076. harmoniz_genes_check = check_gene_symbol_NIH(PRF_gene_symbols = gene_names_full,
  2077. PRF_ref_NIH_expanded = Homo_Sapiens_Gene_info_NIH_expanded,
  2078. PRF_replace_NA_with_old = TRUE)
  2079. harmoniz_genes_check_dict = harmoniz_genes_check$Suggested.Symbol
  2080. names(harmoniz_genes_check_dict) = harmoniz_genes_check$x
  2083. combined_df_no_covar = do.call(rbind, data_list_no_covar)
  2084. combined_df_wth_covar = do.call(rbind, data_list_with_covar)
  2085. combined_df_cell_types = do.call(rbind, data_list_cell_types)
  2087. combined_df_no_covar$Corrected_symbol = mapply(function(x,y){
  2088. combined_symbols = c(x, y)
  2089. combined_symbols = combined_symbols[!is.na(combined_symbols)]
  2091. if (length(combined_symbols)<1){
  2092. return(NA)
  2093. }
  2095. combined_symbols = stri_split_fixed(combined_symbols, pattern = ";")
  2096. combined_symbols = unlist(combined_symbols)
  2097. combined_symbols = combined_symbols[combined_symbols!= ""]
  2098. combined_symbols = combined_symbols[combined_symbols!= " "]
  2099. combined_symbols = combined_symbols[combined_symbols!= "NA"]
  2101. verified_symbols = sapply(combined_symbols, function(x) harmoniz_genes_check_dict[x])
  2102. verified_symbols = unique(verified_symbols)
  2103. verified_symbols = paste0(verified_symbols, collapse = ";")
  2104. return(verified_symbols)
  2105. }, combined_df_no_covar$Gene_symbol, combined_df_no_covar$Gene_symbol_non_hgnc)
  2106. combined_df_wth_covar$Corrected_symbol = mapply(function(x,y){
  2107. combined_symbols = c(x, y)
  2108. combined_symbols = combined_symbols[!is.na(combined_symbols)]
  2110. if (length(combined_symbols)<1){
  2111. return(NA)
  2112. }
  2114. combined_symbols = stri_split_fixed(combined_symbols, pattern = ";")
  2115. combined_symbols = unlist(combined_symbols)
  2116. combined_symbols = combined_symbols[combined_symbols!= ""]
  2117. combined_symbols = combined_symbols[combined_symbols!= " "]
  2118. combined_symbols = combined_symbols[combined_symbols!= "NA"]
  2120. verified_symbols = sapply(combined_symbols, function(x) harmoniz_genes_check_dict[x])
  2121. verified_symbols = unique(verified_symbols)
  2122. verified_symbols = paste0(verified_symbols, collapse = ";")
  2123. return(verified_symbols)
  2124. }, combined_df_wth_covar$Gene_symbol, combined_df_wth_covar$Gene_symbol_non_hgnc)
  2125. combined_df_cell_types$Corrected_symbol = mapply(function(x,y){
  2126. combined_symbols = c(x, y)
  2127. combined_symbols = combined_symbols[!is.na(combined_symbols)]
  2129. if (length(combined_symbols)<1){
  2130. return(NA)
  2131. }
  2133. combined_symbols = stri_split_fixed(combined_symbols, pattern = ";")
  2134. combined_symbols = unlist(combined_symbols)
  2135. combined_symbols = combined_symbols[combined_symbols!= ""]
  2136. combined_symbols = combined_symbols[combined_symbols!= " "]
  2137. combined_symbols = combined_symbols[combined_symbols!= "NA"]
  2139. verified_symbols = sapply(combined_symbols, function(x) harmoniz_genes_check_dict[x])
  2140. verified_symbols = unique(verified_symbols)
  2141. verified_symbols = paste0(verified_symbols, collapse = ";")
  2142. return(verified_symbols)
  2143. }, combined_df_cell_types$Gene_symbol, combined_df_cell_types$Gene_symbol_non_hgnc)
  2145. # Harmonizing symbols for SV dataset
  2146. combined_df_SV = do.call(rbind, data_list_SV)
  2147. dim(combined_df_no_covar)[1] == dim(combined_df_SV)[1] # TRUE -> All genes are already covered in harmoniz_genes_check_dict
  2148. combined_df_SV$Corrected_symbol = mapply(function(x,y){
  2149. combined_symbols = c(x, y)
  2150. combined_symbols = combined_symbols[!is.na(combined_symbols)]
  2152. if (length(combined_symbols)<1){
  2153. return(NA)
  2154. }
  2156. combined_symbols = stri_split_fixed(combined_symbols, pattern = ";")
  2157. combined_symbols = unlist(combined_symbols)
  2158. combined_symbols = combined_symbols[combined_symbols!= ""]
  2159. combined_symbols = combined_symbols[combined_symbols!= " "]
  2160. combined_symbols = combined_symbols[combined_symbols!= "NA"]
  2162. verified_symbols = sapply(combined_symbols, function(x) harmoniz_genes_check_dict[x])
  2163. verified_symbols = unique(verified_symbols)
  2164. verified_symbols = paste0(verified_symbols, collapse = ";")
  2165. return(verified_symbols)
  2166. }, combined_df_SV$Gene_symbol, combined_df_SV$Gene_symbol_non_hgnc)
  2168. ################### Correcting cell type df ###################
  2170. combined_df_cell_types$Cell_type = stri_replace_all_fixed(combined_df_cell_types$Tissue, pattern = "Dorsolateral prefrontal cortex (BA9): ", replacement = "")
  2171. combined_df_cell_types$Study_Cell = paste0(combined_df_cell_types$Study, ": ", combined_df_cell_types$Cell_type)
  2173. table(combined_df_cell_types$Cell_type)
  2174. "
  2175. Ast End ExN InN Mic Mix Oli OPC
  2176. 8510 6145 34875 26001 2093 1305 9025 6806
  2177. "
  2178. # Ast -astrocyte
  2179. # End - endothelial
  2180. # ExN - excitatory neuron
  2181. # InN - inhibitory neuron
  2182. # Mic - microglia
  2183. # Mix - mixed expression profile
  2184. # Oli - oligodendrocytes
  2185. # OPC - oligodendrocyte precursor cells
  2187. combined_df_cell_types_meta = combined_df_cell_types %>%
  2188. group_by(.,Study_Cell, Corrected_symbol) %>%
  2189. summarise(.,avgLog2FC = mean(logFC),
  2190. maxSE = max(SE))
  2191. combined_df_cell_types_meta = ungroup(combined_df_cell_types_meta)
  2192. combined_df_cell_types_meta = as.data.frame(combined_df_cell_types_meta)
  2194. ################### Cohort stats on individual datasets ###################
  2195. length(unique(combined_df_no_covar$Study))
  2196. length(unique(combined_df_wth_covar$Study))
  2197. length(unique(combined_df_SV$Study))
  2199. studies = unique(combined_df_no_covar$Study)
  2201. cohort_stats_vector_no_covariates = vector()
  2202. for (i in 1:length(studies)){
  2204. separator_large = "\n\n"
  2205. if (i == 1){
  2207. string_1 = paste0("Total associations without covariates: ", nrow(combined_df_no_covar))
  2208. string_2 = paste0("Total significant associations (nominal) without covariates: ", nrow(combined_df_no_covar[combined_df_no_covar$P.Value<=0.05,]))
  2209. string_3 = paste0("Average absolute effect size (log2FC) without covariates: ", mean(abs(combined_df_no_covar$logFC)))
  2210. string_4 = paste0("Average absolute effect size (log2FC) (nominal significance) without covariates: ", mean(abs(combined_df_no_covar[combined_df_no_covar$P.Value<=0.05,"logFC"])))
  2211. cohort_stats_vector_no_covariates = c(string_1, string_2, string_3, string_4, separator_large)
  2212. }
  2214. TMP_current_cohort = studies[i]
  2215. TMP_current_dataset = combined_df_no_covar[combined_df_no_covar$Study == TMP_current_cohort,]
  2216. TMP_study_type = unique(TMP_current_dataset$Tissue_type)
  2217. TMP_tech_type = unique(TMP_current_dataset$Technology)
  2219. main_header = paste(c(TMP_current_cohort, TMP_study_type, TMP_tech_type), collapse = "\t")
  2221. aggregated_stats = TMP_current_dataset %>%
  2222. group_by(.,Tissue) %>%
  2223. summarise(
  2224. total_associations = length(ID),
  2225. unique_gene_symbols_p_lt_0_05 = n_distinct(Gene_symbol[P.Value <= 0.05]),
  2226. unique_IDs_p_lt_0_05 = n_distinct(ID[P.Value <= 0.05]),
  2227. unique_gene_symbols_adj_p_lt_0_05 = n_distinct(Gene_symbol[adj.P.Val <= 0.05]),
  2228. unique_IDs_adj_p_lt_0_05 = n_distinct(ID[adj.P.Val <= 0.05]),
  2229. mean_abs_log2FC = mean(abs(logFC)),
  2230. mean_abs_log2FC_nomin_signif = mean(abs(logFC[P.Value <= 0.05]))
  2231. )
  2232. aggregated_stats = as.data.frame(aggregated_stats)
  2233. header = paste(colnames(aggregated_stats), collapse = "\t")
  2234. # char_aggregated_stats = apply(aggregated_stats, 1, function(x) paste(x, collapse = "\t"))
  2236. char_aggregated_stats = knitr::kable(aggregated_stats)
  2238. cohort_stats_vector_no_covariates = c(cohort_stats_vector_no_covariates, main_header, char_aggregated_stats, separator_large)
  2239. if (i == length(studies)){
  2240. cohort_stats_vector_no_covariates = c(cohort_stats_vector_no_covariates, separator_large)
  2241. }
  2244. }
  2245. writeLines(text = cohort_stats_vector_no_covariates, con = "individual_DE_stats_no_covar.txt")
  2247. cohort_stats_vector_with_covariates = vector()
  2248. for (i in 1:length(studies)){
  2250. separator_large = "\n\n"
  2251. if (i == 1){
  2253. string_1 = paste0("Total associations with covariates: ", nrow(combined_df_wth_covar))
  2254. string_2 = paste0("Total significant associations (nominal) with covariates: ", nrow(combined_df_wth_covar[combined_df_wth_covar$P.Value<0.05,]))
  2255. string_3 = paste0("Average absolute effect size (log2FC) with covariates: ", mean(abs(combined_df_wth_covar$logFC)))
  2256. string_4 = paste0("Average absolute effect size (log2FC) (nominal significance) with covariates: ", mean(abs(combined_df_wth_covar[combined_df_wth_covar$P.Value<0.05,"logFC"])))
  2257. cohort_stats_vector_with_covariates = c(string_1, string_2, string_3, string_4, separator_large)
  2258. }
  2260. TMP_current_cohort = studies[i]
  2261. TMP_current_dataset = combined_df_wth_covar[combined_df_wth_covar$Study == TMP_current_cohort,]
  2262. TMP_study_type = unique(TMP_current_dataset$Tissue_type)
  2263. TMP_tech_type = unique(TMP_current_dataset$Technology)
  2265. main_header = paste(c(TMP_current_cohort, TMP_study_type, TMP_tech_type), collapse = "\t")
  2267. aggregated_stats = TMP_current_dataset %>%
  2268. group_by(.,Tissue) %>%
  2269. summarise(
  2270. total_associations = length(ID),
  2271. unique_gene_symbols_p_lt_0_05 = n_distinct(Gene_symbol[P.Value < 0.05]),
  2272. unique_IDs_p_lt_0_05 = n_distinct(ID[P.Value < 0.05]),
  2273. unique_gene_symbols_adj_p_lt_0_05 = n_distinct(Gene_symbol[adj.P.Val < 0.05]),
  2274. unique_IDs_adj_p_lt_0_05 = n_distinct(ID[adj.P.Val < 0.05]),
  2275. mean_abs_log2FC = mean(abs(logFC)),
  2276. mean_abs_log2FC_nomin_signif = mean(abs(logFC[P.Value < 0.05]))
  2277. )
  2278. aggregated_stats = as.data.frame(aggregated_stats)
  2279. header = paste(colnames(aggregated_stats), collapse = "\t")
  2280. # char_aggregated_stats = apply(aggregated_stats, 1, function(x) paste(x, collapse = "\t"))
  2282. char_aggregated_stats = knitr::kable(aggregated_stats)
  2284. cohort_stats_vector_with_covariates = c(cohort_stats_vector_with_covariates, main_header, char_aggregated_stats, separator_large)
  2286. if (i == length(studies)){
  2287. cohort_stats_vector_with_covariates = c(cohort_stats_vector_with_covariates, separator_large)
  2288. }
  2289. }
  2290. writeLines(text = cohort_stats_vector_with_covariates, con = "individual_DE_stats_with_covar.txt")
  2293. cohort_stats_vector_SV = vector()
  2294. for (i in 1:length(studies)){
  2296. separator_large = "\n\n"
  2297. if (i == 1){
  2299. string_1 = paste0("Total associations using surrogate varaibles (SVs): ", nrow(combined_df_SV))
  2300. string_2 = paste0("Total significant associations (nominal) with SVs: ", nrow(combined_df_SV[combined_df_SV$P.Value<0.05,]))
  2301. string_3 = paste0("Average absolute effect size (log2FC) with SVs: ", mean(abs(combined_df_SV$logFC)))
  2302. string_4 = paste0("Average absolute effect size (log2FC) (nominal significance) with SVs: ", mean(abs(combined_df_SV[combined_df_SV$P.Value<0.05,"logFC"])))
  2303. cohort_stats_vector_SV = c(string_1, string_2, string_3, string_4, separator_large)
  2304. }
  2306. TMP_current_cohort = studies[i]
  2307. TMP_current_dataset = combined_df_SV[combined_df_SV$Study == TMP_current_cohort,]
  2308. TMP_study_type = unique(TMP_current_dataset$Tissue_type)
  2309. TMP_tech_type = unique(TMP_current_dataset$Technology)
  2311. main_header = paste(c(TMP_current_cohort, TMP_study_type, TMP_tech_type), collapse = "\t")
  2313. aggregated_stats = TMP_current_dataset %>%
  2314. group_by(.,Tissue) %>%
  2315. summarise(
  2316. total_associations = length(ID),
  2317. unique_gene_symbols_p_lt_0_05 = n_distinct(Gene_symbol[P.Value < 0.05]),
  2318. unique_IDs_p_lt_0_05 = n_distinct(ID[P.Value < 0.05]),
  2319. unique_gene_symbols_adj_p_lt_0_05 = n_distinct(Gene_symbol[adj.P.Val < 0.05]),
  2320. unique_IDs_adj_p_lt_0_05 = n_distinct(ID[adj.P.Val < 0.05]),
  2321. mean_abs_log2FC = mean(abs(logFC)),
  2322. mean_abs_log2FC_nomin_signif = mean(abs(logFC[P.Value < 0.05]))
  2323. )
  2324. aggregated_stats = as.data.frame(aggregated_stats)
  2325. header = paste(colnames(aggregated_stats), collapse = "\t")
  2326. # char_aggregated_stats = apply(aggregated_stats, 1, function(x) paste(x, collapse = "\t"))
  2328. char_aggregated_stats = knitr::kable(aggregated_stats)
  2330. cohort_stats_vector_SV = c(cohort_stats_vector_SV, main_header, char_aggregated_stats, separator_large)
  2332. if (i == length(studies)){
  2333. cohort_stats_vector_SV = c(cohort_stats_vector_SV, separator_large)
  2334. }
  2335. }
  2336. writeLines(text = cohort_stats_vector_SV, con = "individual_DE_stats_SV.txt")
  2338. ################### Cohort inspection for specificity ###################
  2339. table(combined_df_no_covar$Tissue, combined_df_no_covar$Study)
  2341. # GSE101521 -> dorsal lateral prefrontal cortex (Brodmann Area 9) "with an attempt to enrich for gray matter" "samples inevitably also contained a small amount of white matter"
  2342. # GSE102556 -> many regions, no specifics on gray/white
  2343. # GSE144136 -> single cell, "Frozen gray matter samples were dissected from Brodmann area 9 (dlPFC)"
  2344. # GSE199536 -> habenula. "4.4 Habenula In 1872 he described a small mass of grey matter", paper does not specify
  2345. # GSE202537 -> several striatal regions, no specifics on gray/white matter content
  2346. # GSE208338 -> no original paper, GEO: "Total RNA was isolated from ∼100 mg frozen gray matter using 1.0 ml TRIzol reagent"
  2347. # GSE213982 -> "Frozen histological grade samples of gray and white matter were dissected from the dlPFC (Brodmann Area 9)" -> mixed matters
  2348. # GSE243356 -> "The temporal pole, Brodmann Area 20, gray and white matter were dissected for the assays" -> mix of matters
  2349. # Since the dissected brain gyri fold in three dimensions, it is impossible to exclude all white matter while still taking the full thickness (~ 3 mm) of the cortical"
  2350. # GSE247998 -> blood
  2351. # GSE248260 -> ventral white matter (BA 47) -> gray matter was not provided
  2352. # GSE5388 -> DLPFC (BA9) no specifics on gray/white matter content
  2353. # GSE5389 -> OFC (BA11) no specifics on gray/white matter content
  2354. # GSE66937 -> no original paper, no specifics on gray/white matter content, has technical replicates!
  2355. # GSE92538_U133A -> predominantly gray matter, has samples on 2 arrays
  2356. # GSE92538_U133_PLUS2 -> predominantly gray matter, has samples on 2 arrays
  2358. ################### GSE92538 inspection ###################
  2360. GSE92538_U133A = read.csv("Data_preprocessing_analysis/GSE92538_U133A_results/GSE92538_pheno_curated.csv")
  2361. rownames(GSE92538_U133A) = GSE92538_U133A$X
  2362. GSE92538_U133A$X = NULL
  2363. nrow(GSE92538_U133A) # 97
  2365. GSE92538_U133_PLUS2 = read.csv("Data_preprocessing_analysis/GSE92538_U133_PLUS2_results/GSE92538_pheno_curated.csv")
  2366. rownames(GSE92538_U133_PLUS2) = GSE92538_U133_PLUS2$X
  2367. GSE92538_U133_PLUS2$X = NULL
  2368. nrow(GSE92538_U133_PLUS2) # 75
  2369. any(duplicated(GSE92538_U133A$SUBJECT_ID)) # FALSE
  2370. any(duplicated(GSE92538_U133_PLUS2$SUBJECT_ID)) # FALSE
  2372. length(intersect(GSE92538_U133A$SUBJECT_ID, GSE92538_U133_PLUS2$SUBJECT_ID))
  2373. # 0 samples are duplicated
  2376. ################### RRA (all brain, no covar, test) ###################
  2377. # This is a non-paramentic approach to aggregate ranked lists
  2378. combined_df_no_covar_RRA_all = combined_df_no_covar
  2379. combined_df_no_covar_RRA_all = combined_df_no_covar_RRA_all[combined_df_no_covar_RRA_all$Tissue_type == "Brain",]
  2380. table(combined_df_no_covar_RRA_all$Tissue, combined_df_no_covar_RRA_all$Study)
  2381. # GSE102556 many tissues
  2382. # GSE202537 many tissues
  2383. # GSE66937 many tissues
  2385. combined_df_no_covar_RRA_all_1 = combined_df_no_covar_RRA_all[combined_df_no_covar_RRA_all$Study %!in% c("GSE102556", "GSE202537", "GSE66937"), ]
  2386. combined_df_no_covar_RRA_all_2 = combined_df_no_covar_RRA_all[combined_df_no_covar_RRA_all$Study %in% c("GSE102556", "GSE202537", "GSE66937"), ]
  2387. # Caudate https://pubmed.ncbi.nlm.nih.gov/39164232/
  2388. # Nac https://pubmed.ncbi.nlm.nih.gov/38894648/ https://pubmed.ncbi.nlm.nih.gov/38965529/ https://pubmed.ncbi.nlm.nih.gov/39167467/
  2389. # GSE102556 -> Dorsolateral prefrontal cortex (dlPFC; BA8/9)
  2390. # GSE202537 -> Nac 22652 (more probes analyzed)
  2391. # GSE66937 -> prefrontal cortex
  2392. combined_df_no_covar_RRA_all_2 = combined_df_no_covar_RRA_all_2[combined_df_no_covar_RRA_all_2$Tissue %in% c("Dorsolateral prefrontal cortex (dlPFC; BA8/9)",
  2393. "Nac",
  2394. "prefrontal cortex"),]
  2395. table(combined_df_no_covar_RRA_all_2$Tissue)
  2396. "
  2397. Dorsolateral prefrontal cortex (dlPFC; BA8/9) Nac prefrontal cortex
  2398. 21365 22652 70523
  2399. "
  2400. combined_df_no_covar_RRA_all = rbind(combined_df_no_covar_RRA_all_1, combined_df_no_covar_RRA_all_2)
  2401. table(combined_df_no_covar_RRA_all$Tissue, combined_df_no_covar_RRA_all$Study)
  2402. combined_df_no_covar_RRA_all_list = lapply(unique(combined_df_no_covar_RRA_all$Study), function(x){
  2403. df = combined_df_no_covar_RRA_all[combined_df_no_covar_RRA_all$Study == x,]
  2404. return(df)
  2405. })
  2407. combined_df_no_covar_RRA_all_list_reduced = lapply(combined_df_no_covar_RRA_all_list, function(x){
  2409. orig_df = x
  2411. x = x[!stri_detect_fixed(x$Corrected_symbol, pattern = ";"),]
  2412. x = x[!is.na(x$Corrected_symbol),]
  2413. x = x[x$Corrected_symbol != "",]
  2414. x = x %>%
  2415. dplyr::group_by(., Corrected_symbol)%>%
  2416. dplyr::summarise(.,
  2417. avgLog2FC = mean(logFC),
  2418. maxP = max(P.Value)) # Most conservative
  2420. x = dplyr::ungroup(x)
  2421. x$Tissue = unique(orig_df$Tissue)
  2422. x$Tissue_type = unique(orig_df$Tissue_type)
  2423. x$Technology = unique(orig_df$Technology)
  2424. x$Study = unique(orig_df$Study)
  2425. x = as.data.frame(x)
  2426. x = dplyr::arrange(x, maxP)
  2427. return(x)
  2428. })
  2429. combined_df_no_covar_RRA_all_list_reduced = do.call(rbind, combined_df_no_covar_RRA_all_list_reduced)
  2430. total_number = length(unique(combined_df_no_covar_RRA_all_list_reduced$Corrected_symbol))
  2432. combined_df_no_covar_RRA_all_list_reduced_positive = combined_df_no_covar_RRA_all_list_reduced[combined_df_no_covar_RRA_all_list_reduced$avgLog2FC > 0,]
  2433. combined_df_no_covar_RRA_all_list_reduced_negative = combined_df_no_covar_RRA_all_list_reduced[combined_df_no_covar_RRA_all_list_reduced$avgLog2FC < 0,]
  2436. list_sizes = combined_df_no_covar_RRA_all_list_reduced_positive %>%
  2437. group_by(., Study) %>%
  2438. summarize(N=length(Corrected_symbol)) %>%
  2439. ungroup(.)
  2440. mean(list_sizes$N)
  2442. list_sizes = combined_df_no_covar_RRA_all_list_reduced_negative %>%
  2443. group_by(., Study) %>%
  2444. summarize(N=length(Corrected_symbol)) %>%
  2445. ungroup(.)
  2446. mean(list_sizes$N)
  2448. studies = unique(combined_df_no_covar_RRA_all_list_reduced_positive$Study)
  2449. no_covar_RRA_gene_lists_positive = lapply(studies, function(x){
  2450. tmp_df = combined_df_no_covar_RRA_all_list_reduced_positive[combined_df_no_covar_RRA_all_list_reduced_positive$Study == x,]
  2451. tmp_df = dplyr::arrange(tmp_df, maxP)
  2452. ranked_list = tmp_df$Corrected_symbol
  2453. return(ranked_list)
  2454. })
  2455. no_covar_RRA_gene_lists_positive_top = lapply(studies, function(x){
  2456. tmp_df = combined_df_no_covar_RRA_all_list_reduced_positive[combined_df_no_covar_RRA_all_list_reduced_positive$Study == x,]
  2457. tmp_df = dplyr::arrange(tmp_df, maxP)
  2458. tmp_df = tmp_df[tmp_df$maxP<0.05,]
  2459. ranked_list = tmp_df$Corrected_symbol
  2460. return(ranked_list)
  2461. })
  2463. studies = unique(combined_df_no_covar_RRA_all_list_reduced_negative$Study)
  2464. no_covar_RRA_gene_lists_negative = lapply(studies, function(x){
  2465. tmp_df = combined_df_no_covar_RRA_all_list_reduced_negative[combined_df_no_covar_RRA_all_list_reduced_negative$Study == x,]
  2466. tmp_df = dplyr::arrange(tmp_df, maxP)
  2467. ranked_list = tmp_df$Corrected_symbol
  2468. return(ranked_list)
  2469. })
  2470. no_covar_RRA_gene_lists_negative_top = lapply(studies, function(x){
  2471. tmp_df = combined_df_no_covar_RRA_all_list_reduced_negative[combined_df_no_covar_RRA_all_list_reduced_negative$Study == x,]
  2472. tmp_df = dplyr::arrange(tmp_df, maxP)
  2473. tmp_df = tmp_df[tmp_df$maxP<0.05,]
  2474. ranked_list = tmp_df$Corrected_symbol
  2475. return(ranked_list)
  2476. })
  2478. RRA_no_covar_all_brain_positive = aggregateRanks(glist = no_covar_RRA_gene_lists_positive, full = TRUE)
  2479. RRA_no_covar_all_brain_positive$P.adjusted = p.adjust(RRA_no_covar_all_brain_positive$Score, method = "fdr")
  2480. RRA_no_covar_all_brain_positive$N.cohorts = sapply(RRA_no_covar_all_brain_positive$Name, function(x){
  2481. internal_count = sapply(no_covar_RRA_gene_lists_positive, function(sublist) x %in% sublist)
  2482. internal_count = unlist(internal_count)
  2483. internal_count = as.numeric(internal_count)
  2484. internal_count = sum(internal_count)
  2485. return(internal_count)
  2486. })
  2488. RRA_no_covar_all_brain_negative = aggregateRanks(glist = no_covar_RRA_gene_lists_negative, full = TRUE)
  2489. RRA_no_covar_all_brain_negative$P.adjusted = p.adjust(RRA_no_covar_all_brain_negative$Score, method = "fdr")
  2490. RRA_no_covar_all_brain_negative$N.cohorts = sapply(RRA_no_covar_all_brain_negative$Name, function(x){
  2491. internal_count = sapply(no_covar_RRA_gene_lists_positive, function(sublist) x %in% sublist)
  2492. internal_count = unlist(internal_count)
  2493. internal_count = as.numeric(internal_count)
  2494. internal_count = sum(internal_count)
  2495. return(internal_count)
  2496. })
  2498. # Subset interesting genes
  2499. candidates = c("P2RY12", "AC159540.14", "CX3CR1", "DEPP1", "FOS", "GPR34", "KLF4", "PMP2", "SOX9")
  2500. RRA_no_covar_all_brain_positive[RRA_no_covar_all_brain_positive$Name %in% candidates,]
  2501. RRA_no_covar_all_brain_negative[RRA_no_covar_all_brain_negative$Name %in% candidates,]
  2503. # ---> Add RRA in the meta-calculation as columns
  2505. add_RRA_to_meta_df = function(meta_df, initial_study_df_list){
  2507. # Reduce df for RRA (need p-values)
  2508. print("Reducing df")
  2509. initial_study_df_list_reduced = lapply(initial_study_df_list, function(x){
  2510. orig_df = x
  2511. x = x[!stri_detect_fixed(x$Corrected_symbol, pattern = ";"),]
  2512. x = x[!is.na(x$Corrected_symbol),]
  2513. x = x[x$Corrected_symbol != "",]
  2514. x = x %>%
  2515. dplyr::group_by(., Corrected_symbol)%>%
  2516. dplyr::summarise(.,
  2517. avgLog2FC = mean(logFC),
  2518. maxP = max(P.Value)) # Most conservative
  2519. x = dplyr::ungroup(x)
  2520. x$Tissue = unique(orig_df$Tissue)
  2521. x$Tissue_type = unique(orig_df$Tissue_type)
  2522. x$Technology = unique(orig_df$Technology)
  2523. x$Study = unique(orig_df$Study)
  2524. x = as.data.frame(x)
  2525. x = dplyr::arrange(x, maxP)
  2526. return(x)
  2527. })
  2529. initial_study_df_list_reduced = do.call(rbind, initial_study_df_list_reduced)
  2530. RRA_all_list_reduced_positive = initial_study_df_list_reduced[initial_study_df_list_reduced$avgLog2FC > 0,]
  2531. RRA_all_list_reduced_negative = initial_study_df_list_reduced[initial_study_df_list_reduced$avgLog2FC < 0,]
  2533. # Prepare ranked lists
  2534. studies = unique(RRA_all_list_reduced_positive$Study)
  2535. RRA_gene_lists_positive = lapply(studies, function(x){
  2536. tmp_df = RRA_all_list_reduced_positive[RRA_all_list_reduced_positive$Study == x,]
  2537. tmp_df = dplyr::arrange(tmp_df, maxP)
  2538. ranked_list = tmp_df$Corrected_symbol
  2539. return(ranked_list)
  2540. })
  2541. RRA_gene_lists_negative = lapply(studies, function(x){
  2542. tmp_df = RRA_all_list_reduced_negative[RRA_all_list_reduced_negative$Study == x,]
  2543. tmp_df = dplyr::arrange(tmp_df, maxP)
  2544. ranked_list = tmp_df$Corrected_symbol
  2545. return(ranked_list)
  2546. })
  2548. # Run for positive log2FC
  2549. print("Running RRA")
  2550. RRA_all_brain_positive = aggregateRanks(glist = RRA_gene_lists_positive, full = TRUE)
  2551. RRA_all_brain_positive$P.adjusted.RRA = p.adjust(RRA_all_brain_positive$Score, method = "fdr")
  2552. RRA_all_brain_positive$N.cohorts.RRA = sapply(RRA_all_brain_positive$Name, function(x){
  2553. internal_count = sapply(RRA_gene_lists_positive, function(sublist) x %in% sublist)
  2554. internal_count = unlist(internal_count)
  2555. internal_count = as.numeric(internal_count)
  2556. internal_count = sum(internal_count)
  2557. return(internal_count)
  2558. })
  2559. colnames(RRA_all_brain_positive) = c("Name_RRA", "P_RRA", "FDR_RRA", "N_cohorts_RRA")
  2561. # Run for negative log2FC
  2562. RRA_all_brain_negative = aggregateRanks(glist = RRA_gene_lists_negative, full = TRUE)
  2563. RRA_all_brain_negative$P.adjusted.RRA = p.adjust(RRA_all_brain_negative$Score, method = "fdr")
  2564. RRA_all_brain_negative$N.cohorts.RRA = sapply(RRA_all_brain_negative$Name, function(x){
  2565. internal_count = sapply(RRA_gene_lists_negative, function(sublist) x %in% sublist)
  2566. internal_count = unlist(internal_count)
  2567. internal_count = as.numeric(internal_count)
  2568. internal_count = sum(internal_count)
  2569. return(internal_count)
  2570. })
  2571. colnames(RRA_all_brain_negative) = c("Name_RRA", "P_RRA", "FDR_RRA", "N_cohorts_RRA")
  2573. # Attach stats to meta_df
  2574. print("Preparing attachment")
  2575. attachment = mclapply(meta_df$gene, function(x){
  2577. selected_subset = meta_df[meta_df$gene == x,]
  2578. selected_subset = selected_subset[1,] # always ensure 1st row
  2580. meta_LFc = selected_subset$meta_LFc[1]
  2581. cohorts_up = selected_subset$cohorts_up[1]
  2582. cohorts_down = selected_subset$cohorts_down[1]
  2585. if (is.na(meta_LFc)){
  2586. # Number of cohorts less than 5
  2587. # Consider cohort counts only
  2589. if (cohorts_up >= cohorts_down){
  2590. # Positive analysis
  2591. P_RRA = RRA_all_brain_positive[RRA_all_brain_positive$Name_RRA == x, "P_RRA"]
  2592. FDR_RRA = RRA_all_brain_positive[RRA_all_brain_positive$Name_RRA == x, "FDR_RRA"]
  2593. N_cohorts_RRA = RRA_all_brain_positive[RRA_all_brain_positive$Name_RRA == x, "N_cohorts_RRA"]
  2594. } else {
  2595. # Negative analysis
  2596. P_RRA = RRA_all_brain_negative[RRA_all_brain_negative$Name_RRA == x, "P_RRA"]
  2597. FDR_RRA = RRA_all_brain_negative[RRA_all_brain_negative$Name_RRA == x, "FDR_RRA"]
  2598. N_cohorts_RRA = RRA_all_brain_negative[RRA_all_brain_negative$Name_RRA == x, "N_cohorts_RRA"]
  2599. }
  2600. } else {
  2602. if (meta_LFc >= 0){
  2603. # Positive analysis
  2604. P_RRA = RRA_all_brain_positive[RRA_all_brain_positive$Name_RRA == x, "P_RRA"]
  2605. FDR_RRA = RRA_all_brain_positive[RRA_all_brain_positive$Name_RRA == x, "FDR_RRA"]
  2606. N_cohorts_RRA = RRA_all_brain_positive[RRA_all_brain_positive$Name_RRA == x, "N_cohorts_RRA"]
  2607. } else {
  2608. # Negative analysis
  2609. P_RRA = RRA_all_brain_negative[RRA_all_brain_negative$Name_RRA == x, "P_RRA"]
  2610. FDR_RRA = RRA_all_brain_negative[RRA_all_brain_negative$Name_RRA == x, "FDR_RRA"]
  2611. N_cohorts_RRA = RRA_all_brain_negative[RRA_all_brain_negative$Name_RRA == x, "N_cohorts_RRA"]
  2612. }
  2613. }
  2615. att_df_row = data.frame(
  2616. P_RRA = P_RRA,
  2617. FDR_RRA = FDR_RRA,
  2618. N_cohorts_RRA = N_cohorts_RRA
  2619. )
  2621. return(att_df_row)
  2622. }, mc.cores = 9)
  2624. attachment = do.call(rbind, attachment)
  2625. meta_df = cbind(meta_df, attachment)
  2626. return(meta_df)
  2627. }
  2629. ################### Meta-functions ###################
  2631. run_meta_SJ = function(initial_study_list, blood_df_to_compare){
  2634. # Reducing df for meta-analysis
  2635. print("Aggregating Log2FC and SE per gene per study")
  2637. initial_study_list_reduced = lapply(initial_study_list, function(x){
  2639. orig_df = x
  2641. x = x[!stri_detect_fixed(x$Corrected_symbol, pattern = ";"),]
  2642. x = x[!is.na(x$Corrected_symbol),]
  2643. x = x[x$Corrected_symbol != "",]
  2644. x = x %>%
  2645. dplyr::group_by(., Corrected_symbol)%>%
  2646. dplyr::summarise(.,
  2647. avgLog2FC = mean(logFC),
  2648. maxSE = max(SE))
  2650. x = dplyr::ungroup(x)
  2651. x$Tissue = unique(orig_df$Tissue)
  2652. x$Tissue_type = unique(orig_df$Tissue_type)
  2653. x$Technology = unique(orig_df$Technology)
  2654. x$Study = unique(orig_df$Study)
  2655. x = as.data.frame(x)
  2656. return(x)
  2657. })
  2659. initial_study_list_reduced = do.call(rbind, initial_study_list_reduced)
  2660. unique_genes_meta = unique(initial_study_list_reduced$Corrected_symbol)
  2662. print("Total genes: ")
  2663. print(length(unique_genes_meta))
  2665. print("Running meta-analysis per gene")
  2667. meta_analysis_list = mclapply(unique_genes_meta, function(x){
  2670. # meta df
  2671. tmp_df_genes = initial_study_list_reduced[initial_study_list_reduced$Corrected_symbol == x,]
  2672. tmp_cohorts_up = length(tmp_df_genes$avgLog2FC[tmp_df_genes$avgLog2FC > 0])
  2673. tmp_cohorts_down = length(tmp_df_genes$avgLog2FC[tmp_df_genes$avgLog2FC < 0])
  2674. cohorts_total = nrow(tmp_df_genes)
  2676. # blood df
  2677. tmp_df_blood = blood_df_to_compare[blood_df_to_compare$Corrected_symbol == x, ]
  2679. if (nrow(tmp_df_blood) > 0){
  2681. tmp_blood_lfc = mean(tmp_df_blood$logFC)
  2683. if (all(tmp_df_blood$logFC < 0)){
  2684. tmp_lfc_blood_dir = "all negative"
  2685. } else {
  2686. tmp_lfc_blood_dir = "all positive"
  2687. }
  2689. if (all(tmp_df_blood$P.Value < 0.05)){
  2690. tmp_blood_signif = "significant"
  2691. } else {
  2692. tmp_blood_signif = "mixed/non-significant"
  2693. }
  2695. } else {
  2696. tmp_blood_lfc = NA
  2697. tmp_lfc_blood_dir = "Not detected/analyzed"
  2698. tmp_blood_signif = "Not detected/analyzed"
  2699. }
  2701. if (cohorts_total < 5){
  2702. tmp_output_df = data.frame(
  2703. gene = x,
  2704. detected_cohorts = cohorts_total,
  2705. cohorts_up = tmp_cohorts_up,
  2706. cohorts_down = tmp_cohorts_down,
  2707. commentary = "Number of cohorts < 5: gene is excluded",
  2708. meta_LFc = NA,
  2709. meta_se = NA,
  2710. meta_pval = NA,
  2711. tau2 = NA,
  2712. I2 = NA,
  2713. H2 = NA,
  2714. Q = NA,
  2715. Q.p = NA,
  2716. mean_blood_lfc = tmp_blood_lfc,
  2717. blood_dir = tmp_lfc_blood_dir,
  2718. blood_signif = tmp_blood_signif,
  2719. matching_with_blood = NA)
  2720. return(tmp_output_df)
  2721. }
  2723. tmp_meta_model = rma.uni(yi = avgLog2FC,
  2724. vi = maxSE^2,
  2725. data = tmp_df_genes,
  2726. method = "SJ",
  2727. weighted = TRUE)
  2729. logFCcompar = c(tmp_blood_lfc, tmp_meta_model$b)
  2731. if (any(is.na(logFCcompar))){
  2732. matching = NA
  2733. } else if (all(logFCcompar > 0)){
  2734. matching = "Match"
  2735. } else if (all(logFCcompar < 0)){
  2736. matching = "Match"
  2737. } else {
  2738. matching = "Mismatch"
  2739. }
  2741. tmp_output_df = data.frame(
  2742. gene = x,
  2743. detected_cohorts = cohorts_total,
  2744. cohorts_up = tmp_cohorts_up,
  2745. cohorts_down = tmp_cohorts_down,
  2746. commentary = "Number of cohorts >= 5: gene is analyzed",
  2747. meta_LFc = tmp_meta_model$b,
  2748. meta_se = tmp_meta_model$se,
  2749. meta_pval = tmp_meta_model$pval,
  2750. tau2 = tmp_meta_model$tau2,
  2751. I2 = tmp_meta_model$I2,
  2752. H2 = tmp_meta_model$H2,
  2753. Q = tmp_meta_model$QE,
  2754. Q.p = tmp_meta_model$QEp,
  2755. mean_blood_lfc = tmp_blood_lfc,
  2756. blood_dir = tmp_lfc_blood_dir,
  2757. blood_signif = tmp_blood_signif,
  2758. matching_with_blood = matching)
  2760. return(tmp_output_df)
  2762. }, mc.cores = 9)
  2764. meta_analysis_df = do.call(rbind, meta_analysis_list)
  2765. rownames(meta_analysis_df) = NULL
  2766. rm(list = ls(pattern = "tmp_"))
  2767. meta_analysis_df = dplyr::arrange(meta_analysis_df, meta_pval)
  2768. meta_analysis_df$meta_FDR = p.adjust(meta_analysis_df$meta_pval, method = "fdr")
  2769. meta_analysis_df = meta_analysis_df[,c("gene",
  2770. "detected_cohorts",
  2771. "cohorts_up",
  2772. "cohorts_down",
  2773. "commentary",
  2774. "meta_LFc",
  2775. "meta_se",
  2776. "meta_pval",
  2777. "meta_FDR",
  2778. "tau2",
  2779. "I2",
  2780. "H2" ,
  2781. "Q",
  2782. "Q.p",
  2783. "mean_blood_lfc",
  2784. "blood_dir",
  2785. "blood_signif",
  2786. "matching_with_blood")]
  2788. out_list = list()
  2789. out_list[[1]] = meta_analysis_df
  2790. out_list[[2]] = initial_study_list_reduced
  2791. names(out_list) = c("meta_analysis_df", "initial_study_list_reduced")
  2793. return(out_list)
  2794. }
  2796. run_single_gene_meta_SJ = function(gene,
  2797. initial_study_list){
  2800. # Reducing df for meta-analysis
  2801. print("Aggregating Log2FC and SE per gene per study")
  2803. initial_study_list_reduced = lapply(initial_study_list, function(x){
  2805. orig_df = x
  2807. x = x[!stri_detect_fixed(x$Corrected_symbol, pattern = ";"),]
  2808. x = x[!is.na(x$Corrected_symbol),]
  2809. x = x[x$Corrected_symbol != "",]
  2810. x = x %>%
  2811. dplyr::group_by(., Corrected_symbol)%>%
  2812. dplyr::summarise(.,
  2813. avgLog2FC = mean(logFC),
  2814. maxSE = max(SE))
  2816. x = dplyr::ungroup(x)
  2817. x$Tissue = unique(orig_df$Tissue)
  2818. x$Tissue_type = unique(orig_df$Tissue_type)
  2819. x$Technology = unique(orig_df$Technology)
  2820. x$Study = unique(orig_df$Study)
  2821. x = as.data.frame(x)
  2822. return(x)
  2823. })
  2825. initial_study_list_reduced = do.call(rbind, initial_study_list_reduced)
  2826. unique_genes_meta = unique(initial_study_list_reduced$Corrected_symbol)
  2828. print("Running meta-analysis for gene")
  2831. # meta df
  2832. tmp_df_genes = initial_study_list_reduced[initial_study_list_reduced$Corrected_symbol == gene,]
  2833. tmp_cohorts_up = length(tmp_df_genes$avgLog2FC[tmp_df_genes$avgLog2FC > 0])
  2834. tmp_cohorts_down = length(tmp_df_genes$avgLog2FC[tmp_df_genes$avgLog2FC < 0])
  2835. cohorts_total = nrow(tmp_df_genes)
  2838. if (cohorts_total < 5){
  2839. print("N cohorts < 5")
  2840. } else {
  2842. tmp_meta_model = rma.uni(yi = avgLog2FC,
  2843. vi = maxSE^2,
  2844. data = tmp_df_genes,
  2845. method = "SJ",
  2846. weighted = TRUE)
  2848. print(summary(tmp_meta_model))
  2849. }
  2850. }
  2852. run_meta_REML = function(initial_study_list, blood_df_to_compare){
  2855. # Reducing df for meta-analysis
  2856. print("Aggregating Log2FC and SE per gene per study")
  2858. initial_study_list_reduced = lapply(initial_study_list, function(x){
  2860. orig_df = x
  2862. x = x[!stri_detect_fixed(x$Corrected_symbol, pattern = ";"),]
  2863. x = x[!is.na(x$Corrected_symbol),]
  2864. x = x[x$Corrected_symbol != "",]
  2865. x = x %>%
  2866. dplyr::group_by(., Corrected_symbol)%>%
  2867. dplyr::summarise(.,
  2868. avgLog2FC = mean(logFC),
  2869. maxSE = max(SE))
  2871. x = dplyr::ungroup(x)
  2872. x$Tissue = unique(orig_df$Tissue)
  2873. x$Tissue_type = unique(orig_df$Tissue_type)
  2874. x$Technology = unique(orig_df$Technology)
  2875. x$Study = unique(orig_df$Study)
  2876. x = as.data.frame(x)
  2877. return(x)
  2878. })
  2880. initial_study_list_reduced = do.call(rbind, initial_study_list_reduced)
  2881. unique_genes_meta = unique(initial_study_list_reduced$Corrected_symbol)
  2883. print("Running meta-analysis per gene")
  2885. meta_analysis_list = mclapply(unique_genes_meta, function(x){
  2888. # meta df
  2889. tmp_df_genes = initial_study_list_reduced[initial_study_list_reduced$Corrected_symbol == x,]
  2890. tmp_cohorts_up = length(tmp_df_genes$avgLog2FC[tmp_df_genes$avgLog2FC > 0])
  2891. tmp_cohorts_down = length(tmp_df_genes$avgLog2FC[tmp_df_genes$avgLog2FC < 0])
  2892. cohorts_total = nrow(tmp_df_genes)
  2894. # blood df
  2895. tmp_df_blood = blood_df_to_compare[blood_df_to_compare$Corrected_symbol == x, ]
  2897. if (nrow(tmp_df_blood) > 0){
  2899. tmp_blood_lfc = mean(tmp_df_blood$logFC)
  2901. if (all(tmp_df_blood$logFC < 0)){
  2902. tmp_lfc_blood_dir = "all negative"
  2903. } else {
  2904. tmp_lfc_blood_dir = "all positive"
  2905. }
  2907. if (all(tmp_df_blood$P.Value < 0.05)){
  2908. tmp_blood_signif = "significant"
  2909. } else {
  2910. tmp_blood_signif = "mixed/non-significant"
  2911. }
  2913. } else {
  2914. tmp_blood_lfc = NA
  2915. tmp_lfc_blood_dir = "Not detected/analyzed"
  2916. tmp_blood_signif = "Not detected/analyzed"
  2917. }
  2919. if (cohorts_total < 5){
  2920. tmp_output_df = data.frame(
  2921. gene = x,
  2922. detected_cohorts = cohorts_total,
  2923. cohorts_up = tmp_cohorts_up,
  2924. cohorts_down = tmp_cohorts_down,
  2925. commentary = "Number of cohorts < 5: gene is excluded",
  2926. meta_LFc = NA,
  2927. meta_se = NA,
  2928. meta_pval = NA,
  2929. tau2 = NA,
  2930. I2 = NA,
  2931. H2 = NA,
  2932. Q = NA,
  2933. Q.p = NA,
  2934. mean_blood_lfc = tmp_blood_lfc,
  2935. blood_dir = tmp_lfc_blood_dir,
  2936. blood_signif = tmp_blood_signif,
  2937. matching_with_blood = NA)
  2938. return(tmp_output_df)
  2939. }
  2941. tmp_meta_model = rma.uni(yi = avgLog2FC,
  2942. vi = maxSE^2,
  2943. data = tmp_df_genes,
  2944. method = "REML",
  2945. control=list(stepadj=0.5, maxiter=10000),
  2946. test="knha", weighted = TRUE)
  2948. logFCcompar = c(tmp_blood_lfc, tmp_meta_model$b)
  2950. if (any(is.na(logFCcompar))){
  2951. matching = NA
  2952. } else if (all(logFCcompar > 0)){
  2953. matching = "Match"
  2954. } else if (all(logFCcompar < 0)){
  2955. matching = "Match"
  2956. } else {
  2957. matching = "Mismatch"
  2958. }
  2960. tmp_output_df = data.frame(
  2961. gene = x,
  2962. detected_cohorts = cohorts_total,
  2963. cohorts_up = tmp_cohorts_up,
  2964. cohorts_down = tmp_cohorts_down,
  2965. commentary = "Number of cohorts >= 5: gene is analyzed",
  2966. meta_LFc = tmp_meta_model$b,
  2967. meta_se = tmp_meta_model$se,
  2968. meta_pval = tmp_meta_model$pval,
  2969. tau2 = tmp_meta_model$tau2,
  2970. I2 = tmp_meta_model$I2,
  2971. H2 = tmp_meta_model$H2,
  2972. Q = tmp_meta_model$QE,
  2973. Q.p = tmp_meta_model$QEp,
  2974. mean_blood_lfc = tmp_blood_lfc,
  2975. blood_dir = tmp_lfc_blood_dir,
  2976. blood_signif = tmp_blood_signif,
  2977. matching_with_blood = matching)
  2979. return(tmp_output_df)
  2981. }, mc.cores = 9)
  2983. meta_analysis_df = do.call(rbind, meta_analysis_list)
  2984. rownames(meta_analysis_df) = NULL
  2985. rm(list = ls(pattern = "tmp_"))
  2986. meta_analysis_df = dplyr::arrange(meta_analysis_df, meta_pval)
  2987. meta_analysis_df$meta_FDR = p.adjust(meta_analysis_df$meta_pval, method = "fdr")
  2988. meta_analysis_df = meta_analysis_df[,c("gene",
  2989. "detected_cohorts",
  2990. "cohorts_up",
  2991. "cohorts_down",
  2992. "commentary",
  2993. "meta_LFc",
  2994. "meta_se",
  2995. "meta_pval",
  2996. "meta_FDR",
  2997. "tau2",
  2998. "I2",
  2999. "H2" ,
  3000. "Q",
  3001. "Q.p",
  3002. "mean_blood_lfc",
  3003. "blood_dir",
  3004. "blood_signif",
  3005. "matching_with_blood")]
  3007. out_list = list()
  3008. out_list[[1]] = meta_analysis_df
  3009. out_list[[2]] = initial_study_list_reduced
  3010. names(out_list) = c("meta_analysis_df", "initial_study_list_reduced")
  3012. return(out_list)
  3013. }
  3015. ################### Meta (no covar) all brain data ###################
  3016. combined_df_no_covar_meta_all = combined_df_no_covar
  3017. combined_df_no_covar_meta_all = combined_df_no_covar_meta_all[combined_df_no_covar_meta_all$Tissue_type == "Brain",]
  3018. blood_df_no_covar = combined_df_no_covar[combined_df_no_covar$Tissue_type == "Blood",]
  3020. table(combined_df_no_covar_meta_all$Tissue, combined_df_no_covar_meta_all$Study)
  3021. # GSE102556 many tissues
  3022. # GSE202537 many tissues
  3023. # GSE66937 many tissues
  3025. combined_df_no_covar_meta_all_1 = combined_df_no_covar_meta_all[combined_df_no_covar_meta_all$Study %!in% c("GSE102556", "GSE202537", "GSE66937"), ]
  3026. combined_df_no_covar_meta_all_2 = combined_df_no_covar_meta_all[combined_df_no_covar_meta_all$Study %in% c("GSE102556", "GSE202537", "GSE66937"), ]
  3027. # Caudate https://pubmed.ncbi.nlm.nih.gov/39164232/
  3028. # Nac https://pubmed.ncbi.nlm.nih.gov/38894648/ https://pubmed.ncbi.nlm.nih.gov/38965529/ https://pubmed.ncbi.nlm.nih.gov/39167467/
  3029. # GSE102556 -> Dorsolateral prefrontal cortex (dlPFC; BA8/9)
  3030. # GSE202537 -> Nac 22652 (more probes analyzed)
  3031. # GSE66937 -> prefrontal cortex
  3032. combined_df_no_covar_meta_all_2 = combined_df_no_covar_meta_all_2[combined_df_no_covar_meta_all_2$Tissue %in% c("Dorsolateral prefrontal cortex (dlPFC; BA8/9)",
  3033. "Nac",
  3034. "prefrontal cortex"),]
  3035. table(combined_df_no_covar_meta_all_2$Tissue)
  3036. "
  3037. Dorsolateral prefrontal cortex (dlPFC; BA8/9) Nac prefrontal cortex
  3038. 21365 22652 70523
  3039. "
  3040. combined_df_no_covar_meta_full = rbind(combined_df_no_covar_meta_all_1, combined_df_no_covar_meta_all_2)
  3041. table(combined_df_no_covar_meta_full$Tissue, combined_df_no_covar_meta_full$Study)
  3042. combined_df_no_covar_meta_full_list = lapply(unique(combined_df_no_covar_meta_full$Study), function(x){
  3043. df = combined_df_no_covar_meta_full[combined_df_no_covar_meta_full$Study == x,]
  3044. return(df)
  3045. })
  3047. # Run meta-analysis
  3048. analysis_container_no_covar_all_brain = run_meta_SJ(initial_study_list = combined_df_no_covar_meta_full_list,
  3049. blood_df_to_compare = blood_df_no_covar)
  3050. meta_no_covar_all_brain = analysis_container_no_covar_all_brain[["meta_analysis_df"]]
  3051. combined_df_no_covar_meta_full_list_reduced = analysis_container_no_covar_all_brain[["initial_study_list_reduced"]]
  3053. meta_no_covar_all_brain = add_RRA_to_meta_df(meta_no_covar_all_brain, combined_df_no_covar_meta_full_list)
  3054. non_na_substet = meta_no_covar_all_brain[!is.na(meta_no_covar_all_brain$meta_LFc), ]
  3055. cor.test(non_na_substet$meta_pval, non_na_substet$P_RRA, method = "spearman") # 0.1461325
  3056. meta_no_covar_all_brain_signif = meta_no_covar_all_brain[meta_no_covar_all_brain$meta_pval < 0.05,]
  3057. meta_no_covar_all_brain_signif = meta_no_covar_all_brain_signif[!is.na(meta_no_covar_all_brain_signif$meta_pval),]
  3059. nrow(meta_no_covar_all_brain_signif)
  3060. # P2RY12
  3061. # https://pubmed.ncbi.nlm.nih.gov/37396924/
  3062. # https://www.ingentaconnect.com/content/ben/cnr/2022/00000019/00000003/art00004
  3063. # https://pubmed.ncbi.nlm.nih.gov/31739114/
  3064. # https://www.ncbi.nlm.nih.gov/pmc/articles/PMC8908792/
  3065. # https://www.mdpi.com/1422-0067/25/11/5750
  3066. # https://psychiatryonline.org/doi/10.1176/appi.ajp.2021.22010026
  3068. plot_forest_simple_gene(gene_name = "P2RY12", meta_df = combined_df_cell_types_meta, study_column_char = "Study_Cell")
  3069. plot_forest_meta_gene(gene_name = "P2RY12", meta_df = combined_df_no_covar_meta_full_list_reduced)
  3071. # making plots
  3072. dir.create("forest_plots")
  3073. folder = "forest_plots/meta_no_covar_all_brain"
  3074. dir.create(folder)
  3076. for (i in 1:nrow(meta_no_covar_all_brain_signif)){
  3079. gene = meta_no_covar_all_brain_signif$gene[i]
  3080. plot_path = paste0(folder, "/", gene, "_forest.pdf")
  3082. pdf(file = plot_path, width = 10, height = 7)
  3083. print(plot_forest_meta_gene(gene_name = gene, meta_df = combined_df_no_covar_meta_full_list_reduced))
  3084. dev.off()
  3086. }
  3088. # Heterogeneity summary (28644 were excluded from meta due to Ncohorts <5)
  3089. summary(meta_no_covar_all_brain$I2)
  3090. "
  3091. Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
  3092. 0.059 34.431 49.216 48.750 63.528 98.581 28644
  3093. "
  3095. summary(meta_no_covar_all_brain_signif$I2)
  3096. "
  3097. Min. 1st Qu. Median Mean 3rd Qu. Max.
  3098. 0.8334 20.4189 32.9392 34.0485 45.3801 78.5286
  3099. "
  3101. ################### Meta (with covar) all brain data ###################
  3102. combined_df_with_covar_meta_all = combined_df_wth_covar
  3103. combined_df_with_covar_meta_all = combined_df_with_covar_meta_all[combined_df_with_covar_meta_all$Tissue_type == "Brain",]
  3104. blood_df_with_covar = combined_df_wth_covar[combined_df_wth_covar$Tissue_type == "Blood",]
  3106. table(combined_df_with_covar_meta_all$Tissue, combined_df_with_covar_meta_all$Study)
  3108. # GSE102556 many tissues
  3109. # GSE202537 many tissues
  3110. # GSE66937 many tissues
  3112. combined_df_with_covar_meta_all_1 = combined_df_with_covar_meta_all[combined_df_with_covar_meta_all$Study %!in% c("GSE102556", "GSE202537", "GSE66937"), ]
  3113. combined_df_with_covar_meta_all_2 = combined_df_with_covar_meta_all[combined_df_with_covar_meta_all$Study %in% c("GSE102556", "GSE202537", "GSE66937"), ]
  3114. combined_df_with_covar_meta_all_2 = combined_df_with_covar_meta_all_2[combined_df_with_covar_meta_all_2$Tissue %in% c("Dorsolateral prefrontal cortex (dlPFC; BA8/9)",
  3115. "Nac",
  3116. "prefrontal cortex"),]
  3117. table(combined_df_with_covar_meta_all_2$Tissue)
  3118. "
  3119. Dorsolateral prefrontal cortex (dlPFC; BA8/9) Nac prefrontal cortex
  3120. 23331 23919 70523
  3121. "
  3123. combined_df_with_covar_meta_full = rbind(combined_df_with_covar_meta_all_1, combined_df_with_covar_meta_all_2)
  3124. table(combined_df_with_covar_meta_full$Tissue, combined_df_with_covar_meta_full$Study)
  3125. combined_df_with_covar_meta_full_list = lapply(unique(combined_df_with_covar_meta_full$Study), function(x){
  3126. df = combined_df_with_covar_meta_full[combined_df_with_covar_meta_full$Study == x,]
  3127. return(df)
  3128. })
  3130. # Run meta-analysis
  3131. analysis_container_with_covar_all_brain = run_meta_SJ(initial_study_list = combined_df_with_covar_meta_full_list,
  3132. blood_df_to_compare = blood_df_with_covar)
  3136. meta_with_covar_all_brain = analysis_container_with_covar_all_brain[["meta_analysis_df"]]
  3137. combined_df_with_covar_meta_full_list_reduced = analysis_container_with_covar_all_brain[["initial_study_list_reduced"]]
  3139. meta_with_covar_all_brain = add_RRA_to_meta_df(meta_with_covar_all_brain, combined_df_with_covar_meta_full_list)
  3140. meta_with_covar_all_brain_signif = meta_with_covar_all_brain[meta_with_covar_all_brain$meta_pval < 0.05,]
  3141. meta_with_covar_all_brain_signif = meta_with_covar_all_brain_signif[!is.na(meta_with_covar_all_brain_signif$meta_pval),]
  3142. nrow(meta_with_covar_all_brain_signif)
  3144. plot_forest_meta_gene(gene_name = "RB1", meta_df = combined_df_with_covar_meta_full_list_reduced)
  3145. plot_forest_meta_gene(gene_name = "P2RY12", meta_df = combined_df_with_covar_meta_full_list_reduced)
  3147. # Overlap full analysis
  3148. overlapping_genes_meta_full = intersect(meta_no_covar_all_brain_signif$gene, meta_with_covar_all_brain_signif$gene)
  3149. overlapping_genes_meta_full = unique(overlapping_genes_meta_full)
  3150. length(overlapping_genes_meta_full) # 57 genes overlap
  3151. overlapping_subset_allbrain_ncv = meta_no_covar_all_brain_signif[meta_no_covar_all_brain_signif$gene %in% overlapping_genes_meta_full, ]
  3152. overlapping_subset_allbrain_ncv = dplyr::arrange(overlapping_subset_allbrain_ncv, gene)
  3153. overlapping_subset_allbrain_cv = meta_with_covar_all_brain_signif[meta_with_covar_all_brain_signif$gene %in% overlapping_genes_meta_full, ]
  3154. overlapping_subset_allbrain_cv = dplyr::arrange(overlapping_subset_allbrain_cv, gene)
  3156. overlapping_subset_allbrain_stat = cbind(
  3157. overlapping_subset_allbrain_ncv$gene,
  3158. overlapping_subset_allbrain_ncv$meta_LFc,
  3159. overlapping_subset_allbrain_cv$gene,
  3160. overlapping_subset_allbrain_cv$meta_LFc
  3161. )
  3162. overlapping_subset_allbrain_stat = as.data.frame(overlapping_subset_allbrain_stat)
  3163. overlapping_subset_allbrain_stat$V2 = as.numeric(overlapping_subset_allbrain_stat$V2)
  3164. overlapping_subset_allbrain_stat$V4 = as.numeric(overlapping_subset_allbrain_stat$V4)
  3165. overlapping_subset_allbrain_stat$mean = mapply(function(x,y) mean(x,y),overlapping_subset_allbrain_stat$V2, overlapping_subset_allbrain_stat$V4)
  3167. # making plots
  3168. dir.create("forest_plots")
  3169. folder = "forest_plots/meta_with_covar_all_brain"
  3170. dir.create(folder)
  3172. for (i in 1:nrow(meta_with_covar_all_brain_signif)){
  3175. gene = meta_with_covar_all_brain_signif$gene[i]
  3176. plot_path = paste0(folder, "/", gene, "_forest.pdf")
  3178. pdf(file = plot_path, width = 10, height = 7)
  3179. print(plot_forest_meta_gene(gene_name = gene, meta_df = combined_df_with_covar_meta_full_list_reduced))
  3180. dev.off()
  3182. }
  3184. # Heterogeneity estimates
  3185. summary(meta_with_covar_all_brain$I2)
  3187. " Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
  3188. 0.001 25.237 39.893 41.247 56.093 99.471 29942
  3189. "
  3190. ################### Meta cortical regions (no covar) ###################
  3191. combined_df_no_covar_meta_cortical = combined_df_no_covar
  3192. combined_df_no_covar_meta_cortical = combined_df_no_covar_meta_cortical[combined_df_no_covar_meta_cortical$Tissue_type == "Brain",]
  3193. blood_df_no_covar = combined_df_no_covar[combined_df_no_covar$Tissue_type == "Blood",]
  3195. table(combined_df_no_covar_meta_cortical$Tissue, combined_df_no_covar_meta_cortical$Study)
  3196. # Tissues to include
  3197. # Cingulate gyrus 25 (Cg25) is also cortical but Dorsolateral prefrontal cortex (dlPFC; BA8/9) was used instead
  3198. combined_df_no_covar_meta_cortical = combined_df_no_covar_meta_cortical[combined_df_no_covar_meta_cortical$Tissue %in%
  3199. c("DLPFC",
  3200. "dorsal lateral prefrontal cortex (Brodmann Area 9)",
  3201. "Dorsolateral Prefrontal Cortex",
  3202. "Dorsolateral prefrontal cortex (BA9)",
  3203. "Dorsolateral prefrontal cortex (dlPFC; BA8/9)",
  3204. "Orbitofrontal cortex",
  3205. "postmortem dorsolateral prefrontal cortex (DLPFC)",
  3206. "prefrontal cortex",
  3207. "temporal cortex (BA20 and BA36)"),]
  3208. table(combined_df_no_covar_meta_cortical$Tissue, combined_df_no_covar_meta_cortical$Study)
  3210. "
  3211. GSE101521 GSE102556 GSE144136 GSE208338 GSE213982 GSE243356 GSE5388 GSE5389 GSE66937 GSE92538_U133_PLUS2 GSE92538_U133A
  3212. DLPFC 0 0 0 0 0 0 22283 0 0 0 0
  3213. dorsal lateral prefrontal cortex (Brodmann Area 9) 18677 0 0 0 0 0 0 0 0 0 0
  3214. Dorsolateral Prefrontal Cortex 0 0 0 0 0 0 0 0 0 54675 22283
  3215. Dorsolateral prefrontal cortex (BA9) 0 0 19271 0 19191 0 0 0 0 0 0
  3216. Dorsolateral prefrontal cortex (dlPFC; BA8/9) 0 21365 0 0 0 0 0 0 0 0 0
  3217. Orbitofrontal cortex 0 0 0 0 0 0 0 22283 0 0 0
  3218. postmortem dorsolateral prefrontal cortex (DLPFC) 0 0 0 750011 0 0 0 0 0 0 0
  3219. prefrontal cortex 0 0 0 0 0 0 0 0 70523 0 0
  3220. temporal cortex (BA20 and BA36) 0 0 0 0 0 20955 0 0 0 0 0
  3222. "
  3223. combined_df_no_covar_meta_cortical_list = lapply(unique(combined_df_no_covar_meta_cortical$Study), function(x){
  3224. df = combined_df_no_covar_meta_cortical[combined_df_no_covar_meta_cortical$Study == x,]
  3225. return(df)
  3226. })
  3228. # Run meta-analysis
  3229. analysis_container_no_covar_cortical = run_meta_SJ(initial_study_list = combined_df_no_covar_meta_cortical_list,
  3230. blood_df_to_compare = blood_df_no_covar)
  3234. meta_no_covar_cortical = analysis_container_no_covar_cortical[["meta_analysis_df"]]
  3235. combined_df_no_covar_meta_cortical_list_reduced = analysis_container_no_covar_cortical[["initial_study_list_reduced"]]
  3238. meta_no_covar_cortical = add_RRA_to_meta_df(meta_no_covar_cortical, combined_df_no_covar_meta_cortical_list)
  3239. meta_no_covar_cortical_signif = meta_no_covar_cortical[meta_no_covar_cortical$meta_pval < 0.05,]
  3240. meta_no_covar_cortical_signif = meta_no_covar_cortical_signif[!is.na(meta_no_covar_cortical_signif$meta_pval),]
  3241. meta_no_covar_cortical[meta_no_covar_cortical$gene == "P2RY12",]
  3242. nrow(meta_no_covar_cortical_signif)
  3244. # making plots
  3245. dir.create("forest_plots")
  3246. folder = "forest_plots/meta_no_covar_cortical"
  3247. dir.create(folder)
  3249. for (i in 1:nrow(meta_no_covar_cortical_signif)){
  3252. gene = meta_no_covar_cortical_signif$gene[i]
  3253. plot_path = paste0(folder, "/", gene, "_forest.pdf")
  3255. pdf(file = plot_path, width = 10, height = 7)
  3256. print(plot_forest_meta_gene(gene_name = gene, meta_df = combined_df_no_covar_meta_cortical_list_reduced))
  3257. dev.off()
  3259. }
  3261. # Heterogeneity summary (31325 were excluded from meta due to Ncohorts <5)
  3262. summary(meta_no_covar_cortical$I2)
  3263. "
  3264. Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
  3265. 0.021 31.841 48.919 48.347 65.182 98.987 31325
  3266. "
  3267. summary(meta_no_covar_cortical_signif$I2)
  3268. "
  3269. Min. 1st Qu. Median Mean 3rd Qu. Max.
  3270. 1.215 22.009 34.357 35.388 48.552 83.960
  3271. "
  3273. ################### Meta cortical regions (with covar) ###################
  3274. combined_df_with_covar_meta_cortical = combined_df_wth_covar
  3275. combined_df_with_covar_meta_cortical = combined_df_with_covar_meta_cortical[combined_df_with_covar_meta_cortical$Tissue_type == "Brain",]
  3276. blood_df_with_covar = combined_df_wth_covar[combined_df_wth_covar$Tissue_type == "Blood",]
  3278. # Tissues to include
  3279. combined_df_with_covar_meta_cortical = combined_df_with_covar_meta_cortical[combined_df_with_covar_meta_cortical$Tissue %in%
  3280. c("DLPFC",
  3281. "dorsal lateral prefrontal cortex (Brodmann Area 9)",
  3282. "Dorsolateral Prefrontal Cortex",
  3283. "Dorsolateral prefrontal cortex (BA9)",
  3284. "Dorsolateral prefrontal cortex (dlPFC; BA8/9)",
  3285. "Orbitofrontal cortex",
  3286. "postmortem dorsolateral prefrontal cortex (DLPFC)",
  3287. "prefrontal cortex",
  3288. "temporal cortex (BA20 and BA36)"),]
  3289. table(combined_df_with_covar_meta_cortical$Tissue, combined_df_with_covar_meta_cortical$Study)
  3291. "
  3292. GSE101521 GSE102556 GSE144136 GSE208338 GSE213982 GSE243356 GSE5388 GSE5389 GSE66937 GSE92538_U133_PLUS2 GSE92538_U133A
  3293. DLPFC 0 0 0 0 0 0 22283 0 0 0 0
  3294. dorsal lateral prefrontal cortex (Brodmann Area 9) 21197 0 0 0 0 0 0 0 0 0 0
  3295. Dorsolateral Prefrontal Cortex 0 0 0 0 0 0 0 0 0 54675 22283
  3296. Dorsolateral prefrontal cortex (BA9) 0 0 22389 0 24182 0 0 0 0 0 0
  3297. Dorsolateral prefrontal cortex (dlPFC; BA8/9) 0 23331 0 0 0 0 0 0 0 0 0
  3298. Orbitofrontal cortex 0 0 0 0 0 0 0 22283 0 0 0
  3299. postmortem dorsolateral prefrontal cortex (DLPFC) 0 0 0 750011 0 0 0 0 0 0 0
  3300. prefrontal cortex 0 0 0 0 0 0 0 0 70523 0 0
  3301. temporal cortex (BA20 and BA36) 0 0 0 0 0 20955 0 0 0 0 0
  3303. "
  3305. combined_df_with_covar_meta_cortical_list = lapply(unique(combined_df_with_covar_meta_cortical$Study), function(x){
  3306. df = combined_df_with_covar_meta_cortical[combined_df_with_covar_meta_cortical$Study == x,]
  3307. return(df)
  3308. })
  3310. # Run meta-analysis
  3311. analysis_container_with_covar_cortical = run_meta_SJ(initial_study_list = combined_df_with_covar_meta_cortical_list,
  3312. blood_df_to_compare = blood_df_with_covar)
  3316. meta_with_covar_cortical = analysis_container_with_covar_cortical[["meta_analysis_df"]]
  3317. combined_df_with_covar_meta_cortical_list_reduced = analysis_container_with_covar_cortical[["initial_study_list_reduced"]]
  3320. meta_with_covar_cortical = add_RRA_to_meta_df(meta_with_covar_cortical, combined_df_with_covar_meta_cortical_list)
  3321. meta_with_covar_cortical_signif = meta_with_covar_cortical[meta_with_covar_cortical$meta_pval < 0.05,]
  3322. meta_with_covar_cortical_signif = meta_with_covar_cortical_signif[!is.na(meta_with_covar_cortical_signif$meta_pval),]
  3323. nrow(meta_with_covar_cortical_signif) # 205
  3324. meta_with_covar_cortical[meta_with_covar_cortical$gene == "P2RY12",]
  3326. # making plots
  3327. dir.create("forest_plots")
  3328. folder = "forest_plots/meta_with_covar_cortical"
  3329. dir.create(folder)
  3331. for (i in 1:nrow(meta_with_covar_cortical_signif)){
  3334. gene = meta_with_covar_cortical_signif$gene[i]
  3335. plot_path = paste0(folder, "/", gene, "_forest.pdf")
  3337. pdf(file = plot_path, width = 10, height = 7)
  3338. print(plot_forest_meta_gene(gene_name = gene, meta_df = combined_df_with_covar_meta_cortical_list_reduced))
  3339. dev.off()
  3341. }
  3343. # Heterogeneity estimates
  3344. summary(meta_with_covar_cortical$I2)
  3346. " Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
  3347. 0.03 22.37 37.90 39.75 55.92 99.57 33360
  3348. "
  3351. ################### Meta prefrontal regions (no covar) ###################
  3352. combined_df_no_covar_meta_prefrontal = combined_df_no_covar
  3353. combined_df_no_covar_meta_prefrontal = combined_df_no_covar_meta_prefrontal[combined_df_no_covar_meta_prefrontal$Tissue_type == "Brain",]
  3354. blood_df_no_covar = combined_df_no_covar[combined_df_no_covar$Tissue_type == "Blood",]
  3356. table(combined_df_no_covar_meta_prefrontal$Tissue, combined_df_no_covar_meta_prefrontal$Study)
  3357. # Tissues to include
  3358. combined_df_no_covar_meta_prefrontal = combined_df_no_covar_meta_prefrontal[combined_df_no_covar_meta_prefrontal$Tissue %in%
  3359. c("DLPFC",
  3360. "dorsal lateral prefrontal cortex (Brodmann Area 9)",
  3361. "Dorsolateral Prefrontal Cortex",
  3362. "Dorsolateral prefrontal cortex (BA9)",
  3363. "Dorsolateral prefrontal cortex (dlPFC; BA8/9)",
  3364. "postmortem dorsolateral prefrontal cortex (DLPFC)",
  3365. "prefrontal cortex"),]
  3366. table(combined_df_no_covar_meta_prefrontal$Tissue, combined_df_no_covar_meta_prefrontal$Study)
  3368. "
  3369. GSE101521 GSE102556 GSE144136 GSE208338 GSE213982 GSE5388 GSE66937 GSE92538_U133_PLUS2 GSE92538_U133A
  3370. DLPFC 0 0 0 0 0 22283 0 0 0
  3371. dorsal lateral prefrontal cortex (Brodmann Area 9) 18677 0 0 0 0 0 0 0 0
  3372. Dorsolateral Prefrontal Cortex 0 0 0 0 0 0 0 54675 22283
  3373. Dorsolateral prefrontal cortex (BA9) 0 0 19271 0 19191 0 0 0 0
  3374. Dorsolateral prefrontal cortex (dlPFC; BA8/9) 0 21365 0 0 0 0 0 0 0
  3375. postmortem dorsolateral prefrontal cortex (DLPFC) 0 0 0 750011 0 0 0 0 0
  3376. prefrontal cortex 0 0 0 0 0 0 70523 0 0
  3378. "
  3379. combined_df_no_covar_meta_prefrontal_list = lapply(unique(combined_df_no_covar_meta_prefrontal$Study), function(x){
  3380. df = combined_df_no_covar_meta_prefrontal[combined_df_no_covar_meta_prefrontal$Study == x,]
  3381. return(df)
  3382. })
  3384. # Run meta-analysis
  3385. analysis_container_no_covar_prefrontal = run_meta_SJ(initial_study_list = combined_df_no_covar_meta_prefrontal_list,
  3386. blood_df_to_compare = blood_df_no_covar)
  3390. meta_no_covar_prefrontal = analysis_container_no_covar_prefrontal[["meta_analysis_df"]]
  3391. combined_df_no_covar_meta_prefrontal_list_reduced = analysis_container_no_covar_prefrontal[["initial_study_list_reduced"]]
  3394. meta_no_covar_prefrontal = add_RRA_to_meta_df(meta_no_covar_prefrontal, combined_df_no_covar_meta_prefrontal_list)
  3395. meta_no_covar_prefrontal_signif = meta_no_covar_prefrontal[meta_no_covar_prefrontal$meta_pval < 0.05,]
  3396. meta_no_covar_prefrontal_signif = meta_no_covar_prefrontal_signif[!is.na(meta_no_covar_prefrontal_signif$meta_pval),]
  3397. nrow(meta_no_covar_prefrontal_signif)
  3398. meta_no_covar_prefrontal[meta_no_covar_prefrontal$gene == "P2RY12",]
  3400. # making plots
  3401. dir.create("forest_plots")
  3402. folder = "forest_plots/meta_no_covar_prefrontal"
  3403. dir.create(folder)
  3405. for (i in 1:nrow(meta_no_covar_prefrontal_signif)){
  3408. gene = meta_no_covar_prefrontal_signif$gene[i]
  3409. plot_path = paste0(folder, "/", gene, "_forest.pdf")
  3411. pdf(file = plot_path, width = 10, height = 7)
  3412. print(plot_forest_meta_gene(gene_name = gene, meta_df = combined_df_no_covar_meta_prefrontal_list_reduced))
  3413. dev.off()
  3415. }
  3417. # Heterogeneity summary (33247 were excluded from meta due to Ncohorts <5)
  3418. summary(meta_no_covar_prefrontal$I2)
  3419. "
  3420. Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
  3421. 0.00 30.27 48.05 47.93 65.83 99.22 33247
  3422. "
  3424. summary(meta_no_covar_prefrontal_signif$I2)
  3425. "
  3426. Min. 1st Qu. Median Mean 3rd Qu. Max.
  3427. 0.28 21.44 35.29 37.44 50.67 90.73
  3428. "
  3430. ################### Meta prefrontal regions (with covar) ###################
  3431. combined_df_with_covar_meta_prefrontal = combined_df_wth_covar
  3432. combined_df_with_covar_meta_prefrontal = combined_df_with_covar_meta_prefrontal[combined_df_with_covar_meta_prefrontal$Tissue_type == "Brain",]
  3433. blood_df_with_covar = combined_df_wth_covar[combined_df_wth_covar$Tissue_type == "Blood",]
  3435. # Tissues to include
  3436. combined_df_with_covar_meta_prefrontal = combined_df_with_covar_meta_prefrontal[combined_df_with_covar_meta_prefrontal$Tissue %in%
  3437. c("DLPFC",
  3438. "dorsal lateral prefrontal cortex (Brodmann Area 9)",
  3439. "Dorsolateral Prefrontal Cortex",
  3440. "Dorsolateral prefrontal cortex (BA9)",
  3441. "Dorsolateral prefrontal cortex (dlPFC; BA8/9)",
  3442. "postmortem dorsolateral prefrontal cortex (DLPFC)",
  3443. "prefrontal cortex"),]
  3444. table(combined_df_with_covar_meta_prefrontal$Tissue, combined_df_with_covar_meta_prefrontal$Study)
  3446. "
  3447. GSE101521 GSE102556 GSE144136 GSE208338 GSE213982 GSE5388 GSE66937 GSE92538_U133_PLUS2 GSE92538_U133A
  3448. DLPFC 0 0 0 0 0 22283 0 0 0
  3449. dorsal lateral prefrontal cortex (Brodmann Area 9) 21197 0 0 0 0 0 0 0 0
  3450. Dorsolateral Prefrontal Cortex 0 0 0 0 0 0 0 54675 22283
  3451. Dorsolateral prefrontal cortex (BA9) 0 0 22389 0 24182 0 0 0 0
  3452. Dorsolateral prefrontal cortex (dlPFC; BA8/9) 0 23331 0 0 0 0 0 0 0
  3453. postmortem dorsolateral prefrontal cortex (DLPFC) 0 0 0 750011 0 0 0 0 0
  3454. prefrontal cortex 0 0 0 0 0 0 70523 0 0
  3456. "
  3458. combined_df_with_covar_meta_prefrontal_list = lapply(unique(combined_df_with_covar_meta_prefrontal$Study), function(x){
  3459. df = combined_df_with_covar_meta_prefrontal[combined_df_with_covar_meta_prefrontal$Study == x,]
  3460. return(df)
  3461. })
  3463. # Run meta-analysis
  3464. analysis_container_with_covar_prefrontal = run_meta_SJ(initial_study_list = combined_df_with_covar_meta_prefrontal_list,
  3465. blood_df_to_compare = blood_df_with_covar)
  3469. meta_with_covar_prefrontal = analysis_container_with_covar_prefrontal[["meta_analysis_df"]]
  3470. combined_df_with_covar_meta_prefrontal_list_reduced = analysis_container_with_covar_prefrontal[["initial_study_list_reduced"]]
  3473. meta_with_covar_prefrontal = add_RRA_to_meta_df(meta_with_covar_prefrontal, combined_df_with_covar_meta_prefrontal_list)
  3474. meta_with_covar_prefrontal_signif = meta_with_covar_prefrontal[meta_with_covar_prefrontal$meta_pval < 0.05,]
  3475. meta_with_covar_prefrontal_signif = meta_with_covar_prefrontal_signif[!is.na(meta_with_covar_prefrontal_signif$meta_pval),]
  3476. nrow(meta_with_covar_prefrontal_signif)
  3477. meta_with_covar_prefrontal[meta_with_covar_prefrontal$gene == "P2RY12",]
  3479. # plots
  3480. dir.create("forest_plots")
  3481. folder = "forest_plots/meta_with_covar_prefrontal"
  3482. dir.create(folder)
  3484. for (i in 1:nrow(meta_with_covar_prefrontal_signif)){
  3487. gene = meta_with_covar_prefrontal_signif$gene[i]
  3488. plot_path = paste0(folder, "/", gene, "_forest.pdf")
  3490. pdf(file = plot_path, width = 10, height = 7)
  3491. print(plot_forest_meta_gene(gene_name = gene, meta_df = combined_df_with_covar_meta_prefrontal_list_reduced))
  3492. dev.off()
  3494. }
  3496. # Heterogeneity estimates
  3497. summary(meta_with_covar_prefrontal$I2)
  3499. " Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
  3500. 0.03 19.87 35.90 38.37 55.02 99.65 35343
  3501. "
  3503. ################### Volcano plot for every meta-analysis ###################
  3505. plot_meta_volcano = function(meta_df, figure_path, title){
  3507. meta_df = meta_df[!is.na(meta_df$meta_pval),]
  3508. meta_df$adj.P.Val = p.adjust(meta_df$meta_pval, method = "fdr")
  3510. Pval_treshold = meta_df[meta_df$adj.P.Val < 0.05,]
  3511. Pval_treshold = max(Pval_treshold$meta_pval)
  3512. Pval_treshold = -log10(Pval_treshold)
  3513. if (is.na(Pval_treshold)){
  3514. Pval_treshold = -log10(0.001)
  3515. }
  3516. PREFIX_logFC_threshold = 0.2
  3518. meta_df$is.highlight = sapply(meta_df$meta_LFc, function(x){
  3519. if (x > PREFIX_logFC_threshold){
  3520. x = "Up"
  3521. } else if (x < - PREFIX_logFC_threshold){
  3522. x = "Down"
  3523. } else {
  3524. x = "None"
  3525. }
  3526. })
  3527. meta_df = mutate(meta_df, is_annotate = ifelse(-log10(meta_pval) >= -log10(0.05) & is.highlight != "None", "yes", "no"))
  3529. # Make the plot
  3530. plot = ggplot(meta_df, aes(x=meta_LFc, y=-log10(meta_pval))) +
  3532. # Show all points
  3533. geom_point(aes(color= factor(is.highlight, levels = c("None", "Down", "Up"))), alpha=0.6, size=4) +
  3534. scale_color_manual(values = c("grey", "skyblue", "red"))
  3536. # Add standard scale for plot
  3537. plot = plot + scale_x_continuous(limits = c(-1.5,1.5))
  3538. plot = plot + scale_y_continuous(limits = c(0,5))
  3540. # Add title
  3541. plot = plot + ggtitle(label = title)
  3543. # Add pval line
  3544. if (!is.null(Pval_treshold)){
  3545. plot = plot + geom_hline(yintercept=Pval_treshold, linetype="dashed",
  3546. color = "red", size=0.5)
  3547. }
  3549. # Add pval line 0.05
  3550. plot = plot + geom_hline(yintercept= -log10(0.05), linetype="dashed",
  3551. color = "blue", size=0.5)
  3552. # Add logFC lines
  3553. if (min(meta_df$meta_LFc) < -PREFIX_logFC_threshold){
  3554. plot = plot + geom_vline(xintercept= -PREFIX_logFC_threshold, linetype="dashed",
  3555. color = "grey", size=0.5)
  3556. }
  3558. if (max(meta_df$meta_LFc) > PREFIX_logFC_threshold){
  3559. plot = plot + geom_vline(xintercept= PREFIX_logFC_threshold, linetype="dashed",
  3560. color = "grey", size=0.5)
  3561. }
  3563. # Add label using ggrepel to avoid overlapping
  3564. if (any(meta_df$is_annotate == "yes")){
  3565. plot = plot +
  3566. geom_label_repel(data=subset(meta_df, is_annotate=="yes"), aes(label=gene), size=4, force = 10,
  3567. max.overlaps = 50)
  3568. }
  3570. plot = plot +
  3571. labs(x = "meta-estimated Log2FC", y = "-log10 p-value") +
  3572. # Custom the theme:
  3573. theme(
  3574. legend.position="none",
  3575. panel.border = element_blank(),
  3576. panel.grid.major.x = element_blank(),
  3577. panel.grid.major.y = element_line(size = 0.1, linetype = 2, color = "black"), # Modifying horizontal lines in the plot
  3578. panel.background = element_blank(),
  3579. plot.title = element_text(size = 20, face = "bold", hjust = 0.5),
  3580. axis.text.x = element_text(size = 12),
  3581. axis.text.y = element_text(size = 12),
  3582. axis.title.x = element_text(size = 14, face = "bold"),
  3583. axis.title.y = element_text(size = 14, face = "bold")
  3584. )
  3585. ggsave(file = figure_path, plot = plot, width=2560, height=1440, units = "px", scale = 2)
  3586. }
  3588. # making volcanos in a loop
  3589. joined_list_all_metas = list(
  3590. meta_no_covar_all_brain,
  3591. meta_no_covar_cortical,
  3592. meta_no_covar_prefrontal,
  3593. meta_with_covar_all_brain,
  3594. meta_with_covar_cortical,
  3595. meta_with_covar_prefrontal
  3596. )
  3597. names_for_list_all = c("All brain", "Cortical", "DLPFC", "All brain (covar)", "Cortical (covar)", "DLPFC (covar)")
  3598. dir.create("volcano_plots")
  3600. for (i in 1:6){
  3602. curr_df = joined_list_all_metas[[i]]
  3603. plot_path = paste0("volcano_plots/", names_for_list_all[i], ".png")
  3604. plot_meta_volcano(curr_df, figure_path = plot_path, title = names_for_list_all[i])
  3606. }
  3609. ################### Sensitivity analysis ###################
  3611. ################### REML Meta (no covar) all brain data ###################
  3612. combined_df_no_covar_meta_all = combined_df_no_covar
  3613. combined_df_no_covar_meta_all = combined_df_no_covar_meta_all[combined_df_no_covar_meta_all$Tissue_type == "Brain",]
  3614. blood_df_no_covar = combined_df_no_covar[combined_df_no_covar$Tissue_type == "Blood",]
  3616. table(combined_df_no_covar_meta_all$Tissue, combined_df_no_covar_meta_all$Study)
  3617. # GSE102556 many tissues
  3618. # GSE202537 many tissues
  3619. # GSE66937 many tissues
  3621. combined_df_no_covar_meta_all_1 = combined_df_no_covar_meta_all[combined_df_no_covar_meta_all$Study %!in% c("GSE102556", "GSE202537", "GSE66937"), ]
  3622. combined_df_no_covar_meta_all_2 = combined_df_no_covar_meta_all[combined_df_no_covar_meta_all$Study %in% c("GSE102556", "GSE202537", "GSE66937"), ]
  3623. # Caudate https://pubmed.ncbi.nlm.nih.gov/39164232/
  3624. # Nac https://pubmed.ncbi.nlm.nih.gov/38894648/ https://pubmed.ncbi.nlm.nih.gov/38965529/ https://pubmed.ncbi.nlm.nih.gov/39167467/
  3625. # GSE102556 -> Dorsolateral prefrontal cortex (dlPFC; BA8/9)
  3626. # GSE202537 -> Nac 22652 (more probes analyzed)
  3627. # GSE66937 -> prefrontal cortex
  3628. combined_df_no_covar_meta_all_2 = combined_df_no_covar_meta_all_2[combined_df_no_covar_meta_all_2$Tissue %in% c("Dorsolateral prefrontal cortex (dlPFC; BA8/9)",
  3629. "Nac",
  3630. "prefrontal cortex"),]
  3631. table(combined_df_no_covar_meta_all_2$Tissue)
  3632. "
  3633. Dorsolateral prefrontal cortex (dlPFC; BA8/9) Nac prefrontal cortex
  3634. 21365 22652 70523
  3635. "
  3636. combined_df_no_covar_meta_full = rbind(combined_df_no_covar_meta_all_1, combined_df_no_covar_meta_all_2)
  3637. table(combined_df_no_covar_meta_full$Tissue, combined_df_no_covar_meta_full$Study)
  3638. combined_df_no_covar_meta_full_list = lapply(unique(combined_df_no_covar_meta_full$Study), function(x){
  3639. df = combined_df_no_covar_meta_full[combined_df_no_covar_meta_full$Study == x,]
  3640. return(df)
  3641. })
  3643. # Run meta-analysis
  3644. REML_analysis_container_no_covar_all_brain = run_meta_REML(initial_study_list = combined_df_no_covar_meta_full_list,
  3645. blood_df_to_compare = blood_df_no_covar)
  3646. REML_meta_no_covar_all_brain = REML_analysis_container_no_covar_all_brain[["meta_analysis_df"]]
  3648. REML_meta_no_covar_all_brain = add_RRA_to_meta_df(REML_meta_no_covar_all_brain, combined_df_no_covar_meta_full_list)
  3649. REML_meta_no_covar_all_brain_signif = REML_meta_no_covar_all_brain[REML_meta_no_covar_all_brain$meta_pval < 0.05,]
  3650. REML_meta_no_covar_all_brain_signif = REML_meta_no_covar_all_brain_signif[!is.na(REML_meta_no_covar_all_brain_signif$meta_pval),]
  3651. nrow(REML_meta_no_covar_all_brain_signif)
  3653. # Heterogeneity
  3654. summary(REML_meta_no_covar_all_brain$I2)
  3655. "
  3656. Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
  3657. 0.00 0.00 12.71 20.95 39.31 94.73 28644
  3658. "
  3660. ################### REML Meta (with covar) all brain data ###################
  3661. combined_df_with_covar_meta_all = combined_df_wth_covar
  3662. combined_df_with_covar_meta_all = combined_df_with_covar_meta_all[combined_df_with_covar_meta_all$Tissue_type == "Brain",]
  3663. blood_df_with_covar = combined_df_wth_covar[combined_df_wth_covar$Tissue_type == "Blood",]
  3665. table(combined_df_with_covar_meta_all$Tissue, combined_df_with_covar_meta_all$Study)
  3667. # GSE102556 many tissues
  3668. # GSE202537 many tissues
  3669. # GSE66937 many tissues
  3671. combined_df_with_covar_meta_all_1 = combined_df_with_covar_meta_all[combined_df_with_covar_meta_all$Study %!in% c("GSE102556", "GSE202537", "GSE66937"), ]
  3672. combined_df_with_covar_meta_all_2 = combined_df_with_covar_meta_all[combined_df_with_covar_meta_all$Study %in% c("GSE102556", "GSE202537", "GSE66937"), ]
  3673. combined_df_with_covar_meta_all_2 = combined_df_with_covar_meta_all_2[combined_df_with_covar_meta_all_2$Tissue %in% c("Dorsolateral prefrontal cortex (dlPFC; BA8/9)",
  3674. "Nac",
  3675. "prefrontal cortex"),]
  3676. table(combined_df_with_covar_meta_all_2$Tissue)
  3677. "
  3678. Dorsolateral prefrontal cortex (dlPFC; BA8/9) Nac prefrontal cortex
  3679. 23331 23919 70523
  3680. "
  3682. combined_df_with_covar_meta_full = rbind(combined_df_with_covar_meta_all_1, combined_df_with_covar_meta_all_2)
  3683. table(combined_df_with_covar_meta_full$Tissue, combined_df_with_covar_meta_full$Study)
  3684. combined_df_with_covar_meta_full_list = lapply(unique(combined_df_with_covar_meta_full$Study), function(x){
  3685. df = combined_df_with_covar_meta_full[combined_df_with_covar_meta_full$Study == x,]
  3686. return(df)
  3687. })
  3689. # Run meta-analysis
  3690. REML_analysis_container_with_covar_all_brain = run_meta_REML(initial_study_list = combined_df_with_covar_meta_full_list,
  3691. blood_df_to_compare = blood_df_with_covar)
  3692. REML_meta_with_covar_all_brain = REML_analysis_container_with_covar_all_brain[["meta_analysis_df"]]
  3694. REML_meta_with_covar_all_brain = add_RRA_to_meta_df(REML_meta_with_covar_all_brain, combined_df_with_covar_meta_full_list)
  3695. REML_meta_with_covar_all_brain_signif = REML_meta_with_covar_all_brain[REML_meta_with_covar_all_brain$meta_pval < 0.05,]
  3696. REML_meta_with_covar_all_brain_signif = REML_meta_with_covar_all_brain_signif[!is.na(REML_meta_with_covar_all_brain_signif$meta_pval),]
  3697. nrow(REML_meta_with_covar_all_brain_signif)
  3699. ################### REML Meta cortical regions (no covar) ###################
  3700. combined_df_no_covar_meta_cortical = combined_df_no_covar
  3701. combined_df_no_covar_meta_cortical = combined_df_no_covar_meta_cortical[combined_df_no_covar_meta_cortical$Tissue_type == "Brain",]
  3702. blood_df_no_covar = combined_df_no_covar[combined_df_no_covar$Tissue_type == "Blood",]
  3704. table(combined_df_no_covar_meta_cortical$Tissue, combined_df_no_covar_meta_cortical$Study)
  3705. # Tissues to include
  3706. # Cingulate gyrus 25 (Cg25) is also cortical but Dorsolateral prefrontal cortex (dlPFC; BA8/9) was used instead
  3707. combined_df_no_covar_meta_cortical = combined_df_no_covar_meta_cortical[combined_df_no_covar_meta_cortical$Tissue %in%
  3708. c("DLPFC",
  3709. "dorsal lateral prefrontal cortex (Brodmann Area 9)",
  3710. "Dorsolateral Prefrontal Cortex",
  3711. "Dorsolateral prefrontal cortex (BA9)",
  3712. "Dorsolateral prefrontal cortex (dlPFC; BA8/9)",
  3713. "Orbitofrontal cortex",
  3714. "postmortem dorsolateral prefrontal cortex (DLPFC)",
  3715. "prefrontal cortex",
  3716. "temporal cortex (BA20 and BA36)"),]
  3717. table(combined_df_no_covar_meta_cortical$Tissue, combined_df_no_covar_meta_cortical$Study)
  3719. "
  3720. GSE101521 GSE102556 GSE144136 GSE208338 GSE213982 GSE243356 GSE5388 GSE5389 GSE66937 GSE92538_U133_PLUS2 GSE92538_U133A
  3721. DLPFC 0 0 0 0 0 0 22283 0 0 0 0
  3722. dorsal lateral prefrontal cortex (Brodmann Area 9) 18677 0 0 0 0 0 0 0 0 0 0
  3723. Dorsolateral Prefrontal Cortex 0 0 0 0 0 0 0 0 0 54675 22283
  3724. Dorsolateral prefrontal cortex (BA9) 0 0 19271 0 19191 0 0 0 0 0 0
  3725. Dorsolateral prefrontal cortex (dlPFC; BA8/9) 0 21365 0 0 0 0 0 0 0 0 0
  3726. Orbitofrontal cortex 0 0 0 0 0 0 0 22283 0 0 0
  3727. postmortem dorsolateral prefrontal cortex (DLPFC) 0 0 0 750011 0 0 0 0 0 0 0
  3728. prefrontal cortex 0 0 0 0 0 0 0 0 70523 0 0
  3729. temporal cortex (BA20 and BA36) 0 0 0 0 0 20955 0 0 0 0 0
  3731. "
  3732. combined_df_no_covar_meta_cortical_list = lapply(unique(combined_df_no_covar_meta_cortical$Study), function(x){
  3733. df = combined_df_no_covar_meta_cortical[combined_df_no_covar_meta_cortical$Study == x,]
  3734. return(df)
  3735. })
  3737. # Run meta-analysis
  3738. REML_analysis_container_no_covar_cortical = run_meta_REML(initial_study_list = combined_df_no_covar_meta_cortical_list,
  3739. blood_df_to_compare = blood_df_no_covar)
  3740. REML_meta_no_covar_cortical = REML_analysis_container_no_covar_cortical[["meta_analysis_df"]]
  3742. REML_meta_no_covar_cortical = add_RRA_to_meta_df(REML_meta_no_covar_cortical, combined_df_no_covar_meta_cortical_list)
  3743. REML_meta_no_covar_cortical_signif = REML_meta_no_covar_cortical[REML_meta_no_covar_cortical$meta_pval < 0.05,]
  3744. REML_meta_no_covar_cortical_signif = REML_meta_no_covar_cortical_signif[!is.na(REML_meta_no_covar_cortical_signif$meta_pval),]
  3745. nrow(REML_meta_no_covar_cortical_signif)
  3747. # Heterogeneity
  3748. summary(REML_meta_no_covar_cortical$I2)
  3749. "
  3750. Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
  3751. 0.00 0.00 12.86 22.62 42.70 95.93 31325
  3752. "
  3754. ################### REML Meta cortical regions (with covar) ###################
  3755. combined_df_with_covar_meta_cortical = combined_df_wth_covar
  3756. combined_df_with_covar_meta_cortical = combined_df_with_covar_meta_cortical[combined_df_with_covar_meta_cortical$Tissue_type == "Brain",]
  3757. blood_df_with_covar = combined_df_wth_covar[combined_df_wth_covar$Tissue_type == "Blood",]
  3759. # Tissues to include
  3760. combined_df_with_covar_meta_cortical = combined_df_with_covar_meta_cortical[combined_df_with_covar_meta_cortical$Tissue %in%
  3761. c("DLPFC",
  3762. "dorsal lateral prefrontal cortex (Brodmann Area 9)",
  3763. "Dorsolateral Prefrontal Cortex",
  3764. "Dorsolateral prefrontal cortex (BA9)",
  3765. "Dorsolateral prefrontal cortex (dlPFC; BA8/9)",
  3766. "Orbitofrontal cortex",
  3767. "postmortem dorsolateral prefrontal cortex (DLPFC)",
  3768. "prefrontal cortex",
  3769. "temporal cortex (BA20 and BA36)"),]
  3770. table(combined_df_with_covar_meta_cortical$Tissue, combined_df_with_covar_meta_cortical$Study)
  3772. "
  3773. GSE101521 GSE102556 GSE144136 GSE208338 GSE213982 GSE243356 GSE5388 GSE5389 GSE66937 GSE92538_U133_PLUS2 GSE92538_U133A
  3774. DLPFC 0 0 0 0 0 0 22283 0 0 0 0
  3775. dorsal lateral prefrontal cortex (Brodmann Area 9) 21197 0 0 0 0 0 0 0 0 0 0
  3776. Dorsolateral Prefrontal Cortex 0 0 0 0 0 0 0 0 0 54675 22283
  3777. Dorsolateral prefrontal cortex (BA9) 0 0 22389 0 24182 0 0 0 0 0 0
  3778. Dorsolateral prefrontal cortex (dlPFC; BA8/9) 0 23331 0 0 0 0 0 0 0 0 0
  3779. Orbitofrontal cortex 0 0 0 0 0 0 0 22283 0 0 0
  3780. postmortem dorsolateral prefrontal cortex (DLPFC) 0 0 0 750011 0 0 0 0 0 0 0
  3781. prefrontal cortex 0 0 0 0 0 0 0 0 70523 0 0
  3782. temporal cortex (BA20 and BA36) 0 0 0 0 0 20955 0 0 0 0 0
  3784. "
  3786. combined_df_with_covar_meta_cortical_list = lapply(unique(combined_df_with_covar_meta_cortical$Study), function(x){
  3787. df = combined_df_with_covar_meta_cortical[combined_df_with_covar_meta_cortical$Study == x,]
  3788. return(df)
  3789. })
  3791. # Run meta-analysis
  3792. REML_analysis_container_with_covar_cortical = run_meta_REML(initial_study_list = combined_df_with_covar_meta_cortical_list,
  3793. blood_df_to_compare = blood_df_with_covar)
  3794. REML_meta_with_covar_cortical = REML_analysis_container_with_covar_cortical[["meta_analysis_df"]]
  3796. REML_meta_with_covar_cortical = add_RRA_to_meta_df(REML_meta_with_covar_cortical, combined_df_with_covar_meta_cortical_list)
  3797. REML_meta_with_covar_cortical_signif = REML_meta_with_covar_cortical[REML_meta_with_covar_cortical$meta_pval < 0.05,]
  3798. REML_meta_with_covar_cortical_signif = REML_meta_with_covar_cortical_signif[!is.na(REML_meta_with_covar_cortical_signif$meta_pval),]
  3799. nrow(REML_meta_with_covar_cortical_signif)
  3802. ################### REML Meta prefrontal regions (no covar) ###################
  3803. combined_df_no_covar_meta_prefrontal = combined_df_no_covar
  3804. combined_df_no_covar_meta_prefrontal = combined_df_no_covar_meta_prefrontal[combined_df_no_covar_meta_prefrontal$Tissue_type == "Brain",]
  3805. blood_df_no_covar = combined_df_no_covar[combined_df_no_covar$Tissue_type == "Blood",]
  3807. table(combined_df_no_covar_meta_prefrontal$Tissue, combined_df_no_covar_meta_prefrontal$Study)
  3808. # Tissues to include
  3809. combined_df_no_covar_meta_prefrontal = combined_df_no_covar_meta_prefrontal[combined_df_no_covar_meta_prefrontal$Tissue %in%
  3810. c("DLPFC",
  3811. "dorsal lateral prefrontal cortex (Brodmann Area 9)",
  3812. "Dorsolateral Prefrontal Cortex",
  3813. "Dorsolateral prefrontal cortex (BA9)",
  3814. "Dorsolateral prefrontal cortex (dlPFC; BA8/9)",
  3815. "postmortem dorsolateral prefrontal cortex (DLPFC)",
  3816. "prefrontal cortex"),]
  3817. table(combined_df_no_covar_meta_prefrontal$Tissue, combined_df_no_covar_meta_prefrontal$Study)
  3819. "
  3820. GSE101521 GSE102556 GSE144136 GSE208338 GSE213982 GSE5388 GSE66937 GSE92538_U133_PLUS2 GSE92538_U133A
  3821. DLPFC 0 0 0 0 0 22283 0 0 0
  3822. dorsal lateral prefrontal cortex (Brodmann Area 9) 18677 0 0 0 0 0 0 0 0
  3823. Dorsolateral Prefrontal Cortex 0 0 0 0 0 0 0 54675 22283
  3824. Dorsolateral prefrontal cortex (BA9) 0 0 19271 0 19191 0 0 0 0
  3825. Dorsolateral prefrontal cortex (dlPFC; BA8/9) 0 21365 0 0 0 0 0 0 0
  3826. postmortem dorsolateral prefrontal cortex (DLPFC) 0 0 0 750011 0 0 0 0 0
  3827. prefrontal cortex 0 0 0 0 0 0 70523 0 0
  3829. "
  3830. combined_df_no_covar_meta_prefrontal_list = lapply(unique(combined_df_no_covar_meta_prefrontal$Study), function(x){
  3831. df = combined_df_no_covar_meta_prefrontal[combined_df_no_covar_meta_prefrontal$Study == x,]
  3832. return(df)
  3833. })
  3835. # Run meta-analysis
  3836. REML_analysis_container_no_covar_prefrontal = run_meta_REML(initial_study_list = combined_df_no_covar_meta_prefrontal_list,
  3837. blood_df_to_compare = blood_df_no_covar)
  3838. REML_meta_no_covar_prefrontal = REML_analysis_container_no_covar_prefrontal[["meta_analysis_df"]]
  3840. REML_meta_no_covar_prefrontal = add_RRA_to_meta_df(REML_meta_no_covar_prefrontal, combined_df_no_covar_meta_prefrontal_list)
  3841. REML_meta_no_covar_prefrontal_signif = REML_meta_no_covar_prefrontal[REML_meta_no_covar_prefrontal$meta_pval < 0.05,]
  3842. REML_meta_no_covar_prefrontal_signif = REML_meta_no_covar_prefrontal_signif[!is.na(REML_meta_no_covar_prefrontal_signif$meta_pval),]
  3843. nrow(REML_meta_no_covar_prefrontal_signif)
  3845. # Heterogeneity
  3846. summary(REML_meta_no_covar_prefrontal$I2)
  3847. "
  3848. Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
  3849. 0.00 0.00 14.48 24.13 45.66 97.75 33247
  3850. "
  3853. ################### REML Meta prefrontal regions (with covar) ###################
  3854. combined_df_with_covar_meta_prefrontal = combined_df_wth_covar
  3855. combined_df_with_covar_meta_prefrontal = combined_df_with_covar_meta_prefrontal[combined_df_with_covar_meta_prefrontal$Tissue_type == "Brain",]
  3856. blood_df_with_covar = combined_df_wth_covar[combined_df_wth_covar$Tissue_type == "Blood",]
  3858. # Tissues to include
  3859. combined_df_with_covar_meta_prefrontal = combined_df_with_covar_meta_prefrontal[combined_df_with_covar_meta_prefrontal$Tissue %in%
  3860. c("DLPFC",
  3861. "dorsal lateral prefrontal cortex (Brodmann Area 9)",
  3862. "Dorsolateral Prefrontal Cortex",
  3863. "Dorsolateral prefrontal cortex (BA9)",
  3864. "Dorsolateral prefrontal cortex (dlPFC; BA8/9)",
  3865. "postmortem dorsolateral prefrontal cortex (DLPFC)",
  3866. "prefrontal cortex"),]
  3867. table(combined_df_with_covar_meta_prefrontal$Tissue, combined_df_with_covar_meta_prefrontal$Study)
  3869. "
  3870. GSE101521 GSE102556 GSE144136 GSE208338 GSE213982 GSE5388 GSE66937 GSE92538_U133_PLUS2 GSE92538_U133A
  3871. DLPFC 0 0 0 0 0 22283 0 0 0
  3872. dorsal lateral prefrontal cortex (Brodmann Area 9) 21197 0 0 0 0 0 0 0 0
  3873. Dorsolateral Prefrontal Cortex 0 0 0 0 0 0 0 54675 22283
  3874. Dorsolateral prefrontal cortex (BA9) 0 0 22389 0 24182 0 0 0 0
  3875. Dorsolateral prefrontal cortex (dlPFC; BA8/9) 0 23331 0 0 0 0 0 0 0
  3876. postmortem dorsolateral prefrontal cortex (DLPFC) 0 0 0 750011 0 0 0 0 0
  3877. prefrontal cortex 0 0 0 0 0 0 70523 0 0
  3879. "
  3881. combined_df_with_covar_meta_prefrontal_list = lapply(unique(combined_df_with_covar_meta_prefrontal$Study), function(x){
  3882. df = combined_df_with_covar_meta_prefrontal[combined_df_with_covar_meta_prefrontal$Study == x,]
  3883. return(df)
  3884. })
  3886. # Run meta-analysis
  3887. REML_analysis_container_with_covar_prefrontal = run_meta_REML(initial_study_list = combined_df_with_covar_meta_prefrontal_list,
  3888. blood_df_to_compare = blood_df_with_covar)
  3889. REML_meta_with_covar_prefrontal = REML_analysis_container_with_covar_prefrontal[["meta_analysis_df"]]
  3891. REML_meta_with_covar_prefrontal = add_RRA_to_meta_df(REML_meta_with_covar_prefrontal, combined_df_with_covar_meta_prefrontal_list)
  3892. REML_meta_with_covar_prefrontal_signif = REML_meta_with_covar_prefrontal[REML_meta_with_covar_prefrontal$meta_pval < 0.05,]
  3893. REML_meta_with_covar_prefrontal_signif = REML_meta_with_covar_prefrontal_signif[!is.na(REML_meta_with_covar_prefrontal_signif$meta_pval),]
  3894. nrow(REML_meta_with_covar_prefrontal_signif)
  3896. ################### Saving REML sensitivity ###################
  3897. significant_genes_var = ls(pattern = "REML_meta_")
  3898. significant_genes_var = significant_genes_var[stri_detect_regex(significant_genes_var, pattern = "signif")]
  3899. significant_genes_var = significant_genes_var[!stri_detect_regex(significant_genes_var, pattern = "meta_list")]
  3900. significant_genes_counts = sapply(significant_genes_var, function(x) nrow(get(x)))
  3902. stats_meta = paste0("Signifi gene counts for: ", significant_genes_var)
  3903. stats_meta = paste0(stats_meta, " = ", significant_genes_counts)
  3904. writeLines(stats_meta, "stats_meta_REML.txt")
  3906. # saving files
  3907. REML_meta_list_no_covar_signif = list(
  3908. "REML_meta_no_covar_all_brain_signif" = REML_meta_no_covar_all_brain_signif,
  3909. "REML_meta_no_covar_cortical_signif" = REML_meta_no_covar_cortical_signif,
  3910. "REML_meta_no_covar_prefrontal_signif" = REML_meta_no_covar_prefrontal_signif
  3911. )
  3913. REML_meta_list_with_covar_signif = list(
  3914. "REML_meta_with_covar_all_brain_signif" = REML_meta_with_covar_all_brain_signif,
  3915. "REML_meta_with_covar_cortical_signif" = REML_meta_with_covar_cortical_signif,
  3916. "REML_meta_with_covar_prefrontal_signif" = REML_meta_with_covar_prefrontal_signif
  3917. )
  3919. REML_joined_list_significant = c(REML_meta_list_no_covar_signif, REML_meta_list_with_covar_signif)
  3920. REML_joined_list_significant = lapply(REML_joined_list_significant, function(x){
  3921. x = dplyr::arrange(x, -meta_LFc)
  3922. return(x)
  3923. })
  3925. names_for_list = c(
  3926. "REML All brain",
  3927. "REML Cortical regions",
  3928. "REML DLPFC",
  3929. "REML All brain (covar)",
  3930. "REML Cortical regions (covar)",
  3931. "REML DLPFC (covar)"
  3932. )
  3933. names(REML_joined_list_significant)
  3935. wb = createWorkbook()
  3937. for (i in 1:6) {
  3938. # Add a new worksheet with the sheet name
  3939. addWorksheet(wb, names_for_list[i])
  3941. # Write the data frame to the worksheet
  3942. writeData(wb, sheet = names_for_list[i], REML_joined_list_significant[[i]])
  3944. # Adjust column widths to fit the text
  3945. setColWidths(wb, sheet = names_for_list[i], cols = 1:ncol(REML_joined_list_significant[[i]]), widths = "auto")
  3946. }
  3948. saveWorkbook(wb, "Meta_suicide_significant_genes_REML.xlsx", overwrite = TRUE)
  3950. ################### SV meta (all brain) ###################
  3951. combined_df_SV_meta_all = combined_df_SV
  3952. combined_df_SV_meta_all = combined_df_SV_meta_all[combined_df_SV_meta_all$Tissue_type == "Brain",]
  3953. blood_df_SV = combined_df_SV[combined_df_SV$Tissue_type == "Blood",]
  3955. table(combined_df_SV_meta_all$Tissue, combined_df_SV_meta_all$Study)
  3956. # GSE102556 many tissues
  3957. # GSE202537 many tissues
  3958. # GSE66937 many tissues
  3960. combined_df_SV_meta_all_1 = combined_df_SV_meta_all[combined_df_SV_meta_all$Study %!in% c("GSE102556", "GSE202537", "GSE66937"), ]
  3961. combined_df_SV_meta_all_2 = combined_df_SV_meta_all[combined_df_SV_meta_all$Study %in% c("GSE102556", "GSE202537", "GSE66937"), ]
  3962. # Caudate https://pubmed.ncbi.nlm.nih.gov/39164232/
  3963. # Nac https://pubmed.ncbi.nlm.nih.gov/38894648/ https://pubmed.ncbi.nlm.nih.gov/38965529/ https://pubmed.ncbi.nlm.nih.gov/39167467/
  3964. # GSE102556 -> Dorsolateral prefrontal cortex (dlPFC; BA8/9)
  3965. # GSE202537 -> Nac 22652 (more probes analyzed)
  3966. # GSE66937 -> prefrontal cortex
  3967. combined_df_SV_meta_all_2 = combined_df_SV_meta_all_2[combined_df_SV_meta_all_2$Tissue %in% c("Dorsolateral prefrontal cortex (dlPFC; BA8/9)",
  3968. "Nac",
  3969. "prefrontal cortex"),]
  3970. table(combined_df_SV_meta_all_2$Tissue)
  3971. "
  3972. Dorsolateral prefrontal cortex (dlPFC; BA8/9) Nac prefrontal cortex
  3973. 21365 22652 70523
  3974. "
  3975. combined_df_SV_meta_full = rbind(combined_df_SV_meta_all_1, combined_df_SV_meta_all_2)
  3976. table(combined_df_SV_meta_full$Tissue, combined_df_SV_meta_full$Study)
  3977. combined_df_SV_meta_full_list = lapply(unique(combined_df_SV_meta_full$Study), function(x){
  3978. df = combined_df_SV_meta_full[combined_df_SV_meta_full$Study == x,]
  3979. return(df)
  3980. })
  3982. # Run meta-analysis
  3983. analysis_container_SV_all_brain = run_meta_SJ(initial_study_list = combined_df_SV_meta_full_list,
  3984. blood_df_to_compare = blood_df_SV)
  3985. meta_SV_all_brain = analysis_container_SV_all_brain[["meta_analysis_df"]]
  3986. combined_df_SV_meta_full_list_reduced = analysis_container_SV_all_brain[["initial_study_list_reduced"]]
  3988. meta_SV_all_brain = add_RRA_to_meta_df(meta_SV_all_brain, combined_df_SV_meta_full_list)
  3989. meta_SV_all_brain_signif = meta_SV_all_brain[meta_SV_all_brain$meta_pval < 0.05,]
  3990. meta_SV_all_brain_signif = meta_SV_all_brain_signif[!is.na(meta_SV_all_brain_signif$meta_pval),]
  3991. meta_SV_all_brain_signif = arrange(meta_SV_all_brain_signif, -meta_LFc)
  3992. nrow(meta_SV_all_brain_signif)
  3994. # Heterogeneity
  3995. summary(meta_SV_all_brain$I2)
  3996. "
  3997. Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
  3998. 0.031 25.926 40.445 41.201 55.704 98.730 28644
  3999. "
  4001. # plots
  4002. dir.create("forest_plots")
  4003. folder = "forest_plots/meta_SV_all_brain"
  4004. dir.create(folder)
  4006. for (i in 1:nrow(meta_SV_all_brain_signif)){
  4009. gene = meta_SV_all_brain$gene[i]
  4010. plot_path = paste0(folder, "/", gene, "_forest.pdf")
  4012. pdf(file = plot_path, width = 10, height = 7)
  4013. print(plot_forest_meta_gene(gene_name = gene, meta_df = combined_df_SV_meta_full_list_reduced))
  4014. dev.off()
  4016. }
  4018. ################### SV meta (cortical) ###################
  4019. combined_df_SV_meta_cortical = combined_df_SV
  4020. combined_df_SV_meta_cortical = combined_df_SV_meta_cortical[combined_df_SV_meta_cortical$Tissue_type == "Brain",]
  4021. blood_df_SV = combined_df_SV[combined_df_SV$Tissue_type == "Blood",]
  4023. table(combined_df_SV_meta_cortical$Tissue, combined_df_SV_meta_cortical$Study)
  4025. combined_df_SV_meta_cortical = combined_df_SV_meta_cortical[combined_df_SV_meta_cortical$Tissue %in%
  4026. c("DLPFC",
  4027. "dorsal lateral prefrontal cortex (Brodmann Area 9)",
  4028. "Dorsolateral Prefrontal Cortex",
  4029. "Dorsolateral prefrontal cortex (BA9)",
  4030. "Dorsolateral prefrontal cortex (dlPFC; BA8/9)",
  4031. "Orbitofrontal cortex",
  4032. "postmortem dorsolateral prefrontal cortex (DLPFC)",
  4033. "prefrontal cortex",
  4034. "temporal cortex (BA20 and BA36)"),]
  4035. table(combined_df_SV_meta_cortical$Tissue, combined_df_SV_meta_cortical$Study)
  4037. "
  4038. GSE101521 GSE102556 GSE144136 GSE208338 GSE213982 GSE243356 GSE5388 GSE5389 GSE66937 GSE92538_U133_PLUS2 GSE92538_U133A
  4039. DLPFC 0 0 0 0 0 0 22283 0 0 0 0
  4040. dorsal lateral prefrontal cortex (Brodmann Area 9) 18677 0 0 0 0 0 0 0 0 0 0
  4041. Dorsolateral Prefrontal Cortex 0 0 0 0 0 0 0 0 0 54675 22283
  4042. Dorsolateral prefrontal cortex (BA9) 0 0 19271 0 19191 0 0 0 0 0 0
  4043. Dorsolateral prefrontal cortex (dlPFC; BA8/9) 0 21365 0 0 0 0 0 0 0 0 0
  4044. Orbitofrontal cortex 0 0 0 0 0 0 0 22283 0 0 0
  4045. postmortem dorsolateral prefrontal cortex (DLPFC) 0 0 0 750011 0 0 0 0 0 0 0
  4046. prefrontal cortex 0 0 0 0 0 0 0 0 70523 0 0
  4047. temporal cortex (BA20 and BA36) 0 0 0 0 0 20955 0 0 0 0 0
  4048. "
  4050. combined_df_SV_meta_cortical_list = lapply(unique(combined_df_SV_meta_cortical$Study), function(x){
  4051. df = combined_df_SV_meta_cortical[combined_df_SV_meta_cortical$Study == x,]
  4052. return(df)
  4053. })
  4055. # Run meta-analysis
  4056. analysis_container_SV_cortical = run_meta_SJ(initial_study_list = combined_df_SV_meta_cortical_list,
  4057. blood_df_to_compare = blood_df_SV)
  4058. meta_SV_cortical = analysis_container_SV_cortical[["meta_analysis_df"]]
  4059. combined_df_SV_meta_cortical_list_reduced = analysis_container_SV_cortical[["initial_study_list_reduced"]]
  4061. meta_SV_cortical = add_RRA_to_meta_df(meta_SV_cortical, combined_df_SV_meta_cortical_list)
  4062. meta_SV_cortical_signif = meta_SV_cortical[meta_SV_cortical$meta_pval < 0.05,]
  4063. meta_SV_cortical_signif = meta_SV_cortical_signif[!is.na(meta_SV_cortical_signif$meta_pval),]
  4064. meta_SV_cortical_signif = arrange(meta_SV_cortical_signif, -meta_LFc)
  4065. nrow(meta_SV_cortical_signif)
  4067. summary(meta_SV_cortical$I2)
  4068. "
  4069. Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
  4070. 0.033 22.011 37.720 39.168 54.900 97.410 31325
  4071. "
  4073. dir.create("forest_plots")
  4074. folder = "forest_plots/meta_SV_cortical"
  4075. dir.create(folder)
  4077. for (i in 1:nrow(meta_SV_cortical_signif)){
  4080. gene = meta_SV_cortical$gene[i]
  4081. plot_path = paste0(folder, "/", gene, "_forest.pdf")
  4083. pdf(file = plot_path, width = 10, height = 7)
  4084. print(plot_forest_meta_gene(gene_name = gene, meta_df = combined_df_SV_meta_cortical_list_reduced))
  4085. dev.off()
  4087. }
  4090. ################### SV meta (prefrontal) ###################
  4091. combined_df_SV_meta_prefrontal = combined_df_SV
  4092. combined_df_SV_meta_prefrontal = combined_df_SV_meta_prefrontal[combined_df_SV_meta_prefrontal$Tissue_type == "Brain",]
  4093. blood_df_SV = combined_df_SV[combined_df_SV$Tissue_type == "Blood",]
  4095. table(combined_df_SV_meta_prefrontal$Tissue, combined_df_SV_meta_prefrontal$Study)
  4097. combined_df_SV_meta_prefrontal = combined_df_SV_meta_prefrontal[combined_df_SV_meta_prefrontal$Tissue %in%
  4098. c("DLPFC",
  4099. "dorsal lateral prefrontal cortex (Brodmann Area 9)",
  4100. "Dorsolateral Prefrontal Cortex",
  4101. "Dorsolateral prefrontal cortex (BA9)",
  4102. "Dorsolateral prefrontal cortex (dlPFC; BA8/9)",
  4103. "postmortem dorsolateral prefrontal cortex (DLPFC)",
  4104. "prefrontal cortex"),]
  4105. table(combined_df_SV_meta_prefrontal$Tissue, combined_df_SV_meta_prefrontal$Study)
  4106. "
  4107. GSE101521 GSE102556 GSE144136 GSE208338 GSE213982 GSE5388 GSE66937 GSE92538_U133_PLUS2 GSE92538_U133A
  4108. DLPFC 0 0 0 0 0 22283 0 0 0
  4109. dorsal lateral prefrontal cortex (Brodmann Area 9) 18677 0 0 0 0 0 0 0 0
  4110. Dorsolateral Prefrontal Cortex 0 0 0 0 0 0 0 54675 22283
  4111. Dorsolateral prefrontal cortex (BA9) 0 0 19271 0 19191 0 0 0 0
  4112. Dorsolateral prefrontal cortex (dlPFC; BA8/9) 0 21365 0 0 0 0 0 0 0
  4113. postmortem dorsolateral prefrontal cortex (DLPFC) 0 0 0 750011 0 0 0 0 0
  4114. prefrontal cortex 0 0 0 0 0 0 70523 0 0
  4115. "
  4118. combined_df_SV_meta_prefrontal_list = lapply(unique(combined_df_SV_meta_prefrontal$Study), function(x){
  4119. df = combined_df_SV_meta_prefrontal[combined_df_SV_meta_prefrontal$Study == x,]
  4120. return(df)
  4121. })
  4123. # Run meta-analysis
  4124. analysis_container_SV_prefrontal = run_meta_SJ(initial_study_list = combined_df_SV_meta_prefrontal_list,
  4125. blood_df_to_compare = blood_df_SV)
  4126. meta_SV_prefrontal = analysis_container_SV_prefrontal[["meta_analysis_df"]]
  4127. combined_df_SV_meta_prefrontal_list_reduced = analysis_container_SV_prefrontal[["initial_study_list_reduced"]]
  4129. meta_SV_prefrontal = add_RRA_to_meta_df(meta_SV_prefrontal, combined_df_SV_meta_prefrontal_list)
  4130. meta_SV_prefrontal_signif = meta_SV_prefrontal[meta_SV_prefrontal$meta_pval < 0.05,]
  4131. meta_SV_prefrontal_signif = meta_SV_prefrontal_signif[!is.na(meta_SV_prefrontal_signif$meta_pval),]
  4132. meta_SV_prefrontal_signif = arrange(meta_SV_prefrontal_signif, -meta_LFc)
  4133. nrow(meta_SV_prefrontal_signif)
  4135. summary(meta_SV_prefrontal$I2)
  4136. "
  4137. Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
  4138. 0.02 18.40 33.37 35.54 50.55 98.13 33247
  4139. "
  4141. dir.create("forest_plots")
  4142. folder = "forest_plots/meta_SV_prefrontal"
  4143. dir.create(folder)
  4145. for (i in 1:nrow(meta_SV_prefrontal_signif)){
  4148. gene = meta_SV_prefrontal$gene[i]
  4149. plot_path = paste0(folder, "/", gene, "_forest.pdf")
  4151. pdf(file = plot_path, width = 10, height = 7)
  4152. print(plot_forest_meta_gene(gene_name = gene, meta_df = combined_df_SV_meta_prefrontal_list_reduced))
  4153. dev.off()
  4155. }
  4157. ################### Saving SV sensitivity ###################
  4159. # SV numbers
  4160. # GSE208338 -> 1 suggested
  4161. # GSE5388 -> 5 suggested and 5 used
  4162. # GSE5389 -> 0 suggested
  4163. # GSE66937 -> DLPFC 0 suggested
  4164. # GSE199536 -> 5 suggested and 5 used
  4165. # GSE92538_U133A -> 3 suggested
  4166. # GSE92538_U133_PLUS2 -> 2 suggested
  4167. # GSE102556 -> very large number -> 5 used
  4168. # GSE243356 -> very large number (59) -> 5 used
  4169. # GSE248260 -> very large number (22) -> 5 used
  4170. # GSE247998 -> very large number (98) -> 5 used
  4171. # GSE202537 -> very large number (>60) -> 5 used
  4172. # GSE101521 -> very large number (57) -> 5 used
  4173. # GSE144136 -> very large number (30) -> 5 used
  4174. # GSE213982 -> very large number (34) -> 5 used
  4176. significant_genes_var = ls(pattern = "meta_SV")
  4177. significant_genes_var = significant_genes_var[stri_detect_regex(significant_genes_var, pattern = "signif")]
  4178. significant_genes_var = significant_genes_var[!stri_detect_regex(significant_genes_var, pattern = "meta_list")]
  4179. significant_genes_counts = sapply(significant_genes_var, function(x) nrow(get(x)))
  4181. stats_meta = paste0("Signifi gene counts for: ", significant_genes_var)
  4182. stats_meta = paste0(stats_meta, " = ", significant_genes_counts)
  4183. writeLines(stats_meta, "stats_meta_SV.txt")
  4185. # saving files
  4186. meta_list_SV_signif = list(
  4187. "meta_SV_all_brain_signif" = meta_SV_all_brain_signif,
  4188. "meta_SV_cortical_signif" = meta_SV_cortical_signif,
  4189. "meta_SV_prefrontal_signif" = meta_SV_prefrontal_signif
  4190. )
  4192. names_for_list = c(
  4193. "Meta all brain (SVs)",
  4194. "Meta cortical (SVs)",
  4195. "Meta prefrontal cortex (SVs)"
  4196. )
  4198. wb = createWorkbook()
  4200. for (i in 1:3) {
  4201. # Add a new worksheet with the sheet name
  4202. addWorksheet(wb, names_for_list[i])
  4204. # Write the data frame to the worksheet
  4205. writeData(wb, sheet = names_for_list[i], meta_list_SV_signif[[i]])
  4207. # Adjust column widths to fit the text
  4208. setColWidths(wb, sheet = names_for_list[i], cols = 1:ncol(meta_list_SV_signif[[i]]), widths = "auto")
  4209. }
  4211. saveWorkbook(wb, "Meta_suicide_significant_genes_SV.xlsx", overwrite = TRUE)
  4214. ## Save volcano plots
  4215. # making volcanos in a loop
  4216. joined_list_all_metas_SV = list(
  4217. meta_SV_all_brain,
  4218. meta_SV_cortical,
  4219. meta_SV_prefrontal
  4220. )
  4221. names_for_list_all = c("All brain (SVs)", "Cortical (SVs)", "DLPFC (SVs)")
  4223. for (i in 1:3){
  4225. curr_df = joined_list_all_metas_SV[[i]]
  4226. plot_path = paste0("volcano_plots/", names_for_list_all[i], ".png")
  4227. plot_meta_volcano(curr_df, figure_path = plot_path, title = names_for_list_all[i])
  4229. }
  4231. ################### Venn diagrams for SVs ###################
  4233. meta_list_SV_signif = list(
  4234. meta_SV_all_brain_signif,
  4235. meta_SV_cortical_signif,
  4236. meta_SV_prefrontal_signif
  4237. )
  4239. meta_gene_list_overlaps_SV = lapply(meta_list_SV_signif, function(x){
  4240. x = x$gene
  4241. return(x)
  4242. })
  4244. names(meta_gene_list_overlaps_SV) = c(
  4245. "All brain",
  4246. "Cortical regions",
  4247. "DLPFC"
  4248. )
  4250. make_Venn_digram_list(named_list = meta_gene_list_overlaps_SV, palette = 9, plot_full_path = "Venn_SV.pdf")
  4253. dir.create("analyses_comparison_covar_SV_nothing_venn")
  4254. meta_gene_list_SVs = list(
  4255. meta_no_covar_all_brain_signif,
  4256. meta_SV_all_brain_signif,
  4257. meta_with_covar_all_brain_signif
  4258. )
  4259. meta_gene_list_SVs = lapply(meta_gene_list_SVs, function(x){
  4260. x = x$gene
  4261. return(x)
  4262. })
  4264. names(meta_gene_list_SVs) = c(
  4265. "All brain (no covar)",
  4266. "All brain (SVs)",
  4267. "All brain (covar)"
  4268. )
  4270. make_Venn_digram_list(named_list = meta_gene_list_SVs,
  4271. palette = 5,
  4272. label_text_size = 7,
  4273. plot_full_path = "analyses_comparison_covar_SV_nothing_venn/1_Venn_SVs_vs_other_all.png")
  4276. meta_gene_list_SVs = list(
  4277. meta_no_covar_cortical_signif,
  4278. meta_SV_cortical_signif,
  4279. meta_with_covar_cortical_signif
  4280. )
  4281. meta_gene_list_SVs = lapply(meta_gene_list_SVs, function(x){
  4282. x = x$gene
  4283. return(x)
  4284. })
  4286. names(meta_gene_list_SVs) = c(
  4287. "Cortical (no covar)",
  4288. "Cortical (SVs)",
  4289. "Cortical (covar)"
  4290. )
  4292. make_Venn_digram_list(named_list = meta_gene_list_SVs,
  4293. label_text_size = 7,
  4294. palette = 4, plot_full_path = "analyses_comparison_covar_SV_nothing_venn/2_Venn_SVs_vs_other_cortical.png")
  4296. meta_gene_list_SVs = list(
  4297. meta_no_covar_prefrontal_signif,
  4298. meta_SV_prefrontal_signif,
  4299. meta_with_covar_prefrontal_signif
  4300. )
  4301. meta_gene_list_SVs = lapply(meta_gene_list_SVs, function(x){
  4302. x = x$gene
  4303. return(x)
  4304. })
  4306. names(meta_gene_list_SVs) = c(
  4307. "Prefrontal (no covar)",
  4308. "Prefrontal (SVs)",
  4309. "Prefrontal (covar)"
  4310. )
  4312. make_Venn_digram_list(named_list = meta_gene_list_SVs,
  4313. label_text_size = 7,
  4314. palette = 3, plot_full_path = "analyses_comparison_covar_SV_nothing_venn/3_Venn_SVs_vs_other_prefrontal.png")
  4317. # Define cropping function
  4318. crop_image <- function(img, row_range, col_range) {
  4319. # img: 3D array (height, width, channels)
  4320. img[row_range, col_range, , drop = FALSE]
  4321. }
  4323. images_in_folder = list.files("analyses_comparison_covar_SV_nothing_venn", full.names = TRUE)
  4324. images_in_folder = images_in_folder[!grepl("combined_image", images_in_folder)]
  4326. image_list = lapply(images_in_folder, png::readPNG)
  4328. # Define cropping region
  4329. row_range = 50:950
  4330. col_range = 50:950
  4332. cropped_list = lapply(image_list, crop_image, row_range, col_range)
  4333. image_grobs = lapply(cropped_list, rasterGrob)
  4335. # Combine
  4336. height = length(row_range)
  4337. width = length(col_range)
  4338. combined_file_name = "analyses_comparison_covar_SV_nothing_venn/Venn_DE_comparison_combined_image.png"
  4340. png(filename = combined_file_name, width = width*3, height = height, units = "px")
  4341. grid.arrange(grobs = image_grobs[1:3], ncol = 3)
  4342. dev.off()
  4344. ################### Moderator analysis for RE with SJ ###################
  4345. cohort_moderators_data = openxlsx::read.xlsx("Cohort_moderators.xlsx")
  4346. colnames(cohort_moderators_data)[1] = "Study"
  4347. cohort_moderators_data = cohort_moderators_data[cohort_moderators_data$Study != "GSE66937", ]
  4349. cohort_moderators_data$Platform_binary = factor(cohort_moderators_data$Platform_binary, levels = c( "Array", "RNAseq"))
  4350. cohort_moderators_data$Primary_tissue_group = factor(cohort_moderators_data$Primary_tissue_group,
  4351. levels = c( "Cortical", "Non_cortical")) # Better to stratify
  4354. # Moderators to include
  4355. cohort_moderators_data$Platform_binary
  4356. cohort_moderators_data$Mean_PMI_hours
  4357. cohort_moderators_data$Percent_male
  4358. cohort_moderators_data$Percent_Depr
  4359. cohort_mode

OLINK_PSY_SUICIDE_SCRIPT.R at commit 63e5b8e, under GPL-3.0 · at the source

Overview

Authors: Aleksandr V. Sokolov1, Muataz S. Lafta1, Jussi Jokinen2,3, Helgi B. Schiöth1,4
ORCID iDs: Aleksandr V. Sokolov, Helgi B. Schiöth
  1. Department of Surgical Sciences, Functional Pharmacology and Neuroscience, Uppsala University,Uppsala, Sweden
  2. Department of Clinical Neuroscience/Center for Psychiatric Research, Karolinska Institutet,Stockholm, Sweden
  3. Department of Clinical Sciences/Psychiatry, Umeå University,Umeå, Sweden
  4. Laboratory of Pharmaceutical Pharmacology, Latvian Institute of Organic Synthesis,Riga, Latvia
Institutions: Uppsala University (Sweden); Karolinska Institutet (Sweden); Umeå University (Sweden); Latvijas Organiskās Sintēzes Institūts (Latvia)
Journal: Translational psychiatry, volume 16, issue 1, article 222
Dates: received 14 March 2025; accepted 10 March 2026; published online 31 March 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1038/s41398-026-03978-8 · PMID 41916959 · PMCID PMC13039836 · OpenAlex W7144332257
Open access: gold, a free copy (OpenAlex)
Status: code verified (how the code is found and verified)
Categories: genetics / omics (modality), human (organism), cellular / molecular (subfield)
Methods: Statistics
Keywords: Molecular neuroscience, Biomarkers
MeSH: Brain*, Suicide*, Transcriptome*, Cohort Studies, CX3C Chemokine Receptor 1, Gene Expression Profiling, Humans, RNA, Long Noncoding, Single-Cell Gene Expression Analysis, SOX9 Transcription Factor (* major topic)
Topic: Tryptophan and brain disorders (Biological Psychiatry, Neuroscience), according to OpenAlex
Funding: Uppsala Universitet
Citations: cited by 2 papers (Europe PMC); 123 references in the paper

Abstract

Suicide remains a critical global public health issue, accounting for nearly one million deaths annually and imposing profound societal and economic burdens. Despite its urgency, the lack of diagnostic and predictive biomarkers continues to hinder the development of effective prevention and treatment strategies. This study presents a comprehensive meta-analysis that integrates publicly available postmortem brain transcriptomic datasets and a domestic cohort, encompassing 16 cohorts. The transcriptomic data, sourced from the Gene Expression Omnibus repository, were generated using various techniques, including traditional RNA sequencing, microarray methods, and single-cell RNA sequencing. Differential expression analyses were performed across multiple brain regions, with meta-analyses stratified by cortical regions, the dorsolateral prefrontal cortex (DLPFC), and combined. We further analyzed whether covariates may affect the identified genes. Three meta-analytic approaches were employed, complemented by pathway and cell-set enrichment analyses. The unadjusted meta-analysis consistently identified several genes with altered expression, including upregulated P2RY12, CX3CR1, and GPR34, and downregulated SOX9 and PMP2, all at nominal significance. Additionally, multiple genes encoding long non-coding RNAs (lncRNAs) exhibited nominally altered expression in suicide, including RP5-837J1.4, AC159540.14, DNM1P47, AC004158.2, EEF1A1P30, and RP11-339B21.8. Several alternative strategies to run meta-analysis were performed and moderators were investigated. Cell-type-specific expression deconvolution and meta-analysis identified several genes overlapping with bulk expression meta-analysis, and genes were attributed to neuronal lineages. These findings highlight plausible molecular targets for future validation studies, suggesting the involvement of microglia (P2RY12 and CX3CR1), astrocytes (SOX9), immune responses (GPR34), myelin regulation (PMP2), and epigenetic modulation via lncRNAs. This research advances the understanding of the molecular architecture of suicide and provides a foundation for future studies focused on targeted prevention and therapeutic interventions.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 23 matches between paragraphs and lines of code.

AleksandrVSokolov/suicide_meta

License: GPL-3.0
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: 63e5b8ed9b68270cd9b39a74f6031239e177fc49, 10 April 2026
Languages: R (7), Jupyter (5), Python (4), Shell (3)
Size: 21 files, 19 scripts
Software Heritage: not archived
Found in: “Code availability”
Holds: README, license file, 5 notebooks
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: pandas (8 files), clusterProfiler (6 files), data.table (6 files), ggplot2 (6 files), igraph (6 files), limma (6 files), tidyverse (6 files), metafor (4 files), Seurat (4 files), edgeR (3 files)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
21 files
At the source: github.com/AleksandrVSokolov/suicide_meta

Code availability

The source code for all stages of the work is publicly-available at https://github.com/AleksandrVSokolov/suicide_meta.

Reproduced under the paper's license (CC BY), from the paper cited above.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 19 scripts, each with its path and the digest of its content;
  • 23 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability

This analysis is based on open-access transcriptomic datasets available at GEO.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 4 authors, 2 keywords, 10 MeSH terms, 1 funder, 118 references.

Cite

This paper

Sokolov, A. V., Lafta, M. S., Jokinen, J., & Schiöth, H. B. (2026). Identification of suicide brain transcriptomic signatures using meta-analysis of multiple cohorts. Translational psychiatry, 16(1), 222. https://doi.org/10.1038/s41398-026-03978-8

BibTeX

@article{sokolov2026identification,
author = {Sokolov, Aleksandr V. and Lafta, Muataz S. and Jokinen, Jussi and Schiöth, Helgi B.},
title = {{Identification of suicide brain transcriptomic signatures using meta-analysis of multiple cohorts}},
journal = {Translational psychiatry},
year = {2026},
month = mar,
volume = {16},
number = {1},
pages = {222},
publisher = {Nature Publishing Group},
issn = {2158-3188},
doi = {10.1038/s41398-026-03978-8},
url = {https://doi.org/10.1038/s41398-026-03978-8},
pmid = {41916959},
pmcid = {PMC13039836}
}

RIS

TY - JOUR
AU - Sokolov, Aleksandr V.
AU - Lafta, Muataz S.
AU - Jokinen, Jussi
AU - Schiöth, Helgi B.
TI - Identification of suicide brain transcriptomic signatures using meta-analysis of multiple cohorts
T2 - Translational psychiatry
J2 - Transl Psychiatry
PY - 2026
DA - 2026/03/31
VL - 16
IS - 1
SP - 222
SN - 2158-3188
PB - Nature Publishing Group
DO - 10.1038/s41398-026-03978-8
UR - https://doi.org/10.1038/s41398-026-03978-8
LA - en
ER -

CSL-JSON

{
"id": "10.1038/s41398-026-03978-8",
"type": "article-journal",
"title": "Identification of suicide brain transcriptomic signatures using meta-analysis of multiple cohorts",
"container-title": "Translational psychiatry",
"author": [
{
"family": "Sokolov",
"given": "Aleksandr V."
},
{
"family": "Lafta",
"given": "Muataz S."
},
{
"family": "Jokinen",
"given": "Jussi"
},
{
"family": "Schiöth",
"given": "Helgi B."
}
],
"container-title-short": "Transl Psychiatry",
"volume": "16",
"issue": "1",
"page": "222",
"DOI": "10.1038/s41398-026-03978-8",
"PMID": "41916959",
"PMCID": "PMC13039836",
"ISSN": "2158-3188",
"publisher": "Nature Publishing Group",
"URL": "https://doi.org/10.1038/s41398-026-03978-8",
"language": "en",
"issued": {
"date-parts": [
[
2026,
3,
31
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41467-026-71542-5 [code]
Astrocyte fatty acid metabolism as a driver of risk for major depressive disorder.
Journal: Nature communications
In common: edgeR, limma, igraph, 4 other tools, genetics / omics, cellular / molecular, 6 references
[2] doi:10.1186/s12888-026-08170-2
Gene expression meta-analysis in the prefrontal cortex: unraveling biological underpinnings of suicidal risk.
Journal: BMC psychiatry
In common: genetics / omics, 11 references
[3] doi:10.1038/s41467-025-59115-4 [code]
Large-scale transcriptomic analyses of major depressive disorder reveal convergent dysregulation of synaptic pathways in excitatory neurons.
Journal: Nature communications
In common: edgeR, limma, clusterProfiler, 5 other tools, genetics / omics, cellular / molecular, 5 references
[4] doi:10.1007/s00401-026-03035-0
Expression of GPR34 in microglia remains stable in human Alzheimer's disease.
Journal: Acta neuropathologica
In common: cellular / molecular, 9 references
[5] doi:10.1038/s41467-023-38530-5 [code]
Cell type specific transcriptomic differences in depression show similar patterns between males and females but implicate distinct cell types and genes.
Journal: Nature communications
In common: edgeR, limma, clusterProfiler, 5 other tools, genetics / omics, cellular / molecular, 4 references
[6] doi:10.1038/s41593-026-02367-0 [code]
A reproducible three-dimensional model of human brain tissue to investigate physiological and disease-associated microglia phenotypes.
Journal: Nature neuroscience
In common: edgeR, limma, clusterProfiler, 5 other tools, cellular / molecular, 3 references
[7] doi:10.1038/s41586-025-10047-5 [code]
Developmental convergence and divergence in human stem cell models of autism.
Journal: Nature
In common: metafor, edgeR, limma, 6 other tools, genetics / omics, 1 reference
[8] doi:10.3390/ijms27136068 [code]
Loss of Neuropeptide Y Signaling Accompanies the Neural-to-Mesenchymal Transcriptional Transition in Glioblastoma: A Multi-Scale Transcriptomic Analysis.
Journal: International journal of molecular sciences
In common: metafor, edgeR, limma, 5 other tools, genetics / omics, cellular / molecular, 2 references
[9] doi:10.1126/sciadv.adu0726 [code]
CHCHD2 mutant mice link mitochondrial deficits to PD pathophysiology.
Journal: Science advances
In common: metafor, edgeR, limma, 5 other tools, cellular / molecular, 1 reference
[10] doi:10.1038/s41591-024-03150-z [code]
A brain cell atlas integrating single-cell transcriptomes across human brain regions.
Journal: Nature medicine
In common: edgeR, clusterProfiler, igraph, 5 other tools, genetics / omics, 2 references

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Add the badge to the code's README

OSCR: the badge of this paper

If you are an author or a maintainer of this paper's code, you can add this badge. It says what OSCR knows of the paper (its code, its validated tracing map) and links the code to this page. Paste one of these snippets under the README's title. Keep the two marks around it, so that the badge can be replaced in place and is never added twice. OSCR changes nothing in your repository.

Markdown
<!-- oscr-badge:start -->
[![OSCR](https://openscicode.org/badge/doi/10.1038/s41398-026-03978-8.svg)](https://openscicode.org/badge/doi/10.1038/s41398-026-03978-8)
<!-- oscr-badge:end -->
reStructuredText
.. oscr-badge:start

.. image:: https://openscicode.org/badge/doi/10.1038/s41398-026-03978-8.svg
   :target: https://openscicode.org/badge/doi/10.1038/s41398-026-03978-8
   :alt: OSCR

.. oscr-badge:end
HTML
<!-- oscr-badge:start -->
<a href="https://openscicode.org/badge/doi/10.1038/s41398-026-03978-8"><img src="https://openscicode.org/badge/doi/10.1038/s41398-026-03978-8.svg" alt="OSCR"></a>
<!-- oscr-badge:end -->

On GitHub, open the README, choose Edit (the pencil), paste the snippet under the title, then "Commit changes".

If you choose to, you can have the badge kept up to date for you, with oscr badge --write from the command line or with the GitHub Action in your repository's workflows. Every other form and option (a card, AsciiDoc, LaTeX, shields.io) is on the badge page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: code errors and mismatches between the code and the paper about this paper and its code, from readers and its authors, beside the paper and its tracing map.

Report a mismatch between the code and the paper · Report a code error (signed in)

Reproductions: reports from readers who ran the authors' code and did not get the paper's result, what they ran, with which environment, commit and data.

Report a reproduction failure (signed in)

Versions of its code: the releases of its code tied to a version of this paper (the preprint, the accepted manuscript, the version of record…), each with the tracing map versioned with it.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.