Brain Data Alchemy Project: Meta-Analysis of Re-Analyzed Public Transcriptional Profiling Data in the Gemma Database (v.2024)

#Setting GitHub Credentials
#Megan Hagenauer, July 1, 2024

#Checking on my GitHub settings:
library("usethis")
git_sitrep()

#Example Output:
# ── Git global (user) 
# • Name: 'Megan Hastings Hagenauer'
# • Email: 'hagenaue@umich.edu'
# • Global (user-level) gitignore file: <unset>
#   • Vaccinated: FALSE
# ℹ See `?git_vaccinate` to learn more
# ℹ Defaulting to 'https' Git protocol
# • Default Git protocol: 'https'
# • Default initial branch name: <unset>
#   
#   ── GitHub user 
# • Default GitHub host: 'https://github.com'
# • Personal access token for 'https://github.com': <unset>
#   • To create a personal access token, call `create_github_token()`
# • To store a token for current and future use, call `gitcreds::gitcreds_set()`
# ℹ Read more in the 'Managing Git(Hub) Credentials' article:
#   https://usethis.r-lib.org/articles/articles/git-credentials.html
# 
# ── Active usethis project: '/Users/hagenaue/Documents/Example2024_BrainDataAlchemy' ──
# 
# ── Git local (project) 
# • Name: 'Megan Hastings Hagenauer'
# • Email: 'hagenaue@umich.edu'
# • Default branch: 'main'
# • Current local branch -> remote tracking branch:
#   'main' -> 'origin/main'
# 
# ── GitHub project 
# • Type = 'maybe_ours_or_theirs'
# • Host = 'https://github.com'
# • Config supports a pull request = NA
# • origin = 'hagenaue/Example2024_BrainDataAlchemy'
# • upstream = <not configured>
#   • Desc = 'origin' is a GitHub repo and 'upstream' is either not configured or is not a GitHub repo.
# 
# We may be offline or you may need to configure a GitHub personal access
# token. `gh_token_help()` can help with that.
# 
# Read more about what this GitHub remote configurations means at:
#   'https://happygitwithr.com/common-remote-setups.html'

create_github_token()

gitcreds::gitcreds_set()
# -> Adding new credentials...
# -> Removing credentials from cache...
# -> Done.

#Example of a complex query using proper Gemma.R API syntax:

MyQueryTerms<-"(stress* OR \"chronic social defeat\" ) NOT (\"oxidative stress\" OR \"ER Stress\")"

if(!requireNamespace("devtools", quietly=T)){
  install.packages("devtools")
}

devtools:: install_github("PavlidisLab/gemma.R", force=T)

install.packages("plyr")

library(plyr)

install.packages("tidyverse")

library(tidyverse)

MyQueryTerms<-"SSRI fluoxetine sertraline paroxetine citalopram escitalopram fluvoxamine vilazodone vortioxetine amitriptyline imipramine amoxapine desipramine nortriptyline clomipramine trimipramine protriptyline doxepin maprotiline trazadone nefazodone mirtazapine phenelzine nialamide isocarboxazid hydracarbazine tranylcypromine selegiline venlafaxine desvenlafaxine duloxetine levomilnacipran bupropion duloxetine levomilnacipran ketamine esketamine tianeptine brexanolone"


result_MyQueryTerms_RatsMice<- gemma.R ::get_datasets(query=MyQueryTerms, taxa = c("mouse", "rat")) %>% 
  gemma.R:::get_all_pages() 


#To examine the output of the search:
str(result_MyQueryTerms_RatsMice)
#Classes ‘data.table’ and 'data.frame': 111 obs. of  23 variables:

#You can also peruse the results in Rstudio by going to the Environment tab in the righthand corner and clicking on the object

#If you are struggling with getting an R function to work, you can get help using this code:

help("gemma.R ::get_datasets")
??gemma.R ::get_datasets

#To determine how many datasets in our search results are marked as "troubled" within the Gemma database we can make a table:

table(result_MyQueryTerms_NoFilter$experiment.troubled)
# FALSE 
# 194 
#We only want the non-troubled datasets (so FALSE) - looks good!

#We can filter out any troubled datasets using this code (even if you don't have troubled datasets, run this code anyway to make sure the downstream code works):

result_MyQueryTerms_Filtered<-result_MyQueryTerms_RatsMice %>% filter(experiment.troubled==FALSE)

#Outputting the search results in comma separate variable (.csv) file format
#CSV files are easily read into commonly used spreadsheet programs like Excel or GoogleSheets
write.csv(result_MyQueryTerms_Filtered, "result_MyQueryTerms_Filtered.csv")

#Where did my file go?
getwd()

#Let's find out what brain regions we have:

MyResults<-result_MyQueryTerms_Filtered

#Note: some datasets have data from more than one "organism part"
#That means that if we survey all of the organism parts in our datasets
#we may end up with a vector of organism parts that is of unknown length 
#and potentially longer than our number of datasets

#Creating an empty vector to save our results 
#We can just add on to this using concatenate
#That way we don't need to pre-specify length

OrganismPartAnnotations_All<-vector(mode="character")

#Looping over all of the datasets:

#How many datasets do we have?
length(MyResults$experiment.shortName)
nrow(MyResults)

for(i in c(1:nrow(MyResults))){
  
  ExperimentName<-MyResults$experiment.shortName[i]
  
  ExperimentAnnotations<-get_dataset_annotations(dataset=ExperimentName)
  rm(ExperimentName)
  
  #Making sure that there is "organism part" annotation for the dataset before attempting to extract it:
  
  #Do we have annotation for organism part? Here is a true/false vector
  #ExperimentAnnotations$class.name=="organism part"
  #This is the column of terms for the annotation
  #ExperimentAnnotations$term.name 
  #Example of subsetting for organism part:
  #ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="organism part"]
  
  #How many annotations do we have?
  length(ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="organism part"])
  
  if(length(ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="organism part"])>0){
    
  OrganismPartAnnotations<-ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="organism part"]
  rm(ExperimentAnnotations)
  
  OrganismPartAnnotations_All<-c(OrganismPartAnnotations_All, OrganismPartAnnotations)
  rm(OrganismPartAnnotations)
  
  #if there isn't organism part annotation, we just tell the loop to move on to the next dataset entry:
  
  }else{
    rm(ExperimentAnnotations)
  }
}

#This is what our result looks like:
OrganismPartAnnotations_All

#How many datasets have each organism part?  Let's make a summary table!

table(OrganismPartAnnotations_All)

#That is probably worth saving:

write.csv(table(OrganismPartAnnotations_All), "Table_OrganismPartAnnotations.csv")

#We can search for annotations related to hippocampus using the search_annotations function
#e.g., 
annots<-gemma.R ::search_annotations("hippocamp*")

#To view the results you can click on the annots object in the Environment tab 

#That search provides ontological annotations that include terms that start with hippocamp*
#if we want to determine which of these annotations are actually used in the Gemma database, we can run this code to count their usage:

annot_counts <- annots$value.URI %>% sapply(\(x){
  attributes(get_datasets(uris = x))$totalElements
})
annots$counts = annot_counts

#Then we can filter down to just the ontological terms that are actually used in the Gemma database:
annots_wRecords<-annots %>% filter(counts>0)

#You'll note that one of the frequently used terms is an "organism part" hippocampal formation
#That implies that this annotation is used to characterize the actual tissue used in the experiment
#It has the annotation ontological URI: http://purl.obolibrary.org/obo/UBERON_0002421

#If we use "UBERON_0002421" (hippocampal formation) in our search for datasets
# we will get not just results for hippocampal formation, but results matching its child terms as well

#If we wanted to see the hierarchy of ontological terms encompassed by "hippocampal formation" and/or encompassing "hippocampal formation"
#We can go to EMBL-EBI Ontology Lookup Service:
# https://www.ebi.ac.uk/ols4
#And put in the search box UBERON_0002421

#EMBL-EBI will tell us all of the various terms/concepts encompassed by "hippocampal formation" in the "Related from" section under "part of"
#e.g., layer of hippocampus, subiculum, fornix, Ammon's horn, etc

#Whereas to see if "hippocampal formation" can be broadened to include other useful nearby areas, we can click above it in the hierarchy ("cerebral cortex") to see our options

#We can also double check that those child terms for "hippocampal formation" are actually used within Gemma using this function:
get_child_terms("http://purl.obolibrary.org/obo/UBERON_0002421")

#Re-running the search to narrow our datasets down to a particular tissue or ROI:

#Example: Hippocampal formation:

#I'm going to use the URI ontology term for the hippocampus so that I catch all of datasets annotated with the child terms as well:
result_MyQueryTerms_RatsMice_ROI <- gemma.R ::get_datasets(query=MyQueryTerms, filter = 'allCharacteristics.valueUri in (http://purl.obolibrary.org/obo/UBERON_0002421)', taxa = c("mouse", "rat")) %>% 
  gemma.R:::get_all_pages() 

#Output:
str(result_MyQueryTerms_RatsMice_ROI)
#Classes ‘data.table’ and 'data.frame':	31 obs. of  23 variables:

#Filtering out the "troubled" results again:

result_MyQueryTerms_RatsMice_ROI_Filtered<-result_MyQueryTerms_RatsMice_ROI[result_MyQueryTerms_RatsMice_ROI$experiment.troubled==FALSE,]

#Output:

str(result_MyQueryTerms_RatsMice_ROI_Filtered)
#Classes ‘data.table’ and 'data.frame':	31 obs. of  23 variables:

#Adding additional annotation to this basic output from the getDatasets function
#e.g., info about the organism parts, developmental stage annotation, experimental factors used in the analyses

#To begin with, I renamed my results to something generic to make this code easier to functionalize later...
MyResults<-result_MyQueryTerms_RatsMice_ROI_Filtered

#Let's make some empty vectors that are the same length as the columns in our results
#These empty vectors will be used to store our annotations while we loop through the rows of datasets:
OrganismParts<-vector(mode="character", length=nrow(MyResults))
CellTypes<-vector(mode="character", length=nrow(MyResults))
DevelopmentalStages<-vector(mode="character", length=nrow(MyResults))
Treatments<-vector(mode="character", length=nrow(MyResults))
Diseases<-vector(mode="character", length=nrow(MyResults))
DiseaseModels<-vector(mode="character", length=nrow(MyResults))
Genotypes<-vector(mode="character", length=nrow(MyResults))
Strains<-vector(mode="character", length=nrow(MyResults))
Sex<-vector(mode="character", length=nrow(MyResults))

#I'm going to loop over all of the rows (row number =i) in my results (i.e., dataset metadata)
#And collect all of this annotation information
#And then format it in a way so that it can be added into my simple dataframe of results
#And then outputted and read easily in a spreadsheet program like excel

for(i in c(1:nrow(MyResults))){
  
  #Pulling out the name for the dataset in a row (row number=i):
  ExperimentName<-MyResults$experiment.shortName[i]
  
  #Accessing the annotations for the dataset:
  ExperimentAnnotations<-get_dataset_annotations(dataset=ExperimentName)
  #The number and type of annotations for the datasets is quite variable
  
  rm(ExperimentName)
  
  #Determining whether there is any annotation for organism part:
  
  if(length(ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="organism part"])>0){
    
    #If there is organism part annotation, I'm grabbing it:
    
    Annotations<-ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="organism part"]
    
    #And then collapsing that vector of annotations into a single string 
    #that can be easily stashed in a single cell in a data.frame (or Excel spreadsheet) 
    #This will eventually become part of the the row for that dataset in the results
    # e.g., "annotation 1; annotation 2; annotation 3"
    OrganismParts[i]<-paste(Annotations, collapse="; ")
    rm(Annotations)
    
  #If there isn't any annotation for organism part, we move on to the next type of annotation:
  }else{}
  
  #Now grabbing the annotation for cell type in a similar manner: 
  if(length(ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="cell type"])>0){
    
    Annotations<-ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="cell type"]
    
    CellTypes[i]<-paste(Annotations, collapse="; ")
    rm(Annotations)
    
  }else{ }
  
  #Now grabbing the annotation for developmental stage in a similar manner:
  if(length(ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="developmental stage"])>0){
    
    Annotations<-ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="developmental stage"]
    
    DevelopmentalStages[i]<-paste(Annotations, collapse="; ")
    rm(Annotations)
    
  }else{ }
  
  #Now grabbing the annotation for treatment in a similar manner:
  if(length(ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="treatment"])>0){
    
    Annotations<-ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="treatment"]
    
    Treatments[i]<-paste(Annotations, collapse="; ")
    rm(Annotations)
    
  }else{ }
  
  #Now grabbing the annotation for disease in a similar manner:
  if(length(ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="disease"])>0){
    
    Annotations<-ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="disease"]
    
    Diseases[i]<-paste(Annotations, collapse="; ")
    rm(Annotations)
    
  }else{ }
  
  #Now grabbing the annotation for disease model in a similar manner:
  if(length(ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="Disease model"])>0){
    
    Annotations<-ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="Disease model"]
    
    DiseaseModels[i]<-paste(Annotations, collapse="; ")
    rm(Annotations)
    
  }else{ }
  
  #Now grabbing the annotation for genotype in a similar manner:
  if(length(ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="genotype"])>0){
    
    Annotations<-ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="genotype"]
    
    Genotypes[i]<-paste(Annotations, collapse="; ")
    rm(Annotations)
    
  }else{ }
  
  #Now grabbing the annotation for strain in a similar manner:
  if(length(ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="strain"])>0){
    
    Annotations<-ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="strain"]
    
    Strains[i]<-paste(Annotations, collapse="; ")
    rm(Annotations)
    
  }else{ }
  
  #Now grabbing the annotation for biological sex in a similar manner:
  if(length(ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="biological sex"])>0){
    
    Annotations<-ExperimentAnnotations$term.name[ExperimentAnnotations$class.name=="biological sex"]
    
    Sex[i]<-paste(Annotations, collapse="; ")
    rm(Annotations)
    
  }else{ }
  
  rm(ExperimentAnnotations)
}

