We obtained, through NCBI Gene Expression Omnibus (GEO), two single-cell RNA-seq datasets that we employed to derive gene expression signatures for the major brain cell types (i.e., astrocytes, microglia, neurons and oligodendrocytes), one from human temporal lobe [GSE67835 (Darmanis et al., 2015) – only cells from adult samples were used] and the other from mouse cortex [SRP135960 (Zeisel et al., 2015)]. We used a third single-cell RNA-seq dataset [GEO GSE73721 (Zhang et al., 2016)] to independently validate those signatures, considering only cells from the human cortex (12 astrocytes, 1 neuron, 4 oligodendrocytes, and 2 endothelial cells), for consistency.

The AD analysis was based on the temporal cortex RNA-seq dataset from the AMP-AD Knowledge Portal with accession syn3163039 (Allen et al., 2016). We used the table available therein, containing the pre-processed raw read counts for each gene in each sample, for the downstream analyses. We selected only the samples diagnosed as AD and non-AD with RNA integrity number (RIN) ≥8 (Ferrer et al., 2008). We also discarded a non-AD sample with a very low (<0.40) estimated proportion of neurons (Supplementary Figure 1A and Supplementary Table 1). In total, we used 71 AD and 32 non-AD samples.

We used AD dataset GEO GSE104704 (Nativio et al., 2018) for independent validation, less stringently requiring RIN ≥6 to keep enough samples for analysis. Three non-AD samples with abnormally low (<0.40) estimated proportion of neurons were discarded (Supplementary Figure 1B and Supplementary Table 1), leaving a total of 9 AD and 14 non-AD samples.

We fetched the PD RNA-seq dataset from GEO GSE68719 (Dumitriu et al., 2016) and kept samples with RIN ≥7 and from donors older than 60 years, for a better age match between control and diseased samples and given the reported onset of idiopathic PD at around 65 years of age (Nussbaum and Ellis, 2003). When performing principal component analysis (PCA) of the normalized gene expression data (see sections “Statistical Tests” and “Data Processing” below), we identified two samples (SRR2015728 and SRR2015748) with an outlying behavior (Supplementary Figures 2A,B). When clustering samples based on the correlation between their normalized gene expression profiles (see section “Statistical Tests” below), SRR2015728 and SRR2015748 are again shown to be outliers (Supplementary Figure 2C). Moreover, non-PD samples SRR2015714 and SRR2015728 are also those showing an abnormally low (<0.40) estimated neuronal proportion (Supplementary Figure 1C and Supplementary Table 1). As such, we conservatively discarded those three samples from the dataset, leaving 15 PD and 26 non-PD samples.

For independent validation, we used PD gene expression microarray dataset GEO GSE20168 (Zhang et al., 2005). Since the PD RNA-seq dataset only comprised males, we selected the 10 non-diseased and 8 PD male samples from the microarray dataset. Although RINs were not provided for this dataset, we were able to detect possible RNA degradation by using function AffyRNAdeg from the xps R package (Stratowa, 2020). We found some evidence for the expected neuronal loss in PD brains but differences in the distributions of estimated neuronal proportions between PD and non-PD samples are not statistically significant (Supplementary Figure 1D). Although one non-PD sample had a low (<0.40) estimated neuronal proportion, we decided to keep it due to the small number of samples in the dataset (Supplementary Table S1).

All datasets used are summarized in Table 1.

We performed all statistical analyses in R (programming language for statistics and graphics) (R Development Core Team, 2018), extensively using packages from Bioconductor (repository of R tools for the analysis of high-throughput biological data) (Gentleman et al., 2004). We used t-tests (Kalpić et al., 2011) to compare differences in expression of specific marker genes, as well as differences in age distributions between diseased and non-diseased groups. To compare differences in proportions of neural cell types between diseased and non-diseased brains, we used Wilcoxon-signed-rank tests (Rey and Neuhäuser, 2011), and to compare the neuronal proportion densities between diseased and non-diseased brains, we used the Kolmogorov-Smirnov test (Lopes, 2011). For correlation analysis, we used Pearson’s correlation, unless stated otherwise. We chose Euclidean distance for clustering samples based on gene expression correlation, having used the ComplexHeatmap package (Gu et al., 2016) for the purpose and to generate the associated heatmap in Supplementary Figure 2C.

