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

SearchTerms<-"(dentate gyrus stress*) OR 
(hippocamp* stress*) OR 
(Ammon's Horn stress*)"

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

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

result <- gemma.R ::search_datasets(SearchTerms, taxon = 'mouse',limit = 1) 

print(attributes(result)$totalElements)
#[1] 111

result <- gemma.R ::search_datasets(SearchTerms, taxon = 'rat',limit = 1) 

print(attributes(result)$totalElements)
#[1] 36

#Identifying datasets that fit our search terms in mice:

#First, determine how many datasets from mice fit our search:

result <- gemma.R ::search_datasets(SearchTerms, taxon = 'mouse',limit = 1) 

print(attributes(result)$totalElements)
#[1]  111
#There are 111 datasets that fit this search.

#Gemma.R can only download 100 dataset metadata records at a time, so if we have more than 100 results, we will need to perform the search interatively, with the number of iterations determined by the total number of datasets divided by 100:

SearchIterations<-floor(attributes(result)$totalElements/100)

#We can run the search for the first 100 records using the basic syntax:

TempResult<-gemma.R ::search_datasets(SearchTerms, taxon = 'mouse',limit = 100) 

#Then we run a loop to gather the next sets of 100 records, and add those records into the same dataframe as additional rows:

for(i in c(1:SearchIterations)){
  TempResult<-rbind.data.frame(TempResult,(gemma.R ::search_datasets(SearchTerms, taxon = 'mouse',limit = 100, offset=i*100)))
}

#I recommend renaming the object afterwards to make it easier to keep track of things later:
results_Stress_Hippocampus_Mice<-TempResult

#Cleaning up our workspace of objects that aren't needed anymore:
rm(result, TempResult, SearchIterations)

############

#Identifying datasets with our search terms in rats:

result <- gemma.R ::search_datasets(SearchTerms, taxon = 'rat',limit = 1) 

print(attributes(result)$totalElements)
#[1] 36
#There are less than 100 datasets that fit this search, so we do not need to iterate our search.

TempResult<-gemma.R ::search_datasets(SearchTerms, taxon = 'rat',limit = 100) 

#I renamed the object containing all of the rat datasets to make it easier to keep track of later:
results_Stress_Hippocampus_Rat<-TempResult

#Cleaning up our workspace of objects that aren't needed anymore:
rm(result, TempResult, SearchIterations)

combined<-rbind.data.frame(results_Stress_Hippocampus_Mice, results_Stress_Hippocampus_Rat)

unique_combined<- unique(combined)

#How many datasets are marked by Gemma as troubled?
table(unique_combined$experiment.Troubled)

#Excluding datasets marked as troubled:
unique_combined<-unique_combined[unique_combined$experiment.Troubled==FALSE,]

write.csv(unique_combined, "Metadata_forDatasets_toReview.csv")

install.packages("tidyverse")

library(tidyverse)


install.packages("plyr")

library(plyr)

#This code imports the summarized experiment object for dataset "GSE81672":

SummarizedExperiment_Filtered<-gemma.R::get_dataset_object("GSE81672", type = 'se', filter=TRUE, consolidate="average")

#If a summarized experiment object isn't available for a dataset (in this case dataset "GSE81672"), this is the error code that you will get:

SummarizedExperiment_Filtered<-gemma.R::get_dataset_object("GSE81672", type = 'se', filter=TRUE, consolidate="average")
# Error in `dplyr::mutate()`:
#   ! Can't transform a data frame with `NA` or `""` names.
# Run `rlang::last_trace()` to see where the error occurred.


#The expression data should be able to be imported into R using a different function:

GSE85136_Expression<-gemma.R::get_dataset_expression("GSE85136", filter=TRUE)
#There isn't an option to consolidate the data representing the same gene for this version of reading in data

str(GSE85136_Expression)
#Classes ‘data.table’ and 'data.frame':	23821 obs. of  28 variables:


#Renaming the object so the downstream code doesn't break:
#Note: Some of the downstream dimensions will be different
GSE85136_Expression_Filtered<-GSE85136_Expression


#We still need the metadata for this dataset object:


GSE85136_Samples<-gemma.R::get_dataset_design("GSE85136")

str(GSE85136_Samples)
#'data.frame':	24 obs. of  4 variables:

#Double-checking that it is in the same order as the expression data.frame:
row.names(GSE85136_Samples)
# [1] "ControlMale3"   "ControlFemale2" "GFPFemale1"     "GFPMale1"       "StressFemale2"  "CREMale3"      
# [7] "StressMale1"    "CREFemale1"     "StressFemale1"  "GFPMale3"       "CREFemale2"     "GFPMale2"      
# [13] "StressFemale3"  "StressMale2"    "ControlFemale3" "ControlFemale1" "CREMale2"       "GFPFemale3"    
# [19] "ControlMale2"   "GFPFemale2"     "StressMale3"    "ControlMale1"   "CREMale1"       "CREFemale3"    

colnames(GSE85136_Expression_Filtered)
# [1] "Probe"            "GeneSymbol"       "GeneName"         "NCBIid"           "Control Male 3"   "Control Male 1"  
# [7] "Control Male 2"   "Stress Male 2"    "Stress Male 1"    "Stress Male 3"    "GFP Male 3"       "GFP Male 2"      
# [13] "GFP Male 1"       "CRE Male 3"       "CRE Male 2"       "CRE Male 1"       "Control Female 3" "Control Female 1"
# [19] "Control Female 2" "Stress Female 3"  "Stress Female 2"  "Stress Female 1"  "GFP Female 3"     "GFP Female 2"    
# [25] "GFP Female 1"     "CRE Female 3"     "CRE Female 2"     "CRE Female 1" 

#Interestingly, it is *not* in the same order, and the names don't match.
#To make a version that matches, we would need to remove spaces

ExpressionColnames_NoSpaces<- gsub(" ", "", colnames(GSE85136_Expression_Filtered)[-c(1:4)])

#What the names look like with the spaces removed:
ExpressionColnames_NoSpaces
# [1] "ControlMale3"   "ControlMale1"   "ControlMale2"   "StressMale2"    "StressMale1"    "StressMale3"   
# [7] "GFPMale3"       "GFPMale2"       "GFPMale1"       "CREMale3"       "CREMale2"       "CREMale1"      
# [13] "ControlFemale3" "ControlFemale1" "ControlFemale2" "StressFemale3"  "StressFemale2"  "StressFemale1" 
# [19] "GFPFemale3"     "GFPFemale2"     "GFPFemale1"     "CREFemale3"     "CREFemale2"     "CREFemale1"  

#To get the metadata in the same order, we could use join or merge:

#To do this, we need to load the plyr package:
library(plyr)

#We also need to make the gene expression data column names it's own object (a data frame) with a column name ("x"):  
ExpressionColnames_NoSpaces_toJoin<-data.frame(x=ExpressionColnames_NoSpaces)

#We then need to make a data frame for the metadata what has a column of names that is also named "x":
GSE85136_Samples_toJoin<-data.frame(x=row.names(GSE85136_Samples), GSE85136_Samples)

#Then we can join the metadata to the gene expression data column names dataframe using type="left". That will make it so that the metadata is now in the same order as the gene expression data columns.
GSE85136_Samples_toJoin_inOrder<-join(ExpressionColnames_NoSpaces_toJoin, GSE85136_Samples_toJoin, by="x", type="left")

#We should confirm that the metadata is now in the same order as the gene expression data columns:
GSE85136_Samples_toJoin_inOrder$x
# [1] "ControlMale3"   "ControlMale1"   "ControlMale2"   "StressMale2"    "StressMale1"    "StressMale3"   
# [7] "GFPMale3"       "GFPMale2"       "GFPMale1"       "CREMale3"       "CREMale2"       "CREMale1"      
# [13] "ControlFemale3" "ControlFemale1" "ControlFemale2" "StressFemale3"  "StressFemale2"  "StressFemale1" 
# [19] "GFPFemale3"     "GFPFemale2"     "GFPFemale1"     "CREFemale3"     "CREFemale2"     "CREFemale1"   


#Now we can make a summarized experiment object using the gene expression data and the metadata:
GSE85136_SummarizedExperiment<-SummarizedExperiment(assays=list(counts=GSE85136_Expression_Filtered[,-c(1:4)]), colData=GSE85136_Samples_toJoin_inOrder, rowData=GSE85136_Expression_Filtered[,c(1:4)])

#Read in the .csv file containing the MultiQC report information into R:
GSE81672_LibrarySize<-read.csv("GSE81672_LibrarySize.csv", header=TRUE, stringsAsFactors = FALSE)

#The library size information is annotated with sample information using GSM number. 
head(GSE81672_LibrarySize)
# ExternalID X..Aligned M.Aligned X..Dups X..GC M.Seqs
# 1 GSM2166101     85.00%      20.1                   NA
# 2 GSM2166102     86.10%      33.4                   NA

#As far as I can tell, this identifier is not in the Summarized experiment object
#So we need to read in the experimental design info from the Gemma website, which contains GSM#

#To read in the experimental design file that contains the GSM#s for the samples:

GSE81672_expdesign<-read.delim("13458_GSE81672_expdesign.data.txt", sep="\t", comment.char="#", header=TRUE, stringsAsFactors = FALSE)


#This column has the GSM #
GSE81672_expdesign$ExternalID
# [1] "GSM2166189" "GSM2166195" "GSM2166180"

#This column includes the Name... and a bunch of other stuff:
GSE81672_expdesign$Bioassay
# [1] "GSE81672_Biomat_64___BioAssayId=450411Name=Sample102.NAC_susceptible_ketamine_non_responder"  
# [2] "GSE81672_Biomat_58___BioAssayId=450412Name=Sample108.HIP_susceptible_saline"
#...

#This column includes the Name... and a bunch of other stuff:
GSE81672_expdesign$Bioassay
# [1] "GSE81672_Biomat_64___BioAssayId=450411Name=Sample102.NAC_susceptible_ketamine_non_responder"  
# [2] "GSE81672_Biomat_58___BioAssayId=450412Name=Sample108.HIP_susceptible_saline"
#...

#This code splits these long names up at the "Name=" and then combines things back into a matrix again
temp<-do.call(rbind.data.frame, strsplit(GSE81672_expdesign$Bioassay, "Name="))

#This adds the name information to the expdesign matrix:
GSE81672_expdesign$SampleName<-temp[,2]
GSE81672_expdesign$SampleName

#For this dataset, The formatting is still a little different from the Summarized Experiment object, e.g.
GSE81672_expdesign$SampleName[1]
#"Sample102.NAC_susceptible_ketamine_non_responder"

row.names(colData(SummarizedExperiment_Filtered[[1]]))[1]
#"Sample 18: AMY_resilient_saline"
#One version uses a ., the other uses a ": "

#the ": " version is straight off of Gene Expression Omnibus (GEO), so I'm going to assume that for some reason Gemma doesn't like : in the design data.frame

#This is likely to be a dataset specific issue - I'm not sure how to generalize this code.
#For many datasets, things may be fine at this point.

#So we essentially need a find and replace - in R, these are often done with sub() or gsub()
#https://www.digitalocean.com/community/tutorials/sub-and-gsub-function-r
#https://www.rdocumentation.org/packages/base/versions/3.6.2/topics/grep

#this is a little funky, because we're not replacing simple alphanumeric characters, but punctuation that also has a meaning for coding
temp<-gsub("[.]", ": ", GSE81672_expdesign$SampleName)

#But it turns out that the names also differ because of an extra space. 
#To take care of this we just removed all spaces from both vectors because it is easier:

GSE81672_expdesign$SampleName_toJoin<-gsub(" ", "", temp)

SampleName_toJoin<-gsub(" ", "", row.names(colData(SummarizedExperiment_Filtered[[1]])))

#We can't just add library size as a column to experimental design because the samples aren't in the same order
#We have to use a "join" or "merge" function to align them instead

#To join the design matrix with the library size, we need to have two data.frames that have columns that have the same name:
str(GSE81672_LibrarySize)

#In this data.frame, the GSM# is called Sample.Name, whereas in the exp. design data.frame it is called ExternalID, so I will rename it ExternalID:
colnames(GSE81672_LibrarySize)[1]<-"ExternalID"

#Now we can join by these columns (the function "merge" also works for this):
GSE81672_expdesign_wLibrarySize<-join(GSE81672_expdesign, GSE81672_LibrarySize, by="ExternalID", type="left")


#Now we have to add this information into the Summarized Experiment object
#Which also has a different sample order
#I'm going to grab that sample order information and name it after the column with Sample Names in the experimental design object

SamplesOrderedLikeSummarizedExperiment<-data.frame(SampleName_toJoin=SampleName_toJoin)
                                
GSE81672_expdesign_wLibrarySize_Ordered<-join(SamplesOrderedLikeSummarizedExperiment, GSE81672_expdesign_wLibrarySize, by="SampleName_toJoin", type="left")                   


#Double checking that things are actually in the same order
cbind(row.names(colData(SummarizedExperiment_Filtered[[1]])), GSE81672_expdesign_wLibrarySize_Ordered$SampleName_toJoin)

#Then I can finally add library size to the Summarized Experiment object:
colData(SummarizedExperiment_Filtered[[1]])$LibrarySize<-GSE81672_expdesign_wLibrarySize_Ordered$M.Aligned

write.csv(rowData(SummarizedExperiment_Filtered[[1]]), "Example_Annotation.csv")

write.csv(colData(SummarizedExperiment_Filtered[[1]])[,-1], "Example_MetaData.csv")

write.csv(ExpressionData_Filtered, "ExpressionData_Filtered.csv")

#Pull out the matrix of gene expression data from the Summarized Experiment object:
ExpressionData_Filtered<-assay(SummarizedExperiment_Filtered[[1]])

#To observe the distribution of gene expression values, plot a histogram of the gene expression data:
hist(ExpressionData_Filtered)

#The range for the gene expression data can also be explored by examining the minimum and maximum values:
min(ExpressionData_Filtered)
max(ExpressionData_Filtered)

#The sample data - organism part, treatment, batch, etc:
colData(SummarizedExperiment_Filtered[[1]])
head(colData(SummarizedExperiment_Filtered[[1]]))

#To view the distribution of samples from different tissues or brain regions:

table(SummarizedExperiment_Filtered[[1]]$`organism part`)
#Example output:
# basolateral amygdaloid nuclear complex                      nucleus accumbens 
#           25                                                                                      25 
# prefrontal cortex                   ventral,Ammon's horn 
#            25                                     24

#Other common types of variables that you may want to view the distribution of samples for:

#Phenotype:
table(SummarizedExperiment_Filtered[[1]]$phenotype)

#Treatment:
table(SummarizedExperiment_Filtered[[1]]$treatment)

#For an example, I'm using an experiment with a weird/complex design

#This experiment includes samples from many tissues. I am only interested in the hippocampus ("ventral,Ammon's horn").

#I'm interested in antidepressant effects, and it looks like only the stress susceptible subjects were treated with antidepressants
#So I am going to subset down to just the stress susceptible subjects

#I'm going to keep both antidepressant responders and non-responders 
#because I'm looking at antidepressant effects (not antidepressant effectiveness)
#and most of my other datasets will not screen for antidepressant effectiveness


SampleFilter<-
  SummarizedExperiment_Filtered[[1]]$`organism part`=="ventral,Ammon's horn" &
  (SummarizedExperiment_Filtered[[1]]$phenotype=="susceptible toward,Responder"|
     SummarizedExperiment_Filtered[[1]]$phenotype=="Non-responder,susceptible toward"|
     SummarizedExperiment_Filtered[[1]]$phenotype=="susceptible toward")
   
#Subsetting the data to have only subjects with desired characteristics:
SummarizedExperiment_Subset<-SummarizedExperiment_Filtered[[1]][,SampleFilter]


#Sanity Check: Double-checking that the subsetting worked properly:

str(SummarizedExperiment_Subset)
# class: SummarizedExperiment 
# dim: 22782 16 
#This dataset is much smaller now

table(SummarizedExperiment_Subset$`organism part`)
#ventral,Ammon's horn 
#16 
#Only the hippocampal samples remain.

table(SummarizedExperiment_Subset$phenotype, SummarizedExperiment_Subset$treatment)
#                                   imipramine ketamine reference substance role,saline
# Non-responder,susceptible toward   4        3                               0
# susceptible toward                              0        0                               3
# susceptible toward,Responder           3        3                               0

#Pulling out a matrix of gene expression data for this subset (to use in functions that require matrices)

ExpressionData_Subset<-assay(SummarizedExperiment_Subset)

#Creating a matrix showing how each sample correlates with every other sample:
CorMatrix<-cor(ExpressionData_Subset)

#Writing that matrix out to the working directory to save it:
write.csv(CorMatrix, "CorMatrix.csv")

#Creating a hierarchically clustered heatmap illustrating the sample-sample correlation matrix:
heatmap(CorMatrix)

#Creating a boxplot illustrating the distribution of sample-sample correlations for each sample:
boxplot(CorMatrix)
#Outlier samples should be obvious at this point
#Gemma's algorithm uses the cut-off that the median sample-sample correlation for any particular sample should not fall outside the interquartile range (boxes) for all of the other samples. 
#We double-checked that there weren't pairs or small clusters of samples that fell outside the the interquartile range (boxes) for all of the other samples. 

#If there is an outlier sample, you can remove it using subsetting similar to above by identifying it's column name

#E.g.
OutlierFilter<-colnames(ExpressionData_Subset)!="Sample 15: HIP_susceptible_imipramine_non_responder"

OutlierFilter
#[1] FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE

SummarizedExperiment_Subset_noBad<-SummarizedExperiment_Subset[,OutlierFilter]

#Sanity check: the Summarized experiment object should now be one column smaller

SummarizedExperiment_Subset_noBad
# class: SummarizedExperiment 
# dim: 22782 15 

#Then recreate the ExpressionData Matrix from our new summarized experiment object that has had the outlier samples removed:

ExpressionData_Subset_noBad<-assay(SummarizedExperiment_Subset_noBad)

#If we *don't* have any outliers to remove, we can just rename our object so it works with the downstream code:

SummarizedExperiment_Subset_noBad<-SummarizedExperiment_Subset

#Then recreate the ExpressionData Matrix from our new summarized experiment object:

ExpressionData_Subset_noBad<-assay(SummarizedExperiment_Subset_noBad)

#We will need to calculate the sd for each of our experimental groups of interest.
#For example, for this dataset, I might be interested in the experimental groups defined by treatment.

#This function calculates the sd for each treatment group for a particular gene (row of data):
tapply(ExpressionData_Subset_noBad[1,], SummarizedExperiment_Subset_noBad$treatment, sd)
# imipramine                        ketamine reference substance role,saline 
# 0.6547312                       0.9273463                       0.7106250

#We want the minimum sd for all of our groups to not be zero, e.g.:
min(tapply(ExpressionData_Subset_noBad[1,], SummarizedExperiment_Subset_noBad$treatment, sd))!=0
#[1] TRUE

#... and we need to know that for all rows.
GenesToFilter<-apply(ExpressionData_Subset_noBad, 1, function(y) (min(tapply(y, SummarizedExperiment_Subset_noBad$treatment, sd))!=0))
head(GenesToFilter)

#This tells us how many genes (rows) we'll end up keeping:
sum(GenesToFilter)
#[1] 21686

SummarizedExperiment_Subset_noBad_Filtered<-SummarizedExperiment_Subset_noBad[GenesToFilter,]
SummarizedExperiment_Subset_noBad_Filtered

#Remaking the expression set from the filtered summarized experiment object:
ExpressionData_Subset_noBad_Filtered<-assay(SummarizedExperiment_Subset_noBad_Filtered)
str(ExpressionData_Subset_noBad_Filtered)

table(SummarizedExperiment_Subset_noBad_Filtered$block)

# Device=11V6WR1:Run=206:Flowcell=C3HBLACXX:Lane=1    Device=11V6WR1:Run=206:Flowcell=C3HBLACXX:Lane=3 
# 1                                                   3 
# Device=HWI-D00147:Run=59:Flowcell=C2GWMACXX:Lane=2  Device=HWI-D00147:Run=59:Flowcell=C2GWMACXX:Lane=4 
# 3                                                   2 
# Device=HWI-D00147:Run=59:Flowcell=C2GWMACXX:Lane=5  Device=HWI-D00147:Run=59:Flowcell=C2GWMACXX:Lane=6 
# 1                                                   1 
# Device=HWI-D00147:Run=59:Flowcell=C2GWMACXX:Lane=7  Device=HWI-D00147:Run=59:Flowcell=C2GWMACXX:Lane=8 
# 1                                                   1 
# Device=HWI-ST1147:Run=158:Flowcell=C34A7ACXX:Lane=3 Device=HWI-ST1147:Run=158:Flowcell=C34A7ACXX:Lane=7 
# 1                                                   1 
# Device=HWI-ST276:Run=345:Flowcell=C3H70ACXX:Lane=1 
# 1 

#This function breaks apart these "blocks" into specific batch-related variables (device, run, flow cell, lane):
strsplit(SummarizedExperiment_Subset_noBad_Filtered$block, ":")

#To make that into an easier-to-use data.frame
BatchVariables<-do.call(rbind.data.frame, strsplit(SummarizedExperiment_Subset_noBad_Filtered$block, ":")) 

str(BatchVariables)

#I'm going to rename these (note - this will be dataset specific):
colnames(BatchVariables)<-c("Device", "Run", "FlowCell", "Lane")

table(SummarizedExperiment_Subset_noBad_Filtered$treatment, BatchVariables$Device)
#                                   Device=11V6WR1 Device=HWI-D00147 Device=HWI-ST1147 Device=HWI-ST276
# imipramine                                   2                 4                 1                0
# ketamine                                     1                 4                 1                0
# reference substance role,saline              1                 1                 0                1

table(SummarizedExperiment_Subset_noBad_Filtered$treatment, BatchVariables$Run)
#                                   Run=158 Run=206 Run=345 Run=59
# imipramine                            1       2       0      4
# ketamine                              1       1       0      4
# reference substance role,saline       0       1       1      1

#Looks like run and device may be redundant
table(BatchVariables$Device, BatchVariables$Run)
#                     Run=158 Run=206 Run=345 Run=59
# Device=11V6WR1          0       4       0      0
# Device=HWI-D00147       0       0       0      9
# Device=HWI-ST1147       2       0       0      0
# Device=HWI-ST276        0       0       1      0

table(SummarizedExperiment_Subset_noBad_Filtered$treatment, BatchVariables$FlowCell)
#Looks also potentially redundant

table(BatchVariables$FlowCell, BatchVariables$Run)
#                       Run=158 Run=206 Run=345 Run=59
# Flowcell=C2GWMACXX       0       0       0      9
# Flowcell=C34A7ACXX       2       0       0      0
# Flowcell=C3H70ACXX       0       0       1      0
# Flowcell=C3HBLACXX       0       4       0      0
#yep

table(SummarizedExperiment_Subset_noBad_Filtered$treatment, BatchVariables$Lane)
#                                   Lane=1 Lane=2 Lane=3 Lane=4 Lane=5 Lane=6 Lane=7 Lane=8
# imipramine                           1      1      1      1      0      1      1      1
# ketamine                             0      1      2      1      1      0      1      0
# reference substance role,saline      1      1      1      0      0      0      0      0
#Too many Lanes to take into consideration

#Run PCA on the transposed gene expression matrix
#In this case, ExpressionData_Subset_noBad_Filtered is the object containing the gene expression matrix
pca_output<-prcomp(t(ExpressionData_Subset_noBad_Filtered), scale=TRUE)

#Output a scree plot for the PCA results:

#This plot illustrates the proportion of variance explained by each principal component (PC):
png("10 PCA Scree Plot1.png")
plot(summary(pca_output)$importance[2,]~(c(1:length(summary(pca_output)$importance[2,]))), main="Variance Explained by Each Principal Component", xlab="PC #", ylab="Proportion of Variance Explained", col=2)
dev.off()

#Extract out the first four principal components and make them named R objects for graphing purposes:
PC1<-pca_output$x[,1]
PC2<-pca_output$x[,2]
PC3<-pca_output$x[,3]
PC4<-pca_output$x[,4]

#Extract the eigenvectors for PC1-4 and make them a named R object
PCeigenvectors<-pca_output$rotation[ ,c(1:4)]

#Add the gene names (rowData) from the gene expression matrix to the object with the eigenvectors
PCeigenvectors2<-cbind(PCeigenvectors, rowData(SummarizedExperiment_Subset_noBad_Filtered))

#Write this information out as a .csv file:
write.csv(PCeigenvectors2, "PCeigenvectors.csv")

#Output a scatterplot illustrating the relationship between Principal components 1 & 2 (PC1 & PC2):
png("10 PC1 vs PC2.png")
plot(PC1~PC2, main="Principal Components Analysis")
dev.off()

#Examining the relationship between PC1 and a categorical variable using a boxplot
#In this case, the plotting code is referring back to the summarized experiment object containing the gene expression matrix used in the PCA analysis


boxplot(PC1~SummarizedExperiment_Subset_noBad_Filtered$treatment, las=2, xlab="")

#The relationship between the PC and continuous numeric variables can be plotted as a scatterplot using the "plot" function instead of "boxplot" 

#Renaming conditions
#Some of the levels for my variables are a mouthful. 
#I'm going to rename them:

#The names of the conditions (levels) for the variable with the number of subjects in each condition:
table(SummarizedExperiment_Subset_noBad_Filtered$treatment)
#imipramine               ketamine reference substance role,saline 
#7                               6                              11 

#I'm going to rename the control condition because it is a mouthful:
SummarizedExperiment_Subset_noBad_Filtered$treatment[SummarizedExperiment_Subset_noBad_Filtered$treatment=="reference substance role,saline"]<-"saline"

#The new names:
table(SummarizedExperiment_Subset_noBad_Filtered$treatment)
# imipramine   ketamine     saline 
# 7          6         11

#Renaming the phenotype conditions as well:
SummarizedExperiment_Subset_noBad_Filtered$phenotype[SummarizedExperiment_Subset_noBad_Filtered$phenotype=="Non-responder,susceptible toward"]<-"non-responder"
SummarizedExperiment_Subset_noBad_Filtered$phenotype[SummarizedExperiment_Subset_noBad_Filtered$phenotype=="susceptible toward,Responder"]<-"responder"
SummarizedExperiment_Subset_noBad_Filtered$phenotype[SummarizedExperiment_Subset_noBad_Filtered$phenotype=="susceptible toward"]<-"untreated"

#The new names:
table(SummarizedExperiment_Subset_noBad_Filtered$phenotype)
# non-responder     responder     untreated 
# 7             6             3

#Setting up our categorical variables so that they have intuitive reference groups:

#First, we need to define the categorical variable as a factor:
#I tend to make this a new variable in the Summarized experiment object instead of overwriting the old categorical variable:
SummarizedExperiment_Subset_noBad_Filtered$treatment_factor<-as.factor(SummarizedExperiment_Subset_noBad_Filtered$treatment)

#R automatically assigns a reference level for this new factor that may or may not make intuitive sense
#When viewing the levels for a factor, the reference variable is whatever is listed first:
levels(SummarizedExperiment_Subset_noBad_Filtered$treatment_factor)
#[1] "imipramine" "ketamine"   "saline" 

#We want saline to be our control (reference) group:
SummarizedExperiment_Subset_noBad_Filtered$treatment_factor<-relevel(SummarizedExperiment_Subset_noBad_Filtered$treatment_factor, ref="saline")

#Plotting data for a particular gene

#Getting the expression data for a particular gene (in this case, parvalbumin):
Pvalb<-assay(SummarizedExperiment_Subset_noBad_Filtered)[rowData(SummarizedExperiment_Subset_noBad_Filtered)$GeneSymbol=="Pvalb",]

#Plotting a boxplot of the expression of Pvalb Log2 CPM across antidepressant groups:
boxplot(Pvalb~SummarizedExperiment_Subset_noBad_Filtered$treatment, xlab="", ylab="Log2 CPM", main="Pvalb", las=2)

#We can add jittered data points to our graph so that we can see the values for the individual samples in each group.
#To do this, we use the function stripchart and the exact same y~x formula that we used for the boxplot
#The parameter pch determines the size of the data points.
#The parameter "add" places the points on top of the boxplot that we already created.
stripchart(Pvalb~SummarizedExperiment_Subset_noBad_Filtered$treatment, pch = 19, method = "jitter", jitter = 0.2, vertical = TRUE, add=TRUE)

#Example 1: Creating a statistical model matrix that includes experimental treatment (which is a factor, with the reference level defined as saline) and the technical covariate of Library Size:
design <- model.matrix(~treatment_factor+LibrarySize, data=colData(SummarizedExperiment_Subset_noBad_Filtered))

#Example 2: Creating a statistical model matrix that includes experimental treatment (which is a factor, with the reference level defined as saline), the technical covariate of Library Size, and their interaction:
design2 <- model.matrix(~treatment_factor+LibrarySize + treatment_factor*LibrarySize, data=colData(SummarizedExperiment_Subset_noBad_Filtered))

#Running a basic linear regression to double check our differential expression model:

#We have to tell the linear model (lm) function to not add an intercept because the model matrix already defines the intercept
summary.lm(lm(Pvalb~0+design)) 

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

BiocManager::install("limma")

library(limma)

#Applying the differential expression model to all genes using limma:

fit <- lmFit(ExpressionData_Subset_noBad_Filtered, design)

#Adding an eBayes correction to help reduce the influence of outliers/small sample size on estimates
efit <- eBayes(fit, trend=TRUE)

#Applying the differential expression model to all genes using limma:

fit <- lmFit(ExpressionData_Subset_noBad_Filtered, design)

#Adding an eBayes correction to help reduce the influence of outliers/small sample size on estimates
efit <- eBayes(fit)

#Applying a false discovery rate (FDR) correction and writing out the results as a .txt file:
write.fit(efit, adjust="BH", file="Limma_results.txt")

#Writing out the gene annotation for the results as a .csv file:
write.csv(rowData(SummarizedExperiment_Subset_noBad_Filtered), "Annotation_LimmaResults.csv")

#Example code: Reading in a file containing the differential expression results for dataset GSE92718
#Note: This code will only work if the working directory is set to the folder that contains this file of differential expression results
#Note: Depending on file format, you may need to use read.delim instead of read.csv:
TempResultsJoined<-read.csv("Annotation_LimmaResults_GSE92718_all.csv", header=TRUE, stringsAsFactors = FALSE)

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

FilteringDEResults_GoodAnnotation(TempResultsJoined)
# [1] "# of rows with missing NCBI annotation:"
# [1] NA
# [1] "# of rows with missing Gene Symbol annotation:"
# [1] 5
# [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] 16867
# [1] "Outputted object: TempResultsJoined_NoNA_NoMultimapped"

#This code will be dataset specific
#We need to extract the Log2FC ("Coef") and T-statistic ("t.") columns for the statistical contrasts relevant to our meta-analysis and place them into their own matrix

TempResultsJoined_NoNA_NoMultimapped_FoldChanges<-cbind(TempResultsJoined_NoNA_NoMultimapped$Coef.interaction.stress.x.female)

TempResultsJoined_NoNA_NoMultimapped_Tstats<-cbind(TempResultsJoined_NoNA_NoMultimapped$t.interaction.stress.x.female)


#Making the row names for the Log2FC and Tstat matrices the Entrez ID gene annotation:
row.names(TempResultsJoined_NoNA_NoMultimapped_FoldChanges)<-TempResultsJoined_NoNA_NoMultimapped$NCBIid

row.names(TempResultsJoined_NoNA_NoMultimapped_Tstats)<-TempResultsJoined_NoNA_NoMultimapped$NCBIid

#Let's rename our columns to something nicer describing the effect of interest:
#Note - we later discovered that this name needs to include the dataset identifier (GSEID#) for later joining and plotting purposes

ComparisonsOfInterest<-c("GSE85136_StressEffects_InteractionWSex")

colnames(TempResultsJoined_NoNA_NoMultimapped_FoldChanges)<-ComparisonsOfInterest
colnames(TempResultsJoined_NoNA_NoMultimapped_Tstats)<-ComparisonsOfInterest

ExtractingDEResults<-function(GSE_ID, TempResultsJoined_NoNA_NoMultimapped_FoldChanges, TempResultsJoined_NoNA_NoMultimapped_Tstats){
  
  #We calculate the standard error by dividing the log2FC by the tstat
  TempResultsJoined_NoNA_NoMultimapped_SE<-TempResultsJoined_NoNA_NoMultimapped_FoldChanges/TempResultsJoined_NoNA_NoMultimapped_Tstats
  str(TempResultsJoined_NoNA_NoMultimapped_SE)
  
  #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.
  
  TempResultsJoined_NoNA_NoMultimapped_SV<-(TempResultsJoined_NoNA_NoMultimapped_SE)^2
  str(TempResultsJoined_NoNA_NoMultimapped_SV)
  
  TempMasterResults<-list(Log2FC=TempResultsJoined_NoNA_NoMultimapped_FoldChanges, Tstat=TempResultsJoined_NoNA_NoMultimapped_Tstats, SE=TempResultsJoined_NoNA_NoMultimapped_SE, SV=TempResultsJoined_NoNA_NoMultimapped_SV)
  
  assign(paste("DEResults", GSE_ID, sep="_"), TempMasterResults, envir = as.environment(1))
  
  print(paste("Output: Named DEResults", GSE_ID, sep="_"))
  
  rm(TempMasterResults, TempResultsJoined_NoNA_NoMultimapped_SV, TempResultsJoined_NoNA_NoMultimapped_SE, TempResultsJoined_NoNA_NoMultimapped_FoldChanges, TempResultsJoined_NoNA_NoMultimapped_Tstats)
  
}

#Example usage for the extracting results function:

ExtractingDEResults(GSE_ID="GSE92718", TempResultsJoined_NoNA_NoMultimapped_FoldChanges, TempResultsJoined_NoNA_NoMultimapped_Tstats)
# num [1:16862, 1] 0.1655 0.0943 0.1069 0.1443 0.0382 ...
# - attr(*, "dimnames")=List of 2
# ..$ : chr [1:16862] "13654" "215418" "13656" "19252" ...
# ..$ : chr "GSE92718_CuffOperation_vs_Ctrl"
# num [1:16862, 1] 0.02739 0.00889 0.01142 0.02083 0.00146 ...
# - attr(*, "dimnames")=List of 2
# ..$ : chr [1:16862] "13654" "215418" "13656" "19252" ...
# ..$ : chr "GSE92718_CuffOperation_vs_Ctrl"
# [1] "Output: Named DEResults_GSE92718"


#################

#Cleaning up our environment before moving on to the next dataset:

rm(TempResultsJoined, TempResultsJoined_NoNA_NoMultimapped, TempResultsJoined_NoNA_NoMultimapped_FoldChanges, TempResultsJoined_NoNA_NoMultimapped_Tstats, ComparisonsOfInterest)

str(DEResults_GSE92718)
# List of 4
# $ Log2FC: num [1:16862, 1] 1.386 0.668 0.72 0.974 0.245 ...
# ..- attr(*, "dimnames")=List of 2
# .. ..$ : chr [1:16862] "13654" "215418" "13656" "19252" ...
# .. ..$ : chr "GSE92718_CuffOperation_vs_Ctrl"
# $ Tstat : num [1:16862, 1] 8.37 7.09 6.74 6.75 6.41 ...
# ..- attr(*, "dimnames")=List of 2
# .. ..$ : chr [1:16862] "13654" "215418" "13656" "19252" ...
# .. ..$ : chr "GSE92718_CuffOperation_vs_Ctrl"
# $ SE    : num [1:16862, 1] 0.1655 0.0943 0.1069 0.1443 0.0382 ...
# ..- attr(*, "dimnames")=List of 2
# .. ..$ : chr [1:16862] "13654" "215418" "13656" "19252" ...
# .. ..$ : chr "GSE92718_CuffOperation_vs_Ctrl"
# $ SV    : num [1:16862, 1] 0.02739 0.00889 0.01142 0.02083 0.00146 ...
# ..- attr(*, "dimnames")=List of 2
# .. ..$ : chr [1:16862] "13654" "215418" "13656" "19252" ...
# .. ..$ : chr "GSE92718_CuffOperation_vs_Ctrl"

#Sanity Check:

#Calculating the SE from the Log2FC and Tstat for the first gene(row)
1.386/8.37
#[1] 0.1655914

#Calculating the SV for the first gene
0.1655914^2
#[1] 0.02742051
#Looks good (the decimal differences are probably just due to rounding)

#Calculating the SE from the Log2FC and Tstat for the second gene(row)
0.668/7.09
#[1] 0.09421721

#Calculating the SV for the second gene
0.09421721^2
#[1] 0.008876883
#Looks good (the decimal differences are probably just due to rounding)

AligningRatDatasets<-function(ListOfRatDEResults){
  
  Rat_MetaAnalysis_FoldChange_Dfs<-list()
  
  for(i in c(1:length(ListOfRatDEResults))){
    Rat_MetaAnalysis_FoldChange_Dfs[[i]]<-data.frame(Rat_EntrezGene.ID=row.names(ListOfRatDEResults[[i]][[1]]),ListOfRatDEResults[[i]][[1]], stringsAsFactors=FALSE)
  }
  
  print("Rat_MetaAnalysis_FoldChange_Dfs:")
  print(str(Rat_MetaAnalysis_FoldChange_Dfs))
  
  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)
  
  print("Rat_MetaAnalysis_FoldChanges:")
  print(str(Rat_MetaAnalysis_FoldChanges))
  
  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))
  
  rm(Rat_MetaAnalysis_SV_Dfs, Rat_MetaAnalysis_FoldChange_Dfs)
}