#Adding all of those vectors of annotation to my data.frame of results:
MyResults_Annotated<-cbind.data.frame(MyResults, 
                                OrganismParts,
                                CellTypes,
                                DevelopmentalStages,
                                Treatments,
                                Diseases,
                                DiseaseModels,
                                Genotypes,
                                Strains,
                                Sex)

#Let's add some empty columns for taking inclusion/exclusion notes:

ManipulationUnrelatedToTopic<-vector(mode="character", length=nrow(MyResults))
IncorrectDevelopmentalStage<-vector(mode="character", length=nrow(MyResults))
NotBulkDissection_ParticularCellTypeOrSubRegion<-vector(mode="character", length=nrow(MyResults))
NotFullTranscriptome_ChipSeq_TRAP_miRNA<-vector(mode="character", length=nrow(MyResults))
MetadataIssues_MissingInfo_NoPub_Retracted_Duplicated<-vector(mode="character", length=nrow(MyResults))

Excluded<-vector(mode="character", length=nrow(MyResults))
WhyExcluded<-vector(mode="character", length=nrow(MyResults))

MyResults_Annotated<-cbind.data.frame(MyResults_Annotated, ManipulationUnrelatedToTopic, IncorrectDevelopmentalStage, NotBulkDissection_ParticularCellTypeOrSubRegion, NotFullTranscriptome_ChipSeq_TRAP_miRNA, MetadataIssues_MissingInfo_NoPub_Retracted_Duplicated, Excluded, WhyExcluded)

#And then write out the results so that we can snoop through them in a spreadsheet program like Excel:
write.csv(MyResults_Annotated, "MyResults_Annotated.csv")

#We can download the processed expression data for any particular dataset using this code:
Expression<-gemma.R::get_dataset_processed_expression("GSE81672")
#In this case "GSE81672" is the dataset id
#For reference, this example is an RNA-Seq dataset

#Structure of the new Expression object:
str(Expression)
#Classes ‘data.table’ and 'data.frame':	35205 obs. of  103 variables:
#The first four columns are row metadata: Probe, GeneSymbol, GeneName, NCBIid
#The rest of the columns are gene expression values for each subject

#If you want to visualize the distribution of gene expression data for the entire study, you will need to grab all of the columns that aren't row metadata (i.e., not rows 1-4) and force them into the format of a numeric matrix first: 

ExpressionMatrix<-as.matrix(Expression[,-c(1:4)])

#You can visualize the distribution of the gene expression data for the dataset using a histogram
#We can make this pretty by adding a title and x-axis label
#You can also change the color and scaling
hist(ExpressionMatrix, main="Histogram", xlab="Log2 Expression", col="red", cex.axis=1.3, cex.lab=1.3)

#The range of the x-axis is the range of the log2 gene expression values
#The large spike on the left side of the histogram ("floor effect") are all of the genes that aren't truly expressed or have too low of expression to be measurable

#You can save the histogram using "export" in the Plots window

#We can also pull out numeric values summarizing the distribution, e.g.:
min(ExpressionMatrix, na.rm=TRUE)
#[1] -5.8601
median(ExpressionMatrix, na.rm=TRUE)
#[1] -2.1651
max(ExpressionMatrix, na.rm=TRUE)
#[1] 12.312

#Or to get more of an overview:
summary(ExpressionMatrix)

MyDatasets_Screened<-read.csv("MyDatasets_Screened.csv", stringsAsFactors = FALSE, header=TRUE)

#The format of the new object:

str(MyDatasets_Screened)
# 'data.frame':	3 obs. of  24 variables:
#   $ X                         : int  1 2 3
# $ experiment.shortName      : chr  "GSE126678" "GSE181285" "GSE205325"
# $ experiment.name           : chr  "Enduring and sex-specific changes in hippocampal gene expression after subchronic immune challenge" "Expression data from the hippocampus of LPS induced depression  mice model treated with Luteolin" "Modulation of behavioral and hippocampal transcriptomic responses in rat prolonged chronic unpredictable stress"| __truncated__
# $ experiment.ID             : int  15127 21341 24923
# $ experiment.description    : chr  " The goals of this study include examining long-lasting changes in hippocampal gene expression in males and in "| __truncated__ " Gene expression profiling reveals a potential role of Luteolin in LPS induced depression model LPS depression "| __truncated__ " Here, we examine behavioral and brain transcriptomic (RNA-seq) responses in rat prolonged chronic unpredictabl"| __truncated__
# $ experiment.troubled       : logi  FALSE FALSE FALSE
# $ experiment.accession      : chr  "GSE126678" "GSE181285" "GSE205325"
# $ experiment.database       : chr  "GEO" "GEO" "GEO"
# $ experiment.URI            : chr  "https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE126678" "https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE181285" "https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE205325"
# $ experiment.sampleCount    : int  24 16 21
# $ experiment.lastUpdated    : chr  "2023-12-16 09:57:45.253" "2023-12-17 13:14:26.142" "2023-12-18 00:16:39.838"
# $ experiment.batchEffectText: chr  "NO_BATCH_INFO" "NO_BATCH_INFO" "NO_BATCH_INFO"
# $ experiment.batchCorrected : logi  FALSE FALSE FALSE
# $ experiment.batchConfound  : int  0 0 0
# $ experiment.batchEffect    : int  0 0 0
# $ experiment.rawData        : int  1 1 1
# $ geeq.qScore               : num  0.282 0.283 0.28
# $ geeq.sScore               : num  1 0.812 0.75
# $ taxon.name                : chr  "mouse" "mouse" "rat"
# $ taxon.scientific          : chr  "Mus musculus" "Mus musculus" "Rattus norvegicus"
# $ taxon.ID                  : int  2 2 3
# $ taxon.NCBI                : int  10090 10090 10116
# $ taxon.database.name       : chr  "mm10" "mm10" "rn6"
# $ taxon.database.ID         : int  81 81 86


#Pull out the dataset ids (GSE#s) for your datasets:

ExperimentIDs<-MyDatasets_Screened$experiment.shortName

#Format of new ExperimentIDs object:

ExperimentIDs
#[1] "GSE126678" "GSE181285" "GSE205325"

#Alternatively, you can just enter them by hand as a vector:
#ExperimentIDs<-c("GSE126678", "GSE181285", "GSE205325")

GettingResultSetInfoForDatasets<-function(ExperimentIDs){
  
  #Making an empty data.frame to store results:
  
  ResultSets_toScreen<-data.frame(ExperimentID="NA",ResultSetIDs="NA", ContrastIDs="NA", ExperimentIDs="NA", FactorCategory="NA", ExperimentalFactors="NA", BaselineFactors="NA", Subsetted=FALSE, SubsetBy="NA")
  
  str(ResultSets_toScreen)
  # 'data.frame':	1 obs. of  9 variables:
  # $ ExperimentID       : chr "NA"
  # $ ResultSetIDs       : chr "NA"
  # $ ContrastIDs        : chr "NA"
  # $ ExperimentIDs      : chr "NA"
  # $ FactorCategory     : chr "NA"
  # $ ExperimentalFactors: chr "NA"
  # $ BaselineFactors    : chr "NA"
  # $ Subsetted          : logi FALSE
  # $ SubsetBy           : chr "NA"
  
  #We will then loop over each of the datasets:
  
  for(i in c(1:length(ExperimentIDs))){
    
    #For each dataset, we will use Gemma's API to access the experimental design info:
    Design<-gemma.R::get_dataset_differential_expression_analyses(ExperimentIDs[i])
    
    if(nrow(Design)>0){
      #Next, we'll make some empty vectors to store the experimental factor and baseline factor information for each result id for the dataset:
      ExperimentalFactors<-vector(mode="character", length(Design$result.ID))
      BaselineFactors<-vector(mode="character", length(Design$result.ID))
      
      #We will then loop over each of the result ids for the dataset:
      for(j in c(1:length(Design$result.ID))){
        
        #And grab the vector of experimental factors associated with that result id
        ExperimentalFactorVector<-Design$experimental.factors[[j]]$summary
        #And collapse that info down to a single entry that will fit in our data.frame
        ExperimentalFactors[j]<-paste(ExperimentalFactorVector, collapse="; ")
        
        #And then grab the vector of baseline/control/reference values associated with that result id
        BaselineFactorVector<-Design$baseline.factors[[j]]$summary
        #And collapse that info down to a single entry that will fit in our data.frame
        BaselineFactors[j]<-paste(BaselineFactorVector, collapse="; ")
      }
      
      #Some of the datasets are subsetted for the differential expression analyses
      #We will make an empty vector to store subset information for each result id
      SubsetBy<-vector(mode="character", length(Design$result.ID))
      
      #Then we will determine whether the dataset is subsetted:
      if(Design$isSubset[1]==TRUE){
        
        #If it is subsetted, we will loop over each result id for the dataset
        for (j in c(1:length(Design$result.ID))){
          
          #And grab the vector of subsetting information
          SubsetByVector<-Design$subsetFactor[[j]]$summary
          
          #And then collapse that information down to a single entry that will fit in our dataframe
          SubsetBy[j]<-paste(SubsetByVector, collapse="; ")
        }  
        
        #if the dataset wasn't subsetted for the differential expression analysis:
      }else{
        #We'll just make a vector of NA values to put in the "Subsetted by" column
        SubsetBy<-rep(NA, length((Design$result.ID)))
      }
      
      #Then we combine all of the information for all of the result sets for the dataset into a dataframe
      ResultSets_ForExperiment<-cbind.data.frame(ExperimentID=rep(ExperimentIDs[i],length(Design$result.ID)),ResultSetIDs=Design$result.ID, ContrastIDs=Design$contrast.ID, ExperimentIDs=Design$experiment.ID, FactorCategory=Design$factor.category, ExperimentalFactors, BaselineFactors, Subsetted=Design$isSubset, SubsetBy)
      
      #And add that information as rows to our data frame including the result set information for all datasets:
      ResultSets_toScreen<-rbind.data.frame(ResultSets_toScreen, ResultSets_ForExperiment)
      
      #Then clean up our space before looping to the next dataset:
      rm(ResultSets_ForExperiment, Design, ExperimentalFactors, BaselineFactors, SubsetBy)
      
    }else{
      rm(Design)
    }
    
  }
  
  #When we're done, we'll want to remove the initial (empty) row in our data.frame:
  ResultSets_toScreen<-ResultSets_toScreen[-1,]
  
  #We can make some empty vectors that we can use to store screening notes:
  Include<-vector(mode="character", length=nrow(ResultSets_toScreen))
  WrongBaseline<-vector(mode="character", length=nrow(ResultSets_toScreen))
  ResultsNotRegionSpecific<-vector(mode="character", length=nrow(ResultSets_toScreen))
  ReAnalyze<-vector(mode="character", length=nrow(ResultSets_toScreen))
  
  #And add them as columns to our dataframe:            
  ResultSets_toScreen<-cbind.data.frame(ResultSets_toScreen, Include, WrongBaseline, ResultsNotRegionSpecific, ReAnalyze)
  
  #And then write everything out as a .csv file that we can easily mark up in a spreadsheet program:
  write.csv(ResultSets_toScreen, "ResultSets_toScreen.csv")
  
  print("The Result Sets for your Datasets have been outputted into ResultSets_toScreen.csv")
  print(str(ResultSets_toScreen))
  
  #And clean up our environment:
  rm(Include, WrongBaseline, ResultsNotRegionSpecific, ReAnalyze)
}

