GWAS

The principal investigators of all cohorts obtained informed consent from all study participants and approval from Institutional Review Boards (IRB) at their respective institution. We obtained GWAS summary statistics on EA from the SSGAC. The results are based on Okbay et al.8, including the UK Biobank. The PGC shared GWAS summary statistics on SZ with us that were reported in Ripke et al.5, but excluded data from our replication sample (GRAS), yielding a total sample size of n = 34,409 cases and n = 45,670 controls.

All cohorts that were part of both studies5,8 were excluded from the meta-analysis on EA, yielding non-overlapping GWAS samples and n EA = 363,502. The original EA results file contained 12,299,530 genetic markers, compared to 17,221,718 in the SZ results file.

We applied the following additional quality control steps:

1. To maximize statistical power, we excluded SNPs that were missing in large parts of the two samples. Specifically, we continued with SNPs that were available in at least 19 out of 50 cohorts in the SZ results5 (the actual N per SNP was not provided in the SZ GWAS summary statistics) and in N > 200,000 in the EA meta-analysis8. This step excluded 3,778,914 and 6,369,138 genetic markers for EA and SZ, respectively. 2. We dropped SNPs that were not available in both GWAS results files. This step restricted our analyses to the set of available genetic markers that passed the quality-control filters in both the EA and the SZ GWAS results, leaving us with 8,403,560 autosomal SNPs. 3. We dropped six SNPs with non-standard alleles (i.e. not A, C, T or G) and two SNPs with mismatched effective alleles. Furthermore, we dropped 163,272 SNPs in the first and the 99th percentile of the distribution of differences in MAF in the two results files. This final step eliminated SNPs that were likely to be affected by coding errors, strand flips or substantial differences in MAF in the EA and SZ samples.

The remaining 8,240,280 autosomal SNPs were used in the proxy-phenotype and prediction analyses.

Proxy-phenotype method