###########

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)
}

#Example usage for rat datasets:

ListOfRatDEResults<-list(DEResults_GSE172133)

AligningRatDatasets(ListOfRatDEResults)


#Example usage for mouse datasets:

ListOfMouseDEResults<-list(DEResults_GSE85136, DEResults_GSE92718)

AligningMouseDatasets(ListOfMouseDEResults)

#Read in the rat vs. mouse gene orthology database from the working directory:
MouseVsRat_NCBI_Entrez<-read.csv("MouseVsRat_NCBI_Entrez_JacksonLab_20230727.csv", header=TRUE, stringsAsFactors = FALSE, row.names=1, colClasses=c("character", "character", "character"))

#Join the rat vs. mouse gene orthology information to the mouse Log2FC results:
Mouse_MetaAnalysis_FoldChanges_wOrthologs<-join(MouseVsRat_NCBI_Entrez, Mouse_MetaAnalysis_FoldChanges, by="Mouse_EntrezGene.ID", type="full")

#Join the rat vs. mouse gene orthology information to the mouse sampling variance results:
Mouse_MetaAnalysis_SV_wOrthologs<-join(MouseVsRat_NCBI_Entrez, Mouse_MetaAnalysis_SV, by="Mouse_EntrezGene.ID", type="full")


#If there are rat datasets, join the mouse results with the rat results (for both the Log2FC and SV)
MetaAnalysis_FoldChanges<-join(Mouse_MetaAnalysis_FoldChanges_wOrthologs, Rat_MetaAnalysis_FoldChanges, by="Rat_EntrezGene.ID", type="full")

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