#This function makes a data.frame that includes all of the result ids and contrast ids for each dataset with their basic metadata in a format (.csv) that is easily readable in a spreadsheet program.

#You can input this function by running the code discussed above to create the function in your R environment. 

#Alternatively, you can download the script for the function from our Github site and save the file in your working directory:
#https://github.com/hagenaue/BrainDataAlchemy/blob/main/MetaAnalysis_GemmaDEResults_2024/Function_GettingResultSetInfoForDatasets.R

#And then source it from your working directory:
source("Function_GettingResultSetInfoForDatasets.R")

#Example function usage:

GettingResultSetInfoForDatasets(ExperimentIDs)

#The results from this functions will be outputted as a .csv file in the working directory named "ResultSets_toScreen.csv"

ResultSet_contrasts<-read.csv("ResultSets_Screened.csv", header=TRUE, stringsAsFactors = FALSE )

DownloadingDEResults<-function(ResultSet_contrasts){
  
  #Some ResultSets have more than one statistical contrast, so they are present more than once in our data frame, e.g.:
  #ResultSet_contrasts$ResultSetIDs
  #[1] 553805 553805 553805 570552 556647
  
  #To pull down the statistical results, we'll only want the unique result set ids:
  UniqueResultSetIDs<-unique(ResultSet_contrasts$ResultSetIDs)
  
  print("These are the result sets that you identified as being of interest")
  print(UniqueResultSetIDs)
  #553805 570552 556647
  
  differentials <- UniqueResultSetIDs %>% lapply(function(x){
    #   # take the first and only element of the output. the function returns a list 
    #   # because single experiments may have multiple resultSets. Here we use the 
    #   # resultSet argument to directly access the results we need
    get_differential_expression_values(resultSet = x)[[1]]
  })
  
  str(differentials)
  #That code worked. Excellent!
  
  # # some datasets might not have all the advertised differential expression results
  # # calculated due to a variety of factors. here we remove the empty differentials
  missing_contrasts <- differentials %>% sapply(nrow) %>% {.==0}
  #[1] FALSE FALSE FALSE
  differentials <<- differentials[!missing_contrasts]
  UniqueResultSetIDs<<-UniqueResultSetIDs[!missing_contrasts]
  
  print("These are the result sets that had differential expression results:")
  print(UniqueResultSetIDs)
  
  print("Your differential expression results for each of your result sets are stored in the object named differentials. This object is structured as a list of data frames. Each element in the list represetns a result set, with the data frame containing the differential expression results")
  
  #Within any particular Result Set, there are likely to be some contrasts that we want and others that we don't want
  #For example, a result set might contain a variety of stress interventions
  #And maybe we only want the acute stress contrast results
  
  #We already identified which statistical contrasts we wanted during our screening:
  #This is the object with the specific contrast ids that we want:
  #ResultSet_contrasts$ContrastIDs
  
  #Which will be these columns within the listed dataframes of differential expression results:
  
  print("These are the columns for the effect sizes for our statistical contrasts of interest (Log(2) Fold Changes")
  Contrasts_Log2FC<<-paste("contrast_", ResultSet_contrasts$ContrastIDs, "_log2fc", sep="")
  
  print(Contrasts_Log2FC)
  #[1] "contrast_151618_log2fc" "contrast_151617_log2fc" "contrast_151619_log2fc"
  #[4] "contrast_186753_log2fc" "contrast_204289_log2fc"
  
  print("these are the columns for the T-statistics for our statistical contrasts of interest - we will use that information to derive the sampling variances")
  
  Contrasts_Tstat<<-paste("contrast_", ResultSet_contrasts$ContrastIDs, "_tstat", sep="")
  
  print(Contrasts_Tstat)
  # [1] "contrast_151618_tstat" "contrast_151617_tstat" "contrast_151619_tstat"
  # [4] "contrast_186753_tstat" "contrast_204289_tstat"  
  
}

#You can input this function by running the code discussed above to create the function in your R environment. 

#Alternatively, you can download the script for the function from our Github site and save the file in your working directory:
#https://github.com/hagenaue/BrainDataAlchemy/blob/main/MetaAnalysis_GemmaDEResults_2024/Function_DownloadingDEResults.R

#And then source it from your working directory:
source("Function_DownloadingDEResults.R")

#Example usage for the function:

DownloadingDEResults(ResultSet_contrasts)

#Example output:

# [1] "These are the result sets that you identified as being of interest"
# [1] 553805 570552 556647
# List of 3
# $ :Classes ‘data.table’ and 'data.frame':	21693 obs. of  19 variables:
#   ..$ Probe                      : int [1:21693] 20344 20558 110454 20339 17067 68774 21946 21366 76905 16190 ...
# ..$ NCBIid                     : int [1:21693] 20344 20558 110454 20339 17067 68774 21946 21366 76905 16190 ...
# ..$ GeneSymbol                 : chr [1:21693] "Selp" "Slfn4" "Ly6a" "Sele" ...
# ..$ GeneName                   : chr [1:21693] "selectin, platelet" "schlafen 4" "lymphocyte antigen 6 family member A" "selectin, endothelial cell" ...
# ..$ pvalue                     : num [1:21693] 1.70e-10 9.11e-09 5.38e-08 7.73e-08 1.00e-07 ...
# ..$ corrected_pvalue           : num [1:21693] 3.69e-06 9.88e-05 4.00e-04 4.00e-04 4.00e-04 ...
# ..$ rank                       : num [1:21693] 4.61e-05 9.22e-05 1.00e-04 2.00e-04 2.00e-04 ...
# ..$ contrast_151617_coefficient: num [1:21693] 3.102 3.27 0.763 2.868 0.513 ...
# ..$ contrast_151617_log2fc     : num [1:21693] 3.102 3.27 0.763 2.868 0.513 ...
# ..$ contrast_151617_tstat      : num [1:21693] 2.97 3.62 3.86 2.31 4.17 ...
# ..$ contrast_151617_pvalue     : num [1:21693] 0.0077 0.0017 0.001 0.0318 0.0005 ...
# ..$ contrast_151618_coefficient: num [1:21693] 7.01 5.21 1.35 5.86 0.82 ...
# ..$ contrast_151618_log2fc     : num [1:21693] 7.01 5.21 1.35 5.86 0.82 ...
# ..$ contrast_151618_tstat      : num [1:21693] 7.81 6.14 7.02 5.28 6.78 ...
# ..$ contrast_151618_pvalue     : num [1:21693] 1.86e-07 5.64e-06 9.01e-07 3.78e-05 1.47e-06 ...
# ..$ contrast_151619_coefficient: num [1:21693] 1.223 -0.303 -0.243 0.118 -0.214 ...
# ..$ contrast_151619_log2fc     : num [1:21693] 1.223 -0.303 -0.243 0.118 -0.214 ...
# ..$ contrast_151619_tstat      : num [1:21693] 1.0192 -0.2593 -1.1221 0.0785 -1.6535 ...
# ..$ contrast_151619_pvalue     : num [1:21693] 0.32 0.798 0.275 0.938 0.114 ...
# ..- attr(*, ".internal.selfref")=<externalptr> 
#   ..- attr(*, "call")= chr "https://gemma.msl.ubc.ca/rest/v2/resultSets/553805"
# ..- attr(*, "env")=<environment: 0x7f897aec2838> 
#   $ :Classes ‘data.table’ and 'data.frame':	45100 obs. of  11 variables:
#   ..$ Probe                      : chr [1:45100] "1426114_PM_at" "1452267_PM_at" "1441388_PM_at" "1436221_PM_at" ...
# ..$ NCBIid                     : chr [1:45100] "" "224613" "" "100039795" ...
# ..$ GeneSymbol                 : chr [1:45100] "" "Flywch1" "" "Ildr2" ...
# ..$ GeneName                   : chr [1:45100] "" "FLYWCH-type zinc finger 1" "" "immunoglobulin-like domain containing receptor 2" ...
# ..$ pvalue                     : num [1:45100] 3.74e-09 1.57e-08 1.32e-08 8.87e-09 2.37e-08 ...
# ..$ corrected_pvalue           : num [1:45100] 2e-04 2e-04 2e-04 2e-04 2e-04 3e-04 3e-04 3e-04 4e-04 5e-04 ...
# ..$ rank                       : num [1:45100] 2.22e-05 8.87e-05 6.65e-05 4.44e-05 1.00e-04 ...
# ..$ contrast_186753_coefficient: num [1:45100] 0.1587 0.0967 -0.5049 0.7694 0.1039 ...
# ..$ contrast_186753_log2fc     : num [1:45100] 0.1587 0.0967 -0.5049 0.7694 0.1039 ...
# ..$ contrast_186753_tstat      : num [1:45100] 1.65 2.28 -5.4 10.4 2.11 ...
# ..$ contrast_186753_pvalue     : num [1:45100] 1.20e-01 3.83e-02 8.18e-05 3.93e-08 5.26e-02 ...
# ..- attr(*, ".internal.selfref")=<externalptr> 
#   ..- attr(*, "call")= chr "https://gemma.msl.ubc.ca/rest/v2/resultSets/570552"
# ..- attr(*, "env")=<environment: 0x7f897ab59f78> 
#   $ :Classes ‘data.table’ and 'data.frame':	18964 obs. of  11 variables:
#   ..$ Probe                      : int [1:18964] 497918 303604 501110 498386 299052 25636 65166 102548697 100364561 24183 ...
# ..$ NCBIid                     : int [1:18964] 497918 303604 501110 498386 299052 25636 65166 102548697 100364561 24183 ...
# ..$ GeneSymbol                 : chr [1:18964] "Smcr8" "Map3k3" "Gsta6" "Cc2d2a" ...
# ..$ GeneName                   : chr [1:18964] "SMCR8-C9orf72 complex subunit" "mitogen activated protein kinase kinase kinase 3" "glutathione S-transferase alpha 6" "coiled-coil and C2 domain containing 2A" ...
# ..$ pvalue                     : num [1:18964] 0.869 0.885 0.23 0.597 0.518 ...
# ..$ corrected_pvalue           : num [1:18964] 1 1 1 1 1 ...
# ..$ rank                       : num [1:18964] 0.855 0.875 0.191 0.569 0.486 ...
# ..$ contrast_204289_coefficient: num [1:18964] 0.0982 -0.0239 -0.2088 -0.1256 -0.3315 ...
# ..$ contrast_204289_log2fc     : num [1:18964] 0.0982 -0.0239 -0.2088 -0.1256 -0.3315 ...
# ..$ contrast_204289_tstat      : num [1:18964] 0.167 -0.146 -1.242 -0.538 -0.659 ...
# ..$ contrast_204289_pvalue     : num [1:18964] 0.869 0.885 0.23 0.597 0.518 ...
# ..- attr(*, ".internal.selfref")=<externalptr> 
#   ..- attr(*, "call")= chr "https://gemma.msl.ubc.ca/rest/v2/resultSets/556647"
# ..- attr(*, "env")=<environment: 0x7f88feede760> 

#   [1] "These are the result sets that had differential expression results:"
# [1] 553805 570552 556647

# [1] "Your differential expression results for each of your result sets are stored in the object named differentials. This object is structured as a list of data frames. Each element in the list represetns a result set, with the data frame containing the differential expression results"

# [1] "These are the columns for the effect sizes for our statistical contrasts of interest (Log(2) Fold Changes"
# [1] "contrast_151618_log2fc" "contrast_151617_log2fc" "contrast_151619_log2fc"
# [4] "contrast_186753_log2fc" "contrast_204289_log2fc"

# [1] "these are the columns for the T-statistics for our statistical contrasts of interest - we will use that information to derive the sampling variances"
# [1] "contrast_151618_tstat" "contrast_151617_tstat" "contrast_151619_tstat"
# [4] "contrast_186753_tstat" "contrast_204289_tstat"

