Methods
AD: Alzheimer’s disease
All bioinformatics and statistical analysis were performed using R version 4.2.1 software [1].
Study Search and Selection
Available datasets were collected from the Gene Expression Omnibus (GEO) [2] public repository. A systematic search of all published studies in public repositories (2002-2021) was conducted during 2021, following the preferred reporting items for systematic reviews and meta-analyses (PRISMA) guidelines [3]. Keywords employed in the search were “Alzheimer,” “Alzheimer’s Disease”, and “AD”. The following inclusion criteria were applied:
Transcriptomic studies on Homo sapiens.
Control and AD-affected patients included.
Sex, disease/control status, age, and brain region variables registered.
RNA extracted directly from post-mortem brain tissues (no cell lines or cultures).
Brain tissues from either the CT or HP.
Sample size > 3 for case and control groups in both sexes.
Normalized gene expression data of 8 microarray AD datasets (GSE118553, GSE1297, GSE132903, GSE15222, GSE29378, GSE37263, GSE48350, GSE5281, and GSE84422) and the raw counts’ matrix of the GSE125583 RNA-sequencing study from the GEO repository were retrieved. Data downloading was performed using the GEOquery R package [4] for microarray studies and manually for the RNA-sequencing study.
Individual Transcriptomics Analysis
For each selected study, an individual transcriptomics analysis comprised two steps: preprocessing and differential expression analysis. Data preprocessing included the standardization of the terminology for the clinical variables in each study, the homogenization of gene annotation, and exploratory data analysis. For the microarray datasets, the normalization methods performed by the original authors were assessed, log2 transforming data matrices when necessary. All probe sets were annotated to HUGO gene symbols [5] using the annotation provided by each microarray platform. The median of expression values was calculated when dealing with duplicated probe-to-symbol mappings. For the RNA-sequencing dataset, the count matrix was preprocessed using the edgeR R package [6] and transformed using the Voom function included in the limma R package [7]. The exploratory analysis included unsupervised clustering and PCA to detect expression patterns between samples and genes and the presence of batch effects in each study. Differential expression analyses were performed using the limma R package to detect the sex-based differential expression of genes. To achieve this goal, the following comparison was applied:
(AD.female - Control.female) - (AD.male - Control.male)
This comparison allows the detection of genes with a sex-based differential impact on disease (SDID). Genes with a log2 fold change (LFC) greater than zero show either a higher increase or a lesser decrease in expression in female AD patients when comparing the effect of the disease between sexes. On the contrary, genes with an LFC lower than zero have a higher increase or a lesser decrease in expression in male AD patients when comparing the effect of the disease between sexes. To gain a better understanding of the sex-based differential behavior detected with the previous comparison, two additional comparisons were performed: A case vs. control comparison performed only in females (AD.female - Control.female) that evaluates the impact of AD in females (IDF) and another performed only in males (AD.male - Control.male), that informs us about the impact AD in males (IDM). These comparisons were applied separately to all CT and HP samples. The patients’ age was included in the limma linear model as a blocking factor to reduce its impact on the results. P-values were corrected using the Benjamini-Hochberg procedure [8] and considered significant when below a threshold of 0.05.
Gene Expression Meta-Analysis
Differential gene expression results were integrated into a meta-analysis [9] for each brain region (CT and HP). Meta-analyses were implemented with the R package metafor [10] under the DerSimonian & Laird random-effects model [11], considering individual study heterogeneity. This model considers the variability of individual studies by increasing the weights of studies with less variability when computing meta-analysis results. Thus, the most robust functions between studies are highlighted.
P-values, corrected p-values, LFC, LFC standard error (SE), and LFC 95% confidence intervals (CI) were calculated for each evaluated gene. Genes with corrected p-values lower than 0.05 were considered significant, and both funnel and forest plots were computed for each. These representations were evaluated to assess for possible biased results, where the LFC represents the effect size of a gene, and the SE of the LFC serves as a study precision measure [12]. Sensitivity analysis (leave-one-out cross-validation) was conducted for each significant gene to verify alterations in the results due to the inclusion of any study. The Open Targets platform (release 22.09) [13] was used to explore the associations of significant genes with AD.
Sex-based Functional Signature in the CT
Gene meta-analysis of CT data revealed gene sets with a significantly differential expression pattern between male and female AD patients. Several analyses were conducted to identify the functional implications of these differences.
Over-Representation Analysis (ORA) [14] through clusterProfiler and ReactomePA R packages [15, 16] was first used to determine the biological functions and pathways overrepresented in all significant gene subsets. P-values and corrected p-values were calculated for each GO (Gene Ontology) term from the “Biological Processes” GO ontology [17] and each Reactome pathway [18]. Every function and pathway with a corrected p-value lower than 0.05 was labeled as over-represented in each gene set.
Protein-protein interaction (PPI) networks were then calculated using the STRING web tool for each subset of genes [19]. The total number of edges was examined, and PPI enrichment was assessed using the default parameters for each network.
A VIPER analysis [20] with human regulons obtained from the DoRothEA R package [21] was performed to estimate transcription factor activity. Regulons with a confidence level of A, B, C, or D were selected, excluding those with less than 25 genes (n = 217). The p-values were corrected using the Benjamini & Hochberg method. Normalized enrichment scores (NES) were calculated by VIPER as a measure of relative transcription factor activity.
Metafun-AD Web Tool
All data and results generated in the different steps of the meta-analysis are available in the Metafun-AD web tool [22], which is freely accessible to any user and allows the confirmation of the results described in this manuscript and the exploration of other results of interest. The front end was developed using Quarto 1.2 [23] , and the interactive graphics used in this web resource have been implemented with plotly [24].
This easy-to-use resource is divided into eight sections: 1) framework and summary of analysis results in each phase. Then, 2) systematic review conducted to identify studies and, for each of them, the detailed results of the 3) exploratory analysis and differential expression. 4) The gene meta-analysis results from the different meta-analyses. Sections 5-7) provide detailed tables and figures corresponding to the results of the three functional profiling methods (ORA, PPI, and Transcription Factor Activity). Finally, section 8) provides a synthesis of the bioinformatics methods. Through the web, the user can interact with the web tool through graphics and tables and search for specific information for a gene or function.
References
R Core Team (2022) R: A Language and Environment for Statistical Computing
Barrett T, Wilhite SE, Ledoux P, et al (2013) NCBI GEO: archive for functional genomics data sets–update. Nucleic Acids Res 41:D991-995. https://doi.org/10.1093/nar/gks1193
Moher D, Liberati A, Tetzlaff J, et al (2009) Preferred reporting items for systematic reviews and meta-analyses: the PRISMA statement. PLoS Med 6:e1000097. https://doi.org/10.1371/journal.pmed.1000097
Davis S, Meltzer PS (2007) GEOquery: a bridge between the Gene Expression Omnibus (GEO) and BioConductor. Bioinforma Oxf Engl 23:1846–1847. https://doi.org/10.1093/bioinformatics/btm254
Tweedie S, Braschi B, Gray K, et al (2021) Genenames.org: the HGNC and VGNC resources in 2021. Nucleic Acids Res 49:D939–D946. https://doi.org/10.1093/nar/gkaa980
Chen Y, Lun ATL, Smyth GK (2016) From reads to genes to pathways: differential expression analysis of RNA-Seq experiments using Rsubread and the edgeR quasi-likelihood pipeline. F1000Research 5:1438. https://doi.org/10.12688/f1000research.8987.2
Ritchie ME, Phipson B, Wu D, et al (2015) limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 43:e47. https://doi.org/10.1093/nar/gkv007
Benjamini Y, Drai D, Elmer G, et al (2001) Controlling the false discovery rate in behavior genetics research. Behav Brain Res 125:279–284. https://doi.org/10.1016/s0166-4328(01)00297-2
García-García F (2016) Methods of functional enrichment analysis in genomic studies [PhD Thesis]
Viechtbauer W (2010) Conducting Meta-Analyses in R with the metafor Package. J Stat Softw 36:. https://doi.org/10.18637/jss.v036.i03
DerSimonian R, Laird N (1986) Meta-analysis in clinical trials. Control Clin Trials 7:177–188. https://doi.org/10.1016/0197-2456(86)90046-2
Sterne JA, Egger M (2001) Funnel plots for detecting bias in meta-analysis: guidelines on choice of axis. J Clin Epidemiol 54:1046–1055. https://doi.org/10.1016/s0895-4356(01)00377-8
Ochoa D, Hercules A, Carmona M, et al (2021) Open Targets Platform: supporting systematic drug-target identification and prioritisation. Nucleic Acids Res 49:D1302–D1310. https://doi.org/10.1093/nar/gkaa1027
Boyle EI, Weng S, Gollub J, et al (2004) GO::TermFinder–open source software for accessing Gene Ontology information and finding significantly enriched Gene Ontology terms associated with a list of genes. Bioinforma Oxf Engl 20:3710–3715. https://doi.org/10.1093/bioinformatics/bth456
Wu T, Hu E, Xu S, et al (2021) clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innov Camb Mass 2:100141. https://doi.org/10.1016/j.xinn.2021.100141
Yu G, He Q-Y (2016) ReactomePA: an R/Bioconductor package for reactome pathway analysis and visualization. Mol Biosyst 12:477–479. https://doi.org/10.1039/c5mb00663e
Ashburner M, Ball CA, Blake JA, et al (2000) Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet 25:25–29. https://doi.org/10.1038/75556
Gillespie M, Jassal B, Stephan R, et al (2022) The reactome pathway knowledgebase 2022. Nucleic Acids Res 50:D687–D692. https://doi.org/10.1093/nar/gkab1028
Szklarczyk D, Gable AL, Lyon D, et al (2019) STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res 47:D607–D613. https://doi.org/10.1093/nar/gky1131
Alvarez MJ, Shen Y, Giorgi FM, et al (2016) Functional characterization of somatic mutations in cancer using network-based inference of protein activity. Nat Genet 48:838–847. https://doi.org/10.1038/ng.3593
Garcia-Alonso L, Holland CH, Ibrahim MM, et al (2019) Benchmark and integration of resources for the estimation of human transcription factor activities. Genome Res 29:1363–1375. https://doi.org/10.1101/gr.240663.118
López-Cerdán A, Andreu Z, Hidalgo MR, et al (2023) An integrated approach to identifying sex-specific genes, transcription factors, and pathways relevant to Alzheimer’s disease
Allaire JJ, Teague C, Scheidegger C, et al (2022) Quarto
Sievert C (2020) Interactive Web-Based Data Visualization with R, plotly, and shiny. Chapman and Hall/CRC