#If there *aren't* any rat datasets, just rename the mouse Log2FC and SV results in a manner so that the downstream code will still work:

MetaAnalysis_FoldChanges<-Mouse_MetaAnalysis_FoldChanges_wOrthologs

MetaAnalysis_SV<-Mouse_MetaAnalysis_SV_wOrthologs


#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="_")

RunBasicMetaAnalysis<-function(NumberOfComparisons, CutOffForNAs, MetaAnalysis_FoldChanges, MetaAnalysis_SV){
  
  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))
  
  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 (columns 2-10) for each gene (row):
  for(i in c(1:nrow(MetaAnalysis_FoldChanges_ForMeta))){
    
    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(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)
  }
  
  colnames(metaOutput)<-c("Log2FC_estimate", "SE", "pval", "CI_lb", "CI_ub", "Number_Of_Comparisons")
  row.names(metaOutput)<-MetaAnalysis_FoldChanges_ForMeta[,3]
  
  metaOutput<<-metaOutput
  MetaAnalysis_Annotation<<-MetaAnalysis_FoldChanges_ForMeta[,c(1:3)]
  return(metaOutput)
  return(MetaAnalysis_Annotation)
  
  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 actually 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 
# 11043  4823 13346   875

#Let's try running a meta-analysis using genes that were found in at least 2 sets of differential expression results
#Since there are 3 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=3
CutOffForNAs=2
#I have 3 comparisons
#2 NA is too many

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...