SavingGemmaDEResults_forEachResultSet<-function(differentials, UniqueResultSetIDs, ResultSet_contrasts){
  
  for (i in c(1:length(differentials))){
    
    ThisResultSet<-UniqueResultSetIDs[i]
    
    #Pulling out the dataset name from our other data frame
    #For some reason I can find this in the Gemma ResultSet differential expression output
    #Since some datasets have multiple result sets, we just grab the dataset name from the first entry
    ThisDataSet<-ResultSet_contrasts$ExperimentID[ResultSet_contrasts$ResultSetIDs==ThisResultSet][1] 
    
    #Write out a data frame containing the differential expression output for the result set
    #And name it with the dataset id and result set id:
    write.csv(differentials[[i]], paste("DEResults", ThisDataSet, ThisResultSet, ".csv", sep="_"))
    
    rm(ThisDataSet, ThisResultSet)
  }
}

#You can input this function by running the code discussed above to create the function in your R environment. 

#Alternatively, you can download the script for the function from our Github site and save the file in your working directory:
#https://github.com/hagenaue/BrainDataAlchemy/blob/main/MetaAnalysis_GemmaDEResults_2024/Function_SavingGemmaDEResults_forEachResultSet.R

#And then source it from your working directory:
source("Function_SavingGemmaDEResults_forEachResultSet.R")

#Example usage:

SavingGemmaDEResults_forEachResultSet(differentials, UniqueResultSetIDs, ResultSet_contrasts)

#Here is an example of pulling out the results of the first result set in the differentials object:

DE_Results<-differentials[[1]]

#The structure of the new DE_Results object:
str(DE_Results)
#Classes ‘data.table’ and 'data.frame':	21693 obs. of  19 variables:


#After finishing processing that first result set (steps #9.4-#9.8), we would move on to the second result set.
#To pull out the results for the second result set, we would use this code:

DE_Results<-differentials[[2]]

FilteringDEResults_GoodAnnotation<-function(DE_Results){
  
  print("# of rows in results")
  print(nrow(DE_Results))
  
  print("# of rows with missing NCBI annotation:")
  print(sum(DE_Results$NCBIid==""|DE_Results$NCBIid=="null"))
  
  print("# of rows with NA NCBI annotation:")
  print(sum(is.na(DE_Results$NCBIid)))
  
  print("# of rows with missing Gene Symbol annotation:")
  print(sum(DE_Results$GeneSymbol==""|DE_Results$GeneSymbol=="null"))
  
  print("# of rows mapped to multiple NCBI_IDs:")
  print(length(grep('\\|', DE_Results$NCBIid)))
  
  print("# of rows mapped to multiple Gene Symbols:")
  print(length(grep('\\|', DE_Results$GeneSymbol)))
  
  #I only want the subset of data which contains rows that do not contain an NCBI EntrezID of ""
  DE_Results_NoNA<-DE_Results[(DE_Results$NCBIid==""|DE_Results$NCBIid=="null")==FALSE & is.na(DE_Results$NCBIid)==FALSE,]
  
  #I also only want the subset of data that is annotated with a single gene (not ambiguously mapped to more than one gene)
  if(length(grep('\\|', DE_Results_NoNA$NCBIid))==0){
    DE_Results_GoodAnnotation<<-DE_Results_NoNA
  }else{
    #I only want rows annotated with a single Gene Symbol (no pipe):
    DE_Results_GoodAnnotation<<-DE_Results_NoNA[-(grep('\\|', DE_Results_NoNA$NCBIid)),]
  }
  #I used a double arrow in that conditional to place DE_Results_GoodAnnotation back out in the environment outside the function 
  
  print("# of rows with good annotation")
  print(nrow(DE_Results_GoodAnnotation))
  
  #For record keeping (sometimes useful for troubleshooting later)
  write.csv(DE_Results_GoodAnnotation, "DE_Results_GoodAnnotation.csv")
  
  rm(DE_Results_NoNA, DE_Results)
  
  print("Outputted object: DE_Results_GoodAnnotation")
}

#Cleaning up the results for a single result set

#Reading in the function

#You can input this function by running the code discussed above to create the function in your R environment. 

#Alternatively, you can download the script for the function from our Github site and save the file in your working directory:
#https://github.com/hagenaue/BrainDataAlchemy/blob/main/MetaAnalysis_GemmaDEResults_2024/Function_FilteringDEResults_GoodAnnotation.R

#And then source it from your working directory:
source("Function_FilteringDEResults_GoodAnnotation.R")

#Example of using the function for a dataset:

FilteringDEResults_GoodAnnotation(DE_Results)

# [1] "# of rows in results"
# [1] 21693
# [1] "# of rows with missing NCBI annotation:"
# [1] NA
# [1] "# of rows with NA NCBI annotation:"
# [1] 79
# [1] "# of rows with missing Gene Symbol annotation:"
# [1] 79
# [1] "# of rows mapped to multiple NCBI_IDs:"
# [1] 0
# [1] "# of rows mapped to multiple Gene Symbols:"
# [1] 0
# [1] "# of rows with good annotation"
# [1] 21614
# [1] "Outputted object: DE_Results_GoodAnnotation"

str(DE_Results_GoodAnnotation)
# Classes ‘data.table’ and 'data.frame':	21614 obs. of  19 variables:
#   $ Probe                      : int  20344 20558 110454 20339 17067 68774 21946 21366 76905 16190 ...
# $ NCBIid                     : int  20344 20558 110454 20339 17067 68774 21946 21366 76905 16190 ...
# $ GeneSymbol                 : chr  "Selp" "Slfn4" "Ly6a" "Sele" ...
# $ GeneName                   : chr  "selectin, platelet" "schlafen 4" "lymphocyte antigen 6 family member A" "selectin, endothelial cell" ...
# $ pvalue                     : num  1.70e-10 9.11e-09 5.38e-08 7.73e-08 1.00e-07 ...
# $ corrected_pvalue           : num  3.69e-06 9.88e-05 4.00e-04 4.00e-04 4.00e-04 ...
# $ rank                       : num  4.61e-05 9.22e-05 1.00e-04 2.00e-04 2.00e-04 ...
# $ contrast_151617_coefficient: num  3.102 3.27 0.763 2.868 0.513 ...
# $ contrast_151617_log2fc     : num  3.102 3.27 0.763 2.868 0.513 ...
# $ contrast_151617_tstat      : num  2.97 3.62 3.86 2.31 4.17 ...
# $ contrast_151617_pvalue     : num  0.0077 0.0017 0.001 0.0318 0.0005 ...
# $ contrast_151618_coefficient: num  7.01 5.21 1.35 5.86 0.82 ...
# $ contrast_151618_log2fc     : num  7.01 5.21 1.35 5.86 0.82 ...
# $ contrast_151618_tstat      : num  7.81 6.14 7.02 5.28 6.78 ...
# $ contrast_151618_pvalue     : num  1.86e-07 5.64e-06 9.01e-07 3.78e-05 1.47e-06 ...
# $ contrast_151619_coefficient: num  1.223 -0.303 -0.243 0.118 -0.214 ...
# $ contrast_151619_log2fc     : num  1.223 -0.303 -0.243 0.118 -0.214 ...
# $ contrast_151619_tstat      : num  1.0192 -0.2593 -1.1221 0.0785 -1.6535 ...
# $ contrast_151619_pvalue     : num  0.32 0.798 0.275 0.938 0.114 ...
# - attr(*, ".internal.selfref")=<externalptr> 
#   - attr(*, "call")= chr "https://gemma.msl.ubc.ca/rest/v2/resultSets/553805"
# - attr(*, "env")=<environment: 0x7f89998ffea8> 

#This code includes a function within a function.

#This is the inner function:
GetContrastIDsforResultSet<-function(NamesOfFoldChangeColumns){
  #I split apart the column names:
  ColumnNames_BrokenUp<-strsplit(NamesOfFoldChangeColumns, "_")
  #Put them in a matrix format
  MatrixOfColumnNames_BrokenUp<-matrix(unlist(ColumnNames_BrokenUp), ncol=3,byrow=T)
  #And then grab the contrast ids:
  ContrastIDs_inCurrentDF<-MatrixOfColumnNames_BrokenUp[,2]
  rm(ColumnNames_BrokenUp, MatrixOfColumnNames_BrokenUp)
  return(ContrastIDs_inCurrentDF)
}

#This is the outer function:

ExtractingDEResultsForContrasts<-function(DE_Results_GoodAnnotation, Contrasts_Log2FC, Contrasts_Tstat, ResultSet_contrasts){
  
  print("These are all of the columns in the differential expression results for our current result set:")
  
  print(colnames(DE_Results_GoodAnnotation))
  # [1] "Probe"                       "NCBIid"                     
  # [3] "GeneSymbol"                  "GeneName"                   
  # [5] "pvalue"                      "corrected_pvalue"           
  # [7] "rank"                        "contrast_151617_coefficient"
  # [9] "contrast_151617_log2fc"      "contrast_151617_tstat"      
  # [11] "contrast_151617_pvalue"      "contrast_151618_coefficient"
  # [13] "contrast_151618_log2fc"      "contrast_151618_tstat"      
  # [15] "contrast_151618_pvalue"      "contrast_151619_coefficient"
  # [17] "contrast_151619_log2fc"      "contrast_151619_tstat"      
  # [19] "contrast_151619_pvalue" 
  
  print("These are the names of the Log(2) Fold Change Columns for our statistical contrasts of interest within the differential expression results for this particular result set:")
  
  NamesOfFoldChangeColumns<<-colnames(DE_Results_GoodAnnotation)[colnames(DE_Results_GoodAnnotation)%in%Contrasts_Log2FC]
  
  print(NamesOfFoldChangeColumns)
  #[1] "contrast_151617_log2fc" "contrast_151618_log2fc" "contrast_151619_log2fc"
  
  print("These are the names of the T-statistic Columns for our statistical contrasts of interest within the differential expression results for this particular result set:")
  
  NamesOfTstatColumns<<-colnames(DE_Results_GoodAnnotation)[colnames(DE_Results_GoodAnnotation)%in%Contrasts_Tstat]
  
  print(NamesOfTstatColumns)
  #[1] "contrast_151617_tstat" "contrast_151618_tstat" "contrast_151619_tstat"
  
  #Next we're going to pull out the contrast IDs associated with each result:
  
  ContrastIDs_inCurrentDF<-GetContrastIDsforResultSet(NamesOfFoldChangeColumns)
  
  print("These are the contrast ids for the statistical contrasts of interest within your current result set:")
  print(ContrastIDs_inCurrentDF)
  #[1] "151617" "151618" "151619"
  
  #Next we're going to grab some metadata to go with those statistical contrasts
  
  print("This is the dataset id for the result set and statistical contrasts:")
  Datasets_inCurrentDF<-ResultSet_contrasts$ExperimentID[ResultSet_contrasts$ContrastIDs%in%ContrastIDs_inCurrentDF]
  
  GSE_ID<<-Datasets_inCurrentDF[1]
  
  print(GSE_ID)
  #[1] "GSE126678"
  
  #And I would like the experimental factor information for our statistical contrast:
  Factors_inCurrentDF<-ResultSet_contrasts$ExperimentalFactors[ResultSet_contrasts$ContrastIDs%in%ContrastIDs_inCurrentDF]
  
  #We can combine those to make an interpretable unique identifier for each statistical comparison:
  ComparisonsOfInterest<<-paste(Datasets_inCurrentDF, Factors_inCurrentDF, sep="_" )
  
  print("These are the current names for your statistical contrasts of interest - if they are unwieldy, you may want to change them")
  print(ComparisonsOfInterest)
  
  #cleaning up the workspace:
  rm(Datasets_inCurrentDF, Factors_inCurrentDF, ContrastIDs_inCurrentDF)
  
}

#Read in the function:

#You can input this function by running the code discussed above to create the function in your R environment. 

#Alternatively, you can download the script for the function from our Github site and save the file in your working directory:
#https://github.com/hagenaue/BrainDataAlchemy/blob/main/MetaAnalysis_GemmaDEResults_2024/Function_ExtractingDEResultsForContrasts.R

#And then source it from your working directory:
source("Function_ExtractingDEResultsForContrasts.R")

