OSCR

The Converging Effects of Different Categories of Antidepressants on the Brain: A Systematic Meta-Analysis of Public Transcriptional Profiling Data From the Hippocampus and Cortex.

Code ↔ Paper

22 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 22 matches
  1. [1] § Results › Cortex › Cortical Meta‐Analysis Reveals One Gene Consistently Upregulated by Antidepressants (Atp6v1b2) ↔ revisions/HPC_CTX_RunMetaRegression_fullCode.R, lines 816–886 · score 0.80 · vitamin B6, enzyme encoded, Pyridoxamine, catalyzes, synthesis, terminal
  2. [2] § Methods › Assessment of Meta‐Analysis Result Robustness and Validity ↔ revisions/CTX_RunMetaAnalysis_wHeterogeneity_PubBias_Robustness_Influence.R, lines 94–237 · score 0.77 · leave1out, Publication bias, robustness, regtest, asymmetry, funnel
  3. [3] § Methods › Assessment of Meta‐Analysis Result Robustness and Validity ↔ revisions/HPC_RunMetaAnalysis_wHeterogeneity_PubBias_Robustness_Influence.R, lines 28–171 · score 0.77 · leave1out, Publication bias, robustness, regtest, asymmetry, funnel
  4. [4] § Results › Cortex › Cortical Meta‐Analysis Reveals One Gene Consistently Upregulated by Antidepressants (Atp6v1b2) ↔ revisions/CTX_RunMetaAnalysis_wHeterogeneity_PubBias_Robustness_Influence.R, lines 479–566 · score 0.76 · transporting V1 subunit, ATPase, Atp6v1b2, FDR, gene, meta
  5. [5] § Methods › Result Extraction and Meta‐Analysis ↔ CTX/PFC_Aligning_Datasets.R, lines 48–85 · score 0.67 · Mouse Ortholog Database, Jackson Labs, trimmed, Log2FC, Fold Change, species
  6. [6] § Methods › Result Extraction and Meta‐Analysis ↔ HPC/whole/AligningResults.R, lines 42–79 · score 0.67 · Mouse Ortholog Database, Jackson Labs, trimmed, Log2FC, Fold Change, species
  7. [7] § Results › Hippocampus › Meta‐Analysis Reveals Consistent Differential Expression Across Antidepressant Categories That Was Enriched in Many Cell Types and Pathways Previously Linked to Depression ↔ CTX/First_iteration_Sophia_Meta.R, lines 764–849 · score 0.67 · hierarchically clustered heatmap, mouse gene symbol, Entrez ID, Discovery Rate, expressed gene, fold change
  8. [8] § Methods ↔ HPC/whole/DatasetSearch.R, lines 127–198 · score 0.61 · manipulation unrelated, incorrect developmental stage, metadata, Seq, cell, transcriptional
  9. [9] § Results › Hippocampus › Exploratory Analyses: Traditional and Nontraditional Antidepressants Have Overlapping Effects on the Hippocampus ↔ HPC/Nontraditional/ComparisonTrad_NonTrad.R, lines 126–167 · score 0.61 · Spearman rank correlation, top genes, expressed gene, datapoint, Trad, Log2FC
  10. [10] § Methods › Functional Patterns ↔ HPC/whole/NondirectionalGSEA.R, lines 45–96 · score 0.61 · Brain.GMT, abs, nperm, nondirectional, enriched, min
  11. [11] § Methods › Assessment or Residual Heterogeneity ↔ revisions/CTX_RunMetaAnalysis_wHeterogeneity_PubBias_Robustness_Influence.R, lines 94–237 · score 0.60 · sampling variability, residual heterogeneity, metafor, Tau2, tailed, ratio
  12. [12] § Methods › Assessment or Residual Heterogeneity ↔ revisions/HPC_RunMetaAnalysis_wHeterogeneity_PubBias_Robustness_Influence.R, lines 28–171 · score 0.60 · sampling variability, residual heterogeneity, metafor, Tau2, tailed, ratio
  13. [13] § Methods ↔ CTX/First_iteration_Sophia_Meta.R, lines 270–308 · score 0.60 · manipulation unrelated, incorrect developmental stage, metadata, Seq, cell, transcriptional
  14. [14] § Methods › Initial Filtering ↔ CTX/First_iteration_Sophia_Meta.R, lines 149–188 · score 0.60 · cerebral cortex, organism part, metadata, filtered, cell, treatment
  15. [15] § Results › Hippocampus › Exploratory Analyses: Traditional and Nontraditional Antidepressants Have Overlapping Effects on the Hippocampus ↔ HPC/whole/DatasetSearch.R, lines 1–39 · score 0.60 · norepinephrine reuptake inhibitor, selective serotonin reuptake, Electroconvulsive therapy, SNRI, Tricyclic, SSRI
  16. [16] § Results › Hippocampus › Assessment of Result Robustness and Validity ↔ revisions/HPC_RunMetaAnalysis_wHeterogeneity_PubBias_Robustness_Influence.R, lines 471–531 · score 0.58 · publication bias, Calb1, Egr4, Grm7, Tdo2, Jun
  17. [17] § Methods › Dataset Identification and Search ↔ HPC/Nontraditional/Nontraditional_Meta_Analysis.R, lines 1–52 · score 0.57 · gemma.R, traditional antidepressants, agomelatine, tianeptine, ketamine, TMS
  18. [18] § Methods › Functional Patterns ↔ HPC/traditional only/Updated_Adrienne_Phi_Traditional_fgsea.R, lines 41–80 · score 0.56 · Brain.GMT, nperm, Log2FC, fgsea, enriched, min
  19. [19] § Results › Hippocampus › Meta‐Analysis Reveals Consistent Differential Expression Across Antidepressant Categories That Was Enriched in Many Cell Types and Pathways Previously Linked to Depression ↔ HPC/Nontraditional/Nontraditional_Meta_Analysis.R, lines 622–670 · score 0.55 · hierarchically clustered heatmap, gene symbol, Mm, Rn, fold change, Discovery
  20. [20] § Results › Cortex › Assessment of Result Robustness and Validity ↔ revisions/HPC_RunMetaAnalysis_wHeterogeneity_PubBias_Robustness_Influence.R, lines 471–531 · score 0.51 · weak evidence, publication bias, Egger, Robustness, FDR, gene
  21. [21] § Methods › Reprocessing ↔ CTX/First_iteration_Sophia_Meta.R, lines 351–397 · score 0.51 · transformed twice, GSE84183, GSE118670, variable, Gemma, treatment
  22. [22] § Results › Cortex › Cortical Meta‐Analysis Reveals One Gene Consistently Upregulated by Antidepressants (Atp6v1b2) ↔ revisions/CTX_RunMetaAnalysis_wHeterogeneity_PubBias_Robustness_Influence.R, lines 829–889 · score 0.51 · positive correlation, Atp6v1b2, CTX, Log2FC, Forest, Spearman

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 1,083 lines · 55 KB · no license · 4 matches

  1. #This code document includes the code for a function that is designed to run a basic meta-analysis of Log2FC and sampling variance values using our previously generated objects MetaAnalysis_FoldChanges & MetaAnalysis_SV
  2. #Megan Hagenauer
  3. #Original version: July 25 2024
  4. #In response to reviewers' comments, this function has been updated to include heterogeneity statistics, publication bias statistics, and robustness statistics
  5. #In the process, we caught an error in the code (inclusion of GSE84185, which is labeled ACG in the Gemma metadata but is actually the blood data from the same paper as GSE84183) and fixed it
  6. #Updated version: March 10, 2026
  9. ### This version is for the PFC
  12. ######################
  14. #Grabbing input data and setting the working directory:
  16. setwd("~/Library/CloudStorage/[email hidden]/My Drive/BrainAlchemyProject/ProjectFolders/2024_EvaGeoghegan_Antidepressants_Hippocampus/R_Code_And_Workspaces/Meta-analysis")
  18. #Eva directed me to this PFC workspace, but I think it is not the final one
  19. # load("~/Library/CloudStorage/[email hidden]/My Drive/BrainAlchemyProject/ProjectFolders/2024_EvaGeoghegan_Antidepressants_Hippocampus/R_Code_And_Workspaces/Meta-analysis/04112025_PFC_Meta_Run_Reanalysis_Included.RData")
  21. load("~/Library/CloudStorage/[email hidden]/My Drive/BrainAlchemyProject/ProjectFolders/2024_SophiaEspinoza_Antidepressants_FrontalCortex/ROutput_And_Results/New Results_Reanalysis/April.RData")
  23. load("~/Library/CloudStorage/[email hidden]/My Drive/BrainAlchemyProject/ProjectFolders/2024_EvaGeoghegan_Antidepressants_Hippocampus/ROutput_And_Results/Revisions/PFC/April212025Workspace_16Comparisons_6NAcutoff/Workspace_16comparisons6NAcutoff_CorrectComparisons.RData")
  25. setwd("~/Library/CloudStorage/[email hidden]/My Drive/BrainAlchemyProject/ProjectFolders/2024_EvaGeoghegan_Antidepressants_Hippocampus/ROutput_And_Results/Revisions/PFC/April212025Workspace_16Comparisons_6NAcutoff")
  28. ######################
  30. #This workspace contains an error - GSE84185 is actually the blood data from the same paper as GSE84183
  31. #We eventually decided to cut it out - but I don't see that in the workspace/code
  32. #And the meta-analysis output in our supplemental file also erroneously reports 17 contrasts - so we must have written up the paper using the wrong workspace/code/output
  33. #I don't know where the final workspace/code/output is with GSE84185 removed, so I'm going to fix our input objects here to remove GSE84185
  34. #Since we're fixing the meta-analysis anyway, I'm also going to make the minimum number of contrasts for inclusion in the meta-analysis the same as what we have for the hippocampus (11)
  36. colnames(MetaAnalysis_FoldChanges)
  38. # [1] "Rat_EntrezGene.ID" "Mouse_EntrezGene.ID" "MouseVsRat_EntrezGene.ID"
  39. # [4] "GSE26836_Amitriptyline" "GSE84183_fluoxetine" "GSE118670_Fluoxetine"
  40. # [7] "GSE28644_fluoxetine" "GSE93041_ketamine" "GSE81672_ketamine"
  41. # [10] "GSE81672_imipramine" "GSE150264_imipramine" "GSE84185_fluoxetine"
  42. # [13] "GSE168172_duloxetine" "GSE168172_sertraline" "GSE138802_ketamine"
  43. # [16] "GSE129359_duloxetine" "GSE45229_quetiapine_low_dose" "GSE45229_quetiapine_high_dose"
  44. # [19] "GSE230149_TMS" "GSE253280_MDMA"
  46. MetaAnalysis_FoldChanges<-MetaAnalysis_FoldChanges[,-12]
  48. colnames(MetaAnalysis_FoldChanges)
  49. # [1] "Rat_EntrezGene.ID" "Mouse_EntrezGene.ID" "MouseVsRat_EntrezGene.ID"
  50. # [4] "GSE26836_Amitriptyline" "GSE84183_fluoxetine" "GSE118670_Fluoxetine"
  51. # [7] "GSE28644_fluoxetine" "GSE93041_ketamine" "GSE81672_ketamine"
  52. # [10] "GSE81672_imipramine" "GSE150264_imipramine" "GSE168172_duloxetine"
  53. # [13] "GSE168172_sertraline" "GSE138802_ketamine" "GSE129359_duloxetine"
  54. # [16] "GSE45229_quetiapine_low_dose" "GSE45229_quetiapine_high_dose" "GSE230149_TMS"
  55. # [19] "GSE253280_MDMA"
  57. #correct now
  59. colnames(MetaAnalysis_SV)
  61. MetaAnalysis_SV<-MetaAnalysis_SV[,-12]
  63. colnames(MetaAnalysis_SV)
  64. # [1] "Rat_EntrezGene.ID" "Mouse_EntrezGene.ID" "MouseVsRat_EntrezGene.ID"
  65. # [4] "GSE26836_Amitriptyline" "GSE84183_fluoxetine" "GSE118670_Fluoxetine"
  66. # [7] "GSE28644_fluoxetine" "GSE93041_ketamine" "GSE81672_ketamine"
  67. # [10] "GSE81672_imipramine" "GSE150264_imipramine" "GSE168172_duloxetine"
  68. # [13] "GSE168172_sertraline" "GSE138802_ketamine" "GSE129359_duloxetine"
  69. # [16] "GSE45229_quetiapine_low_dose" "GSE45229_quetiapine_high_dose" "GSE230149_TMS"
  70. # [19] "GSE253280_MDMA"
  72. #Correct now
  74. save.image("~/Library/CloudStorage/[email hidden]/My Drive/BrainAlchemyProject/ProjectFolders/2024_EvaGeoghegan_Antidepressants_Hippocampus/ROutput_And_Results/Revisions/PFC/April212025Workspace_16Comparisons_6NAcutoff/Workspace_16comparisons6NAcutoff_CorrectComparisons.RData")
  76. #For documentation:
  77. write.csv(MetaAnalysis_FoldChanges, "MetaAnalysis_FoldChanges_PFC.csv")
  78. write.csv(MetaAnalysis_FoldChanges_ForMeta, "MetaAnalysis_FoldChanges_ForMeta_PFC.csv")
  79. write.csv(MetaAnalysis_SV, "MetaAnalysis_SV_PFC.csv")
  80. write.csv(MetaAnalysis_SV_ForMeta, "MetaAnalysis_SV_ForMeta_PFC.csv")
  83. ######################
  85. #Installing and loading relevant code packages:
  86. install.packages("metafor")
  87. library(metafor)
  88. library(plyr)
  90. ######################
  92. #Function:
  94. RunBasicMetaAnalysis<-function(NumberOfComparisons, CutOffForNAs, MetaAnalysis_FoldChanges, MetaAnalysis_SV){
  96. #The function first provides information about how many of the statistical contrasts have NA values as their differential expression results for each gene:
  97. MetaAnalysis_FoldChanges_NAsPerRow<-apply(MetaAnalysis_FoldChanges[,-c(1:3)], 1, function(y) sum(is.na(y)))
  99. print("Table of # of NAs per Row (Gene):")
  100. print(table(MetaAnalysis_FoldChanges_NAsPerRow))
  102. #Then any row (gene) that has too many NAs is removed from the analysis:
  103. MetaAnalysis_FoldChanges_ForMeta<<-MetaAnalysis_FoldChanges[MetaAnalysis_FoldChanges_NAsPerRow<CutOffForNAs,]
  104. MetaAnalysis_SV_ForMeta<<-MetaAnalysis_SV[MetaAnalysis_FoldChanges_NAsPerRow<CutOffForNAs,]
  106. print("MetaAnalysis_FoldChanges_ForMeta:")
  107. print(str(MetaAnalysis_FoldChanges_ForMeta))
  109. #I'm going to make an empty matrix to store the results of my meta-analysis:
  110. #This matrix was originally 6 columns
  111. #It was made larger to incorporate heterogeneity statistics (8 stats)
  112. #Then it was made larger to incorporate publication bias (3 stats)
  113. #...And even larger to incorporate robustness stats (4 stats)
  114. metaOutput<-matrix(NA, nrow(MetaAnalysis_FoldChanges_ForMeta), 21)
  116. influence_dfbs<-matrix(NA, nrow(MetaAnalysis_FoldChanges_ForMeta), ncol(MetaAnalysis_FoldChanges_ForMeta[-c(1:3)]))
  117. influence_cookd<-matrix(NA, nrow(MetaAnalysis_FoldChanges_ForMeta), ncol(MetaAnalysis_FoldChanges_ForMeta[-c(1:3)]))
  118. influence_TF<-matrix(NA, nrow(MetaAnalysis_FoldChanges_ForMeta), ncol(MetaAnalysis_FoldChanges_ForMeta[-c(1:3)]))
  120. colnames(influence_dfbs)<-colnames(MetaAnalysis_FoldChanges_ForMeta)[-c(1:3)]
  121. colnames(influence_cookd)<-colnames(MetaAnalysis_FoldChanges_ForMeta)[-c(1:3)]
  122. colnames(influence_TF)<-colnames(MetaAnalysis_FoldChanges_ForMeta)[-c(1:3)]
  124. row.names(influence_dfbs)<-MetaAnalysis_FoldChanges_ForMeta[,3]
  125. row.names(influence_cookd)<-MetaAnalysis_FoldChanges_ForMeta[,3]
  126. row.names(influence_TF)<-MetaAnalysis_FoldChanges_ForMeta[,3]
  128. #And then run a loop that run's a meta-analysis on the differential expression results (i.e., the columns that aren't annotation) for each gene (row):
  129. for(i in c(1:nrow(MetaAnalysis_FoldChanges_ForMeta))){
  131. print(i)
  133. #When pulling out the log2FC values and sampling variances (SV) for each gene, we use the function as.numeric to make sure they are in numeric matrix format because this is the required input format for the meta-analysis function that we will use:
  134. effect<-as.numeric(MetaAnalysis_FoldChanges_ForMeta[i,-c(1:3)])
  135. var<-as.numeric(MetaAnalysis_SV_ForMeta[i,-c(1:3)])
  137. #I added a function tryCatch that double-checks that the meta-analysis function (rma) doesn't produce errors (which breaks the loop):
  138. skip_to_next <- FALSE
  139. tryCatch(TempMeta<-rma(effect, var), error = function(e) {skip_to_next <<- TRUE})
  141. #If everything looks good, we move on to running the meta-analysis using a model that treats the variation in Log2FC across studies as random effects:
  142. if(skip_to_next){}else{
  143. TempMeta<-rma(effect, var)
  144. metaOutput[i, 1]<-TempMeta$b #gives estimate Log2FC
  145. metaOutput[i, 2]<-TempMeta$se #gives standard error
  146. metaOutput[i, 3]<-TempMeta$pval #gives pval
  147. metaOutput[i, 4]<-TempMeta$ci.lb #gives confidence interval lower bound
  148. metaOutput[i, 5]<-TempMeta$ci.ub #gives confidence interval upper bound
  149. metaOutput[i, 6]<-NumberOfComparisons-sum(is.na(effect))#Number of comparisons with data
  150. metaOutput[i, 7]<-TempMeta$k #Metafor output: number of studies (contrasts) included in the analysis - sanity check, should be the same as column 6
  151. metaOutput[i, 8]<-TempMeta$p #Metafor output: number of coefficients in model
  152. metaOutput[i, 9]<-TempMeta$tau2 #estimated amount of (residual) heterogeneity
  153. metaOutput[i, 10]<-TempMeta$se.tau2 #SE of the estimated amount of (residual) heterogeneity
  154. metaOutput[i, 11]<-TempMeta$QE #test statistic of the test for (residual) heterogeneity (Cochran’s Q-test)
  155. metaOutput[i, 12]<-TempMeta$QEp #p-value of the test for (residual) heterogeneity (Cochran’s Q-test)
  156. metaOutput[i, 13]<-TempMeta$I2 #the I 2 statistic, which estimates (in percent) how much of the total variability in the observed effect sizes or outcomes can be attributed to heterogeneity among the true effects
  157. metaOutput[i, 14]<-TempMeta$H2 #the H2 statistic, which estimates the ratio of the total amount of variability in the observed effect sizes or outcomes to the amount of sampling variability
  159. #Testing for evidence of publication bias (funnel plot asymmetry)
  160. skip_to_next <- FALSE
  161. tryCatch(PubBias<-regtest(TempMeta), error = function(e) {skip_to_next <<- TRUE})
  162. if(skip_to_next){}else{
  164. PubBias<-regtest(TempMeta) #The regression test by Egger et al. (1997)
  166. metaOutput[i, 15]<-PubBias$zval #the value of the Egger test statistic
  167. metaOutput[i, 16]<-PubBias$pval #the corresponding Egger test p-value
  168. metaOutput[i, 17]<-PubBias$dfs #the corresponding Egger test degrees of freedom
  169. rm(PubBias)
  170. }
  172. #Testing for robustness (sensitivity of the meta-analysis to the results of individual contrasts)
  173. skip_to_next <- FALSE
  174. tryCatch(Robustness<-leave1out(TempMeta), error = function(e) {skip_to_next <<- TRUE})
  175. if(skip_to_next){}else{
  177. Robustness<-leave1out(TempMeta) #In a leave-one-out analysis, the same model is repeatedly fitted, leaving out one study at a time. By doing so, we can assess how much the results are influenced by each individual study
  179. metaOutput[i, 18]<-min(Robustness$estimate, na.rm=TRUE) #minimum estimated effect in leave-one-out analysis
  180. metaOutput[i, 19]<-max(Robustness$estimate, na.rm=TRUE) #maximum estimated effect in leave-one-out analysis
  181. metaOutput[i, 20]<-min(Robustness$pval, na.rm=TRUE) #minimum pval in leave-one-out analysis
  182. metaOutput[i, 21]<-max(Robustness$pval, na.rm=TRUE) #maximum pval in leave-one-out analysis
  183. rm(Robustness)
  184. }
  186. #Testing for influential datapoints:
  187. skip_to_next <- FALSE
  188. tryCatch(InfluenceStats<-influence(TempMeta), error = function(e) {skip_to_next <<- TRUE})
  189. if(skip_to_next){}else{
  191. InfluenceStats<-influence(TempMeta) #Pulling out some stats for influential datapoints (contrasts)
  193. influence_dfbs[i,InfluenceStats$ids]<-InfluenceStats$dfbs$intrcpt
  195. influence_cookd[i,InfluenceStats$ids]<-InfluenceStats$inf$cook.d
  197. influence_TF[i,InfluenceStats$ids]<-InfluenceStats$is.infl
  199. rm(InfluenceStats)
  200. }
  202. rm(TempMeta)
  203. }
  204. rm(effect, var)
  205. }
  207. #Naming the columns in our output:
  208. colnames(metaOutput)<-c("Log2FC_estimate", "SE", "pval", "CI_lb", "CI_ub", "Number_Of_Comparisons", "Number_of_ Contrasts", "Number_of_Coefficients", "tau2_ResidualHeterogeneity", "SE_tau2_ResidualHeterogeneity", "QE_CochransQ_Teststat", "QEp_CochransQ_pval", "I2_PercentVar_TrueHeterogeneity", "H2_Ratio_EffectHetero_overSamplVar", "PubBias_Egger_Zstat", "PubBias_Egger_pval", "PubBias_Egger_DF", "Leave1Out_Min_Log2FC","Leave1Out_Max_Log2FC", "Leave1Out_Min_Pval","Leave1Out_Max_Pval")
  210. #The row names for our output are the combined mouse-rat entrez ids:
  211. row.names(metaOutput)<-MetaAnalysis_FoldChanges_ForMeta[,3]
  213. #We return this output back into our global environment
  214. metaOutput<<-metaOutput
  215. MetaAnalysis_Annotation<<-MetaAnalysis_FoldChanges_ForMeta[,c(1:3)]
  216. influence_dfbs<<-influence_dfbs
  217. influence_cookd<<-influence_cookd
  218. influence_TF<<-influence_TF
  220. return(metaOutput)
  221. return(MetaAnalysis_Annotation)
  222. return(influence_dfbs)
  223. return(influence_cookd)
  224. return(influence_TF)
  226. #... and provide the user with an update about the newly created object:
  228. print("metaOutput:")
  229. print(str(metaOutput))
  231. print("Top of metaOutput:")
  232. print(head(metaOutput))
  234. print("Bottom of metaOutput")
  235. print(tail(metaOutput))
  237. }
  239. ######################
  242. #Example Usage:
  244. NumberOfComparisons=16
  245. CutOffForNAs=6
  246. #I have 16 comparisons
  247. #5 NA is too many was our original code
  248. #MH: Ideally, this probably should match the HPC - so I changed it to 6 NAs max, with 11 comparisons min
  250. metaOutput<-RunBasicMetaAnalysis(NumberOfComparisons, CutOffForNAs, MetaAnalysis_FoldChanges, MetaAnalysis_SV)
  251. #Note: this function can take a while to run, especially if you have a lot of data
  252. #Plug in your computer, take a break, grab some coffee...
  254. ######################
  256. #Example Output:
  259. str(metaOutput)
  260. # num [1:15583, 1:21] -0.0068 -0.05028 -0.01136 0.0244 -0.00281 ...
  261. # - attr(*, "dimnames")=List of 2
  262. # ..$ : chr [1:15583] "23825_114087" "18585_191569" "66514_246307" "20480_65041" ...
  263. # ..$ : chr [1:21] "Log2FC_estimate" "SE" "pval" "CI_lb" ...
  264. #
  266. head(metaOutput)
  268. # Log2FC_estimate SE pval CI_lb CI_ub Number_Of_Comparisons
  269. # 23825_114087 -0.006797601 0.01234194 0.5817897 -0.03098736 0.01739216 16
  270. # 18585_191569 -0.050282937 0.05406125 0.3523139 -0.15624104 0.05567517 15
  271. # 66514_246307 -0.011362509 0.03010752 0.7058781 -0.07037217 0.04764715 16
  272. # 20480_65041 0.024402883 0.01819127 0.1797708 -0.01125135 0.06005712 16
  273. # 13726_25437 -0.002809407 0.02198028 0.8982955 -0.04588996 0.04027115 16
  274. # 16952_25380 0.104529113 0.09556532 0.2740438 -0.08277548 0.29183371 15
  275. # Number_of_ Contrasts Number_of_Coefficients tau2_ResidualHeterogeneity
  276. # 23825_114087 16 1 0.000000000
  277. # 18585_191569 15 1 0.034622852
  278. # 66514_246307 16 1 0.009398246
  279. # 20480_65041 16 1 0.002190424
  280. # 13726_25437 16 1 0.003170320
  281. # 16952_25380 15 1 0.078255969
  282. # SE_tau2_ResidualHeterogeneity QE_CochransQ_Teststat QEp_CochransQ_pval
  283. # 23825_114087 0.0006861533 6.264427 9.749554e-01
  284. # 18585_191569 0.0164109329 98.513674 9.132359e-15
  285. # 66514_246307 0.0050766920 47.894933 2.644191e-05
  286. # 20480_65041 0.0017510669 27.503192 2.489408e-02
  287. # 13726_25437 0.0026182397 28.621530 1.798692e-02
  288. # 16952_25380 0.0486414187 44.273614 5.347797e-05
  289. # I2_PercentVar_TrueHeterogeneity H2_Ratio_EffectHetero_overSamplVar PubBias_Egger_Zstat
  290. # 23825_114087 0.00000 1.000000 -0.224277068
  291. # 18585_191569 87.79823 8.195535 0.939284909
  292. # 66514_246307 79.17924 4.802898 0.007966891
  293. # 20480_65041 48.60820 1.945836 -0.433221252
  294. # 13726_25437 49.43049 1.977476 0.272109623
  295. # 16952_25380 73.00519 3.704416 -0.361673909
  296. # PubBias_Egger_pval PubBias_Egger_DF Leave1Out_Min_Log2FC Leave1Out_Max_Log2FC
  297. # 23825_114087 0.8225417 NA -0.009474074 -0.002006621
  298. # 18585_191569 0.3475845 NA -0.069023226 -0.012744959
  299. # 66514_246307 0.9936434 NA -0.022676290 0.010608597
  300. # 20480_65041 0.6648540 NA 0.007195653 0.032225433
  301. # 13726_25437 0.7855377 NA -0.011418324 0.005550514
  302. # 16952_25380 0.7175957 NA 0.054765112 0.148668538
  303. # Leave1Out_Min_Pval Leave1Out_Max_Pval
  304. # 23825_114087 0.4506641 0.8760513
  305. # 18585_191569 0.1952691 0.7146329
  306. # 66514_246307 0.4522039 0.9469670
  307. # 20480_65041 0.0858462 0.5809216
  308. # 13726_25437 0.5661581 0.9816990
  309. # 16952_25380 0.1201277 0.5449782
  311. tail(metaOutput)
  313. # Log2FC_estimate SE pval CI_lb CI_ub Number_Of_Comparisons
  314. # 83673_NA -0.001042998 0.01737201 0.95212455 -0.03509150 0.033005508 12
  315. # 93673_NA 0.071301867 0.05154496 0.16657451 -0.02972439 0.172328123 11
  316. # 93739_NA 0.005242300 0.01195665 0.66106573 -0.01819230 0.028676903 12
  317. # 93887_NA -0.030616422 0.02023995 0.13036285 -0.07028600 0.009053157 12
  318. # 97775_NA 0.030957476 0.02119296 0.14408589 -0.01057996 0.072494916 12
  319. # 99870_NA -0.078221628 0.02498038 0.00174021 -0.12718228 -0.029260980 12
  320. # Number_of_ Contrasts Number_of_Coefficients tau2_ResidualHeterogeneity
  321. # 83673_NA 12 1 4.729851e-06
  322. # 93673_NA 11 1 5.775706e-03
  323. # 93739_NA 12 1 1.220826e-05
  324. # 93887_NA 12 1 0.000000e+00
  325. # 97775_NA 12 1 0.000000e+00
  326. # 99870_NA 12 1 1.092589e-03
  327. # SE_tau2_ResidualHeterogeneity QE_CochransQ_Teststat QEp_CochransQ_pval
  328. # 83673_NA 0.0021979446 14.341406 0.21467434
  329. # 93673_NA 0.0099761851 15.998416 0.09967774
  330. # 93739_NA 0.0006108196 6.391316 0.84602120
  331. # 93887_NA 0.0021467157 6.523351 0.83626703
  332. # 97775_NA 0.0024535679 5.736181 0.89037428
  333. # 99870_NA 0.0025942165 15.504391 0.16055004
  334. # I2_PercentVar_TrueHeterogeneity H2_Ratio_EffectHetero_overSamplVar PubBias_Egger_Zstat
  335. # 83673_NA 0.05194076 1.000520 -0.04131045
  336. # 93673_NA 24.24320735 1.320014 2.62835275
  337. # 93739_NA 0.55460484 1.005577 1.12462809
  338. # 93887_NA 0.00000000 1.000000 0.01116656
  339. # 97775_NA 0.00000000 1.000000 -0.24323181
  340. # 99870_NA 15.42112766 1.182328 1.19660635
  341. # PubBias_Egger_pval PubBias_Egger_DF Leave1Out_Min_Log2FC Leave1Out_Max_Log2FC
  342. # 83673_NA 0.967048400 NA -0.0056245805 0.003855966
  343. # 93673_NA 0.008579949 NA 0.0142399379 0.118459201
  344. # 93739_NA 0.260746665 NA 0.0002093189 0.021542383
  345. # 93887_NA 0.991090563 NA -0.0481637033 -0.022645821
  346. # 97775_NA 0.807825822 NA 0.0235604830 0.034605070
  347. # 99870_NA 0.231460026 NA -0.0957609890 -0.060583847
  348. # Leave1Out_Min_Pval Leave1Out_Max_Pval
  349. # 83673_NA 0.7470933934 0.97238365
  350. # 93673_NA 0.0433114715 0.68403637
  351. # 93739_NA 0.2046787053 0.98652188
  352. # 93887_NA 0.0957746159 0.27411386
  353. # 97775_NA 0.1093217519 0.50643807
  354. # 99870_NA 0.0008481506 0.00964769
  357. write.csv(metaOutput, "metaOutput_wHeterogeneityPubBiasRobustMeasures.csv")
  358. write.csv(MetaAnalysis_Annotation, "MetaAnalysis_Annotation_for_metaOutput_wHeterogeneityPubBiasRobustMeasures.csv")
  360. colnames(metaOutput)
  362. write.csv(influence_dfbs, "influence_dfbs.csv")
  363. write.csv(influence_cookd, "influence_cookd.csv" )
  364. write.csv(influence_TF, "influence_TF.csv")
  366. #############################
  368. #I should probably run an FDR correction for all of these pvalues
  369. #That function needs to be adapted too
  372. if (!require("BiocManager", quietly = TRUE))
  373. install.packages("BiocManager")
  374. BiocManager::install("multtest")
  375. library(multtest)
  378. dim(metaOutput)
  379. colnames(metaOutput)
  381. ####################
  383. #Function:
  385. FalseDiscoveryCorrection<-function(metaOutput, HOM_MouseVsRat, MetaAnalysis_Annotation){
  387. #For the meta-analysis pval:
  389. #This calculates the false discovery rate, or q-value, for each of our p-values using the Benjamini-Hochberg procedure:
  390. tempPvalAdjMeta<-mt.rawp2adjp(metaOutput[,3], proc=c("BH"))
  392. #Then we put those results back into the order of our orginal output:
  393. metaPvalAdj<-tempPvalAdjMeta$adjp[order(tempPvalAdjMeta$index),]
  395. #And bind the false discovery rate (FDR) to the rest of the meta-analysis output:
  396. metaOutputFDR<-cbind(metaOutput, metaPvalAdj[,2])
  398. #And name that column FDR:
  399. colnames(metaOutputFDR)[22]<-"FDR"
  401. rm(tempPvalAdjMeta, metaPvalAdj)
  403. #For the QEp:
  405. #This calculates the false discovery rate, or q-value, for each of our p-values using the Benjamini-Hochberg procedure:
  406. tempPvalAdjMeta<-mt.rawp2adjp(metaOutput[,12], proc=c("BH"))
  408. #Then we put those results back into the order of our orginal output:
  409. metaPvalAdj<-tempPvalAdjMeta$adjp[order(tempPvalAdjMeta$index),]
  411. #And bind the false discovery rate (FDR) to the rest of the meta-analysis output:
  412. metaOutputFDR<-cbind(metaOutputFDR, metaPvalAdj[,2])
  414. #And name that column FDR:
  415. colnames(metaOutputFDR)[23]<-"QEp_FDR"
  417. rm(tempPvalAdjMeta, metaPvalAdj)
  420. #For the Egger pval
  422. #This calculates the false discovery rate, or q-value, for each of our p-values using the Benjamini-Hochberg procedure:
  423. tempPvalAdjMeta<-mt.rawp2adjp(metaOutput[,16], proc=c("BH"))
  425. #Then we put those results back into the order of our orginal output:
  426. metaPvalAdj<-tempPvalAdjMeta$adjp[order(tempPvalAdjMeta$index),]
  428. #And bind the false discovery rate (FDR) to the rest of the meta-analysis output:
  429. metaOutputFDR<-cbind(metaOutputFDR, metaPvalAdj[,2])
  431. #And name that column FDR:
  432. colnames(metaOutputFDR)[24]<-"Egger_FDR"
  434. rm(tempPvalAdjMeta, metaPvalAdj)
  437. #These results are returned to our global environment:
  438. metaOutputFDR<<-metaOutputFDR
  440. #We let the user know the basic structure of the meta-analysis output with FDR added to it (just to make sure everything still looks good)
  441. print("metaOutputFDR:")
  442. print(str(metaOutputFDR))
  444. #Then we make a dataframe that adds the annotation to that output:
  445. TempDF<-cbind.data.frame(metaOutputFDR, MetaAnalysis_Annotation)
  446. #And then adds even more detailed gene annotation:
  448. #First the detailed annotation for the mouse genes:
  449. TempDF2<-join(TempDF, HOM_MouseVsRat[,c(4:5,9:11)], by="Mouse_EntrezGene.ID", type="left", match="first")
  451. #Next the annnotation for the rat genes:
  452. TempDF3<-join(TempDF2, HOM_MouseVsRat[,c(15:16,20:22)], by="Rat_EntrezGene.ID", type="left", match="first")
  454. #This is renamed and returned to our global environment:
  455. metaOutputFDR_annotated<-TempDF3
  456. metaOutputFDR_annotated<<-metaOutputFDR_annotated
  458. #And written out into our working directory:
  459. write.csv(metaOutputFDR_annotated, "metaOutputFDR_annotated.csv")
  461. #Then we make a version of the output in order by p-value:
  462. metaOutputFDR_OrderbyPval<<-metaOutputFDR_annotated[order(metaOutputFDR_annotated[,5]),]
  464. #Let's write out a version of the output in order by p-value:
  465. write.csv(metaOutputFDR_OrderbyPval, "metaOutputFDR_orderedByPval.csv")
  467. #And give the user some information about their results:
  469. print("Do we have any genes that are statistically significant following traditional false discovery rate correction (FDR<0.05)?")
  470. print(sum(metaOutputFDR_annotated[,24]<0.05, na.rm=TRUE))
  472. print("What are the top results?")
  473. print(head(metaOutputFDR_annotated[order(metaOutputFDR_annotated[,5]),]))
  476. }
  479. FalseDiscoveryCorrection(metaOutput, HOM_MouseVsRat, MetaAnalysis_Annotation)
  481. # [1] "metaOutputFDR:"
  482. # num [1:15583, 1:24] -0.0068 -0.05028 -0.01136 0.0244 -0.00281 ...
  483. # - attr(*, "dimnames")=List of 2
  484. # ..$ : chr [1:15583] "23825_114087" "18585_191569" "66514_246307" "20480_65041" ...
  485. # ..$ : chr [1:24] "Log2FC_estimate" "SE" "pval" "CI_lb" ...
  486. # NULL
  487. # [1] "Do we have any genes that are statistically significant following traditional false discovery rate correction (FDR<0.05)?"
  488. # [1] 1
  489. # [1] "What are the top results?"
  490. # Rat_EntrezGene.ID Mouse_EntrezGene.ID Log2FC_estimate SE pval
  491. # 11966_117596 117596 11966 0.04871522 0.009256685 1.419497e-07
  492. # 26557_29547 29547 26557 -0.08048271 0.018976741 2.224035e-05
  493. # 16768_297596 297596 16768 0.10684491 0.025237295 2.299678e-05
  494. # 15275_25058 25058 15275 0.04066582 0.010019813 4.938024e-05
  495. # 227394_432363 432363 227394 -0.11618422 0.028947089 5.978050e-05
  496. # 22143_500929 500929 22143 0.04172231 0.010399663 6.023441e-05
  497. # CI_lb CI_ub Number_Of_Comparisons Number_of_ Contrasts Number_of_Coefficients
  498. # 11966_117596 0.03057245 0.06685799 16 16 1
  499. # 26557_29547 -0.11767644 -0.04328898 15 15 1
  500. # 16768_297596 0.05738072 0.15630910 15 15 1
  501. # 15275_25058 0.02102735 0.06030429 16 16 1
  502. # 227394_432363 -0.17291947 -0.05944897 14 14 1
  503. # 22143_500929 0.02133935 0.06210527 14 14 1
  504. # tau2_ResidualHeterogeneity SE_tau2_ResidualHeterogeneity QE_CochransQ_Teststat
  505. # 11966_117596 0.000000e+00 0.0003758395 16.04169
  506. # 26557_29547 0.000000e+00 0.0016113693 13.46161
  507. # 16768_297596 1.890744e-06 0.0028468385 17.85922
  508. # 15275_25058 0.000000e+00 0.0004066250 11.46787
  509. # 227394_432363 5.744005e-07 0.0037060975 13.90855
  510. # 22143_500929 0.000000e+00 0.0004729750 10.60494
  511. # QEp_CochransQ_pval I2_PercentVar_TrueHeterogeneity H2_Ratio_EffectHetero_overSamplVar
  512. # 11966_117596 0.3792863 0.000000000 1.000000
  513. # 26557_29547 0.4905413 0.000000000 1.000000
  514. # 16768_297596 0.2132679 0.015739554 1.000157
  515. # 15275_25058 0.7187780 0.000000000 1.000000
  516. # 227394_432363 0.3803157 0.003850019 1.000039
  517. # 22143_500929 0.6438748 0.000000000 1.000000
  518. # PubBias_Egger_Zstat PubBias_Egger_pval PubBias_Egger_DF Leave1Out_Min_Log2FC
  519. # 11966_117596 -1.1337063 0.2569178 NA 0.04384071
  520. # 26557_29547 -0.4287161 0.6681299 NA -0.09112955
  521. # 16768_297596 -0.6453592 0.5186945 NA 0.09577661
  522. # 15275_25058 -1.0429272 0.2969820 NA 0.03655320
  523. # 227394_432363 -1.1726924 0.2409192 NA -0.15494844
  524. # 22143_500929 0.4669214 0.6405561 NA 0.03308265
  525. # Leave1Out_Max_Log2FC Leave1Out_Min_Pval Leave1Out_Max_Pval FDR QEp_FDR
  526. # 11966_117596 0.05219320 4.748323e-08 0.0001265629 0.002212002 0.5064626
  527. # 26557_29547 -0.06942375 4.387497e-06 0.0028695781 0.106984945 0.6102590
  528. # 16768_297596 0.11851389 5.949228e-06 0.0200083873 0.106984945 0.3276177
  529. # 15275_25058 0.04417552 1.887013e-05 0.0012571081 0.106984945 0.8043718
  530. # 227394_432363 -0.10685521 5.471650e-05 0.0002522097 0.106984945 0.5074024
  531. # 22143_500929 0.04470365 3.396796e-05 0.0046424510 0.106984945 0.7447889
  532. # Egger_FDR MouseVsRat_EntrezGene.ID Mouse_Symbol Mouse_Genetic.Location
  533. # 11966_117596 0.9818284 11966_117596 Atp6v1b2 Chr8 cM
  534. # 26557_29547 1.0000000 26557_29547 Homer2 Chr7 cM
  535. # 16768_297596 1.0000000 16768_297596 Lag3 Chr6 cM
  536. # 15275_25058 0.9838184 15275_25058 Hk1 Chr10 cM
  537. # 227394_432363 0.9818284 227394_432363 Slco4c1 Chr1 cM
  538. # 22143_500929 1.0000000 22143_500929 Tuba1b Chr15 cM
  539. # Mouse_Genome.Coordinates..mouse..GRCm39.human..GRCh38.
  540. # 11966_117596 Chr8:69541388-69566370(+)
  541. # 26557_29547 Chr7:81250229-81356673(-)
  542. # 16768_297596 Chr6:124881324-124888668(-)
  543. # 15275_25058 Chr10:62104634-62215687(-)
  544. # 227394_432363 Chr1:96744918-96800027(-)
  545. # 22143_500929 Chr15:98829310-98832271(-)
  546. # Mouse_Name Rat_Symbol
  547. # 11966_117596 ATPase, H+ transporting, lysosomal V1 subunit B2 Atp6v1b2
  548. # 26557_29547 homer scaffolding protein 2 Homer2
  549. # 16768_297596 lymphocyte-activation gene 3 Lag3
  550. # 15275_25058 hexokinase 1 Hk1
  551. # 227394_432363 solute carrier organic anion transporter family, member 4C1 Slco4c1
  552. # 22143_500929 tubulin, alpha 1B Tuba1b
  553. # Rat_Genetic.Location Rat_Genome.Coordinates..mouse..GRCm39.human..GRCh38.
  554. # 11966_117596 Chr16 p14 NA
  555. # 26557_29547 Chr1 q31 NA
  556. # 16768_297596 Chr4 q42 NA
  557. # 15275_25058 Chr20 q11 NA
  558. # 227394_432363 Chr9 q36 NA
  559. # 22143_500929 Chr7 NA
  560. # Rat_Name
  561. # 11966_117596 ATPase H+ transporting V1 subunit B2
  562. # 26557_29547 homer scaffold protein 2
  563. # 16768_297596 lymphocyte activating 3
  564. # 15275_25058 hexokinase 1
  565. # 227394_432363 solute carrier organic anion transporter family, member 4C1
  566. # 22143_500929 tubulin, alpha 1B
  568. #################
  570. #Peeking at the results:
  572. colnames(metaOutputFDR_OrderbyPval)
  574. max(metaOutputFDR_OrderbyPval$Number_Of_Comparisons, na.rm=TRUE)
  575. #[1] 16
  577. min(metaOutputFDR_OrderbyPval$Number_Of_Comparisons, na.rm=TRUE)
  578. #[1] 11
  580. length(metaOutputFDR_OrderbyPval$Log2FC_estimate)
  581. #[1] 15583
  582. sum(is.na(metaOutputFDR_OrderbyPval$Log2FC_estimate)==FALSE)
  583. #[1] 15454
  585. sum(metaOutputFDR_OrderbyPval$FDR<0.05, na.rm=TRUE)
  586. #[1] 1
  588. sum(metaOutputFDR_OrderbyPval$FDR<0.10, na.rm=TRUE)
  589. #[1] 1
  591. metaOutputFDR_OrderbyPval$Mouse_Symbol[which(metaOutputFDR_OrderbyPval$FDR<0.05)]
  592. #[1] "Atp6v1b2"
  595. #histogram of I2
  596. hist(metaOutputFDR_OrderbyPval$I2_PercentVar_TrueHeterogeneity, breaks=100)
  597. #huge spike at 0
  599. summary(metaOutputFDR_OrderbyPval$I2_PercentVar_TrueHeterogeneity)
  600. # Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
  601. # 0.000 0.303 38.836 39.203 68.887 99.590 129
  603. #For comparison with HPC results:
  605. summary(metaOutputFDR_OrderbyPval$I2_PercentVar_TrueHeterogeneity[c(1:58)])
  606. # Min. 1st Qu. Median Mean 3rd Qu. Max.
  607. # 0.0000 0.0000 0.0234 9.9152 9.3123 82.1796
  609. summary(metaOutputFDR_OrderbyPval$I2_PercentVar_TrueHeterogeneity[-c(1:58)])
  610. # Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
  611. # 0.0000 0.4152 39.0024 39.3136 68.9891 99.5897 129
  614. #histogram of QEp_CochransQ_pval
  615. hist(metaOutputFDR_OrderbyPval$QEp_CochransQ_pval, breaks=100)
  616. #Large spike by 0 - so many genes have significant heterogeneity (not surprising)
  618. sum(metaOutputFDR_OrderbyPval$QEp_FDR<0.05, na.rm=TRUE)
  619. #[1] 6516
  621. sum(metaOutputFDR_OrderbyPval$QEp_FDR[metaOutputFDR_OrderbyPval$FDR<0.05]<0.05, na.rm=TRUE)
  622. #[1] 0
  624. sum(metaOutputFDR_OrderbyPval$QEp_CochransQ_pval[metaOutputFDR_OrderbyPval$FDR<0.05]<0.05, na.rm=TRUE)
  625. #[1] 0
  626. #but not our 1 sig gene lol
  628. #For comparison with HPC:
  629. sum(metaOutputFDR_OrderbyPval$QEp_CochransQ_pval[c(1:58)]<0.05, na.rm=TRUE)
  630. #[1] 7
  631. #Seven of the top 58 genes have nominally significant heterogeneity
  632. sum(metaOutputFDR_OrderbyPval$QEp_FDR[c(1:58)]<0.05, na.rm=TRUE)
  633. #[1] 4
  634. #four of the top 58 genes have heterogeneity that is sig at FDR<0.05
  636. #histogram of Egger_pval
  637. hist(metaOutputFDR_OrderbyPval$PubBias_Egger_pval, breaks=100)
  638. #Very small little spike at 0
  640. sum(metaOutputFDR_OrderbyPval$Egger_FDR<0.05, na.rm=TRUE)
  641. #[1] 7
  643. metaOutputFDR_OrderbyPval$Mouse_Symbol[which(metaOutputFDR_OrderbyPval$Egger_FDR<0.05)]
  644. #[1] "Pim1" "Gpr149" "Vgll2" "Hspa1a" "Cartpt" "Sox1" "Scg2"
  646. sum(metaOutputFDR_OrderbyPval$PubBias_Egger_pval<0.05 & metaOutputFDR_OrderbyPval$FDR<0.05, na.rm=TRUE)
  647. #[1] 0
  649. sum(metaOutputFDR_OrderbyPval$Egger_FDR<0.05 & metaOutputFDR_OrderbyPval$FDR<0.05, na.rm=TRUE)
  650. #[1] 0
  653. library(ggplot2)
  655. TempDF<-data.frame(PubBias=metaOutputFDR_OrderbyPval$Egger_FDR<0.05, Top58=c(rep(TRUE,58), rep(FALSE, (15583-58))), I2=metaOutputFDR_OrderbyPval$I2_PercentVar_TrueHeterogeneity)
  657. pdf("ViolinPlot_I2_vs_PubBias.pdf", height=5, width=4)
  658. ggplot(TempDF, aes(x=PubBias, y=I2)) +
  659. geom_violin(width=1.4) +
  660. geom_boxplot(width=0.1, color="grey", alpha=0.2)
  661. dev.off()
  663. pdf("ViolinPlot_I2_vs_Top58.pdf", height=5, width=4)
  664. ggplot(TempDF, aes(x=Top58, y=I2)) +
  665. geom_violin(width=1.4) +
  666. geom_boxplot(width=0.1, color="grey", alpha=0.2)
  667. dev.off()
  670. sum(metaOutputFDR_OrderbyPval$Leave1Out_Max_Pval>0.05 & metaOutputFDR_OrderbyPval$FDR<0.05, na.rm=TRUE)
  671. #[1] 0
  673. metaOutputFDR_OrderbyPval$Leave1Out_Max_Pval[which(metaOutputFDR_OrderbyPval$FDR<0.05)]
  674. #[1] 0.0001265629
  676. max(metaOutputFDR_OrderbyPval$pval[metaOutputFDR_OrderbyPval$FDR<0.05], na.rm=TRUE)
  677. #[1] 1.419497e-07
  679. sum(metaOutputFDR_OrderbyPval$Leave1Out_Min_Pval<1.419497e-07, na.rm=TRUE)
  680. #[1] 2
  682. metaOutputFDR_OrderbyPval$Mouse_Symbol[which(metaOutputFDR_OrderbyPval$Leave1Out_Min_Pval<1.419497e-07 & metaOutputFDR_OrderbyPval$FDR>0.05)]
  683. #[1] "Gapdhs"
  685. #How many weren't just borderline to begin with?
  686. metaOutputFDR_OrderbyPval$Mouse_Symbol[which(metaOutputFDR_OrderbyPval$Leave1Out_Min_Pval<1.419497e-07 & metaOutputFDR_OrderbyPval$FDR>0.10)]
  687. #[1] "Gapdhs"
  690. #... and how many of those aren't driven by a single study/contrast?
  691. metaOutputFDR_OrderbyPval$Mouse_Symbol[which(metaOutputFDR_OrderbyPval$Leave1Out_Min_Pval<1.419497e-07 & metaOutputFDR_OrderbyPval$Leave1Out_Max_Pval<0.05 & metaOutputFDR_OrderbyPval$FDR>0.10)]
  692. #[1] "Gapdhs"
  694. #Interesting.
  697. #Looking at the influence stats:
  699. str(influence_cookd)
  700. head(influence_cookd)
  702. str(influence_dfbs)
  703. head(influence_dfbs)
  705. str(influence_TF)
  706. head(influence_TF)
  708. boxplot(influence_cookd)
  709. boxplot(influence_dfbs)
  711. Influence_TF_Total<-apply(influence_TF, 2, function(y) sum(y, na.rm=TRUE))
  712. Influence_TF_Total
  713. write.csv(Influence_TF_Total, "Influence_TF_Total.csv")
  714. #The two extreme studies are GSE84183_fluoxetine and GSE28644_fluoxetine
  716. Influence_CooksD_MoreThan1<-apply(influence_cookd, 2, function(y) sum(y>1, na.rm=TRUE))
  717. Influence_CooksD_MoreThan1
  718. write.csv(Influence_CooksD_MoreThan1, "Influence_CooksD_MoreThan1.csv")
  720. Influence_Dfbs_MoreThan1<-apply(influence_dfbs, 2, function(y) sum(abs(y)>1, na.rm=TRUE))
  721. Influence_Dfbs_MoreThan1
  722. write.csv(Influence_Dfbs_MoreThan1, "Influence_Dfbs_MoreThan1.csv")
  725. colnames(MetaAnalysis_FoldChanges_ForMeta[,-c(1:3)])
  727. Covariates<-data.frame(row.names=colnames(MetaAnalysis_FoldChanges_ForMeta[,-c(1:3)]), ADType=rep(0, 16), Dissection=rep(0, 16), Platform=rep(0, 16), DepressionModel=rep(0, 16), SampleSize=rep(0, 16))
  729. #I think I better fill this in outside of R
  730. write.csv(Covariates, "PFC_Covariates.csv")
  732. #coding:
  733. # ADType: 0.5=NonTrad, -0.5=Trad
  734. # Dissection_ACG: 0.5=ACg, -0.5=PFC, Parietal, Cortex, 0=FC (could have ACg in the mix)
  735. # Dissection_PFC: 0.5=PFC, -0.5=ACg, Parietal, Cortex, 0=FC (could have PFC in the mix)
  736. # Platform: 0.5=microarray, -0.5=RNA-Seq
  737. # DepressionModel: 0.5=depression model included, -0.5=no depression model
  738. # Sample size: treatment + control for that contrast
  740. Covariates<-read.csv("PFC_Covariates.csv", header=TRUE, stringsAsFactors = FALSE)
  741. str(Covariates)
  743. Covariates_AndInfluence<-cbind.data.frame(Covariates, Influence_TF_Total, Influence_CooksD_MoreThan1, Influence_Dfbs_MoreThan1)
  745. cor(as.matrix(Covariates_AndInfluence[,-1]), method="spearman")
  747. write.csv(cor(as.matrix(Covariates_AndInfluence[,-1]), method="spearman"), "CorMatrix_Covariates_spearman.csv")
  749. cor(as.matrix(Covariates_AndInfluence[,-1]), method="pearson")
  751. write.csv(cor(as.matrix(Covariates_AndInfluence[,-1]), method="pearson"), "CorMatrix_Covariates_pearson.csv")
  754. str(MetaAnalysis_FoldChanges_ForMeta[,-c(1:3)])
  756. str(as.matrix(MetaAnalysis_FoldChanges_ForMeta[,-c(1:3)]))
  759. CorMatrix_Log2FC<-cor(as.matrix(MetaAnalysis_FoldChanges_ForMeta[,-c(1:3)]), method="spearman", use="pairwise.complete.obs")
  761. write.csv(CorMatrix_Log2FC, "CorMatrix_Log2FC.csv")
  763. heatmap(CorMatrix_Log2FC)
  764. #col = colorRampPalette(c("blue", "white", "red"))(20)
  766. pdf("Heatmap_CorMatrix_Log2FC_Spearman.pdf", width=8, height=8)
  767. heatmap(CorMatrix_Log2FC, margins = c(20, 20))
  768. dev.off()
  770. pdf("Scatterplot_InfluenceTF_vs_SampleSize.pdf", height=5, width=4)
  771. plot(Influence_TF_Total~SampleSize, data=Covariates_AndInfluence, xlab="Sample Size for Contrast", ylab="# of Genes for which Contrast is Deemed Influential (Outlier)")
  772. dev.off()
  774. pdf("Scatterplot_Influence_CooksD_MoreThan1_vs_SampleSize.pdf", height=5, width=4)
  775. plot(Influence_CooksD_MoreThan1~SampleSize, data=Covariates_AndInfluence, xlab="Sample Size for Contrast", ylab="# of Genes for which Contrast has Cook's D>1 (Outlier)")
  776. dev.off()
  778. pdf("Scatterplot_Influence_Dfbs_MoreThan1_vs_SampleSize.pdf", height=5, width=4)
  779. plot(Influence_Dfbs_MoreThan1~SampleSize, data=Covariates_AndInfluence, xlab="Sample Size for Contrast", ylab="# of Genes for which Contrast has |DFBeta|>1 (Outlier)")
  780. dev.off()
  783. ################
  785. #Redoing the heatmap of the top genes
  788. gene_top_50 <- rownames(metaOutputFDR_OrderbyPval)[1:50]
  789. names_genes_top_50<-paste("Mm", metaOutputFDR_OrderbyPval$Mouse_EntrezGene.ID[1:50], metaOutputFDR_OrderbyPval$Mouse_Symbol[1:50], ";", "Rn", metaOutputFDR_OrderbyPval$Rat_EntrezGene.ID[1:50], metaOutputFDR_OrderbyPval$Rat_Symbol[1:50], sep=" ")
  791. #Rows_Interest<-MetaAnalysis_FoldChanges_ForMeta$MouseVsRat_EntrezGene.ID%in%gene_top_50
  792. #Log2FC_Subsetted <- MetaAnalysis_FoldChanges_ForMeta[Rows_Interest, ]
  793. #Not in the original order:
  794. #cbind(Log2FC_Subsetted$MouseVsRat_EntrezGene.ID, gene_top_50)
  795. #row.names(Log2FC_Subsetted_Matrix)<-Log2FC_Subsetted$MouseVsRat_EntrezGene.ID
  797. Log2FC_Subsetted<-join(data.frame(MouseVsRat_EntrezGene.ID=gene_top_50), MetaAnalysis_FoldChanges_ForMeta, by="MouseVsRat_EntrezGene.ID", type="left")
  799. str(Log2FC_Subsetted)
  801. Log2FC_Subsetted_Matrix <- as.matrix(Log2FC_Subsetted[,-c(1:3)])
  803. row.names(Log2FC_Subsetted_Matrix)<-names_genes_top_50
  805. str(Log2FC_Subsetted_Matrix)
  807. library(pheatmap)
  808. library(dichromat)
  811. pdf("Heatmap_TopMetaGenes_50_color_scale_narrow.pdf",
  812. height = 11, width = 8.5)
  814. pheatmap(Log2FC_Subsetted_Matrix,
  815. color = colorRampPalette(c("#2166ac", "white", "#b2182b"))(100),
  816. scale = "none",
  817. breaks = seq(-1.5, 1.5, length.out = 101),
  818. cluster_rows = TRUE,
  819. cluster_cols = TRUE,
  820. fontsize_row = 8,
  821. fontsize_col = 8,
  822. width = 8.5,
  823. height = 11,
  824. border_color = NA)
  826. dev.off()
  829. ##############
  831. #Redoing the forest plot for Atp6v1b2:
  832. #Looks like we don't need to re-do it because it didn't include GSE84185 anyway
  834. ##############
  836. #Comparing HPC and PFC results:
  838. setwd("~/Library/CloudStorage/[email hidden]/My Drive/BrainAlchemyProject/ProjectFolders/2024_EvaGeoghegan_Antidepressants_Hippocampus/ROutput_And_Results/Revisions")
  840. list.files()
  842. metaOutputFDR_HPC<-read.csv("metaOutputFDR_orderedByPval_wHeterogeneityPubBiasRobustness.csv", header=TRUE, stringsAsFactors = FALSE)
  844. colnames(metaOutputFDR_HPC)[4:27]<-paste("HPC", colnames(metaOutputFDR_HPC)[4:27], sep="_")
  846. setwd("~/Library/CloudStorage/[email hidden]/My Drive/BrainAlchemyProject/ProjectFolders/2024_EvaGeoghegan_Antidepressants_Hippocampus/ROutput_And_Results/Revisions/PFC/April212025Workspace_16Comparisons_6NAcutoff")
  848. metaOutputFDR_PFC<-metaOutputFDR_OrderbyPval
  850. colnames(metaOutputFDR_PFC)[3:26]<-paste("PFC", colnames(metaOutputFDR_PFC)[3:26], sep="_")
  852. library(plyr)
  854. metaOutputFDR_HPC_PFC<-join(metaOutputFDR_HPC, metaOutputFDR_PFC, by="MouseVsRat_EntrezGene.ID", type="full")
  856. str(metaOutputFDR_HPC_PFC)
  857. #'data.frame': 17488 obs. of 60 variables
  859. pdf("Scatterplot_MetaAnalysisLog2FC_HPC_vs_PFC.pdf", height=5, width=4)
  860. plot(metaOutputFDR_HPC_PFC$PFC_Log2FC_estimate~metaOutputFDR_HPC_PFC$HPC_Log2FC_estimate, xlab="HPC: Antidepressant Log2FC", ylab="CTX: Antidepressant Log2FC")
  861. TrendLine<-lm(PFC_Log2FC_estimate~HPC_Log2FC_estimate, data=metaOutputFDR_HPC_PFC[is.na(metaOutputFDR_HPC_PFC$PFC_Log2FC_estimate)==FALSE & is.na(metaOutputFDR_HPC_PFC$HPC_Log2FC_estimate)==FALSE,])
  862. abline(TrendLine, col=2, lwd=3)
  863. dev.off()
  864. #Slight positive correlation
  866. summary.lm(TrendLine)
  868. # Call:
  869. # lm(formula = PFC_Log2FC_estimate ~ HPC_Log2FC_estimate, data = metaOutputFDR_HPC_PFC[is.na(metaOutputFDR_HPC_PFC$PFC_Log2FC_estimate) ==
  870. # FALSE & is.na(metaOutputFDR_HPC_PFC$HPC_Log2FC_estimate) ==
  871. # FALSE, ])
  872. #
  873. # Residuals:
  874. # Min 1Q Median 3Q Max
  875. # -0.36828 -0.01725 -0.00006 0.01740 0.37711
  876. #
  877. # Coefficients:
  878. # Estimate Std. Error t value Pr(>|t|)
  879. # (Intercept) 0.0036720 0.0003074 11.95 <2e-16 ***
  880. # HPC_Log2FC_estimate 0.2856445 0.0064575 44.23 <2e-16 ***
  881. # ---
  882. # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  883. #
  884. # Residual standard error: 0.03691 on 14434 degrees of freedom
  885. # Multiple R-squared: 0.1194, Adjusted R-squared: 0.1193
  886. # F-statistic: 1957 on 1 and 14434 DF, p-value: < 2.2e-16
  889. cor.test(metaOutputFDR_HPC_PFC$PFC_Log2FC_estimate, metaOutputFDR_HPC_PFC$HPC_Log2FC_estimate, method="spearman", use="pairwise.complete.obs")
  891. # Spearman's rank correlation rho
  892. #
  893. # data: metaOutputFDR_HPC_PFC$PFC_Log2FC_estimate and metaOutputFDR_HPC_PFC$HPC_Log2FC_estimate
  894. # S = 3.5589e+11, p-value < 2.2e-16
  895. # alternative hypothesis: true rho is not equal to 0
  896. # sample estimates:
  897. # rho
  898. # 0.2902152
  901. pdf("Scatterplot_MetaAnalysisLog2FC_HPC_vs_PFC_HPC_FDR10.pdf", height=5, width=4)
  902. plot(PFC_Log2FC_estimate~HPC_Log2FC_estimate,data=metaOutputFDR_HPC_PFC[metaOutputFDR_HPC_PFC$HPC_FDR<0.10,], xlab="HPC: Antidepressant Log2FC", ylab="CTX: Antidepressant Log2FC", pch=16, col="darkgrey")
  903. points(PFC_Log2FC_estimate~HPC_Log2FC_estimate,data=metaOutputFDR_HPC_PFC[metaOutputFDR_HPC_PFC$HPC_FDR<0.10,])
  904. TrendLine<-lm(PFC_Log2FC_estimate~HPC_Log2FC_estimate, data=metaOutputFDR_HPC_PFC[is.na(metaOutputFDR_HPC_PFC$PFC_Log2FC_estimate)==FALSE & is.na(metaOutputFDR_HPC_PFC$HPC_Log2FC_estimate)==FALSE & metaOutputFDR_HPC_PFC$HPC_FDR<0.10,])
  905. abline(TrendLine, col=2, lwd=3)
  906. dev.off()
  908. summary.lm(TrendLine)
  910. # Call:
  911. # lm(formula = PFC_Log2FC_estimate ~ HPC_Log2FC_estimate, data = metaOutputFDR_HPC_PFC[is.na(metaOutputFDR_HPC_PFC$PFC_Log2FC_estimate) ==
  912. # FALSE & is.na(metaOutputFDR_HPC_PFC$HPC_Log2FC_estimate) ==
  913. # FALSE & metaOutputFDR_HPC_PFC$HPC_FDR < 0.1, ])
  914. #
  915. # Residuals:
  916. # Min 1Q Median 3Q Max
  917. # -0.162829 -0.020470 -0.000796 0.022155 0.114044
  918. #
  919. # Coefficients:
  920. # Estimate Std. Error t value Pr(>|t|)
  921. # (Intercept) 0.004396 0.003608 1.218 0.225597
  922. # HPC_Log2FC_estimate 0.138171 0.036366 3.799 0.000233 ***
  923. # ---
  924. # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  925. #
  926. # Residual standard error: 0.03908 on 116 degrees of freedom
  927. # Multiple R-squared: 0.1107, Adjusted R-squared: 0.103
  928. # F-statistic: 14.44 on 1 and 116 DF, p-value: 0.0002325
  930. cor.test(metaOutputFDR_HPC_PFC$PFC_Log2FC_estimate[metaOutputFDR_HPC_PFC$HPC_FDR<0.10], metaOutputFDR_HPC_PFC$HPC_Log2FC_estimate[metaOutputFDR_HPC_PFC$HPC_FDR<0.10], method="spearman", use="pairwise.complete.obs")
  932. # Spearman's rank correlation rho
  933. #
  934. # data: metaOutputFDR_HPC_PFC$PFC_Log2FC_estimate[metaOutputFDR_HPC_PFC$HPC_FDR < 0.1] and metaOutputFDR_HPC_PFC$HPC_Log2FC_estimate[metaOutputFDR_HPC_PFC$HPC_FDR < 0.1]
  935. # S = 179286, p-value = 0.0001412
  936. # alternative hypothesis: true rho is not equal to 0
  937. # sample estimates:
  938. # rho
  939. # 0.345239
  942. pdf("Scatterplot_MetaAnalysisLog2FC_HPC_vs_PFC_HPC_FDR05.pdf", height=5, width=4)
  943. plot(PFC_Log2FC_estimate~HPC_Log2FC_estimate,data=metaOutputFDR_HPC_PFC[metaOutputFDR_HPC_PFC$HPC_FDR<0.05,], xlab="HPC: Antidepressant Log2FC", ylab="CTX: Antidepressant Log2FC", pch=16, col="darkgrey")
  944. points(PFC_Log2FC_estimate~HPC_Log2FC_estimate,data=metaOutputFDR_HPC_PFC[metaOutputFDR_HPC_PFC$HPC_FDR<0.05,])
  945. TrendLine<-lm(PFC_Log2FC_estimate~HPC_Log2FC_estimate, data=metaOutputFDR_HPC_PFC[is.na(metaOutputFDR_HPC_PFC$PFC_Log2FC_estimate)==FALSE & is.na(metaOutputFDR_HPC_PFC$HPC_Log2FC_estimate)==FALSE & metaOutputFDR_HPC_PFC$HPC_FDR<0.05,])
  946. abline(TrendLine, col=2, lwd=3)
  947. dev.off()
  949. summary.lm(TrendLine)
  951. # Call:
  952. # lm(formula = PFC_Log2FC_estimate ~ HPC_Log2FC_estimate, data = metaOutputFDR_HPC_PFC[is.na(metaOutputFDR_HPC_PFC$PFC_Log2FC_estimate) ==
  953. # FALSE & is.na(metaOutputFDR_HPC_PFC$HPC_Log2FC_estimate) ==
  954. # FALSE & metaOutputFDR_HPC_PFC$HPC_FDR < 0.05, ])
  955. #
  956. # Residuals:
  957. # Min 1Q Median 3Q Max
  958. # -0.162364 -0.024904 0.002058 0.028847 0.114328
  959. #
  960. # Coefficients:
  961. # Estimate Std. Error t value Pr(>|t|)
  962. # (Intercept) 0.002533 0.006423 0.394 0.6949
  963. # HPC_Log2FC_estimate 0.127907 0.057740 2.215 0.0313 *
  964. # ---
  965. # Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
  966. #
  967. # Residual standard error: 0.04541 on 50 degrees of freedom
  968. # Multiple R-squared: 0.08937, Adjusted R-squared: 0.07116
  969. # F-statistic: 4.907 on 1 and 50 DF, p-value: 0.03133
  971. cor.test(metaOutputFDR_HPC_PFC$PFC_Log2FC_estimate[metaOutputFDR_HPC_PFC$HPC_FDR<0.05], metaOutputFDR_HPC_PFC$HPC_Log2FC_estimate[metaOutputFDR_HPC_PFC$HPC_FDR<0.05], method="spearman", use="pairwise.complete.obs")
  972. #
  973. # Spearman's rank correlation rho
  974. #
  975. # data: metaOutputFDR_HPC_PFC$PFC_Log2FC_estimate[metaOutputFDR_HPC_PFC$HPC_FDR < 0.05] and metaOutputFDR_HPC_PFC$HPC_Log2FC_estimate[metaOutputFDR_HPC_PFC$HPC_FDR < 0.05]
  976. # S = 17282, p-value = 0.06059
  977. # alternative hypothesis: true rho is not equal to 0
  978. # sample estimates:
  979. # rho
  980. # 0.2622727
  982. metaOutputFDR_HPC_PFC$Mouse_Symbol[which(metaOutputFDR_HPC_PFC$HPC_FDR<0.05 & metaOutputFDR_HPC_PFC$PFC_pval<0.05)]
  983. #[1] "Tlr9" "Ccdc160" "Tnni1" "Stxbp5" "Dusp9" "Dmrtb1" "Insrr"
  986. ##############
  989. save.image("~/Library/CloudStorage/[email hidden]/My Drive/BrainAlchemyProject/ProjectFolders/2024_EvaGeoghegan_Antidepressants_Hippocampus/ROutput_And_Results/Revisions/PFC/April212025Workspace_16Comparisons_6NAcutoff/Workspace_16comparisons6NAcutoff_CorrectComparisons.RData")
  991. load("~/Library/CloudStorage/[email hidden]/My Drive/BrainAlchemyProject/ProjectFolders/2024_EvaGeoghegan_Antidepressants_Hippocampus/ROutput_And_Results/Revisions/PFC/April212025Workspace_16Comparisons_6NAcutoff/Workspace_16comparisons6NAcutoff_CorrectComparisons.RData")
  994. ##########
  997. #Re-running fGSEA - code copied from Brain.GMT github and tweaked to fit this analysis
  999. if (!require("BiocManager", quietly = TRUE))
  1000. install.packages("BiocManager")
  1002. BiocManager::install("fgsea", force = TRUE)
  1004. library(fgsea)
  1006. #This analysis assumes a differential expression (DE) output file structure similar to that produced by the Limma or EdgeR pipelines
  1007. #Rows=all genes included in the DE analysis, columns=gene annotation and DE statistical output
  1008. #At least one of the annotation columns must be official gene symbol
  1009. #At least one of the columns of differential statistics must include DE effect size (e.g., Log2 Fold Change)
  1011. #Read in the full DE results for a condition from the working directory
  1012. #Replace "DEResults.csv" in the code with your file name
  1013. DEResults<-metaOutputFDR_OrderbyPval
  1015. #Remove rows of DE results that are missing gene symbol annotation or effect size information
  1016. #Replace $gene_symbol in the code with the column name containing gene symbols in your DE output
  1017. #Replace $Log2FC in the code with the column name containing effect sizes in your DE output
  1018. DEResults_noNA<-DEResults[is.na(DEResults$Mouse_Symbol)==FALSE & is.na(DEResults$Log2FC_estimate)==FALSE,]
  1020. #The analysis only works if there is one effect size (e.g., log2 fold change or Log2FC) per gene symbol.
  1021. #One way to deal with multiple effect sizes mapping to the same gene (e.g., multiple transcripts or probes) is to average them:
  1022. #Replace $Log2FC in the code with the column name containing effect sizes in your DE output
  1023. #Replace $gene_symbol in the code with the column name containing gene symbols in your DE output
  1024. DEResults_Log2FC_forGSEA<-tapply(X=DEResults_noNA$Log2FC_estimate, INDEX=DEResults_noNA$Mouse_Symbol, FUN=mean)
  1025. names(DEResults_Log2FC_forGSEA)<-names(table(DEResults_noNA$Mouse_Symbol))
  1027. #The effect sizes should be ordered from smallest to largest:
  1028. #Replace $Log2FC in the code with the column name containing effect sizes in your DE output
  1029. DEResults_Log2FC_forGSEA_Ranked<-DEResults_Log2FC_forGSEA[order(DEResults_Log2FC_forGSEA)]
  1031. str(DEResults_Log2FC_forGSEA_Ranked)
  1032. # num [1:14791(1d)] -0.405 -0.369 -0.365 -0.362 -0.334 ...
  1033. # - attr(*, "dimnames")=List of 1
  1034. # ..$ : chr [1:14791] "Cav3" "Slc22a7" "Gpr101" "Slc9a4" ...
  1036. #Read in Brain.GMT for your species of interest (this example uses rat)
  1037. #If you get a warning about an incomplete line in the .gmt file, just ignore it
  1038. setwd("~/Library/CloudStorage/[email hidden]/My Drive/BrainAlchemyProject/ProjectFolders/2024_EvaGeoghegan_Antidepressants_Hippocampus/R_Code_And_Workspaces/Revisions_Code")
  1039. BrainGMT<-gmtPathways("BrainGMTv2_wGO_MouseOrthologs.gmt.txt")
  1042. setwd("~/Library/CloudStorage/[email hidden]/My Drive/BrainAlchemyProject/ProjectFolders/2024_EvaGeoghegan_Antidepressants_Hippocampus/ROutput_And_Results/Revisions/PFC/April212025Workspace_16Comparisons_6NAcutoff")
  1044. #Run fast fGSEA on your ranked, averaged effect sizes:
  1045. #This code should be compatible with updated fgsea packages - if you have an updated package, this code will run as fgseaSimple()
  1046. GSEA_Results<-fgsea(BrainGMT, DEResults_Log2FC_forGSEA_Ranked, nperm=10000, minSize = 10, maxSize = 1000)
  1048. #Pull out the names for the genes that are driving the enrichment of differential expression in each gene set:
  1049. GSEA_Results$leadingEdge<-vapply(GSEA_Results$leadingEdge, paste, collapse= ",", character(1L))
  1051. #Write out the results:
  1052. write.csv(GSEA_Results, "GSEA_Results.csv")
  1054. #You can easily view these results in Excel
  1055. # Sort by p-value
  1056. # padj: false discovery rate (FDR) corrected p-value. This value is normally used to set the threshold for significance (FDR<0.05)
  1057. # ES & NES: Enrichment Score and Normalized Enrichment Score for each gene set.
  1058. # Positive ES & NES values mean that the gene set is enriched with upregulation in response to your variable of interest
  1059. # Negative ES & NES values mean that the gene set is enriched with downregulation in response to your variable of interest
  1061. # Other aspects of the output can be deciphered by referencing the original GSEA publication: Subramanian et al. 2005
  1062. # https://www.pnas.org/doi/10.1073/pnas.0506580102
  1065. #Non-directional version:
  1067. DEResults_Log2FC_forGSEA_Ranked_NonDirectional<-DEResults_Log2FC_forGSEA[order(abs(DEResults_Log2FC_forGSEA))]
  1069. str(DEResults_Log2FC_forGSEA_Ranked_NonDirectional)
  1070. # num [1:14791(1d)] -1.74e-06 3.56e-06 4.11e-06 6.08e-06 -1.07e-05 ...
  1071. # - attr(*, "dimnames")=List of 1
  1072. # ..$ : chr [1:14791] "Areg" "Uprt" "Sbds" "Ep300" ...
  1074. #Run fast fGSEA on your ranked, averaged effect sizes:
  1075. #This code should be compatible with updated fgsea packages - if you have an updated package, this code will run as fgseaSimple()
  1076. GSEA_Results_NonDirectional<-fgsea(BrainGMT, DEResults_Log2FC_forGSEA_Ranked_NonDirectional, nperm=10000, minSize = 10, maxSize = 1000)
  1078. #Pull out the names for the genes that are driving the enrichment of differential expression in each gene set:
  1079. GSEA_Results_NonDirectional$leadingEdge<-vapply(GSEA_Results_NonDirectional$leadingEdge, paste, collapse= ",", character(1L))
  1081. #Write out the results:
  1082. write.csv(GSEA_Results_NonDirectional, "GSEA_Results_NonDirectional.csv")