# [1] "Table of # of NAs per Row (Gene):"
# MetaAnalysis_FoldChanges_NAsPerRow
# 0     1     2     3 
# 11043  4823 13346   875 
# [1] "MetaAnalysis_FoldChanges_ForMeta:"
# 'data.frame':	15866 obs. of  6 variables:
# $ Rat_EntrezGene.ID            : chr  "498097" "114521" "360576" "24628" ...
# $ Mouse_EntrezGene.ID          : chr  "68980" "21665" "237858" "18591" ...
# $ MouseVsRat_EntrezGene.ID     : chr  "68980_498097" "21665_114521" "237858_360576" "18591_24628" ...
# $ StressEffects_InteractionWSex: num  0.0831 -0.0504 -0.7557 -0.0782 0.161 ...
# $ CuffOperation_vs_Ctrl        : num  0.01656 0.00964 -0.0427 0.01155 NA ...
# $ ChronicConstriction_vs_Ctrl  : num  -0.0608 0.0631 -0.4516 0.4982 1.4642 ...
# NULL

#Example Output:

str(metaOutput)
# num [1:15866, 1:6] 0.012 0.0143 -0.1497 0.157 1.3608 ...
# - attr(*, "dimnames")=List of 2
# ..$ : chr [1:15866] "68980_498097" "21665_114521" "237858_360576" "18591_24628" ...
# ..$ : chr [1:6] "Log2FC_estimate" "SE" "pval" "CI_lb" ...