#Example function usage:
ExtractingDEResultsForContrasts(DE_Results_GoodAnnotation, Contrasts_Log2FC, Contrasts_Tstat, ResultSet_contrasts)
  
# [1] "These are all of the columns in the differential expression results for our current result set:"
# [1] "Probe"                       "NCBIid"                     
# [3] "GeneSymbol"                  "GeneName"                   
# [5] "pvalue"                      "corrected_pvalue"           
# [7] "rank"                        "contrast_151617_coefficient"
# [9] "contrast_151617_log2fc"      "contrast_151617_tstat"      
# [11] "contrast_151617_pvalue"      "contrast_151618_coefficient"
# [13] "contrast_151618_log2fc"      "contrast_151618_tstat"      
# [15] "contrast_151618_pvalue"      "contrast_151619_coefficient"
# [17] "contrast_151619_log2fc"      "contrast_151619_tstat"      
# [19] "contrast_151619_pvalue"     
# [1] "These are the names of the Log(2) Fold Change Columns for our statistical contrasts of interest within the differential expression results for this particular result set:"
# [1] "contrast_151617_log2fc" "contrast_151618_log2fc" "contrast_151619_log2fc"
# [1] "These are the names of the T-statistic Columns for our statistical contrasts of interest within the differential expression results for this particular result set:"
# [1] "contrast_151617_tstat" "contrast_151618_tstat" "contrast_151619_tstat"
# [1] "These are the contrast ids for the statistical contrasts of interest within your current result set:"
# [1] "151617" "151618" "151619"
# [1] "This is the dataset id for the result set and statistical contrasts:"
# [1] "GSE126678"
# [1] "These are the current names for your statistical contrasts of interest - if they are unwieldy, you may want to change them"
# [1] "GSE126678_lipopolysaccharide has modifier Acute immune challenge , vehicle; lipopolysaccharide has modifier Acute immune challenge , vehicle"                              
# [2] "GSE126678_lipopolysaccharide has modifier Long-term subchronic immune challenge + acute immune challenge"                                                                  
# [3] "GSE126678_lipopolysaccharide has modifier long-term subchronic immune challenge , vehicle; lipopolysaccharide has modifier long-term subchronic immune challenge , vehicle"


#Those names are super unwieldy. I'm going to rename them:
#Note: Make sure the name still includes the GEO dataset id number (GSE...)
ComparisonsOfInterest<-c("GSE126678_LPS_Acute", "GSE126678_LPS_SubchronicPlusAcute", "GSE126678_Subchronic")

CollapsingDEResults_OneResultPerGene<-function(GSE_ID, DE_Results_GoodAnnotation, ComparisonsOfInterest, NamesOfFoldChangeColumns, NamesOfTstatColumns){
  
  print("Double check that the vectors containing the fold change and tstat column names contain the same order as the group comparisons of interest - otherwise this function won't work properly!  If the order matches, proceed:")
  
  print("# of rows with unique NCBI IDs:")
  print(length(unique(DE_Results_GoodAnnotation$NCBIid)))
  
  print("# of rows with unique Gene Symbols:")
  print(length(unique(DE_Results_GoodAnnotation$GeneSymbol)))
  
  #This makes a folder named after the dataset ID to store your results
  dir.create(paste("./", GSE_ID, sep=""))
  
  #And then sets the working directory to be that folder:
  setwd(paste("./", GSE_ID, sep=""))
  
  #Calculating the average Log2FC and T-stat for each gene:
  
  #Creating an empty list to store our results from each of our statistical contrasts:
  DE_Results_GoodAnnotation_FoldChange_Average<-list()
  DE_Results_GoodAnnotation_Tstat_Average<-list()
  DE_Results_GoodAnnotation_SE_Average<-list()
  
  #For each of the columns containing fold change and t-stat information for our statistical contrasts of interest:
  
  for(i in c(1:length(NamesOfFoldChangeColumns))){
    
    #Grab our Log2FC column of interest:
    FoldChangeColumn<-select(DE_Results_GoodAnnotation, NamesOfFoldChangeColumns[i])
    
    #Grab our Tstat column of interest:
    TstatColumn<-select(DE_Results_GoodAnnotation, NamesOfTstatColumns[i])
    
    #Calculate the SE:
    DE_Results_GoodAnnotation_SE<-FoldChangeColumn[[1]]/TstatColumn[[1]]
    
    #Calculate the average Log2FC per gene:
    DE_Results_GoodAnnotation_FoldChange_Average[[i]]<-tapply(FoldChangeColumn[[1]], DE_Results_GoodAnnotation$NCBIid, mean)
    
    #Calculate the average Tstat per gene:
    DE_Results_GoodAnnotation_Tstat_Average[[i]]<-tapply(TstatColumn[[1]], DE_Results_GoodAnnotation$NCBIid, mean)
    
    #Calculate the average SE per gene:
    DE_Results_GoodAnnotation_SE_Average[[i]]<-tapply(DE_Results_GoodAnnotation_SE, DE_Results_GoodAnnotation$NCBIid, mean)
    
  }
  
  #Currently we have all of this averaged Log2FC information stored in a list, with one entry for each statistical contrast
  #Let's convert this to a data frame 
  DE_Results_GoodAnnotation_FoldChange_AveragedByGene<-do.call(cbind, DE_Results_GoodAnnotation_FoldChange_Average)
  
  print("Dimensions of Fold Change matrix, averaged by gene symbol:")
  print(dim(DE_Results_GoodAnnotation_FoldChange_AveragedByGene))
  
  #Name the columns in the dataframe in a manner that describes the dataset and factors for each statistic contrast
  colnames(DE_Results_GoodAnnotation_FoldChange_AveragedByGene)<-ComparisonsOfInterest
  
  #and then write it out to save it:
  write.csv(DE_Results_GoodAnnotation_FoldChange_AveragedByGene, "DE_Results_GoodAnnotation_FoldChange_AveragedByGene.csv")
  
  #And do the same for the T-stats:
  DE_Results_GoodAnnotation_Tstat_AveragedByGene<-do.call(cbind, DE_Results_GoodAnnotation_Tstat_Average)
  
  colnames(DE_Results_GoodAnnotation_Tstat_AveragedByGene)<-ComparisonsOfInterest
  
  write.csv(DE_Results_GoodAnnotation_Tstat_AveragedByGene, "DE_Results_GoodAnnotation_Tstat_AveragedByGene.csv")
  
  #And the same for the SE:
  DE_Results_GoodAnnotation_SE_AveragedByGene<-do.call(cbind, DE_Results_GoodAnnotation_SE_Average)
  
  colnames(DE_Results_GoodAnnotation_SE_AveragedByGene)<-ComparisonsOfInterest
  
  write.csv(DE_Results_GoodAnnotation_SE_AveragedByGene, "DE_Results_GoodAnnotation_SE_AveragedByGene.csv")
  
  #For running our meta-analysis, we are actually going to need the sampling variance instead of the standard error
  #The sampling variance is just the standard error squared.
  
  DE_Results_GoodAnnotation_SV<-(DE_Results_GoodAnnotation_SE_AveragedByGene)^2
  
  #Writing that information out to store it:
  write.csv(DE_Results_GoodAnnotation_SV, "DE_Results_GoodAnnotation_SV.csv")
  
  #Compiling all of our results into a single, big object 
  TempMasterResults<-list(Log2FC=DE_Results_GoodAnnotation_FoldChange_AveragedByGene, Tstat=DE_Results_GoodAnnotation_Tstat_AveragedByGene, SE=DE_Results_GoodAnnotation_SE_AveragedByGene, SV=DE_Results_GoodAnnotation_SV)
  
  #And then sending that object out into our environment outside the function
  #With a name that we can read and understand (DEResults for a particular GSEID)
  assign(paste("DEResults", GSE_ID, sep="_"), TempMasterResults, envir = as.environment(1))
  
  #Letting us know what that name is:
  print(paste("Output: Named DEResults", GSE_ID, sep="_"))
  
  #And then cleaning up our environment of temporary results:
  rm(TempMasterResults, DE_Results_GoodAnnotation, DE_Results_GoodAnnotation_SV, DE_Results_GoodAnnotation_SE, DE_Results_GoodAnnotation_FoldChange_AveragedByGene, DE_Results_GoodAnnotation_FoldChange_Average, DE_Results_GoodAnnotation_Tstat_AveragedByGene, DE_Results_GoodAnnotation_Tstat_Average, DE_Results_GoodAnnotation_SE_Average, FoldChangeColumn, TstatColumn, GSE_ID, ComparisonsOfInterest, NamesOfFoldChangeColumns, NamesOfTstatColumns)
  
  #This sets the working directory back to the parent directory
  setwd("../")
  
}

#Reading in the function:

#You can input this function by running the code discussed above to create the function in your R environment. 

#Alternatively, you can download the script for the function from our Github site and save the file in your working directory:
#https://github.com/hagenaue/BrainDataAlchemy/blob/main/MetaAnalysis_GemmaDEResults_2024/Function_CollapsingDEResults_OneResultPerGene.R

#And then source it from your working directory:
source("Function_CollapsingDEResults_OneResultPerGene.R")

#Notes about parameters for function CollapsingDEResults_OneResultPerGene()
#GSE_ID is a string indicating the name of the Gemma dataset
#DE_Results_GoodAnnotation is the data frame outputted by our previous function
#ComparisonsOfInterest is a character vector containing the names of the group comparisons of interest within this dataset. Important: These group comparisons should be listed in exactly the same order as the order that you provide the column names for their associated Fold Change and Tstat output.
#NamesOfFoldChangeColumns is a vector containing the names of the columns of DE_Results_GoodAnnotation containing the FoldChange results for your comparisons of interes, in the same order as the ComparisonsOfInterest vector
#NamesOfTstatColumns is a vector containing the names of the columns of DE_Results_GoodAnnotation containing the Tstat results for your comparisons of interes, in the same order as the ComparisonsOfInterest vector

#Example usage:
CollapsingDEResults_OneResultPerGene(GSE_ID, DE_Results_GoodAnnotation, ComparisonsOfInterest, NamesOfFoldChangeColumns, NamesOfTstatColumns)

# [1] "Double check that the vectors containing the fold change and tstat column names contain the same order as the group comparisons of interest - otherwise this function won't work properly!  If the order matches, proceed:"
# [1] "# of rows with unique NCBI IDs:"
# [1] 21614
# [1] "# of rows with unique Gene Symbols:"
# [1] 21614
# [1] "Dimensions of Fold Change matrix, averaged by gene symbol:"
# [1] 21614     3
# [1] "Output: Named DEResults_GSE126678"

#A function for aligning all of our mouse differential expression results from different datasets into a single data frame for Log2FCs and sampling variances (SVs):

#This function works the same way as the rat alignment function:

AligningMouseDatasets<-function(ListOfMouseDEResults){
  
  Mouse_MetaAnalysis_FoldChange_Dfs<-list()
  
  for(i in c(1:length(ListOfMouseDEResults))){
    Mouse_MetaAnalysis_FoldChange_Dfs[[i]]<-data.frame(Mouse_EntrezGene.ID=row.names(ListOfMouseDEResults[[i]][[1]]),ListOfMouseDEResults[[i]][[1]], stringsAsFactors=FALSE)
  }
  
  print("Mouse_MetaAnalysis_FoldChange_Dfs:")
  print(str(Mouse_MetaAnalysis_FoldChange_Dfs))
  
  Mouse_MetaAnalysis_FoldChanges<<-join_all(Mouse_MetaAnalysis_FoldChange_Dfs, by="Mouse_EntrezGene.ID", type="full")
  #This function could be join_all (if there are more than 2 datasets) or merge/merge_all (if the plyr package isn't working)
  
  print("Mouse_MetaAnalysis_FoldChanges:")
  print(str(Mouse_MetaAnalysis_FoldChanges))
  
  Mouse_MetaAnalysis_SV_Dfs<-list()
  
  for(i in c(1:length(ListOfMouseDEResults))){
    Mouse_MetaAnalysis_SV_Dfs[[i]]<-data.frame(Mouse_EntrezGene.ID=row.names(ListOfMouseDEResults[[i]][[4]]),ListOfMouseDEResults[[i]][[4]], stringsAsFactors=FALSE)
  }
  
  print("Mouse_MetaAnalysis_SV_Dfs:")
  print(str(Mouse_MetaAnalysis_SV_Dfs))
  
  Mouse_MetaAnalysis_SV<<-join_all(Mouse_MetaAnalysis_SV_Dfs, by="Mouse_EntrezGene.ID", type="full")
  #This function could be join_all (if there are more than 2 datasets) or merge/merge_all (if the plyr package isn't working)
  
  print("Mouse_MetaAnalysis_SV:")
  print(str(Mouse_MetaAnalysis_SV))
  
  rm(Mouse_MetaAnalysis_SV_Dfs, Mouse_MetaAnalysis_FoldChange_Dfs)
}