CTX_RunMetaAnalysis_wHeterogeneity_PubBias_Robustness_Influence.R at commit bff5902, no license · at the source

Overview

Authors: Eva M Geoghegan1, Megan H Hagenauer2, Erin Hernandez3, Sophia Espinoza4, Elizabeth I Flandreau5, Phi T Nguyen1, Adrienne N Santiago1, Mubashshir Ra’eed Bhuiyan2, Sophie Mensch2, Stanley J Watson Jr2, Huda Akil2, René Hen1
ORCID iDs: Eva M Geoghegan, Megan H Hagenauer, Erin Hernandez, Sophia Espinoza, Elizabeth I Flandreau, Phi T Nguyen, Adrienne N Santiago, Mubashshir Ra’eed Bhuiyan, Stanley J Watson Jr, Huda Akil
  1. Columbia University, New York, New York, USA
  2. University of Michigan, Ann Arbor, Michigan, USA
  3. University of Chicago, Chicago, Illinois, USA
  4. Michigan State University, East Lansing, Michigan, USA
  5. Grand Valley State University, Allendale, Michigan, USA
Institutions: Columbia University (United States); University of Michigan (United States); University of Chicago (United States); Michigan State University (United States); Grand Valley State University (United States)
Journal: Journal of neurochemistry, volume 170, issue 7, article e70502
Dates: received 21 January 2026; accepted 12 June 2026; published online 3 July 2026; in print July 2026
Type: Review · Language: English
License: CC BY
Identifiers: DOI 10.1111/jnc.70502 · PMID 42396602 · PMCID PMC13329748 · OpenAlex W4409728392
Open access: hybrid, a free copy (OpenAlex)
Preprint: doi:10.1101/2025.04.21.648805
Status: code verified
Categories: genetics / omics (modality), mouse (organism), rat (organism), depression (population), cellular / molecular (subfield)
Methods: Statistics, Connectivity, fMRI & imaging, Physiology & signal measures
Keywords: antidepressant, hippocampus, meta‐analysis, RNA‐seq
MeSH: Antidepressive Agents*, Cerebral Cortex*, Gene Expression Profiling*, Hippocampus*, Animals, Mice, Rats (* major topic)
Topic: Nuclear Receptors and Signaling (Cellular and Molecular Neuroscience, Neuroscience), according to OpenAlex
Funding: Hope for Depression Research Foundation; Grinnell College Center for Careers, Life, and Service; National Institute on Drug Abuse (NIDA U01 DA043098); The Pritzker Neuropsychiatric Disorders Research Foundation; University of Michigan Undergraduate Research Opportunities Program; NIDA NIH HHS (NIDA U01 DA043098)
Citations: not cited yet (Europe PMC); 166 references in the paper

