Study cohort description

We used data on a total of 172,270 individuals of Japanese and East Asian ancestry. Of these, data on 165,098 individuals were obtained from BBJ, which has enrolled 200,000 participants to date. BBJ is a multi-institutional hospital-based genome cohort that collected participants affected with at least one of 47 diseases20. We excluded (1) individuals with low genotyping call rates (<98%), (2) closely related individuals (PI_HAT0.125 by PLINK, v.1.90b4.4; https://www.cog-genomics.org/plink/) and (3) outliers from the Japanese cluster based on principal component analysis (PCA) using PLINK2 (v.2.00a2.3 and v.2.00a3; https://www.cog-genomics.org/plink/2.0/) with samples of the 1000 Genomes Projects. Further, we separated the BBJ individuals into two Japanese clusters22,27 the mainland cluster (n=156,151) and Ryukyu cluster (n=8,947), by visual inspection based on the PCA plot (Supplementary Fig. 1). All the participants provided written informed consent approved from ethics committees of RIKEN Center for Integrative Medical Sciences, and the Institute of Medical Sciences, the University of Tokyo.

The Japanese subjects in replication cohorts were collected from three Japanese population-based cohorts (the Nagahama cohort study, JBIC and the Osaka University healthy cohort). The Nagahama cohort study is a community-based cohort in Nagahama city, Shiga prefecture, Japan. The study recruited healthy individuals between the ages of 30 and 74 (ref. 46). JBIC consists of EpsteinBarr virus-transformed B lymphoblast cell lines of unrelated Japanese individuals47. Osaka University healthy cohort is a volunteer-based cohort study recruited from the Osaka University Graduate School of Medicine, the University of Tokyo and the University of Tsukuba48. For each cohort, we also excluded individuals with a low genotyping call rate, a high heterozygosity rate, closely related individuals (PI_HAT0.125) and PCA outliers from EAS populations28,48,49. In addition, we extracted the EAS individuals from UKB. UKB is a population-based cohort that recruited approximately 500,000 individuals between 40 and 69 years of age from across the United Kingdom50. We obtained EAS individuals from unrelated UKB individuals based on PCA visualization combined with the 1000 Genomes Projects (Supplementary Fig. 2). Finally, we included 16,119 individuals in the replication study (n=8,947 from BBJ Ryukyu, n=1,275 from Osaka University healthy cohort, n=2,945 from the Nagahama cohort study, n=1,110 from JBIC and n=1,842 from UKB EAS). This study was approved by the ethical committee of Osaka University Graduate School of Medicine.

BBJ collected baseline clinical information and dietary and activity habits information through interviews and reviews of medical records using a standardized questionnaire. We selected 81 traits (57 anthropometric traits and biomarkers, 11 dietary habits, six behavioural traits, six diseases and one dummy; Supplementary Tables 24). We used these data from participants above the age of 18, and drinking and smoking traits from those above the age of 20. We normalized each anthropometric trait and biomarker traits by applying rank-based inverse normal transformation as previously reported (Supplementary Table 8)51,52,53. For each dietary habit, the participants were asked to clarify the frequency of consumption on a four-point scale, and we assigned the corresponding values to their responses as previously described26, where almost every day=7, 34 days per week=3.5, 12 days per week=1.5 and rarely=0. Behavioural traits included ever versus never drinking and ever versus never smoking54 as binary traits, and the frequency of four PAs (light-PA, gymnastics, walking and sports). For each PA, participants were also asked for the frequency and the length of time per week on a seven-point scale, and we quantified the activity by converting the responses to total minutes of activity time per week (minweek1), where 30 (15) minday1=210 (105), <30 (15) minday1=140 (70), three to four times a week for 30 (15) min=105 (52.5), three to four times a week for <30 (15) min=70 (35), one to two times a week for 30 (15) min=45 (22.5), one to two times a week for <30 (15) min=30 (15) and rarely=0 (the number in parentheses indicates gymnastics time).

For disease phenotypes, cases with myocardial infarction, stable angina and unstable angina were reclassified as CAD. We then selected six diseases from the target disease of BBJ (T2D, dyslipidaemia, cataract, CAD, arrhythmia and ischaemic stroke), where the number of cases exceeded 10,000 individuals55.

In addition, we set a dummy phenotype as a negative control. We generated a phenotype with heritability (h2=0.5) from 10,000 causal variants randomly sampled from BBJ GWAS data using GCTA GWAS simulation56. The phenotype followed the model yj=gj+ej, where gj=i(Wiji) and Wij=(xij2pi)[2pi(1pi)]1/2, where xij is the genotype for the ith causal variant of the jth individual, pi is the allele frequency of the ith causal variant within a population and ej is the residual effect generated from a normal distribution with mean 0 and variance Var(gj)(1h2)/h2. i is the effect size of the ith causal variant generated from a normal distribution with mean 0 and variance 1 (ref. 57). The values were normalized by applying a rank-based inverse normal transformation.

The BBJ GWAS data were genotyped using the Illumina HumanOmniExpressExome BeadChip or a combination of the Illumina HumanOmniExpress and HumanExome BeadChips. The quality control of the genotypes was described elsewhere51. In brief, we excluded variants satisfying the following criteria: (1) call rate <99%, (2) P value for HWE<1.0106, (3) number of heterozygotes <5 and (4) a concordance rate <99.5% or a non-reference concordance rate between the GWAS array and whole genome sequencing. The genotype data were phased by Eagle (v.2; https://alkesgroup.broadinstitute.org/Eagle/), and imputed with the 1000 Genomes Project Phase3 (v.5) and BBJ1K using Minimac3 software (v.2.0.1; https://genome.sph.umich.edu/wiki/Minimac3). After imputation, we excluded variants with an imputation quality of R-square (Rsq)<0.7 and those with a minor allele frequency (MAF)<1%.

As for the other Japanese datasets, JBIC was genotyped using Illumina HumanCoreExome Beadchip. As stringent quality control filters, we excluded the variants that satisfied (1) call rate<0.99, (2) MAF<1% and (3) HWE P<1.0107 (ref. 47). Osaka University healthy cohort was genotyped using Illumina Infinium Asian Screening Array. We excluded the variants that satisfied (1) call rate<0.99, (2) minor allele count<5 and (3) HWE P<1.0105 (ref. 48). The Nagahama cohort study was genotyped using six genotype arrays. We then selected two platforms (Illumina Human610-Quad Beadchip and Illumina HumanOmni2.5-4v1 Beadchip) with a large number of samples. For each of the two datasets, we excluded variants with (1) call rate<0.98, (2) MAF<1% and (3) HWE P<1.0106 (ref. 28). Genotype data were phased by Shapeit (v.2; https://mathgen.stats.ox.ac.uk/genetics_software/shapeit/shapeit.html) or Eagle, and imputed with the reference panel from the 1000 Genomes Project Phase3 (v.5) and BBJ1K using Mimimac3. After imputation, we excluded variants with an imputation quality of Rsq<0.7 and MAF<1%.

The UKB project was genotyped using either Applied Biosystems UK BiLEVE Axiom Array or Applied Biosystems UKB Axiom Array. The genotypes were imputed using the Haplotype Reference Consortium, UK10K and the 1000 Genomes Phase 3 reference panel by IMPUTE4. The detailed characteristics of the cohort and genotypephenotype data were described elsewhere50. We extracted EAS individuals and excluded variants with INFO score 0.8 and MAF1%.

As independent external reference GWASs or genotype data of Japanese ancestry were not publicly available, we adopted a tenfold LOGO meta-analysis to maintain both the accuracy of the GWAS statistics and the statistical power in PGS21. We first randomly split the BBJ mainland samples into the 10 target subsets. GWAS was performed on 81 complex traits for samples excluding the target subset using GCTA-fastGWA (v.1.93.3beta2; https://cnsgenomics.com/software/gcta/#Overview) as a MLM approach with 7,401,847 autosomal variants23,24. For GCTA-fastGWA, we computed a sparse genetic relationship matrix (GRM) for BBJ participants (n=182,961) using slightly LD-pruning variants (LD-pruning parameters in PLINK: indep-pairwise 1000 100 0.9, and MAF1%, sparse GRM parameter: make-bK-sparse 0.05). Regarding the 57 anthropometric traits and biomarkers, the 11 dietary traits, the four PA traits and the two binary traits in the behavioural traits, we fitted age, age-squared, sex, the top 20 PCs and 47 disease status as covariates. For the six diseases, we also fitted age, age-squared, sex and the top 20 PCs as covariates. We also performed GWAS using GCTA-fastGWA for all individuals in the BBJ mainland cluster to apply to other Japanese or EAS datasets. LD score regression (LDSC, v.1.0.0; https://github.com/bulik/ldsc) was applied to the summary statistics of the whole-sample GWAS to estimate potential population stratification. We adopted the HapMap3 SNPs, excluding the human leukocyte antigen region, using precomputed LD scores from 1KG EAS downloaded from the LDSC software website (Supplementary Table 5)58.

To estimate phenotypic variances explained by imputed data for some of the traits, we applied GREML-LDMS using GCTA (v.1.93.2beta; https://cnsgenomics.com/software/gcta/#Overview)57. We created the GRM using all variants for BBJ mainland samples. We estimated LD scores using default parameters in GCTA, and stratified SNPs into LD score quartiles. Next, we divided the SNPs within each LD score quartile into six MAF groups (MAF<5%, 5%MAF<10%, 10%MAF<20%, 20%MAF<30%, 30%MAF<40%, 40%MAF) and generated 24 GRMs. We calculated the phenotypic variance for each GRM and summed them to derive the total phenotypic variance (Supplementary Table 7). In the calculations, we randomly sampled 50,000 unrelated individuals (GRM<0.05) randomly downsampled from BBJ mainland individuals to avoid computational burden and used the same normalized values for quantitative traits and covariates for binary traits as used in the GWAS analysis.

To derive PGSs of individuals in each of the target subsets, we applied PRS-CS (https://github.com/getian107/PRScs) to construct PGSs that included genome-wide HapMap3 variants. PRS-CS is one of the beta shrinkage methods, which applies a Bayesian regression framework to identify posterior variant effect sizes based on continuous shrinkage before using both GWAS summary data and the external LD reference panel25. When the training sample size was large enough and the casecontrol imbalance was small, the automated optimization model (PRS-CS-auto) had the same precision as the grid model59,60. Therefore, for each of the target folds, we estimated the posterior mean effects of SNPs from the MLM-GWAS summary data of all training samples using PRS-CS-auto with the precomputed HapMap3 SNP LD reference panel from 1KG EAS downloaded from the PRS-CS website. We calculated PGSodd and PGSeven of individuals within the target subset using the estimated posterior effect of SNPs by PLINK2 score function. We normalized the calculated PGSs for each trait in each target subset to compare the effect sizes across the phenotypes.

We quantified the trait variance explained by the derived PGSs in individuals within one withheld subgroup. Each trait was modelled as a combination of PGS and all covariates. The null hypothesis used the same model without the PGS term. We calculated the adjusted R2 for quantitative traits and the Nagelkerkes R2 for binary traits (Supplementary Table 5).

For GPD estimation, we performed PCA of even and odd number chromosomes for each of the target subsets. We then estimated GPD using a linear regression method following the formula based on the original study18:

$${mathrm{PGS}_{rm{odd}}}approx theta _{mathrm{{even}}_{mathrm{to}}_{mathrm{odd}}}{mathrm{PGS}_{rm{even}}} + 20{mathrm{PCs}_{rm{even}}}$$

$${mathrm{PGS}_{rm{even}}}approx theta _{mathrm{{odd}}_{mathrm{to}}_{mathrm{even}}}{mathrm{PGS}_{rm{odd}}} + 20{mathrm{PCs}_{rm{odd}}}$$

where PGS is the scaled polygenic score, PCs are the results of the PCA and is the estimate of GPD. We further meta-analysed the GPD estimate from each of the ten subsets using the fixed effect method using metafor (v.1.9-9; http://www.metafor-project.org/doku.php/metafor) implemented in R (v.3.4.0; https://www.r-project.org/). We also estimated the GPD for the other Japanese and EAS datasets using the summary results of the whole BBJ sample GWASs by PRS-CS-auto. Finally, we performed a meta-analysis on the GPD estimates from the BBJ and other Japanese and EAS datasets by the fixed effect method using metafor. We estimated the P value of meta-analysed GPD using the Wald test.

To assess the robustness of our analysis to the chosen grouping of chromosomes, we altered the combinations of chromosomes such that the number of SNPs was the same in the two groups: (1) first half and second half; chromosomes 18 versus chromosomes 922, and (2) pseudo-random chromosomes; chromosomes 1, 3, 5, 6, 9, 10, 13, 14, 17 and 18 versus chromosomes 2, 4, 7, 8, 11, 12, 15, 16, 19, 20, 21 and 22. Using the two alternative combinations, we estimated the GPD for each cohort and meta-analysed the results.

We also calculated the theoretical GPD derived in the original study18. The theoretical value () followed the formula,

$$theta = frac{{rho f_0}}{{2 - rho (2 - f_0)}}left[ {1 + frac{{M(1 - rho )}}{{nh_{mathrm{{eq}}}^2}}left{ {1 + frac{{rho f_0}}{{2(1 - rho )}}} right}^{ - 3}} right]^{ - 1}$$

where (rho = rh_{mathrm{{eq}}}^2), r is a phenotypic correlation between spouses, (h_{mathrm{eq}}^2) is an equilibrium heritability of the phenotype, (f_0 approx f_{mathrm{{eq}}}/(1 - rho )), (f_{{rm{eq}}} = h_{{rm{snp}}}^2/h_{{rm{eq}}}^2,) (h_{{rm{snp}}}^2) is a SNP-based heritability, M is the number of causal variants and n is the sample size of the GWAS.

We analysed individuals of white British ancestry determined by PCA (n=337,139) from UKB by adopting the tenfold LOGO approach to the six available traits (adult height, BMI, T2D, CAD, duration of light-PA and yoghurt consumption)50. When adult height and BMI were measured multiple times, we adopted the mean value to obtain a single value per participant and normalized the values using the rank-based inverse normal transformation method. Regarding T2D, the case was defined following the ICD-10 codes and probable T2D and possible T2D in a T2D inference algorithm based on Eastwood et al.61. We also defined individuals without any diabetes status as the T2D control based on ICD-10 and the inference algorithm. As for CAD, the case was extracted following ICD-10 codes, surgical procedure recodes, self-reported illness codes and self-reported operation codes based on Fall et al.62. Regarding the duration of light-PA (Data-Field 104920), we extracted the data from instance 0 (n=70,692) and converted the coding to categorical values. Regarding the consumption of yoghurt, we extracted data from instance 0 within consumers of yoghurt/ice cream as binary traits (n=70,692 and Data-Field 102080). From the imputed GWAS data, we excluded the variants that satisfied MAF1% and INFO score 0.8, and fastGWA conducted generalized MLM approaches for nine subset samples with adjustment for age, age-squared, sex, top 20 PCs, ascertainment centre information and batch information as covariates. For the six phenotypes, we estimated the PGSs for odd and even chromosomes by PRS-CS-auto using genome-wide HapMap SNPs and the 1KG EUR LD reference panel, and the GPD was estimated in the same way as described in the Japanese study. We further meta-analysed the GPD estimate from each of the ten subgroups by the fixed effect method using metafor.

Further information on research design is available in the Nature Research Reporting Summary linked to this article.

See the original post:
Genetic footprints of assortative mating in the Japanese population - Nature.com