Principal component analysis (PCA), enabling the identification of the linear combinations of variables that contribute the most to data variance (Ringnér, 2008), was implemented through the singular value decomposition (SVD) algorithm provided by the PCA function from R package FactoMineR (Lê et al., 2008).

Where applicable and not indicated otherwise, p-values were corrected for multiple testing using Benjamin-Hochberg’s False Discovery Rate (FDR).

For all the RNA-Seq datasets with no pre-processed data available (Table 1), we aligned the reads against the human transcriptome [hg38 Gencode annotation (Frankish et al., 2019)] with Kallisto (Bray et al., 2016) using the default parameters.

For both single-cell datasets, we performed state-of-the-art procedures for quality assessment (Lun et al., 2016b), such as checking for library size discrepancies between cells, the number of expressed genes per cell and the proportion of reads aligning to mitochondrial genes (Lun et al., 2016b; Supplementary Figures 3A,C). We removed low-quality cells that presented a median absolute deviation (MAD) <−3 for the library size, MAD < −3 for the number of expressed genes or a MAD > 3 for the proportion of mitochondrial reads. Additionally, we kept for downstream analysis only genes whose log₁₀[average read counts per million (CPM)] > 0 (Supplementary Figures 3B,D) and whose variance in expression was significantly associated with the biological component (i.e., the cell type) as assessed through the usage of the decomposeVar function from the scran R package (Lun et al., 2016b). Briefly, the variance in expression for each gene was decomposed into their biological and technical components. The technical component is estimated by fitting the mean-dependent trend of the variance. The biological component of the variance is then calculated by subtracting the technical component from the overall variance (Lun et al., 2016b). This last step avoids prioritizing genes whose expression is highly variable due to technical factors such as sampling noise during RNA capture and library preparation (Lun et al., 2016b).

Furthermore, the t-Distributed Stochastic Neighbor Embedding (tSNE) (van der Maaten and Hinton, 2008) plot of human single-cell (Darmanis) gene expression shows a few cells not clustered together with those of their respective annotated type (Supplementary Figure 4A). Moreover, all of them appear to have been misclassified also based on single-cell trajectories [i.e., cells’ ordering according to their inferred biological state (Qiu et al., 2017b)] obtained with the monocle package (Trapnell et al., 2014; Qiu et al., 2017a, b; Supplementary Figure 4B) or the nearest shrunken centroid classification, implemented in R package pamr (Tibshirani et al., 2002; Supplementary Figures 4C,D). Therefore, they were discarded from our analysis (Supplementary Figure 4E). No potentially misclassified cells were detected in the mouse dataset (Supplementary Figure 4F).

After the filtering steps mentioned above, summed expression values across pools of cells were deconvolved in cell-based factors for normalization of the Darmanis and the mouse single-cell gene expression datasets (Lun et al., 2016a). All bulk RNA-seq datasets were quantile-normalized using the voom function, from the limma R package (Ritchie et al., 2015). The rma function from the affy R package (Gautier et al., 2004) was used to normalize and summarize the PD microarray dataset.

Moreover, we used the ComBat function from the sva package to correct for batch effects. This function requires possible technical effects to be encoded as categorical variables (Leek et al., 2012). Thus, for the AD MayoClinic dataset, RIN was defined as high if >8.5 and low if ≤8.5, in the AD Nativio dataset high if >7.3 and low if ≤7.3, and in the PD Dumitriu dataset it was defined as high if >7.8 and low if ≤7.8. For the PD Zhang dataset, the RNA degradation slope, derived from the average intensities per relative 5′–3′ position of probes in their target transcripts across probe sets (Draghici, 2012), was used and defined as low if ≤5 and high if >5.

We quantified gene expression from RNA-seq data in counts per million (CPM) and kept only genes with an average CPM higher of 10/L, where L is the minimum library size in million reads (Chen et al., 2016), in at least N samples, where N is the smallest sample size in our analyses (Supplementary Figure 5). For the microarray dataset, gene expression was quantified by normalized intensities.