Abstract

Depression can be treated with traditional pharmaceuticals targeting monoaminergic function, nontraditional drug classes and neuromodulatory interventions. To identify mechanisms of action shared across clinically‐effective antidepressant treatment categories, we performed two systematic meta‐analyses of public transcriptional profiling data from adult laboratory rodents (rats, mice). The outcome variable was gene expression, measured by microarray or RNA‐Seq from bulk‐dissected tissue from two depression‐related brain regions (hippocampus, cortex). Relevant datasets were identified in the Gemma database of curated, reprocessed transcriptional profiling data using predefined search terms and inclusion/exclusion criteria (hippocampus: June 24, 2024, cortex: July 10, 2024). Differential expression results were extracted for all genes, minimizing bias. For each gene, a random effects meta‐analysis model was fit to antidepressant vs. control effect sizes (Log2 Fold Changes) from each study for each brain region, with follow‐up analyses exploring sources of effect heterogeneity. For the hippocampus, 15 relevant studies were identified, containing 22 antidepressant vs. control group comparisons (collective n = 313 samples), with approximately half representing traditional versus nontraditional antidepressants. Of 16 439 analyzed genes, 58 were consistently differentially expressed (False Discovery Rate (FDR) < 0.05) following treatment. Antidepressant effects were enriched in the dentate gyrus and in gene sets related to stress regulation, brain growth and plasticity, vasculature and glia, and immune function. Comparisons with single nucleus RNA‐Seq confirmed effects on specific hippocampal cell types, including potential rejuvenation of dentate granule neurons. For the cortex, 13 studies were identified, containing 16 antidepressant vs. control group comparisons (collective n = 233 samples). Of 15 583 analyzed genes, only one was consistently differentially expressed (FDR < 0.05: Atp6v1b2), but overall expression patterns moderately resembled the hippocampus. These genes and pathways showing consistent differential expression across treatment categories may be promising targets for novel therapies. Future work should explore relevance to human clinical populations and potential heterogeneity introduced by sex and subregion.

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 22 matches between paragraphs and lines of code.