head(metaOutput)
# Log2FC_estimate         SE       pval       CI_lb      CI_ub
# 68980_498097       0.01203600 0.05083559 0.81284048 -0.08759993 0.11167192
# 21665_114521       0.01426912 0.03385763 0.67343034 -0.05209062 0.08062887
# 237858_360576     -0.14968156 0.17385300 0.38925664 -0.49042717 0.19106405
# 18591_24628        0.15695379 0.17760171 0.37683643 -0.19113918 0.50504675
# 229323_688737      1.36081008 0.45349519 0.00269346  0.47197584 2.24964432
# 240549_309227      0.04938634 0.14315287 0.73010174 -0.23118814 0.32996081
# Number_Of_Comparisons
# 68980_498097                      3
# 21665_114521                      3
# 237858_360576                     3
# 18591_24628                       3
# 229323_688737                     2
# 240549_309227                     2                    

tail(metaOutput)
# Log2FC_estimate         SE       pval       CI_lb      CI_ub
# 97423_NA      0.08335277 0.22999011 0.71703901 -0.36741957 0.53412510
# 97775_NA      0.19686310 0.17990823 0.27384896 -0.15575056 0.54947676
# 98303_NA      0.06781550 0.04093082 0.09755343 -0.01240744 0.14803843
# 98736_NA      0.03495936 0.04096462 0.39343552 -0.04532983 0.11524854
# 98870_NA      0.01442392 0.03195141 0.65167762 -0.04819970 0.07704753
# 99100_NA      0.11964317 0.11449473 0.29603839 -0.10476237 0.34404872
# Number_Of_Comparisons
# 97423_NA                     2
# 97775_NA                     2
# 98303_NA                     2
# 98736_NA                     2
# 98870_NA                     2
# 99100_NA                     2