We employed CIBERSORTx (Newman et al., 2019) to infer both human and mouse gene expression signatures for each of the major brain cell types (Supplementary Tables 2, 3) and subsequently used them to estimate the cellular composition of brain samples from their bulk transcriptomes.

We followed three different approaches to assess the accuracy of human and mouse CIBERSORTx-derived signatures in correctly identifying the major cell types in human brain samples:

We used CIBERSORTx deconvolution (Newman et al., 2019), relying on human and mouse gene expression signatures for the major brain cell types derived as described above, to estimate the cellular composition of all AD, PD, and non-diseased brain samples from their bulk transcriptomes. Moreover, as CIBERSORTx options, we enabled batch normalization, disabled quantile normalization and used 100 permutations for significance analysis. Following CIBERSORTx’s user guidelines, the B-mode batch normalization was chosen to perform deconvolution when using the human signature, since the single cell data used to derive it were generated with SMART-seq2 (Picelli et al., 2013), and the S-mode batch normalization when using the mouse signature, since it is tailored for signatures derived from data generated with the 10x Genomics Chromium platform, as was the case (Newman et al., 2019).

We performed differential gene expression using the limma (Ritchie et al., 2015) and edgeR (Dai et al., 2014) packages.

For each coefficient in the linear model, the magnitude of differences in gene expression was measured in log₂ fold-change and their significance was given by the FDR-adjusted p-value of the moderated t-statistic (an ordinary t-statistic with its standard errors moderated across genes), along with the empirical Bayes statistic (B statistic - log-odds ratio of a gene being differentially expressed) (Smyth, 2004). Moreover, we also used the moderated t-statistic to assess the differential gene expression coherence between different datasets.

We linearly modeled gene expression in the AD datasets according to the following:

Here GE_x is the expression of gene x; Disease is the sample’s centered disease status; RIN is the categorized sample’s RNA Integrity Number (1 for high and 0 for low); Neuronal proportion is given by the sample’s estimated proportion of neurons centered; Age is the age of the sample’s donor in years; Sex is the biological sex of the sample’s donor (1 for male and 0 for female); Interaction is the interaction between Disease and the Neuronal proportion effects, given by the product of the two and interpretable as the differential effect of the loss of neurons between AD and non-diseased samples or, equivalently, the part of AD effect that is dependent of the sample’s neuronal contents; βs are the unknown coefficients, to be estimated from fitting that linear model to the gene expression data, for each of the aforementioned variables hypothesized to impact gene expression; ε states the error of the model, that is the remaining variance not explained by the model. Disease and Neuronal proportion were centered to diminish the correlation between their associated estimated coefficients, thereby using a model more consistent with the purpose of estimating independent effects (Afshartous and Preston, 2011). We thus shifted the “prediction center” (i.e., the virtual reference) to the average sample (Afshartous and Preston, 2011) by turning the variables’ means to 0 through the usage of the scale function from the built-in R package base (R Development Core Team, 2018), with the scale argument turned to “false.”

Likewise, we modeled gene expression in the PD Dumitriu dataset as following:

Unknown batch corresponds to a batch effect of unknown source detected by PCA (Supplementary Figure S8) that was thereby adjusted for (Supplementary Figure 12C).

For validation with the independent PD microarray dataset (Zhang), we used the following linear model:

RNA degradation is given by its slope grouping for each sample (1 for high and 0 for low).

We considered a gene differentially expressed if FDR < 0.05, except for the Zhang PD microarray dataset, where we considered FDR < 0.11. This arbitrary cut-off was used to “rescue” a reasonable number genes for further analyses, given the small sample size of the Zhang dataset and the consequent lower statistical power of the associated differential expression analysis. This arbitrary looseness in specificity is dealt with by subsequent filtering (v. section on permutation analyses below).

Genes already reported to play a role in AD and PD were gathered from the DisGeNET database (Piñero et al., 2019). Only genes with a human gene-disease association (GDA) score >0.1 and an evidence index ≥0.9 (180 genes for AD and 112 genes for PD) were considered as such in our analyses.