evageoghegan/AntidepressantMetaAnalysis

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: bff5902031faf7c12a6501eeb067496f881205da, 5 May 2026
Languages: R (46)
Size: 53 files, 46 scripts
Software Heritage: not archived
Found in: the text, “Methods”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: tidyverse (17 files), metafor (12 files), clusterProfiler (6 files), circlize (2 files), ComplexHeatmap (2 files), ggplot2 (2 files), pheatmap (2 files)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
47 files
At the source: github.com/evageoghegan/AntidepressantMetaAnalysis

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;
  • 46 scripts, each with its path and the digest of its content;
  • 22 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

Datasets cited

Data Availability Statement

This paper centers on secondary data analysis using publicly available datasets (see Table 1 for full list of accession numbers).

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 2, 28 September 2026

  • Publisher: n/a → Wiley
  • Authors: added Eva M Geoghegan (0009-0008-3799-6944); Erin Hernandez (0009-0004-4468-7921); Sophia Espinoza (0009-0002-8517-6496); Elizabeth I Flandreau (0000-0003-3634-1694); Phi T Nguyen (0000-0001-5598-9213); Adrienne N Santiago (0000-0002-5022-1675); Stanley J Watson Jr (0000-0003-4980-5523); Huda Akil (0000-0003-0623-1056); removed Eva M Geoghegan; Erin Hernandez; Sophia Espinoza; Elizabeth I Flandreau; Phi T Nguyen; Adrienne N Santiago; Stanley J Watson Jr; Huda Akil

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 12 authors, 4 keywords, 7 MeSH terms, 6 funders, 156 references.