#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:

HOM_MouseVsRat <- read.csv("HOM_MouseVsRat.csv", header = TRUE, row.names = 1)
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){
  
  tempPvalAdjMeta<-mt.rawp2adjp(metaOutput[,3], proc=c("BH"))
  
  metaPvalAdj<-tempPvalAdjMeta$adjp[order(tempPvalAdjMeta$index),]
  
  metaOutputFDR<-cbind(metaOutput, metaPvalAdj[,2])
  
  colnames(metaOutputFDR)[7]<-"FDR"
  
  metaOutputFDR<<-metaOutputFDR
  
  print("metaOutputFDR:")
  print(str(metaOutputFDR))
  
  write.csv(metaOutputFDR, "metaOutputFDR.csv")
  
  #a version of the output in order by p-value:
  metaOutputFDR_OrderbyPval<<-metaOutputFDR[order(metaOutputFDR[,3]),]
  
  #Let's write out a version of the output in order by p-value:
  write.csv(metaOutputFDR_OrderbyPval, "metaOutputFDR_orderedByPval_wHDRFData.csv")
  
  print("Do we have any genes that are statistically significant following false discovery rate correction?")
  print(sum(metaOutputFDR[,7]<0.10, na.rm=TRUE))
  
  print("What are the top results?")
  print(head(metaOutputFDR[order(metaOutputFDR[,3]),]))
  
  rm(tempPvalAdjMeta, metaPvalAdj)
  
}