#Reading in the function:

#You can input this function by running the code discussed above to create the function in your R environment. 

#Alternatively, you can download the script for the function from our Github site and save the file in your working directory:
#https://github.com/hagenaue/BrainDataAlchemy/blob/main/MetaAnalysis_GemmaDEResults_2024/Function_AligningDEResults.R

#And then source it from your working directory:
source("Function_AligningDEResults.R")

#Aligning the mouse datasets with each other:

#Example Usage:
#In this example, there are differential expression results from two mouse datasets (GSE126678 and GSE181285):

ListOfMouseDEResults<-list(DEResults_GSE126678, DEResults_GSE181285)


AligningMouseDatasets(ListOfMouseDEResults)

# [1] "Mouse_MetaAnalysis_FoldChange_Dfs:"
# List of 2
# $ :'data.frame':	21614 obs. of  4 variables:
#   ..$ Mouse_EntrezGene.ID              : chr [1:21614] "11287" "11298" "11302" "11303" ...
# ..$ GSE126678_LPS_Acute              : num [1:21614] 1.9397 0.0805 0.0595 0.0306 0.276 ...
# ..$ GSE126678_LPS_SubchronicPlusAcute: num [1:21614] 1.2967 -0.0472 -0.1459 0.1367 1.5651 ...
# ..$ GSE126678_LPS_Subchronic         : num [1:21614] 0.0582 0.203 -0.1144 0.1361 -0.0051 ...
# $ :'data.frame':	18563 obs. of  2 variables:
#   ..$ Mouse_EntrezGene.ID: chr [1:18563] "100008567" "100009600" "100012" "100017" ...
# ..$ GSE181285_LPS_Acute: num [1:18563] 0.0198 0.0225 0.0641 -0.0049 -0.0588 ...
# NULL
# [1] "Mouse_MetaAnalysis_FoldChanges:"
# 'data.frame':	24287 obs. of  5 variables:
#   $ Mouse_EntrezGene.ID              : chr  "11287" "11298" "11302" "11303" ...
# $ GSE126678_LPS_Acute              : num  1.9397 0.0805 0.0595 0.0306 0.276 ...
# $ GSE126678_LPS_SubchronicPlusAcute: num  1.2967 -0.0472 -0.1459 0.1367 1.5651 ...
# $ GSE126678_LPS_Subchronic         : num  0.0582 0.203 -0.1144 0.1361 -0.0051 ...
# $ GSE181285_LPS_Acute              : num  -0.042 -0.0368 -0.0534 0.1067 -0.5258 ...
# NULL
# [1] "Mouse_MetaAnalysis_SV_Dfs:"
# List of 2
# $ :'data.frame':	21614 obs. of  4 variables:
#   ..$ Mouse_EntrezGene.ID              : chr [1:21614] "11287" "11298" "11302" "11303" ...
# ..$ GSE126678_LPS_Acute              : num [1:21614] 0.62127 0.14737 0.00437 0.01624 1.34369 ...
# ..$ GSE126678_LPS_SubchronicPlusAcute: num [1:21614] 0.66434 0.14559 0.00438 0.01567 0.98359 ...
# ..$ GSE126678_LPS_Subchronic         : num [1:21614] 0.84004 0.14211 0.00439 0.01592 1.40671 ...
# $ :'data.frame':	18563 obs. of  2 variables:
#   ..$ Mouse_EntrezGene.ID: chr [1:18563] "100008567" "100009600" "100012" "100017" ...
# ..$ GSE181285_LPS_Acute: num [1:18563] 0.11456 0.01612 0.00329 0.00487 0.00719 ...
# NULL
# [1] "Mouse_MetaAnalysis_SV:"
# 'data.frame':	24287 obs. of  5 variables:
#   $ Mouse_EntrezGene.ID              : chr  "11287" "11298" "11302" "11303" ...
# $ GSE126678_LPS_Acute              : num  0.62127 0.14737 0.00437 0.01624 1.34369 ...
# $ GSE126678_LPS_SubchronicPlusAcute: num  0.66434 0.14559 0.00438 0.01567 0.98359 ...
# $ GSE126678_LPS_Subchronic         : num  0.84004 0.14211 0.00439 0.01592 1.40671 ...
# $ GSE181285_LPS_Acute              : num  0.00391 0.02738 0.00601 0.0101 0.03332 ...
# NULL

#A function for aligning all of our rat differential expression results from different datasets into a single data frame for Log2FCs and sampling variances (SVs):

AligningRatDatasets<-function(ListOfRatDEResults){
  
  #Making an empty list to hold our results:
  Rat_MetaAnalysis_FoldChange_Dfs<-list()
  
  #Looping over all of the rat differential expression results:
  for(i in c(1:length(ListOfRatDEResults))){
    
    #Placing each of the log2FC results for each dataset into a single list
    #Each element in the list is formatted so that the rownames are Rat Entrez Gene ID and then there are columns containing the Log2FC for the differential expression results:
    Rat_MetaAnalysis_FoldChange_Dfs[[i]]<-data.frame(Rat_EntrezGene.ID=row.names(ListOfRatDEResults[[i]][[1]]),ListOfRatDEResults[[i]][[1]], stringsAsFactors=FALSE)
  }
  
  #Letting the user know the structure of the set of Log2FC dataframes that we are starting out with:
  print("Rat_MetaAnalysis_FoldChange_Dfs:")
  print(str(Rat_MetaAnalysis_FoldChange_Dfs))
  
  #Running "join all" on the list to align all of the results by Entrez Gene ID and make them a single data frame:
  Rat_MetaAnalysis_FoldChanges<<-join_all(Rat_MetaAnalysis_FoldChange_Dfs, by="Rat_EntrezGene.ID", type="full")
  #This function could be join_all (if there are more than 2 datasets) or merge/merge_all (if the plyr package isn't working)
  
  #Letting the user know the structure of the dataframe that we just created:
  print("Rat_MetaAnalysis_FoldChanges:")
  print(str(Rat_MetaAnalysis_FoldChanges))
  
  #Doing the same steps for the sampling variances:
  
  Rat_MetaAnalysis_SV_Dfs<-list()
  
  for(i in c(1:length(ListOfRatDEResults))){
    Rat_MetaAnalysis_SV_Dfs[[i]]<-data.frame(Rat_EntrezGene.ID=row.names(ListOfRatDEResults[[i]][[4]]),ListOfRatDEResults[[i]][[4]], stringsAsFactors=FALSE)
  }
  
  print("Rat_MetaAnalysis_SV_Dfs:")
  print(str(Rat_MetaAnalysis_SV_Dfs))
  
  Rat_MetaAnalysis_SV<<-join_all(Rat_MetaAnalysis_SV_Dfs, by="Rat_EntrezGene.ID", type="full")
  #This function could be join_all (if there are more than 2 datasets) or merge/merge_all (if the plyr package isn't working)
  
  print("Rat_MetaAnalysis_SV:")
  print(str(Rat_MetaAnalysis_SV))
  
  #Cleaning up our environment to remove unneeded objects:
  rm(Rat_MetaAnalysis_SV_Dfs, Rat_MetaAnalysis_FoldChange_Dfs)
}

#Reading in the function:

#You can input this function by running the code discussed above to create the function in your R environment. 

#Alternatively, you can download the script for the function from our Github site and save the file in your working directory:
#https://github.com/hagenaue/BrainDataAlchemy/blob/main/MetaAnalysis_GemmaDEResults_2024/Function_AligningDEResults.R

#And then source it from your working directory:
source("Function_AligningDEResults.R")

#Example Usage:

#This is an example in which there are only differential expression results from one rat dataset (GSE205325):

ListOfRatDEResults<-list(DEResults_GSE205325)

AligningRatDatasets(ListOfRatDEResults)
# [1] "Rat_MetaAnalysis_FoldChange_Dfs:"
# List of 1
# $ :'data.frame':	17196 obs. of  2 variables:
#   ..$ Rat_EntrezGene.ID    : chr [1:17196] "24153" "24157" "24158" "24159" ...
# ..$ GSE205325_LPS_Chronic: num [1:17196] 0.3485 0.0288 -0.1887 -0.2126 0.1235 ...
# NULL
# [1] "Rat_MetaAnalysis_FoldChanges:"
# 'data.frame':	17196 obs. of  2 variables:
#   $ Rat_EntrezGene.ID    : chr  "24153" "24157" "24158" "24159" ...
# $ GSE205325_LPS_Chronic: num  0.3485 0.0288 -0.1887 -0.2126 0.1235 ...
# NULL
# [1] "Rat_MetaAnalysis_SV_Dfs:"
# List of 1
# $ :'data.frame':	17196 obs. of  2 variables:
#   ..$ Rat_EntrezGene.ID    : chr [1:17196] "24153" "24157" "24158" "24159" ...
# ..$ GSE205325_LPS_Chronic: num [1:17196] 0.0745 0.0229 0.0383 0.0257 0.0135 ...
# NULL
# [1] "Rat_MetaAnalysis_SV:"
# 'data.frame':	17196 obs. of  2 variables:
#   $ Rat_EntrezGene.ID    : chr  "24153" "24157" "24158" "24159" ...
# $ GSE205325_LPS_Chronic: num  0.0745 0.0229 0.0383 0.0257 0.0135 ...
# NULL

#We have a ortholog database that we downloaded from Jackson Labs on April 25, 2024
#The reformatted version can be downloaded from our github site:
#https://github.com/hagenaue/BrainDataAlchemy/blob/main/MetaAnalysis_GemmaDEResults_2024/MouseVsRat_NCBI_Entrez_JacksonLab_20240425.csv 
#Make sure it is saved in your working directory 


MouseVsRat_NCBI_Entrez<-read.csv("MouseVsRat_NCBI_Entrez_JacksonLab_20240425.csv", header=TRUE, stringsAsFactors = FALSE, row.names=1, colClasses=c("character", "character", "character"))

#We want to join this ortholog database to our mouse results (Log2FC and SV):

Mouse_MetaAnalysis_FoldChanges_wOrthologs<-join(MouseVsRat_NCBI_Entrez, Mouse_MetaAnalysis_FoldChanges, by="Mouse_EntrezGene.ID", type="full")

#Structure of the new object:
str(Mouse_MetaAnalysis_FoldChanges_wOrthologs)
#'data.frame':	25288 obs. of  6 variables:

Mouse_MetaAnalysis_SV_wOrthologs<-join(MouseVsRat_NCBI_Entrez, Mouse_MetaAnalysis_SV, by="Mouse_EntrezGene.ID", type="full")

#Structure of the new object:
str(Mouse_MetaAnalysis_SV_wOrthologs)
#'data.frame':	25288 obs. of  6 variables:

#*If there are rat datasets*, we then want to join our mouse Log2FC and SV results to the rat results using the ortholog information:

MetaAnalysis_FoldChanges<-join(Mouse_MetaAnalysis_FoldChanges_wOrthologs, Rat_MetaAnalysis_FoldChanges, by="Rat_EntrezGene.ID", type="full")

#Structure of the new object:
str(MetaAnalysis_FoldChanges)
#'data.frame':	28101 obs. of  7 variables:

MetaAnalysis_SV<-join(Mouse_MetaAnalysis_SV_wOrthologs, Rat_MetaAnalysis_SV, by="Rat_EntrezGene.ID", type="full")

#Structure of the new object:
str(MetaAnalysis_SV)
#'data.frame':	28101 obs. of  7 variables:

#*If there aren't any rat datasets*, we just rename the dataframes so that our downstream code works:
MetaAnalysis_FoldChanges<-Mouse_MetaAnalysis_FoldChanges_wOrthologs
str(MetaAnalysis_FoldChanges)

MetaAnalysis_SV<-Mouse_MetaAnalysis_SV_wOrthologs
str(MetaAnalysis_SV)