Cite

This paper

Geoghegan, E. M., Hagenauer, M. H., Hernandez, E., Espinoza, S., Flandreau, E. I., Nguyen, P. T., Santiago, A. N., Bhuiyan, M. R., Mensch, S., Watson, S. J., Akil, H., & Hen, R. (2026). The Converging Effects of Different Categories of Antidepressants on the Brain: A Systematic Meta-Analysis of Public Transcriptional Profiling Data From the Hippocampus and Cortex. Journal of neurochemistry, 170(7), e70502. https://doi.org/10.1111/jnc.70502

BibTeX

@article{geoghegan2026converging,
author = {Geoghegan, Eva M and Hagenauer, Megan H and Hernandez, Erin and Espinoza, Sophia and Flandreau, Elizabeth I and Nguyen, Phi T and Santiago, Adrienne N and Bhuiyan, Mubashshir Ra’eed and Mensch, Sophie and Watson, Stanley J and Akil, Huda and Hen, René},
title = {{The Converging Effects of Different Categories of Antidepressants on the Brain: A Systematic Meta-Analysis of Public Transcriptional Profiling Data From the Hippocampus and Cortex}},
journal = {Journal of neurochemistry},
year = {2026},
month = jul,
volume = {170},
number = {7},
pages = {e70502},
publisher = {Wiley},
issn = {0022-3042},
doi = {10.1111/jnc.70502},
url = {https://doi.org/10.1111/jnc.70502},
pmid = {42396602},
pmcid = {PMC13329748}
}

RIS

TY - JOUR
AU - Geoghegan, Eva M
AU - Hagenauer, Megan H
AU - Hernandez, Erin
AU - Espinoza, Sophia
AU - Flandreau, Elizabeth I
AU - Nguyen, Phi T
AU - Santiago, Adrienne N
AU - Bhuiyan, Mubashshir Ra’eed
AU - Mensch, Sophie
AU - Watson, Stanley J
AU - Akil, Huda
AU - Hen, René
TI - The Converging Effects of Different Categories of Antidepressants on the Brain: A Systematic Meta-Analysis of Public Transcriptional Profiling Data From the Hippocampus and Cortex
T2 - Journal of neurochemistry
J2 - J Neurochem
PY - 2026
DA - 2026/07/01
VL - 170
IS - 7
SP - e70502
SN - 0022-3042
PB - Wiley
DO - 10.1111/jnc.70502
UR - https://doi.org/10.1111/jnc.70502
LA - en
ER -

CSL-JSON

{
"id": "10.1111/jnc.70502",
"type": "article-journal",
"title": "The Converging Effects of Different Categories of Antidepressants on the Brain: A Systematic Meta-Analysis of Public Transcriptional Profiling Data From the Hippocampus and Cortex",
"container-title": "Journal of neurochemistry",
"author": [
{
"family": "Geoghegan",
"given": "Eva M"
},
{
"family": "Hagenauer",
"given": "Megan H"
},
{
"family": "Hernandez",
"given": "Erin"
},
{
"family": "Espinoza",
"given": "Sophia"
},
{
"family": "Flandreau",
"given": "Elizabeth I"
},
{
"family": "Nguyen",
"given": "Phi T"
},
{
"family": "Santiago",
"given": "Adrienne N"
},
{
"family": "Bhuiyan",
"given": "Mubashshir Ra’eed"
},
{
"family": "Mensch",
"given": "Sophie"
},
{
"family": "Watson",
"given": "Stanley J"
},
{
"family": "Akil",
"given": "Huda"
},
{
"family": "Hen",
"given": "René"
}
],
"container-title-short": "J Neurochem",
"volume": "170",
"issue": "7",
"page": "e70502",
"DOI": "10.1111/jnc.70502",
"PMID": "42396602",
"PMCID": "PMC13329748",
"ISSN": "0022-3042",
"publisher": "Wiley",
"URL": "https://doi.org/10.1111/jnc.70502",
"language": "en",
"issued": {
"date-parts": [
[
2026,
7,
1
]
]
}
}

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-71360-9 [code]
Perinatal brain developmental transition revealed by transcriptomic and proteomic analyses of Bama miniature pigs.
Journal: Nature communications
In common: circlize, clusterProfiler, ComplexHeatmap, 3 other tools, genetics / omics, 2 references
[2] 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: circlize, clusterProfiler, ComplexHeatmap, 3 other tools, cellular / molecular, 2 references
[3] doi:10.1038/s41380-026-03497-4 [code]
Transcriptome-informed brain cartography of polygenic risk and association with brain structure in major psychiatric disorders.
Journal: Molecular psychiatry
In common: circlize, clusterProfiler, ComplexHeatmap, 3 other tools, depression, genetics / omics, cellular / molecular, 1 reference
[4] doi:10.1038/s41593-026-02384-z [code]
cGAS-mediated type I IFN signaling contributes to disease progression in drug-refractory epilepsy.
Journal: Nature neuroscience
In common: circlize, clusterProfiler, ComplexHeatmap, 3 other tools, mouse, cellular / molecular, 1 reference
[5] doi:10.1038/s44318-026-00768-2 [code]
The H3K36me3 methyltransferase SETD2 contributes to PAF1C interactions with RNA Pol II and is required for neuronal differentiation.
Journal: The EMBO journal
In common: circlize, clusterProfiler, ComplexHeatmap, 3 other tools, mouse, 2 references
[6] doi:10.1172/jci.insight.207270 [code]
Progressive hypothalamic neuroinflammation in ovariectomized mice parallels aging-related transcriptomic changes in the female human hypothalamus.
Journal: JCI insight
In common: circlize, clusterProfiler, ComplexHeatmap, 3 other tools, genetics / omics, mouse, cellular / molecular, 1 reference
[7] doi:10.1016/j.isci.2026.116412 [code]
KOLF2.1J iTF-Microglia: A standardized platform to study microglial transcriptional regulatory networks in CNS disease.
Journal: iScience
In common: circlize, clusterProfiler, ComplexHeatmap, 3 other tools, genetics / omics, 1 reference
[8] doi:10.1038/s41467-026-76341-6 [code]
Neonatal inflammation disrupts a temporally restricted postnatal Numb-enriched microglial state in mice.
Journal: Nature communications
In common: circlize, clusterProfiler, ComplexHeatmap, 3 other tools, mouse, 1 reference
[9] doi:10.1038/s41380-026-03629-w [code]
Maternal fasting during early gestation induces epigenetic alterations and schizophrenia-related phenotypes.
Journal: Molecular psychiatry
In common: circlize, clusterProfiler, ComplexHeatmap, 3 other tools, genetics / omics, mouse, cellular / molecular, 1 reference
[10] doi:10.1186/s11689-026-09713-0 [code]
DRP1 mutations associated with EMPF1 encephalopathy perturb the transcriptional profile and maturation of cortical neurons.
Journal: Journal of neurodevelopmental disorders
In common: circlize, clusterProfiler, ComplexHeatmap, 3 other tools, 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.

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: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

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