#Example usage:

FalseDiscoveryCorrection(metaOutput)
# [1] "metaOutputFDR:"
# num [1:11562, 1:7] 0.06022 -0.01543 -0.06567 -0.07498 -0.00186 ...
# - attr(*, "dimnames")=List of 2
# ..$ : chr [1:11562] "A4galt" "Aaas" "Aacs" "Aadac" ...
# ..$ : chr [1:7] "Log2FC_estimate" "SE" "pval" "CI_lb" ...
# NULL
# [1] "Do we have any genes that are statistically significant following false discovery rate correction?"
# [1] 797
# [1] "What are the top results?"
# Log2FC_estimate         SE         pval      CI_lb       CI_ub Number_Of_Comparisons          FDR
# Parm1       -0.1798498 0.02389052 5.149049e-14 -0.2266743 -0.13302524                     7 5.953331e-10
# Lpl          0.2138084 0.03015666 1.341870e-12  0.1547024  0.27291434                     7 7.757348e-09
# Gabrr2       0.1965106 0.02859386 6.309657e-12  0.1404677  0.25255358                     7 2.431742e-08
# Pla2g5       0.1593702 0.02375177 1.948623e-11  0.1128176  0.20592280                     7 5.632494e-08
# Akap8l      -0.1174562 0.01848935 2.116479e-10 -0.1536947 -0.08121775                     7 4.894145e-07
# Dusp1       -0.4415588 0.07101670 5.045755e-10 -0.5807489 -0.30236858                     7 9.723171e-07

MakeForestPlots<-function(EntrezIDAsCharacter){
  
  pdf(paste("ForestPlot_", EntrezIDAsCharacter, ".pdf", sep=""), height=5, width=8)
  
  effect<-as.numeric(MetaAnalysis_FoldChanges_ForMeta[MetaAnalysis_FoldChanges_ForMeta$Mouse_EntrezGene.ID==EntrezIDAsCharacter,-c(1:3)])
  var<-as.numeric(MetaAnalysis_SV_ForMeta[MetaAnalysis_FoldChanges_ForMeta$Mouse_EntrezGene.ID==EntrezIDAsCharacter,-c(1:3)])
  
  forest.rma(rma(effect, var),slab=colnames(MetaAnalysis_FoldChanges_ForMeta)[-c(1:3)],  xlim=c(-3, 3))
  
  mtext(paste(EntrezIDAsCharacter), line=-1.5, cex=2)
  dev.off()
}

#Example Usage:

MakeForestPlots("215418")