Note: The analysis that produced this target list is available to all registered AMP PD users in Terra and can be reused in your work.

Introduction

This analysis uses proteomics data described here.

Differential expression

• Cohorts for DE were created by filtering for visit month = 0 and condition = case or control.

• Tools:

  • R
  • Olink Analyze Package
  • olink_ttest

• Input:

  • TNT_CSF-PPEA-D01.csv
  • TNT_CSF-PPEA-D02.csv
  • TNT_PLA-PPEA-D01.csv
  • TNT_PLA-PPEA-D02.csv

Note: The analysis that produced this target list is available to all registered AMP PD users in Terra and can be reused in your work.

Introduction

This analysis uses proteomics data described here.

Differential expression

• Cohorts for DE were created by filtering for visit month = 0 and visit month = 48 and using the visit months as the groups for the t-test. Case and Control samples were separated into different data frames prior to the t-test.

• Tools:

  • R
  • Olink Analyze Package
  • olink_ttest

• Input:

  • TNT_CSF-PPEA-D01.csv
  • TNT_CSF-PPEA-D02.csv
  • TNT_PLA-PPEA-D01.csv
  • TNT_PLA-PPEA-D02.csv

https://github.com/tgen/ppmi-rnaseq-wb-paper

https://github.com/tgen/ppmi-rnaseq-wb-paper

https://github.com/tgen/ppmi-rnaseq-wb-paper

https://github.com/tgen/ppmi-rnaseq-wb-paper

https://github.com/tgen/ppmi-rnaseq-wb-paper

https://github.com/tgen/ppmi-rnaseq-wb-paper

Note: The analysis that produced this target list is available to all registered AMP PD users in Terra and can be reused in your work.

Introduction

The AMP-PD Target Explorer uses transcriptomic data described here.

3 Studies were included in this analysis:

Methods:

Outliers and covariates

Outlier samples

• Principal components were correlated to metadata and QC metrics. RIN and PCT USABLE BASES (Picard RnaSeqMetrics) were used to set cut offs for outlier detection.
• Tools:
  • R

Covariates

• To help identify covariates, principal components were correlated to metadata and QC metrics. variancePartition was used to validate potential covariates. Finally a neutrophil score was developed to fill the gap in CBC data.
• Tools:
  • R
  • variancePartition
• Input:
  • fullTableRNAseqMeta_TNT.csv created by: Terra Notebook
  • neutAndLymph.csv generated using: Lymphocyte and neutrophil enriched genes, from Protein Atlas Human Proteome
• Output:
  • default covariates: sex + plate + age + neutPerc + lymphPerc

Expression

Differential expression

• Cohorts for DE were created by filtering metadata for case/control status, visit month, and genetic mutation status. DE was performed with and without neutrophil score as a covariate.
• Tools:
  • R
  • edgeR
  • limma
  • Transcriptomics Differential Expression Terra Workspace
• Input:
  • matrix.featureCounts.tsv - raw count data (output from Subread featureCounts)
  • fullTableRNAseqMeta_TNT.csv - table of relevant sample metadata
• Output:
  • _limmaResults.tsv - differential expression table including fold change, p-values, and average expression
  • Plots (volcano and mean difference)