We performed permutation tests to identify genes with consistent differential expression ranking between datasets. For each gene, we multiplied its t-statistic values for the intrinsic disease effect (Disease in the linear models) in each of the two datasets (MayoClinic and Nativio for AD; Dumitriu and Zhang for PD) and compared that product with the distribution of those resulting from 5000 random permutations of the disease status labeling of samples. The proportion of random products more extreme than the empirical one was taken as its False Discovery Rate (FDR).

To assess the similarity between the intrinsic AD and PD Disease effects on gene expression, we compared the aforementioned FDRs. For each disease, when t-statistics for both datasets were positive, we used −log₁₀(FDR), when both negative, we used log10(FDR), and when contradictory (i.e., showing different signs) we set this value to 0. When FDRs were originally zero, we equaled them to 1e-5 (half of the FDR resolution) to avoid infinite values when computing their logarithms. Then, for each gene, we multiplied those scores of AD and PD and compared this product with the distribution yielded by 1 000 000 permutations of randomly shuffled product scores. The proportion of random products more extreme than the empirical one was taken as its FDR (Figure 7A).

We identified KEGG (Kanehisa et al., 2016) pathways dysregulated in AD and PD datasets using the Piano R package (Väremo et al., 2013) to perform gene set enrichment analysis (GSEA) (Mootha et al., 2003; Subramanian et al., 2005), by default on t-statistics, but also on B-statistics of differential gene expression for the AD Disease and Neuronal proportion effects. We also used the AD cell-type marker genes defined by Kelley et al. (2018) as a gene set. For GSEA on genes commonly changed in AD and PD (Figure 7B), we used −log₁₀(FDR) when both AD and PD scores were positive, log₁₀(FDR) when both negative, and zero when signs were contradictory.

We used cTRAP (de Almeida et al., 2020) to compare the changes in gene expression induced by thousands of drugs in human cell lines, available in the Connectivity Map (CMap) (Subramanian et al., 2017), with those in human brains that we have inferred to be related to the intrinsic (i.e., systemic) AD and PD effects. As input for cTRAP, we used the aforementioned scores for the Disease and Neuronal proportion effects, thereby ranking changes that are coherent between the MayoClinic and Nativio datasets for AD and the Dumitriu and Zhang datasets for PD, as well as those coherent between AD and PD. The compounds, in clinical trials or launched, with their perturbation z-scores (Subramanian et al., 2017) exhibiting the 20 most negative and the 20 most positive average (across different cell lines) Spearman’s correlation with the Disease effect scores across common genes, and with an average absolute Spearman’s correlation with the Neuronal proportion effect scores <0.05 (to avoid confounding between effects), were selected for AD (Supplementary Figure 15A) and PD (Supplementary Figure 15B) as the top candidates for reversal or induction of disease-associated gene expression alterations for discussion. Noteworthily, cMap includes data for the same compounds tested with different concentrations and at different time points, as well as run in different plate types (ASG, CPD, HOG, etc.) (Supplementary Tables 16–18).

Most neuronal markers [DCX (Tanapat, 2013), MAP2 (Tanapat, 2013), NEFM (Tanapat, 2013), NEFH (Tanapat, 2013), NEFL (Tanapat, 2013), RBFOX3 (Tanapat, 2013), SYP (Tanapat, 2013)] are significantly downregulated in AD temporal cortex samples from the MayoClinic dataset (Figure 1A). In contrast, all astrocytic [ALDH1L1 (Preston et al., 2019), GFAP (Preston et al., 2019)), SLC1A3 (Preston et al., 2019)] and a few microglial and oligodendrocytic markers [CD40 (Ponomarev et al., 2006), OLIG1 (Ligon et al., 2004) and OLIG2 (Ligon et al., 2004)] are significantly upregulated in AD brains.