#For simplicity's sake, I'm going to replace that Mouse-Rat Entrez annotation
#Because it is missing entries for any genes in the datasets that *don't* have orthologs
MetaAnalysis_FoldChanges$MouseVsRat_EntrezGene.ID<-paste(MetaAnalysis_FoldChanges$Mouse_EntrezGene.ID, MetaAnalysis_FoldChanges$Rat_EntrezGene.ID, sep="_")

MetaAnalysis_SV$MouseVsRat_EntrezGene.ID<-paste(MetaAnalysis_SV$Mouse_EntrezGene.ID, MetaAnalysis_SV$Rat_EntrezGene.ID, sep="_")

#Example scatter plot comparing the effects (Log2FCs) for all of the genes represented in two datasets:
plot(MetaAnalysis_FoldChanges$GSE126678_LPS_Acute~MetaAnalysis_FoldChanges$GSE181285_LPS_Acute)

#Note - many people prefer to plot these relationships using RRHOs (Rank rank hypergeometric overlap plots)
#I like using both.
#The code for the RRHOs is a little complicated, but I'm happy to share if folks are interested.

#Here's code for looking at the correlation of all of our log2FC results with all of our other log2FC results
#This is called a correlation matrix
#This one uses non-parametric (Spearman) correlation coefficients

write.csv(cor(as.matrix(MetaAnalysis_FoldChanges[,-c(1:3)]), use="pairwise.complete.obs", method="spearman"), "CorMatrix_AllLog2FCs.csv")

#An illustration of the correlation matrix using a hierarchically clustered heatmap:
heatmap(cor(as.matrix(MetaAnalysis_FoldChanges[,-c(1:3)]), use="pairwise.complete.obs", method="spearman"))

RunBasicMetaAnalysis<-function(NumberOfComparisons, CutOffForNAs, MetaAnalysis_FoldChanges, MetaAnalysis_SV){
  
  #The function first provides information about how many of the statistical contrasts have NA values as their differential expression results for each gene:
  MetaAnalysis_FoldChanges_NAsPerRow<-apply(MetaAnalysis_FoldChanges[,-c(1:3)], 1, function(y) sum(is.na(y)))
  
  print("Table of # of NAs per Row (Gene):")
  print(table(MetaAnalysis_FoldChanges_NAsPerRow))
  
  #Then any row (gene) that has too many NAs is removed from the analysis:
  MetaAnalysis_FoldChanges_ForMeta<<-MetaAnalysis_FoldChanges[MetaAnalysis_FoldChanges_NAsPerRow<CutOffForNAs,]
  MetaAnalysis_SV_ForMeta<<-MetaAnalysis_SV[MetaAnalysis_FoldChanges_NAsPerRow<CutOffForNAs,]
  
  print("MetaAnalysis_FoldChanges_ForMeta:")
  print(str(MetaAnalysis_FoldChanges_ForMeta))
  
  #I'm going to make an empty matrix to store the results of my meta-analysis:
  metaOutput<-matrix(NA, nrow(MetaAnalysis_FoldChanges_ForMeta), 6)
  
  #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):
  for(i in c(1:nrow(MetaAnalysis_FoldChanges_ForMeta))){
    
    #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:
    effect<-as.numeric(MetaAnalysis_FoldChanges_ForMeta[i,-c(1:3)])
    var<-as.numeric(MetaAnalysis_SV_ForMeta[i,-c(1:3)])
    
    #I added a function tryCatch that double-checks that the meta-analysis function (rma) doesn't produce errors (which breaks the loop):
    skip_to_next <- FALSE
    tryCatch(TempMeta<-rma(effect, var), error = function(e) {skip_to_next <<- TRUE})
    
    #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:
    if(skip_to_next){}else{
      TempMeta<-rma(effect, var)
      metaOutput[i, 1]<-TempMeta$b #gives estimate Log2FC
      metaOutput[i, 2]<-TempMeta$se #gives standard error
      metaOutput[i, 3]<-TempMeta$pval #gives pval
      metaOutput[i, 4]<-TempMeta$ci.lb #gives confidence interval lower bound
      metaOutput[i, 5]<-TempMeta$ci.ub #gives confidence interval upper bound
      metaOutput[i, 6]<-NumberOfComparisons-sum(is.na(effect))#Number of comparisons with data
      rm(TempMeta)
    }
    rm(effect, var)
  }
  
  #Naming the columns in our output:
  colnames(metaOutput)<-c("Log2FC_estimate", "SE", "pval", "CI_lb", "CI_ub", "Number_Of_Comparisons")
  
  #The row names for our output are the combined mouse-rat entrez ids: 
  row.names(metaOutput)<-MetaAnalysis_FoldChanges_ForMeta[,3]
  
  #We return this output back into our global environment
  metaOutput<<-metaOutput
  MetaAnalysis_Annotation<<-MetaAnalysis_FoldChanges_ForMeta[,c(1:3)]
  return(metaOutput)
  return(MetaAnalysis_Annotation)
  
  #... and provide the user with an update about the newly created object:
  
  print("metaOutput:")
  print(str(metaOutput))
  
  print("Top of metaOutput:")
  print(head(metaOutput))
  
  print("Bottom of metaOutput")
  print(tail(metaOutput))
  
}

#We can only run a meta-analysis if there are differential expression results from more than one comparison.
#Since I know that the differential expression results from the same study (dataset) are artificially correlated, I would prefer that there are results from more than one dataset.

#How many genes satisfy this criteria?

#This code caculates the number of NAs in each row:
MetaAnalysis_FoldChanges_NAsPerRow<-apply(MetaAnalysis_FoldChanges[,-c(1:3)], 1, function(y) sum(is.na(y)))


#I'm going to make a histogram of the results because I'm curious to see how they are distributed
hist(MetaAnalysis_FoldChanges_NAsPerRow)

#Or, since there are a limited number of integer answers (0-3), I could make a table of the results:
table(MetaAnalysis_FoldChanges_NAsPerRow)
# 0     1     2     3     4     5 
# 13355  3200  5059   277  5293   917

#For this dataset, I'm going to try running a meta-analysis using genes that were found in at least 4 sets of differential expression results
#Since there are 5 sets of differential expression results, that means that the genes that we are including need to have 1 or fewer NAs in their results
#I set this conservatively, because there are so few studies in this meta-analysis.
#2 NA is too many


#Example Usage:

NumberOfComparisons=5
CutOffForNAs=2
#I have 5 statistical contrasts total (comparisons)
#and 2 NA is too many


#Running a basic meta-analysis:

#This function  is designed to run a basic meta-analysis of Log2FC and sampling variance values using our previously generated objects MetaAnalysis_FoldChanges & MetaAnalysis_SV

#Reading in the function:

#You can input this function by running the code discussed above to create the function in your R environment. 

#Alternatively, you can download the script for the function from our Github site and save the file in your working directory:
#https://github.com/hagenaue/BrainDataAlchemy/blob/main/MetaAnalysis_GemmaDEResults_2024/Function_RunBasicMetaAnalysis.R

#And then source it from your working directory:
source("Function_RunBasicMetaAnalysis.R")                         

#Example usage:

metaOutput<-RunBasicMetaAnalysis(NumberOfComparisons, CutOffForNAs, MetaAnalysis_FoldChanges, MetaAnalysis_SV)
#Note: this function can take a while to run, especially if you have a lot of data  
#Plug in your computer, take a break, grab some coffee...


#Example output:

# [1] "Table of # of NAs per Row (Gene):"
# MetaAnalysis_FoldChanges_NAsPerRow
# 0     1     2     3     4     5 
# 13355  3200  5059   277  5293   917 
# [1] "MetaAnalysis_FoldChanges_ForMeta:"
# 'data.frame':	16555 obs. of  8 variables:
#   $ Rat_EntrezGene.ID                 : chr  "114087" "191569" "246307" "65041" ...
# $ Mouse_EntrezGene.ID               : chr  "23825" "18585" "66514" "20480" ...
# $ MouseVsRat_EntrezGene.ID          : chr  "23825_114087" "18585_191569" "66514_246307" "20480_65041" ...
# $ LPS_SubchronicAndAcute_vs_Vehicle : num  -0.0239 -0.0858 -0.0686 0.0891 0.0376 ...
# $ LPS_Acute_vs_Vehicle              : num  -0.0938 0.2524 -0.109 0.0788 -0.0475 ...
# $ LPS_Subchronic_vs_Vehicle         : num  0.0355 0.2735 0.0415 -0.0116 0.0299 ...
# $ LPS_Acute850ugPerKg_vs_Vehicle    : num  0.1754 0.1651 0.1878 -0.0146 0.103 ...
# $ LPS_Chronic_vs_Vehicle_AllStressed: num  0.223 -0.3824 -0.197 -0.0051 0.0393 ...
# NULL

str(metaOutput)
# num [1:16555, 1:6] 0.0234 0.1315 -0.023 0.0437 0.0769 ...
# - attr(*, "dimnames")=List of 2
# ..$ : chr [1:16555] "23825_114087" "18585_191569" "66514_246307" "20480_65041" ...
# ..$ : chr [1:6] "Log2FC_estimate" "SE" "pval" "CI_lb" ...

head(metaOutput)
# Log2FC_estimate         SE       pval        CI_lb      CI_ub Number_Of_Comparisons
# 23825_114087      0.02335239 0.04876404 0.63202016 -0.072223377 0.11892815                     5
# 18585_191569      0.13152931 0.06603932 0.04640599  0.002094618 0.26096399                     5
# 66514_246307     -0.02304066 0.05315084 0.66465470 -0.127214402 0.08113308                     5
# 20480_65041       0.04374446 0.03359984 0.19294219 -0.022110021 0.10959893                     5
# 13726_25437       0.07691445 0.04407249 0.08095348 -0.009466050 0.16329495                     5
# 16952_25380       0.04465958 0.10253423 0.66315765 -0.156303828 0.24562298                     5                   

tail(metaOutput)
# Log2FC_estimate         SE       pval       CI_lb      CI_ub Number_Of_Comparisons
# 108168734_NA      0.11566278 0.07742137 0.13519162 -0.03608031 0.26740587                     4
# 108168883_NA      0.16641238 0.09078313 0.06679127 -0.01151929 0.34434405                     4
# 108168923_NA      0.07757892 0.16987291 0.64789533 -0.25536587 0.41052370                     4
# 108168987_NA     -0.13375235 0.07407534 0.07097679 -0.27893734 0.01143264                     4
# 108169023_NA     -0.06106329 0.04794113 0.20276476 -0.15502618 0.03289960                     4
# 113002583_NA      0.02982367 0.07843726 0.70377969 -0.12391053 0.18355788                     4

#Installing and loading code package multtest

if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("multtest")                                 
library(multtest)

#We can add some additional gene annotation at this point too:

#Reading in a database containing more detailed gene annotation:
HOM_MouseVsRat <- read.csv("HOM_MouseVsRat_20240425.csv", header = TRUE, row.names = 1)

colnames(HOM_MouseVsRat)
# [1] "DB.Class.Key"                                          
# [2] "Mouse_Common.Organism.Name"                            
# [3] "Mouse_NCBI.Taxon.ID"                                   
# [4] "Mouse_Symbol"                                          
# [5] "Mouse_EntrezGene.ID"                                   
# [6] "Mouse_Mouse.MGI.ID"                                    
# [7] "Mouse_HGNC.ID"                                         
# [8] "Mouse_OMIM.Gene.ID"                                    
# [9] "Mouse_Genetic.Location"                                
# [10] "Mouse_Genome.Coordinates..mouse..GRCm39.human..GRCh38."
# [11] "Mouse_Name"                                            
# [12] "Mouse_Synonyms"                                        
# [13] "Rat_Common.Organism.Name"                              
# [14] "Rat_NCBI.Taxon.ID"                                     
# [15] "Rat_Symbol"                                            
# [16] "Rat_EntrezGene.ID"                                     
# [17] "Rat_Mouse.MGI.ID"                                      
# [18] "Rat_HGNC.ID"                                           
# [19] "Rat_OMIM.Gene.ID"                                      
# [20] "Rat_Genetic.Location"                                  
# [21] "Rat_Genome.Coordinates..mouse..GRCm39.human..GRCh38."  
# [22] "Rat_Name"                                              
# [23] "Rat_Synonyms" 