The data preprocessing and analyses were performed using R version 3.5.3 and SPSS version 25. Between-group analyses for demographic variables were conducted using two-sided one-way analysis of variance for normally distributed continuous data or Pearson's chi-square test for categorical variables. Analysis of covariance was performed when an association between classical AD CSF biomarker and age and/or sex was detected. Adjustment for multiple testing was done using the Bonferroni method. Non-Gaussian distributed data were analyzed using the Kruskal-Wallis Test. For the CSF proteome data, differences in protein abundance between pairs of clinical groups were evaluated using nested linear models. Each individual protein feature was assessed to determine if its addition to a base model containing age and gender contributed to model fit. Multiplicity was taken into account by controlling the False Discovery Rate (FDR) at q ≤ 0.05 based on the number of features analyzed. To identify the best CSF protein combination (CSF panels) that could discriminate the groups of interest while keeping the number of markers to a minimum, binary classification signatures (DLB vs. CON and DLB vs. AD) were constructed using penalized generalized linear modeling (GLM) with an elastic net penalty. Age and sex were included as covariates. The elastic net penalty, which is a linear combination of lasso and ridge penalties, enabled estimation in settings with a high feature to sample ratio. It also performed automatic feature decorrelation and feature selection. Multiple models were compared, considering different values for the elastic-net mixing parameter and the maximum number of proteins that could be selected under each model (up to 21 markers maximum). The optimal penalty parameters were determined based on balanced 10-fold cross-validation of the model likelihood. The predictive performance of all models was assessed using Receiver Operating Characteristic (ROC) curves and Area Under the ROC Curves (AUCs). The model with the highest AUC and lowest number of markers for each classification signature was selected. The fold-based selection proportions for each marker were assessed to identify and select the most promising markers within each model. To reflect the manual selection pressure for the final marker sets, each final logistic signature was subjected to ridge-regularization with a penalty parameter of 0.1. The performance (AUC) was evaluated by internal validation using repeated 5-fold cross-validation with 1000 repeats. External validations were also performed to assess the performance of the final models with the markers of interest in validation cohorts using ROC analysis. Non-parametric correlation analysis was conducted to understand the associations between the proteins within the CSF panels and the classical AD CSF biomarkers or cognitive function (MMSE score). This analysis was performed using the complete discovery cohort without stratifying per diagnostic category and conditioning on age and sex as covariates. Functional enrichment analysis was performed using Metascape, selecting GO Biological Processes as the ontology source. All the CSF proteins optimally analyzed with Olink arrays (a total of 645 protein gene products) were included as the enrichment background. Default parameters were used for the analysis, and terms with a p-value < 0.01, a minimum count of 3, and an enrichment factor > 1.5 were collected and grouped into clusters based on their membership similarities.

The following brief summaries are taken from the Materials and methods section of the original publication.

Data quality control Whole-genome sequence data were available from participants in AMP-PD cohorts. The remainder of samples were genotyped with the Illumina HumanCoreExome array (TPD), Illumina HumanCoreExome-12 v.1.1 or Illumina Infinium HumanCoreExome-24 v.1.1 arrays (OPDC) and the Illumina Infinium Multi-Ethnic Global (MEGA) array (DIGPD). Time-to-event genome-wide survival study and meta-analysis A time-to-event genome-wide survival study (GWSS) was performed in R in each cohort, using the Cox proportional hazards (CPH) function in the survival package, in which time to PDD was regressed against each single nucleotide polymorphism (SNP), with age at diagnosis, sex and first five principal components as covariates. Conditional analysis To understand whether one or more genome-wide significant variants at the same locus were contributing to the signal, we performed conditional analysis on single SNPs using a conditional and joint association analysis approach. We used the GWSS meta-analysis summary statistics and the entire AMP-PD cohort (n = 10 418) as the reference sample for linkage disequilibrium estimation. The reference sample was subjected to the same QC steps as described before. We then used CGTA-COJO software to perform association analysis conditional on SNPs of interest. Colocalization analysis To investigate whether there is an overlap between PDD loci and expression quantitative trait loci (eQTLs), we used the coloc R package. We took a Bayesian inference approach to test the H4 null hypothesis that there is a shared causal variant associated with both progression to PDD and gene expression regulation. Signal interaction between APOE and LRP1B Given the affinity of LRP1B for ApoE-carrying lipoproteins, we conducted a survival analysis based on APOE ε4 allele and LRP1B rs80306347 carrier status to understand whether the effect of LRP1B rs80306347 signal was dependent on APOE. Candidate loci analysis We additionally performed a candidate loci analysis of specific loci or variants of interest in the combined cohorts to increase power (n = 3923), using CPH models adjusted for age at diagnosis, sex, the first five principal components and a cohort term. The regions of interest consisted of genetic variants or loci previously identified in association with cognitive impairment in PD and/or dementia with Lewy bodies: APOE ε4 allele (rs429358), GBA variants E365K (or E326K, rs2230288), T408M (or T369M, rs75548401) and N409S (or N370S, rs76763715), SNCA (rs356219, rs7680557, rs7681440, rs11931074, rs7684318), MAPT H1 haplotype (rs1800547), RIMS2 (rs182987047), TMEM108 (rs138073281) and WWOX (rs8050111). In addition, participants from DIGPD and a subset of individuals from the TPD study were Sanger sequenced for GBA (n = 1793). Genetic risk scores To understand whether there is overlap in the risk of development of PDD and the risk of PD or Alzheimer’s disease, we performed a genetic risk score (GRS) analysis using PLINK v.1.9 software. Association of clinical phenotype and APOE genotype with CSF biomarkers A subset of AMP-PD participants [from the Investigation for New Discovery of Biomarkers (BioFIND) and Parkinson’s Progression Markers Initiative (PPMI) studies] included in the analysis have longitudinal CSF Alzheimer's disease biomarker data available (n = 434). We investigated the association of phenotype (PDD versus PD) and APOE ε4 carrier status with average levels of amyloid beta (Aβ) 42, total tau and tau phosphorylated at threonine 181 (p-Tau181) using unpaired two-sample Wilcoxon rank-sum tests (R stats package, v.4.1.2) at baseline, 12, 24 and 36 months of follow-up. Significance was set at α = 0.05.