CIBERSORTx (Newman et al., 2019), a tool that estimates cell type abundances in tissues from their bulk transcriptomes and machine learning-inferred cell-type-specific gene expression profiles, was used to derive the composition in major cell types (astrocytes, microglia, neurons and oligodendrocytes) of AD brain samples (see section “Materials and Methods”). These estimates (Figure 1B) are concordant with the observations in Figure 1A, including significant increase and decrease, respectively, in the proportions of astrocytes and neurons in AD brain samples. Despite the known differences in gene expression between mouse and human brain cells (Hodge et al., 2019), the same trends can be seen using the mouse signature (Supplementary Figures 9A,B). We also performed principal component analysis (PCA) on normalized gene expression in the MayoClinic brain samples. The neuronal composition, along with the disease effect, is correlated with the first principal component (PC1), i.e., that retaining the most data variance (Figure 1C – rho = −0.88; p < 2.2e-16). Sex also shows a strong association with PC1 (Supplementary Figure 9E) but there is no significant difference in age or neuronal proportion between female and male individuals (Supplementary Figure 9F).

We linearly modeled gene expression in the MayoClinic brain samples as a function of technical (RIN) and biological variables, such as neuronal proportion (reflecting neuronal loss), systemic AD, Age, Sex and interaction between neuronal proportion and AD (Supplementary Table 6, Figure 2, and Supplementary Figure 10A). We were thereby able to discriminate genes whose expression is significantly systemically affected by AD (Figures 2A,B) from those essentially showing a strong association with neuronal loss (Figures 2C,D). For instance, LIAS, CTB-171A8.1, COX18, and ETV4 exemplify genes that show a strong intrinsic AD effect, independent of neuronal proportion (Figure 2B). Moreover, genes that have previously been reported as playing a role in AD in the DisGeNet database, namely CDK5 (Liu et al., 2016), CDK5R1 (Moncini et al., 2017), FERMT2 (Shulman et al., 2014), and HSD17B10 (Marques et al., 2009), were actually found to be associated with neuronal loss rather than with the disease component (Figure 2D). Additionally, with the Interaction effect we were also able to detect genes, such as PNPLA5 and PTPN20A (Supplementary Figure 10C), whose expression was differentially altered with cellular composition in AD brains compared with non-AD samples. It is worth noting that AD genes specific of a cell type, as defined by Kelley et al. (2018), are also more related with the neuronal composition effect (p_GSEA = 0.0001) than with the disease effect (p_GSEA = 0.6) (Figures 2A,C). Additionally, most up-regulated AD-specific genes seem to be related with cell survival and immune pathways, whereas down-regulated ones with oxidative phosphorylation and Parkinson’s disease (Supplementary Figure 10C).

In order to validate those results, we used the independent AD RNA-seq dataset (lateral temporal lobe) (Nativio et al., 2018), herein named Nativio (Table 1), that we found to better match the larger MayoClinic dataset (temporal cortex) (Allen et al., 2016) in terms of brain area. Although the Nativio dataset did not present significant differences in cellular composition between AD and non-AD samples (Supplementary Figure 11A), its samples were from significantly younger donors (p = 2.7e-10) and the neuronal proportion of its AD samples is significantly different (p = 0.033) from MayoClinic AD samples’ (Supplementary Figure 11C), we found consistency in AD-associated gene expression changes between the MayoClinic and the Nativio datasets (Figure 3 and Supplementary Tables 7–9).

We analyzed the PD datasets following similar approaches to those used on AD transcriptomes. In the RNA-seq PD dataset (Dumitriu), only two neural (DCX and MAP2) and two microglial (CD40 and ITGAM) markers showed significant alterations between PD and non-PD samples (Figure 4A), concordantly with no significant differences in cellular composition as estimated by CIBERSORTx (Figure 4B) using both the human and the mouse signatures (Supplementary Figures 12A,B). However, the neuronal composition of samples, along with the disease effect, has a significant association with PC1 (Figure 4C – rho = −0.66; p = 3.2e-6) of gene expression. Age also shows a relationship with PC1 (Supplementary Figures 12C,E) but no significant age difference exists between PD and non-PD sample donors (Supplementary Figure 12F).