Look-up: We conducted our proxy-phenotype analyses following a pre-registered analysis plan (https://osf.io/dnhfk/), using the 8,240,280 autosomal SNPs that passed quality control. We selected 10−5 as the default P value threshold to identify EA-associated SNPs prior to carrying out the proxy-phenotype analyses (Supplementary Note 2).

To select approximately independent SNPs from the EA GWAS results, we applied the clumping procedure in PLINK version 1.944,45 using r2 > 0.1 and 1,000,000 kb as the clumping parameters and the 1000 Genomes phase 1 version 3 European reference panel46 to estimate LD among SNPs. This algorithm assigns the SNP with the smallest P value as the lead SNP in its ‘clump’. All SNPs in the vicinity of 1,000,000 kb around the lead SNP that are correlated with it at r2 > 0.1 are assigned to this clump. The next clump is formed around the SNP with the next smallest P value, consisting of SNPs that have not been already assigned to the first clump. This process is iterated until no SNPs remain with P < 10−5, leading to 506 approximately independent EA-associated lead SNPs. 108 of the 506 EA-associated lead SNPs are genome-wide significant (P < 5 × 10−8).

We looked up the SZ GWAS results for these 506 EA-associated lead SNPs. Results for all 506 SNPs are reported in Supplementary Data 2 and Fig. 2.

In order to investigate the novelty of the findings, we extracted all the SNPs in LD with these 21 SNPs at r2 ≥ 0.1 with a maximum distance of 1000 kb using the 1000 Genomes phase 1 European reference panel.

Bayesian credibility of results: We probed the credibility of our proxy-phenotype association results using a heuristic Bayesian calculation following Rietveld et al. (Supplementary Information pp. 13–15)25. We focus on the 21 EA-associated lead SNPs that are also associated with SZ after Bonferroni correction.

Bayes’ rule implies that the probability that an association is true given that we observe significance is given by

$$\begin{array}{*{20}{l}} {P\left( {H_1|t > t_{\alpha /2}} \right)} \hfill & = \hfill & {\frac{{P\left( {t > t_{\alpha /2}|H_1} \right)P\left( {H_1} \right)}}{{P\left( {t > t_{\alpha /2}|H_1} \right)P\left( {H_1} \right) + P\left( {t > t_{\alpha /2}|H_0} \right)P\left( {H_0} \right)}}} \hfill \\ {} \hfill & = \hfill & {\frac{{({\mathrm{power}})(\pi )}}{{\left( {\mathrm{power}} \right)\left( \pi \right) + (\alpha )(1 - \pi )}}} \hfill \end{array}$$

‘Power’, as well as the significance test, are two-sided, π is the prior belief that the SNP is truly associated and α is the significance threshold used for testing (in our case, α = \(\frac{{0.05}}{{506}}\) = 9.88 × 10−5).

To calculate power for each SNP, we computed the winner’s curse corrected OR using the procedure described in Rietveld et al. (Supplementary Information pp. 7–13)47 for the α threshold of 9.88 × 10−5. Because the actual sample size per SNP is not reported in the SZ GWAS summary statistics, we furthermore assumed that each SNP was available in the entire sample of 34,409 cases and 45,670 controls (i.e. the PGC results from Ripke et al.5 excluding the GRAS data collection).

An important question is which prior beliefs are reasonable starting points for these Bayesian calculations. For an arbitrarily chosen SNP, the most conservative reasonable prior would assume that each truly associated SNP has the same effect size as the strongest effect size that was actually observed in the data. If one divides the SNP-based heritability of the trait by that effect size in R2 units, one obtains a lower bound for the number of SNPs that can be assumed to be truly associated. To aid this line of thinking, we converted the winner’s curse corrected OR of our 21 SNPs into R2 using

$$R^2 = \left( {\frac{d}{{\sqrt {d^2 + a} }}} \right)^2$$

where d is Cohen’s d, which is calculated as

$$d = {\mathrm{ln}}(\rm{Odds})\frac{{\sqrt 3 }}{\pi }$$

and a is a correction factor that adjusts for the MAF of the SNP. This correction factor is calculated as

$$a = \frac{{\left( {n_1 + n_2} \right)^2}}{{n_1n_2}}$$

where n 1 = N × MAF and n 2 = N × (1 − MAF)48.

The largest effect size in R2 that we observe in our results is rs4378243 with 0.044%. The SNP-based heritability of SZ is ≈21%49. Thus, if all causal SZ SNPs would have an effect of R2 = 0.044%, we would expect that ≈500 truly causal loci exist. The chance of finding any one of them by chance from a set of ≈500,000 independent loci in the human genome is ≈0.1%. (Our pruning algorithm of SNPs that passed QC leads to only 223,065 independent loci. Thus, assuming 500,000 independent loci in these calculations is conservative.) However, in reality most truly associated loci for SZ will surely have smaller effects than that. Thus, a prior belief of ≈0.1% is certainly too conservative.

Furthermore, the SNPs we investigate are not arbitrary but selected based on their association with another, genetically related cognitive trait (EA) in a very large, independent sample. Thus, a prior belief of 1 or 5% that these SNPs are also associated with SZ is probably more reasonable. As an upper bound, we assume that 10% of all loci are causal. Thus, the chance to pick any one of them by chance would be 10%.

Table 1 displays the winner’s curse corrected effect size of the 21 EA-associated lead SNPs that are also associated with SZ after Bonferroni correction. It also shows the posterior probability that these SNPs are truly associated with SZ given our results for prior beliefs ranging from 0.1, 1, 5 to 10%. Thirteen of these SNPs have posterior probabilities of being true positives of >50% for even the most conservative prior. For a more realistic prior belief of 5%, all 21 SNPs are likely or almost certain to be true positives.

Sign concordance: We compared the signs of the beta coefficients of the 506 EA lead SNPs (P EA < 10−5) with the beta coefficients for SZ. If the signs were aligned, we assigned a ‘1’ to the SNP and ‘0’ otherwise. By chance, sign concordance is expected to be 50%. We tested if the observed sign concordance is different from 50% using the binomial probability test50. 263 of the 506 SNPs have the same sign (52%, P = 0.40, two-sided).

Sign concordance is 58% (P = 0.10, two-sided) in the set of 132 EA lead SNPs that are also nominally significant for SZ (P EA < 1 × 10−5 and P SZ < 0.05).

Finally, for the 21 SNPs that passed Bonferroni correction for SZ (P EA < 1 × 10−5 and P SZ < 9.88 × 10−5), sign concordance is 62% (P = 0.38, twosided).

Raw enrichment factor (not corrected for LD score of SNPs): Because EA and SZ are highly polygenic, we tested for enrichment by taking the actual distribution of P values in the GWAS result files into account.

Due to the polygenic architecture of both traits, it is expected to find some EA-associated SNPs that are also associated with SZ just by chance even if both traits are genetically independent. Under this null hypothesis, the expected number of EA-associated lead SNPs that are also significantly associated with SZ is

$$E_{H_0}\left[ {N_{S,{\rm{EA}} \to {\rm{SZ}}}} \right] = N_{T,{\rm{EA}}} \times \tau _{P_{{\rm{EA}}}} \times \tau _{P_{{\rm{SZ}}}}$$

where N T,EA is the total number of independent lead SNPs in the EA GWAS results, and \(\tau _{P_{{\rm{EA}}}}\) and \(\tau _{P_{\rm{{SZ}}}}\) are the shares of SNPs in N T,EA that have P values for EA and SZ below a certain threshold, respectively.

We define the raw enrichment factor as

$$N_{S,{\rm{EA}} \to {\rm{SZ}}}/E\left[ {N_{S,{\rm{EA}} \to {\rm{SZ}}}} \right]$$

where N S,EA→SZ is the observed independent number of SNPs that pass both the P value thresholds P EA and P SZ .

We obtained N T,EA by applying the clumping procedure described above (PPM) without a P value threshold for EA, leading to 222,289 independent EA lead SNPs in our merged GWAS results file. For P EA < 10−5, we found 506 SNP (\(\tau _{P_{\rm{EA}}}\) = \(\frac{{506}}{{222,289}}\) = 0.2276%).

The Bonferroni threshold for testing 506 independent hypothesis is P SZ < \(\frac{{0.05}}{{506}}\) = 9.88 × 10−5. There are 341 independent SNPs in the SZ results that pass this threshold, thus \(\tau _{P_{\rm{SZ}}}\) = \(\frac{{341}}{{222,289}}\) = 0.1534%. Therefore, we expect [N S,EA→SZ ] = 222,289 × 0.2276% × 0.1534% = 0.776 (i.e. less than one) SNP to be jointly associated with both traits under the hull hypothesis of no genetic overlap. At these P value thresholds, we actually observe N S,EA→SZ = 21 SNPs, implying a raw enrichment factor of \(\frac{{21}}{{0.776}}\) = 27.

For P SZ < 0.05, we found 17,935 SNP (\(\tau _{P_{\rm{SZ}}}\) = \(\frac{{17,935}}{{222,289}}\) = 8.068%). Thus, [N S,EA→SZ ] = 222,289 × 0.2276% × 8.068% = 41. At this more liberal P value threshold, we actually observe N S,EA→SZ = 132 SNPs, implying a raw enrichment factor of \(\frac{{132}}{{41}}\) = 3.23.

Raw enrichment P value (not corrected for LD score of SNPs): Following Okbay et al.26, we performed a non-parametric test of joint enrichment that probes whether the EA lead SNPs are more strongly associated with SZ than randomly chosen sets of SNPs with MAF within one percentage point of the lead SNP. To perform our test, we randomly drew ten matched SNPs for each of the 506 EA lead SNPs with P EA < 10−5.

We then ranked the 506 × 10 randomly matched SNPs and the original 506 lead EA SNPs by P value and conducted a Mann–Whitney test51 of the null hypothesis that the P value distribution of the 506 EA lead SNPs are drawn from the same distribution as the 506 × 10 randomly matched SNPs. We reject the null hypothesis with P = 6.872 × 10−10 (Z = 6.169, two-sided). As a negative control test, we also calculated the raw enrichment P value of the first randomly drawn, MAF-matched set of SNPs against the remaining nine sets, yielding P = 0.17.

Repeating this raw enrichment test for the subset of 21 EA-associated SNPs that remained significantly associated with SZ after Bonferroni correction (threshold P SZ < \(\frac{{0.05}}{{506}}\) = 9.88 × 10−5) yields P = 5.44 × 10−14 (Z = 7.521, two-sided). The negative control test based on the raw enrichment P value of the first randomly drawn, MAF-matched set of SNPs against the remaining nine sets yields P = 0.34.

GWAS catalogue look-up

In order to investigate the novelty of the 21 SNP associations that were found significant for SZ after Bonferroni correction, reported in Table 1, we performed a look-up in the GWAS catalogue52 (revision 2016-08-25, downloaded on 2016-08-29, https://www.ebi.ac.uk/gwas/api/search/downloads/full) with the SNPs and all their ‘LD partners’ (i.e. all SNPs with an r2 > 0.5 within a 250 kb window). The LD partners were extracted with PLINK44 using a version of the 1000G reference panel specifically harmonized to combine 1000G phase 1 and phase 3 imputed data53, and the reference panel has been described previously26. The result of the GWAS catalogue look-up is reported in Supplementary Data 3.

Prediction of future GWAS loci for SZ

To identify LD partners and to clump our GWAS results, we used a threshold of r2 > 0.1 and a 1,000,000 kb window in the 1000 Genomes phase 1 version 3 European reference panel. Our SZ summary statistics contained 51,721 approximately independent SNPs with P SZ < 0.05. We identified 21,430 SNPs in LD with the 50 novel SNPs reported in ref.28 and 54,425 SNPs in LD with the 128 genome-wide significant loci that were previously reported5. We removed SNPs in LD with the previously GWAS hits from our analyses because those SNPs could (by definition) not be identified as novel. The remaining set of 51,528 approximately independent SNPs with P SZ < 0.05 in our SZ GWAS results contained one proxy for each of the 50 novel SNPs in ref.28. After removing SNPs in LD with previous GWAS hits, 110 SNPs with P SZ < 0.05 also exhibited P EA < 10−5 in the independent EA GWAS sample. Of those 110 SNPs, six were identified as novel SZ loci in the most recent GWAS dataset expansion28. Using Fisher’s exact test, we rejected the null hypothesis that the proportion of novel SNPs (6/110 vs 50/51528) is equal in the two sets (P = 2.4 × 10−9, two-sided). Furthermore, as a robustness check, we performed the analysis again by excluding the SNPs with MAF ≤ 0.1 and found similar results (P = 1.2 × 10−6). Thus, we conclude that conditioning GWAS results on SZ with independent GWAS evidence on EA significantly outperforms pure chance in predicting GWAS results on SZ from even larger samples.

Pleiotropy between EA and SZ

To explore if the loci identified by our PPM may have direct pleiotropic effects on EA and SZ, we applied genetic fine mapping using full GWAS results for both traits. Our procedure was as follows:

First, the SZ and EA GWAS results were merged into a single file and aligned such that the reference allele is identical. Ambiguous SNPs or SNPs which may be subject to strands flips were removed. Second, all SNPs within 500 kb upstream and downstream of the 21 significant lead SNPs from the PPM analyses were extracted. The pairwise LD between all SNPs in each window was computed. We then ran PAINTOR 3.030 which estimates the posterior probability of any SNP within a locus to be causal. We applied this procedure for EA and SZ. We then selected a 90% credibility set for EA and SZ, which reflects the broadest possible set of SNPs whose posterior probability covers 90% of the total posterior probability at that locus. We predetermined the maximum number of true causal loci to be 2. For the EA 90% credibility set, we then determined the posterior probability that this set contains the causal locus for SZ and vice versa. As the size of this set fluctuates between <1% of the locus size to ~30% of the total locus size (i.e. a more narrow set can be identified for some loci compared to others), we also computed 80, 65 and 50% credibility sets for EA and SZ (Supplementary Note 3 and Supplementary Data 4), which all have increasingly narrower sets of SNPs (the 50% credibility set is the narrowest set we investigated). Finally, we compute the ratio of the cross trait credibility for the 90% sets and the proportion of SNPs in the locus, which reflects the enrichment of signal over the baseline where each SNP is equally credible. We classify the probability of a locus being pleiotropic as low, medium, or high if the posterior probability of both the EA set on SZ and the SZ set on EA are <15, 15–45 or >45% respectively.

Biological annotations

To gain insights into possible biological pathways that are indicated by the PPM results, we applied DEPICT8,54 using a false discovery rate threshold of ≤0.05. DEPICT is a data-driven integrative method that uses reconstituted gene sets based on massive numbers of experiments measuring gene expression to (1) prioritize genes and gene sets and (2) identify tissues and cell types wereprioritised genes are highly expressed. The input for our analyses (DEPICT version 1 release 194) were the 132 EA lead SNPs that are also nominally associated with SZ.

For these 132 EA lead SNPs, we also used DEPICT to determine the enrichment of expression in particular tissues and cell types by testing whether the genes overlapping the GWAS loci are highly expressed in any of 209 Medical Subject Heading (MeSH) annotations.

To identify independent biological groupings, we computed the pairwise Pearson correlations of all significant gene sets using the ‘network_plot.py’ script provided with DEPICT. Next, we used the Affinity Propagation method on the Pearson distance matrix for clustering55. The Affinity Propagation method automatically chooses an exemplar for each cluster.

Furthermore, we prioritize genes using DEPICT. Any particular locus centred on a SNP may contain multiple genes. One straightforward approach is to nominate a gene that is closest to the SNP. But this approach does not consider if the expression of the gene is likely to be altered or regulated by the causal site in the locus. Therefore, we used DEPICT to map genes to associated loci, which prioritize important genes that share similar annotations in bioinformatic databases.

Significant reconstituted gene sets, tissues, cell types and prioritized genes identified by DEPICT are described in Supplementary Note 4.

LD-aware enrichment of PPM results across different traits

For SNP i in trait j, the expected χ2 statistic can be calculated as

$$E\left[ {Z_{ij}^2} \right] = \left( {N_j \times h^2_j \times {\rm{LDscore}}_i{\mathrm{/}}M} \right) + (1 + {\rm{Na}})_j$$

where N is the sample size of the target trait j, h2 is the heritability of trait j, LDscore i = \(\mathop {\sum }\limits_{k = 1}^M {\kern 1pt} r_{ik}^2\) for SNP i is calculated using HapMap3 SNPs from European-ancestry, M is the number of SNPs included in the calculation of the LD score (n = 1,173,569 SNPs), \(r_{jk}^2\) is the squared correlation between SNPs j and k in the HapMap3 reference panel and 1 + Na is the LD score regression intercept for trait j.

To determine whether a particular realization is significantly larger than expected (and thus the ratio χ2 observed /χ2 expected is significantly greater than one), we test each particular observed Z statistic (the square root of the χ2) for SNP j against a normal distribution with variance: (N j × h2 j × LDscore i /M) + (1 + Na) j .

We used precomputed LD scores available from the LDSC software38. As recommended by Bulik-Sullivan et al.38, we restricted our analysis to HapMap3 SNPs (using the merge-alleles flag) because these seem to be well-imputed in most studies. Out of 132 SNPs with P EA < 1 × 10−5and P SZ < 0.05, only 30 SNPs are directly present in HapMap3 SNP list (Supplementary Data 9). Therefore, we extracted proxy SNPs with r2 > 0.8 and a maximum distance of 500 kb to our missing EA lead SNPs and chose the one with the highest r2 as a proxy. After this step, we could include 105 (out of 132) SNPs in our analyses. For each of these 105 SNPs, we observed the Z-statistics in the publicly available GWAS results of the traits. Z-statistics were converted into χ2 statistics by squaring them. The LD score corrected enrichment per SNP for each trait is the ratio of the observed to the expected χ2.

Furthermore, since the SNPs considered for enrichment are independent, we can use Fisher’s method to combine the enrichment P values per SNP into a single P value per trait. The latter P value reflects excess enrichment for the set of SNPs beyond what is expected if these SNPs are part of the infinitesimal genetic contribution to the trait in question.

The results are shown in Supplementary Fig. 5 (and in Supplementary Data 10).

Our LD-aware enrichment test has two limitations. First, LD score regression assumes that allele frequency (AF) does not correlate with effect size, an assumption which has been empirically shown to be violated for low-frequency alleles56. Second, our test assumes the absence of selection on the trait. Variation in AF and the degree of negative selection could explain excess signal in low LD SNPs57. However, our raw enrichment P value is robust to this because it takes the AF of the candidate SNPs explicitly into account.

Replication of PPM results in the GRAS data collection

The GRAS data collection complies with the Helsinki Declaration and were approved by the Ethics Committees of the Universities of Göttingen and Greifswald, Germany or of collaborating centres. All subjects and/or their authorized legal representatives gave written informed consent.

We showed in our pre-registered analysis plan that our replication sample (GRAS) is not large enough to replicate individual SNPs (https://osf.io/dnhfk/). Instead, we decided at the outset to attempt replication of the proxy-phenotype analysis results using a PGS that consists of the >80 most strongly associated, independent SNPs. The set that best meets this criterion are the 132 independent EA lead SNPs that are also nominally associated with SZ (P SZ < 0.05). PGS for this set of 132 candidate SNPs were constructed using either the β coefficient estimates of the EA or the SZ GWAS meta-analysis, resulting in two different scores (named EA_132 and SZ_132).

In addition, we also constructed PGS for EA, SZ, BIP and neuroticism in the GRAS data collection using all available SNPs as control variables for multivariate prediction analyses (named EA_all, SZ_all, BIP_all and Neuro_all). Technical details are described below.

Polygenic score calculations in the GRAS data collection

PGS were calculated using PLINK version 1.944,45. We calculated eight different scores, which are described below. Supplementary Fig. 6 shows the distribution of the MAFs of all SNPs that were included in these PGS.

SZ scores: We received the GWAS summary statistics for SZ from the PGC excluding the data from our replication sample (GRAS). We constructed a PGS using the 132 EA lead SNPs (P EA < 10−5) that are also nominally associated with SZ (P SZ < 0.05). This score (SZ_132) is used for replication of the proxy-phenotype analyses.

Furthermore, we constructed a PGS using all 8,240,280 SNPs that survived quality control (SZ_all). Next, we applied the clumping procedure using r² > 0.1 and 1000 kb as the clumping parameters and the 1000 Genomes phase 1 version 3 European reference panel to estimate LD among SNPs, eventually leaving a set of 349,357 SNPs ready for profile scoring. The vast majority of the SNPs we included in the SZ PGS (79%) have MAF > 1% and MAF < 99% (Supplementary Fig. 6).

For secondary analyses on the prediction of SZ symptoms, we constructed two PGS using the 349,357 SNPs that have concordant (+ and +; or – and –) or discordant signs (+ and –; or – and +) for EA and SZ, respectively. This resulted in 174,734 and 174,623 independent SNPs with concordant or discordant, respectively. We used these approximately independent SNPs for profile scoring and call the resulting PGS Concordant and Discordant. We note that this approach of constructing the Concordant and Discordant scores is based on a very conservative LD-pruning algorithm because we first LD-pruned all SNPs that survived quality control for both phenotypes and then sort them based on sign concordance. Thus, this approach prunes for LD both within and across the Concordant and Discordant scores.

As an alternative, we also used a less conservative approach that only prunes for LD within scores. Specifically, this alternative approach first sorts all SNPs that pass quality control for both phenotypes (8,240,280 SNPs) based on sign concordance (4,147,926 SNPs with concordant and 4,092,354 SNPs with discordant signs) and then LD-prunes the resulting set, yielding 260,441 and 261,062 independent SNPs with concordant or discordant, respectively. We call the PGS resulting from this approach Concordant_more_SNPs and Discordant_more_SNPs.

To enable tests for the sensitivity of our results to the MAF distribution of SNPs included in our scores, we also constructed PGS using SNPs with MAF > 1% (i.e. dropping SNPs with allele frequency ≤0.01 or ≥0.99) and with MAF > 10% (i.e. dropping SNPs with allele frequency ≤0.10 or ≥0.90), respectively, from the 8,240,280 SNPs that survived quality control. Next, we used r2 > 0.1 and 1000 kb as the clumping parameters and the 1000 Genomes phase 1 version 3 European reference panel to estimate LD among SNPs, eventually leaving a set of 302,150 approximately independent lead-SNPs with MAF > 1% and 108,075 SNPs with MAF > 10% ready for profile scoring. The resulting scores are called SZ_all_1% and SZ_all_10%, respectively. We also constructed PGS from this set of SNPs that took the sign concordance between EA and SZ into account and call these PGS Concordant_1% (151,265 SNPs), Discordant_1% (150,885 SNPs) and Concordant_10% (54,287 SNPs), Discordant_10% (53,788 SNPs), respectively.

All scores were calculated using the –score function in PLINK using the natural log of the OR of the SNPs for SZ as effect sizes.

Educational attainment scores: Beta coefficients for the EA GWAS were approximated using \(\widehat {\beta _j}\) = \(\frac{{z_j}}{{\sqrt {N_j \ast 2 \ast {\rm{MAF}}_j \ast (1 - {\rm{MAF}}_j)} }}\), see Rietveld et al.47 for the derivation. Using these betas, we constructed a PGS using the 132 EA lead SNPs (P EA < 1 × 10−5) that are also nominally associated with SZ (P SZ < 0.05). The resulting score is called EA_132.

Furthermore, we constructed a PGS using all 8,240,280 SNPs that survived quality control. Next, we applied the clumping procedure using r2 > 0.1 and 1000 kb as the clumping parameters and the 1000 Genomes phase 1 version 3 European reference panel to estimate LD among SNPs, eventually leaving a set of 348,429 SNPs ready for profile scoring (Supplementary Fig. 6). The resulting score is called EA_all. Eighty percent of the SNPs we included in the EA PGS have MAF > 1% & MAF < 99% (Supplementary Fig. 6).

We also constructed PGS using SNPs with MAF > 1% i.e. dropping SNPs with allele frequency ≤0.01 or ≥0.99) and with MAF > 10% (i.e. dropping SNPs with allele frequency ≤0.10 or ≥0.90), respectively, from the 8,240,280 SNPs that survived quality control. Next, we used r2 > 0.1 and 1000 kb as the clumping parameters and the 1000 Genomes phase 1 version 3 European reference panel to estimate LD among SNPs, eventually leaving a set of 306,977 approximately independent lead-SNPs with MAF > 1% and 106,607 SNPs with MAF > 10% ready for profile scoring. The resulting scores are called EA_all_1% and EA_all_10%, respectively.

BIP score: We obtained GWAS summary statistics on BIP from the PGC58. We used the LD-pruned GWAS summary from PGC (‘pgc.bip.clump.2012–04.txt’) with a set of 108,834 LD-pruned SNPs ready for profile scoring. PGS for the GRAS data collection were calculated by the application of the –score function in PLINK using the natural log of the OR. The resulting score is called BIP_all.

Neuroticism score: We obtained GWAS summary statistics on Neuroticism from the SSGAC. The results are based on the analyses reported in Okbay et al.26 containing 6,524,432 variants. We applied the clumping procedure using r2 > 0.1 and 1000 kb as the clumping parameters and the 1000 Genomes phase 1 version 3 European reference panel to estimate LD among SNPs, eventually leaving a set of 232,483 SNPs ready for profile scoring (Supplementary Fig. 6). PGS for the GRAS data collection were calculated by the application of –score function in PLINK using the Neuroticism beta values. The resulting score is called Neuro_all.

Note that our replication sample (GRAS) was not included in the GWAS summary statistics of any of these traits.

Polygenic score correlations: We calculated Pearson correlations between all PGS that we constructed in the GRAS data collection (SZ_all, SZ_132, EA_all, EA_132, Concordant, Discordant, BIP_all and Neuro_all). Results for SZ patients and healthy controls together are reported in Supplementary Data 12a. We found very similar results among the SZ cases (Supplementary Data 12b) and healthy controls when we analysed them separately from each other (Supplementary Data 12c). These results were used to inform the correct multiple regression model specification for the polygenic prediction analyses (Supplementary Note 6 and 8).

GWIS

A GWIS infers genome-wide summary statistics for a (non-linear) function of phenotypes for which GWAS summary statistics are available41. Here, in particular, we wish to infer for each SNP the effect on SZ, conditioned upon its effect on BIP. One possible approximation involves a GWIS of the following linear regression function:

$${\rm{SZ}} = \beta \ast {\rm{BIP}} + e$$

where the parameter β is estimated from the genetic covariance between SZ and BIP and the genetic variance in BIP as β = \(\frac{{{\rm{cov}}_g\left( {\rm{SZ,BIP}} \right)}}{{{\rm{var}}_g(\rm{BIP})}}\). The residual (e) is actually our trait of intrest, for which we use the term SZ (min BIP). Using GWIS we infer the genome-wide summary statistics for SZ (min BIP) given the most recent PGC GWAS results for SZ (omitting the GRAS data collection)5 and BIP59. The effect size with respect to SZ (min BIP) for a single SNP is computed as:

$${\rm{eff}}_{\rm{sz}} - \beta \ast {\rm{eff}}_{\rm{BIP}} = {\rm{eff}}_e$$

The standard error for each SNP effect is approximated using the delta method and accounts for the possible effect of sample overlap between the SZ and BIP GWAS.

As data input, we used the GWAS results on SZ (excluding the GRAS data collection). GWAS results for BIP59 (6990 cases; 4820 controls) were obtained from the website of the PGC (https://www.med.unc.edu/pgc/files/resultfiles/pgc.cross.bip.zip).

Using the same method and data, we also ‘purged’ the genetic association results for BIP of their overlap with SZ, obtaining ‘unique’ BIP (min SZ) results.

Code availability

Source code for GWIS and LD-aware enrichment analyses are available at https://github.com/MichelNivard/EA_SZ.

Data availability

The GWAS summary statistics that support the findings of this study are available on the website of the SSGAC: http://www.thessgac.org/#!data/kuzq8. The GRAS data collection is not publicly available due to data protection laws in Germany that strictly safeguard the privacy ofstudy participants. To request access, contact the study’s principal investigator Prof. Dr. Hannelore Ehrenreich (ehrenreich@em.mpg.de)