“We performed the largest meta-GWAS of PD to date, involving the analysis of 7.8M SNPs in 37.7K cases, 18.6K UK Biobank proxy-cases (having a first degree relative with PD), and 1.4M controls. We carried out a meta-analysis of this GWAS data to nominate novel loci. We then evaluated heritable risk estimates and predictive models using this data. We also utilized large gene expression and methylation resources to examine possible functional consequences as well as tissue, cell type and biological pathway enrichments for the identified risk factors. Additionally we examined shared genetic risk between PD and other phenotypes of interest via genetic correlations followed by Mendelian randomization.”

Our objective was to nominate the most probable genes to be involved in PD from each GWAS locus based on the most recent PD GWAS (see Figure 1 for the study protocol). To do so, we first defined all the genes and SNPs that are within these loci (see below) and used a machine learning approach to nominate the top genes in each locus. Based on the previous literature and consensus between authors, we identified seven genes from well-established loci associated with PD that can be considered the likeliest driving genes of their respective loci (GBA1, LRRK2, SNCA, GCH1, MAPT, TMEM175, VPS13C). We then acquired data for multiple features, including different distance measures from top SNPs, different QTLs, expression in relevant tissues and cell types and predictions of variant consequences (78 features out of 284 were used after removal of redundant features, Supplementary Table 1). Using the seven well-established PD genes, which were labeled as positive, and 212 genes in the same loci that received negative labels (i.e. not likely to drive the association with PD, since the PD-driving gene is already well established), we trained a machine learning model. This model enabled us to generate a prediction score for each gene within each locus, assessing their potential involvement in PD. The gene with the highest score in each locus is the nominated gene to be associated with PD. We then performed multiple post hoc analyses to further validate and explore our results: burden tests for rare variants in the top-scoring genes, pathway enrichment and pathway PRS analyses, differential expression analyses and structural analyses for candidate coding variants.

To assess PD risk, summary statistics from Chang et al. PD GWAS meta-analysis involving 26,035 PD cases and 403,190 controls of European ancestry were used as the reference dataset for the primary analysis to define risk allele weights. Individual-level genotyping data not included in Chang et al. [7] and from the last GWAS meta-analysis [17] was then randomly divided as the training and testing datasets. The training dataset used to construct the PRS consisted of 7218 PD cases and 9424 controls, while the testing dataset to validate the results consisted of 5429 PD cases and 5814 controls, all of European ancestry. A polygenic effect score (PES) was generated to estimate polygenic risk for each of the 2199 gene sets representative of biological pathways and then tested for association with PD. PES was calculated based on the weighted allele dose as implemented in PRSice2 (v2.1.1). The sequence kernel association test-optimal (SKAT-O) was implemented using default parameters in RVTESTS [35] to determine the difference in the aggregate burden of rare coding genetic variants (minor allele count≥ 3) between PD cases and controls for the nominated gene-sets by PRS. ANNOVAR was used for variant annotation. Baseline peri-diagnostic RNA sequencing data derived from the blood for 1612 PD patients and 1042 healthy subjects available from the Parkinson Progression Marker Initiative (PPMI) was used to construct a network of expression communities based on a graph model with Louvain clusters. This cleaned and normalized data was downloaded from the Accelerating Medicines Partnership for Parkinson’s disease (AMP-PD) on March 1st, 2020. Scikit-learn’s extraTreeClassifer option was used to extract coding gene features for inclusion in the network builds that are likely to contribute to classifying cases versus controls under default settings in the feature selection phase, leaving 8.3 k protein-coding genes for candidate networks. Following this feature extraction phase, controls were excluded, and case-only correlations were calculated for all remaining gene features. Next, this correlation structure was converted to a graph object using NetworkX. Subsequently, the Louvain algorithm was employed to build network communities within this graph object derived from the selected feature set. Finally, pathway enrichment analysis within expression communities was performed to further dissect its biological function using the function g:GOSt from g:ProfleR. The significance of each pathway was tested by hypergeometric tests with Bonferroni correction to calculate the error rate of each network. Single-cell RNA sequencing data [25] based on a total of 9970 cells obtained from several mouse brain regions (neocortex, hippocampus, hypothalamus, striatum, and midbrain) was used to explore cell types associated with PD risk. Linear regression adjusted by the number of SNPs included in the PRS was performed to assess the trend of increased PRS R2 per decile of cell-type expression specificity. Two-sample SMR was applied to explore the enrichment of cis eQTLs within the 46 gene-sets nominated by our large scale PRS analysis. The methodology can be interpreted as an analysis to test if the effect size of genetic variants influencing PD risk is mediated by gene expression or methylation to prioritize genes underlying these gene-sets for follow-up functional studies. Additionally, we studied expression patterns in blood from the largest eQTL meta-analysis so far. The number of genes tested per gene-set were Bonferroni corrected, and a Chisq test was applied to assess whether the proportion of QTLs per geneset was significantly higher than expected by chance.