Gene expression was modeled for each gene as a function of the technical variables (RIN and unknown confounder), neuronal proportion (i.e., neuronal loss), intrinsic PD and Age (Figure 5A, Supplementary Figure 13A and Supplementary Table 10). No PD-neuronal loss interaction effect was considered because no significant differences in cell type proportions between PD and non-PD samples were detected (Figure 4B). We were thereby able to discriminate genes with a strong disease effect (Figure 5B) from those essentially altered by neuronal loss (Figure 5D). In fact, according to our analysis, genes reported as playing a role in PD [ABL1 (Mahul-Mellier et al., 2014), COMT (Jiménez-Jiménez et al., 2014), GRK5 (Liu et al., 2010), and APT1A3 (Haq et al., 2019)] were found associated with neuronal loss rather than the disease itself (Figures 5A,C).

Although RNA-seq provides a more precise quantification of gene expression than microarrays (Wang et al., 2009), we could not find any other independent PD RNA-seq dataset matching, in terms of brain region, the Dumitriu study and therefore resorted to the Zhang microarray study. This independent dataset did not present any significant cellular composition alteration between PD and non-PD samples (Supplementary Figure 14A) either. Additionally, we found no significant differences in neuronal proportion estimates or age between samples from the Dumitriu and Zhang studies (Supplementary Figure 14C). We find gene expression changes that are consistent between the two analyzed PD datasets (Figures 6A,B and Supplementary Tables 10–13), including for LRRC40 and ABCB6 (Figure 6C). However, from the selected PD candidates as examples shown in Figure 5, only ADAMTS2 and ADCYAP1 were profiled in the Zhang dataset and did not recapitulate the changes observed in the Dumitriu dataset (Figure 6B).

Alzheimer’s diseases and PD-associated gene expression changes in human brains are very correlated (Figure 7A and Supplementary Tables 14, 15), suggesting commonalities in the molecular mechanisms underlying both diseases. Although some neuronal markers are amongst the genes commonly altered by AD and PD, most of them are not, indicating effectivity in decoupling the neuronal loss effect on gene expression (Figure 7A). Genes consistently up-regulated in both diseases are linked to Wnt signaling (basal cell carcinoma pathway) and NF-KB signaling (acute myeloid leukemia) (Figure 7B). Indeed, the genes driving the basal cell carcinoma pathway are FZD9, FZD7, FZD2, DVL1, and AXIN1, all playing a role in the Wnt signaling pathway, which has already been linked to AD and PD (Harvey and Marchetti, 2014). The genes contributing the most to the acute myeloid leukemia pathway (RAF1, RELA, and IKBKB) are related with NF-KB signaling, a process already known to also play a role in AD and PD (Mattson and Meffert, 2006). Genes consistently down-regulated in both diseases are linked essentially to cell metabolism (Figure 7B). Although the magnitude of disease-induced changes in gene expression is generally modest (as expected, as samples from the same type of tissue are being compared), reassuringly they are overall quite independent from the neuronal composition of the brain samples (Figure 7C).

We used cTRAP (de Almeida et al., 2020) to identify drugs that, when delivered to human cell lines, cause similar (correlated) or opposite (anti-correlated) gene expression changes to those we observed as intrinsically associated to AD and PD (Figures 8A,C). Interestingly, gene expression perturbations induced by drugs known to be used in the clinic to treat AD (Siavelis et al., 2016) (donepezil, galantamine, memantine, and rivastigmine) and PD (Zahoor et al., 2018) (amantadine, bromocriptine, cabergoline, carbidopa, entacapone, levodopa, lisuride, and selegiline) were not amongst the most correlated with those by the respective target diseases (Figures 8A,B and Supplementary Tables 16–18). Siavelis et al. (2016) had followed a similar approach, although they did not decouple the neuronal loss effect, having identified 27 drugs linked to the AD phenotype (Figure 8A). Chloroquine and scriptaid seem promising drug candidates for AD since their known targets are indeed overexpressed in AD and vary very little with neuronal loss (Figure 8B). Scriptaid also seems promising for PD therapeutics for similar reasons (Figure 8D). Additionally, gene expression changes upon metaraminol administration showed up as being the most anti-correlated with those commonly induced by AD and PD (Figure 9), being metaraminol therefore a potential candidate drug to be tested for repurposing. Gene expression changes upon wortmannin administration are, in a dose dependent manner, the most correlated with those commonly induced by AD and PD (Figure 9).