#Renaming the columns so that we can easily join the annotation to our meta-analysis results:
HOM_MouseVsRat$Mouse_EntrezGene.ID <- as.character(HOM_MouseVsRat$Mouse_EntrezGene.ID)

HOM_MouseVsRat$Rat_EntrezGene.ID <- as.character(HOM_MouseVsRat$Rat_EntrezGene.ID)
 

FalseDiscoveryCorrection<-function(metaOutput, HOM_MouseVsRat, MetaAnalysis_Annotation){
  
  #This calculates the false discovery rate, or q-value, for each of our p-values using the Benjamini-Hochberg procedure:
  tempPvalAdjMeta<-mt.rawp2adjp(metaOutput[,3], proc=c("BH"))
  
  #Then we put those results back into the order of our orginal output:
  metaPvalAdj<-tempPvalAdjMeta$adjp[order(tempPvalAdjMeta$index),]
  
  #And bind the false discovery rate (FDR) to the rest of the meta-analysis output:
  metaOutputFDR<-cbind(metaOutput, metaPvalAdj[,2])
  
  #And name that column FDR:
  colnames(metaOutputFDR)[7]<-"FDR"
  
  #These results are returned to our global environment:
  metaOutputFDR<<-metaOutputFDR
  
  #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)
  print("metaOutputFDR:")
  print(str(metaOutputFDR))
  
  #Then we make a dataframe that adds the annotation to that output:
  TempDF<-cbind.data.frame(metaOutputFDR, MetaAnalysis_Annotation)
  #And then adds even more detailed gene annotation:
  
  #First the detailed annotation for the mouse genes:
  TempDF2<-join(TempDF, HOM_MouseVsRat[,c(4:5,9:11)], by="Mouse_EntrezGene.ID", type="left", match="first")
  
  #Next the annnotation for the rat genes:
  TempDF3<-join(TempDF2, HOM_MouseVsRat[,c(15:16,20:22)], by="Rat_EntrezGene.ID", type="left", match="first")
  
  #This is renamed and returned to our global environment:
  metaOutputFDR_annotated<-TempDF3
  metaOutputFDR_annotated<<-metaOutputFDR_annotated
  
  #And written out into our working directory:
  write.csv(metaOutputFDR_annotated, "metaOutputFDR_annotated.csv")
  
  #Then we make a version of the output in order by p-value:
  metaOutputFDR_OrderbyPval<<-metaOutputFDR_annotated[order(metaOutputFDR_annotated[,5]),]
  
  #Let's write out a version of the output in order by p-value:
  write.csv(metaOutputFDR_OrderbyPval, "metaOutputFDR_orderedByPval.csv")
  
  #And give the user some information about their results:
  
  print("Do we have any genes that are statistically significant following loose false discovery rate correction (FDR<0.10)?")
  print(sum(metaOutputFDR_annotated[,9]<0.10, na.rm=TRUE))
  
  print("Do we have any genes that are statistically significant following traditional false discovery rate correction (FDR<0.05)?")
  print(sum(metaOutputFDR_annotated[,9]<0.05, na.rm=TRUE))
  
  print("What are the top results?")
  print(head(metaOutputFDR_annotated[order(metaOutputFDR_annotated[,5]),]))
  
  rm(tempPvalAdjMeta, metaPvalAdj)
  
} 

#This code runs a function that corrects the meta-analysis output to take into account the fact that we are running the statistical calculations thousands of times and therefore have a heightened risk of false discovery (false discovery rate correction) 

#Reading in the function:

#You can input this function by running the code discussed above to create the function in your R environment. 

#Alternatively, you can download the script for the function from our Github site and save the file in your working directory:
#https://github.com/hagenaue/BrainDataAlchemy/blob/main/MetaAnalysis_GemmaDEResults_2024/Function_FalseDiscoveryCorrection.R

#And then source it from your working directory:
source("Function_FalseDiscoveryCorrection.R")

#Example usage:

FalseDiscoveryCorrection(metaOutput, HOM_MouseVsRat, MetaAnalysis_Annotation)

#Example output:

# [1] "metaOutputFDR:"
# num [1:16555, 1:7] 0.0234 0.1315 -0.023 0.0437 0.0769 ...
# - attr(*, "dimnames")=List of 2
# ..$ : chr [1:16555] "23825_114087" "18585_191569" "66514_246307" "20480_65041" ...
# ..$ : chr [1:7] "Log2FC_estimate" "SE" "pval" "CI_lb" ...
# NULL
# [1] "Do we have any genes that are statistically significant following loose false discovery rate correction (FDR<0.10)?"
# [1] 253
# [1] "Do we have any genes that are statistically significant following traditional false discovery rate correction (FDR<0.05)?"
# [1] 136
# [1] "What are the top results?"
# Rat_EntrezGene.ID Mouse_EntrezGene.ID Log2FC_estimate         SE         pval      CI_lb       CI_ub
# 13170_24309               24309               13170      -0.4293013 0.07014204 9.330764e-10 -0.5667771 -0.29182542
# 74772_362645             362645               74772      -0.1313814 0.02330647 1.729185e-08 -0.1770612 -0.08570151
# 192188_282580            282580              192188       0.3562365 0.06392332 2.505735e-08  0.2309491  0.48152389
# 54006_83632               83632               54006      -0.1635032 0.02996583 4.861038e-08 -0.2222351 -0.10477121
# 216565_305556            305556              216565       0.1764604 0.03363489 1.551430e-07  0.1105372  0.24238352
# 268445_360575            360575              268445      -0.1584491 0.03081248 2.712959e-07 -0.2188405 -0.09805779
# Number_Of_Comparisons          FDR MouseVsRat_EntrezGene.ID Mouse_Symbol Mouse_Genetic.Location
# 13170_24309                       5 1.544708e-05              13170_24309          Dbp               Chr7  cM
# 74772_362645                      5 1.382748e-04             74772_362645      Atp13a2               Chr4  cM
# 192188_282580                     5 1.382748e-04            192188_282580        Stab2              Chr10  cM
# 54006_83632                       5 2.011862e-04              54006_83632        Deaf1               Chr7  cM
# 216565_305556                     5 5.136784e-04            216565_305556        Ehbp1              Chr11  cM
# 268445_360575                     5 7.485505e-04            268445_360575     Ankrd13b              Chr11  cM
# Mouse_Genome.Coordinates..mouse..GRCm39.human..GRCh38.                              Mouse_Name Rat_Symbol
# 13170_24309                                Chr7:45354658-45359579(+) D site albumin promoter binding protein        Dbp
# 74772_362645                             Chr4:140714184-140734641(+)                        ATPase type 13A2    Atp13a2
# 192188_282580                             Chr10:86677062-86843889(-)                              stabilin 2      Stab2
# 54006_83632                              Chr7:140877093-140907603(-)             DEAF1, transcription factor      Deaf1
# 216565_305556                             Chr11:21955825-22237086(-)             EH domain binding protein 1      Ehbp1
# 268445_360575                             Chr11:77361311-77380504(-)               ankyrin repeat domain 13b   Ankrd13b
# Rat_Genetic.Location Rat_Genome.Coordinates..mouse..GRCm39.human..GRCh38.                                    Rat_Name
# 13170_24309               Chr1 q22                                                   NA D-box binding PAR bZIP transcription factor
# 74772_362645              Chr5 q36                                                   NA             ATPase cation transporting 13A2
# 192188_282580                 Chr7                                                   NA                                  stabilin 2
# 54006_83632               Chr1 q41                                                   NA                  DEAF1 transcription factor
# 216565_305556            Chr14 q22                                                   NA                 EH domain binding protein 1
# 268445_360575            Chr10 q24                                                   NA                   ankyrin repeat domain 13B

MakeForestPlots<-function(metaOutputFDR_annotated, EntrezIDAsCharacter, species){
  
  #I originally wrote this function using only mouse Entrez IDs as input 
  #but then I realized that the function didn't work for genes that were only found in rats
  #so now the function allows either rat or mouse Entrez ids as input, and includes a conditional (if/else) statement
  
  if(species=="Mouse"){
    
    #This grabs the mouse gene symbol for the EntrezID from our meta-analysis annotation:
    MouseGeneSymbol<-metaOutputFDR_annotated$Mouse_Symbol[which(metaOutputFDR_annotated$Mouse_EntrezGene.ID==EntrezIDAsCharacter)]
    #This grabs the rat gene symbol for the EntrezID from our meta-analysis annotation:
    RatGeneSymbol<-metaOutputFDR_annotated$Rat_Symbol[which(metaOutputFDR_annotated$Mouse_EntrezGene.ID==EntrezIDAsCharacter)]
   
    #This grabs the Log2FC values for the EntrezID
    effect<-as.numeric(MetaAnalysis_FoldChanges_ForMeta[which(MetaAnalysis_FoldChanges_ForMeta$Mouse_EntrezGene.ID==EntrezIDAsCharacter),-c(1:3)])
    
    #This grabs the sampling variance (SV) values for the EntrezID
    var<-as.numeric(MetaAnalysis_SV_ForMeta[which(MetaAnalysis_FoldChanges_ForMeta$Mouse_EntrezGene.ID==EntrezIDAsCharacter),-c(1:3)]) 
    
  }else if(species=="Rat"){
    
    #This set of code does all of the same processes as above, but interpreting the EntrezID as a Rat Entrez ID:
    
    RatGeneSymbol<-metaOutputFDR_annotated$Rat_Symbol[which(metaOutputFDR_annotated$Rat_EntrezGene.ID==EntrezIDAsCharacter)]
    
    MouseGeneSymbol<-metaOutputFDR_annotated$Mouse_Symbol[which(metaOutputFDR_annotated$Rat_EntrezGene.ID==EntrezIDAsCharacter)]
    
    effect<-as.numeric(MetaAnalysis_FoldChanges_ForMeta[which(MetaAnalysis_FoldChanges_ForMeta$Rat_EntrezGene.ID==EntrezIDAsCharacter),-c(1:3)])
    
    var<-as.numeric(MetaAnalysis_SV_ForMeta[which(MetaAnalysis_FoldChanges_ForMeta$Rat_EntrezGene.ID==EntrezIDAsCharacter),-c(1:3)]) 
    
  }else{
    
    print("Please use either 'Mouse' or 'Rat' to indicate whether you are using annotation for mouse or rat genes")
    
  }
    
#This code makes the Forest Plot
  
  #First it opens up a .pdf file to output the plot into:
  #It automatically names that file with the mouse and rat ene symbols
  pdf(paste("ForestPlot_Mouse", MouseGeneSymbol,"Rat", RatGeneSymbol,".pdf", sep="_"), height=5, width=8)
  
  #This code makes the forest plot:
  #Note that the x-axis limits are currently set to -3 to 3
  #This may be too big or too small for visualizing the results for some genes.
  forest.rma(rma(effect, var), slab=colnames(MetaAnalysis_FoldChanges_ForMeta)[-c(1:3)],  xlim=c(-3, 3))
  
  #This code labels the forest plot with the mouse and rat gene symbols:
  mtext(paste("Mouse", MouseGeneSymbol, "Rat", RatGeneSymbol, sep="_"), line=-1.5, cex=2)
  
  #This closes the connection to the .pdf file, finishing the plot
  dev.off()
  
}

#Reading in the function:

#You can input this function by running the code discussed above to create the function in your R environment. 

#Alternatively, you can download the script for the function from our Github site and save the file in your working directory:
#https://github.com/hagenaue/BrainDataAlchemy/blob/main/MetaAnalysis_GemmaDEResults_2024/Function_MakeForestPlots.R

#And then source it from your working directory:
source("Function_MakeForestPlots.R")

#Example Usage:

#Note - this function currently uses Entrez ID (NCBI ID) as it's input
#It needs to be formatted as a character (not an integer) to work
#It can accept either mouse annotation (Entrez ID) or rat annotation (Entrez ID) as its input
                     
#Making a forest plot for gene Dbp (Mouse Entrez ID 13170):
                       
MakeForestPlots(metaOutputFDR_annotated, EntrezIDAsCharacter="13170", species="Mouse")

#Making a forest plot for gene Atp13a2 (Mouse Entrez ID 74772):

MakeForestPlots(metaOutputFDR_annotated, EntrezIDAsCharacter="74772", species="Mouse")

#Making a forest plot for gene Stab2 (Rat Entrez ID 282580):

MakeForestPlots(metaOutputFDR_annotated, EntrezIDAsCharacter="282580", species="Rat")