In an effort to prioritize the top genes within significant gene-sets showing the highest cumulative effect on PD risk, individual gene-based SKAT-O analyses were performed considering a MAF threshold ≤3% and three functional categories (missense, loss of function and CADD score >12). Using this approach, gene-level prioritization is highlighted in Supplementary Table 7.

To assess PD risk, summary statistics from Chang et al. PD GWAS meta-analysis involving 26,035 PD cases and 403,190 controls of European ancestry were used as the reference dataset for the primary analysis to define risk allele weights. Individual-level genotyping data not included in Chang et al. [7] and from the last GWAS meta-analysis [17] was then randomly divided as the training and testing datasets. The training dataset used to construct the PRS consisted of 7218 PD cases and 9424 controls, while the testing dataset to validate the results consisted of 5429 PD cases and 5814 controls, all of European ancestry. A polygenic effect score (PES) was generated to estimate polygenic risk for each of the 2199 gene sets representative of biological pathways and then tested for association with PD. PES was calculated based on the weighted allele dose as implemented in PRSice2 (v2.1.1). The sequence kernel association test-optimal (SKAT-O) was implemented using default parameters in RVTESTS [35] to determine the difference in the aggregate burden of rare coding genetic variants (minor allele count≥ 3) between PD cases and controls for the nominated gene-sets by PRS. ANNOVAR was used for variant annotation. Baseline peri-diagnostic RNA sequencing data derived from the blood for 1612 PD patients and 1042 healthy subjects available from the Parkinson Progression Marker Initiative (PPMI) was used to construct a network of expression communities based on a graph model with Louvain clusters. This cleaned and normalized data was downloaded from the Accelerating Medicines Partnership for Parkinson’s disease (AMP-PD) on March 1st, 2020. Scikit-learn’s extraTreeClassifer option was used to extract coding gene features for inclusion in the network builds that are likely to contribute to classifying cases versus controls under default settings in the feature selection phase, leaving 8.3 k protein-coding genes for candidate networks. Following this feature extraction phase, controls were excluded, and case-only correlations were calculated for all remaining gene features. Next, this correlation structure was converted to a graph object using NetworkX. Subsequently, the Louvain algorithm was employed to build network communities within this graph object derived from the selected feature set. Finally, pathway enrichment analysis within expression communities was performed to further dissect its biological function using the function g:GOSt from g:ProfleR. The significance of each pathway was tested by hypergeometric tests with Bonferroni correction to calculate the error rate of each network. Single-cell RNA sequencing data [25] based on a total of 9970 cells obtained from several mouse brain regions (neocortex, hippocampus, hypothalamus, striatum, and midbrain) was used to explore cell types associated with PD risk. Linear regression adjusted by the number of SNPs included in the PRS was performed to assess the trend of increased PRS R2 per decile of cell-type expression specificity. Two-sample SMR was applied to explore the enrichment of cis eQTLs within the 46 gene-sets nominated by our large scale PRS analysis. The methodology can be interpreted as an analysis to test if the effect size of genetic variants influencing PD risk is mediated by gene expression or methylation to prioritize genes underlying these gene-sets for follow-up functional studies. Additionally, we studied expression patterns in blood from the largest eQTL meta-analysis so far. The number of genes tested per gene-set were Bonferroni corrected, and a Chisq test was applied to assess whether the proportion of QTLs per geneset was significantly higher than expected by chance.

SMR revealed functional genomic associations with eQTLs in 201 genes (Supplementary Table 12, online resource) of which 88 were found to be part of the network communities significantly associated with PD in our transcriptome community map (Supplementary Table 13, online resource).