Ethics approval and consent to participate<\/p>\n
All cohorts involving human participants that were included in this research have Institutional Review Board (IRB) approval following a full review of IRB applications and study consent forms by appropriate ethics committees. Details on each of the participating cohorts and IRB approvals are listed below:<\/p>\n
Hispanic ASD samples, IRB title: The Study of Novel Autism Genes, protocol no. 2012p001018, IRB of record: Mass General Brigham (previously known as PARTNERS), principal investigator on IRB: M. Daly.<\/p>\n
Chinese samples, IRB title: Neuropsychiatric genetics of a Chinese population from Shanghai Province, protocol no. IRB17-1379, IRB of record: Harvard T.H. Chan School of Public Health (HSPH), principal investigator on IRB: B. Neale.<\/p>\n
NeuroGAP-Psychosis:Ethiopia and South African samples, IRB title: Neuropsychiatric Genetics of African Populations \u2013 Psychosis (NeuroGAP-Psychosis), protocol no. IRB17- 0822, IRB of record: Harvard T. H. Chan School of Public Health (HSPH), principal investigator on IRB: K. Koenen.<\/p>\n
Paisa and Genomic Psychiatry Cohort, IRB title: Molecular Profiling of Psychiatric Disease, protocol no. 2014p001342, IRB of record: Mass General Brigham (previously known as Partners Healthcare), principal investigator on IRB: B. Neale.<\/p>\n
Inflammatory Bowel Disease samples, IRB title: The Broad Institute Study of Inflammatory Bowel Disease Genetics, protocol no. 2013P002634, IRB of record: Mass General Brigham (previously known as Partners Healthcare), principal investigator on IRB: R. Xavier.<\/p>\n
All participants who joined the studies listed under the above IRBs consented to the de-identified use of their data in disease research studies. During the consenting process, study staff explain the research goals and the basics of genetic processing activities. Participants are also given the option to opt out of studies if they choose. Samples used in this research and method development are all de-identified and have all been fully consented for inclusion in the study and any resulting publications.<\/p>\n
Wet lab protocolBiological samples<\/p>\n
BGE DNA samples have been extracted from saliva and whole blood specimens outside our laboratory. Samples derived from whole blood are generally preferred. Lower alignment rates of sequencing reads have been observed in saliva samples due to the presence of bacterial species. Furthermore, saliva samples can be more difficult to accurately quantify due to their consistency and heterogeneity. Following DNA extraction, we follow the BGE protocol described previously and shown in Supplementary Fig. 1<\/a>.<\/p>\n PCR-free library preparation<\/p>\n After an initial quantification, DNA is normalized to 50\u2009ng\u2009\u03bcl\u22121 and transferred into a 384-well plate. Normalized DNA is purified using an automated 2.75\u00d7 solid-phase reversible immobilization (SPRI) clean-up with Ampure XP Agencourt beads (Beckman Coulter). Cleaned DNA is then quantified by spectrophotometry (Lunatic, Unchained Labs) and normalized to 25\u2009ng\u2009\u03bcl\u22121. DNA (target of 134\u2009ng input) undergoes a reduced and customized fragmentation\/end-repair\/A-tailing reaction for Illumina-compatible PCR-free library construction with custom NEBNext Ultra II FS DNA Library Preparation Kits (New England Biolabs) using the following conditions: 37\u2009\u00b0C for 42.57\u2009min, 65\u2009\u00b0C for 30\u2009min. Unique, dual-indexed adaptors (NEBNext Unique Dual Index UMI Adaptors, New England Biolabs) are ligated to fragments (20\u2009\u00b0C for 20\u2009min), and libraries undergo two consecutive SPRI size selections (0.5\u00d7 and 0.55\u00d7). Cleaned and size-selected PCR-free libraries are quantified by qPCR (Kapa Library Quantification Kit, Roche), then normalized (to ~0.3\u20132\u2009nM, depending on concentrations) before pooling into a single tube and concentrated (Unagi, Unchained Labs).<\/p>\n The oligonucleotides used are commercially available. From New England Biolabs, we use four plates containing 384 unique adapters (cat. no. E7874L, E7876L, E7878L and E7395L), which also include the forward and reverse primer mix for exome capture.<\/p>\n Exome-captured library preparation<\/p>\n An aliquot from the prenormalized and prepooled PCR free libraries is used as input for PCR amplification (cycles are 98\u2009\u00b0C for 30\u2009s; 12 cycles of 98\u2009\u00b0C for 10\u2009sec and 65\u2009\u00b0C for 75\u2009s; 65\u2009\u00b0C for 5\u2009min) using the NEBNext Ultra II FS Library Preparation Kit and primers from the indexed adaptor kit (New England Biolabs). PCR-amplified libraries are quantified by spectrophotometer (Lunatic), purified with a 1\u00d7 SPRI-cleanup (Ampure) and normalized to 70\u2009ng\u2009\u03bcl\u22121. Samples are then pooled and undergo exome capture (Twist Alliance Clinical Research Exome probes from Twist Biosciences) using the recommended hybridization-capture protocol for xGen Hybridization Capture Core Reagents (Integrated DNA Technologies).<\/p>\n Blending of PCR-free genome- and exome-captured libraries<\/p>\n PCR-free and exome-captured pools are both qPCR-quantified on the same qPCR run. Nanomolar concentrations are taken into account to calculate the appropriate volumes to blend 33% WES with 67% WGS. BGE samples are again qPCR-quantified for sequencer loading calculations.<\/p>\n Sequencing<\/p>\n BGE blended pools with 384 samples containing unique barcodes are sequenced across six lanes of NovaSeqS4 (Illumina) with 2\u00d7 150-bp runs.<\/p>\n Datasets analyzed to evaluate BGE quality<\/p>\n BGE data used in these analyses were generated at the Broad Clinical Lab. Sample cohorts included in the dataset were recruited and submitted from the collaborating institutions for the PUMAS project, which includes cohorts from the GPC, Paisa population and NeuroGAP-Psychosis, as described below, in Table 2<\/a> and in Supplementary Table 1<\/a>.<\/p>\n The GPC is a multi-institutional collaboration led by Rutgers University. The GPC resource includes a National Institute of Mental Health (NIMH)-managed repository of genomic samples at SAMPLED, genotypic and sequence data, and detailed clinical and demographic data for investigations of schizophrenia, bipolar disorder, obsessive\u2013compulsive disorder and coronavirus disease from a variety of ancestries collected in the USA. Within this analysis, we included 4,553 samples with 3,926 passing QC filters.<\/p>\n The Paisa population is a genetic isolate from Colombia that has expanded rapidly following a series of migration-related bottlenecks. They have been the focus of genetics studies in neuropsychiatric disorders in the last decade. Paisa BGE data generated here emerged from a long-standing collaboration between teams at the Universidad de Antioquia, Medellin, Colombia, and University of California, Los Angeles, USA, to recruit samples from this population. Within the analysis, we included 9,007 samples, with 8,200 passing QC filters.<\/p>\n The Stanley Center at the Broad Institute of MIT and Harvard initiated the NeuroGAP-Psychosis project in 201527<\/a>, a collaboration with colleagues at Addis Ababa University in Ethiopia, KEMRI-Wellcome Trust in Kenya, Makerere University in Uganda, Moi University\/Moi Teaching and Referral Hospital in Kenya, the University of Cape Town in South Africa, and the Harvard T. H. Chan School of Public Health in the USA. With initial plans to recruit 35,000 participants (half cases with a diagnosis of schizophrenia or bipolar disorder and half controls), the target was expanded to 39,000 participants via the PUMAS Project awarded by NIMH. Sample recruitment ultimately exceeded the revised target, with >42,000 samples collected across the five NeuroGAP-Psychosis collection sites. Within the analysis, we included 39,886 samples, with 35,138 passing QC filters.<\/p>\n The Australian Genetics of Bipolar Study19<\/a> led by the QIMR team is a national cohort of adults diagnosed with either type I or type II bipolar disorder. Genotypes were obtained using the Isohelix GeneFix GFX-02 2-ml saliva collection devices. The QSkin Sun and Health Study20<\/a> includes individuals from Queensland randomly sampled from the Australian Electoral Roll, with saliva samples for genotyping obtained using the Oragene DNA self-collection kit. Given the differing kits and enrollment strategies used in sample collection, we analyzed these cohorts separately for the imputation concordance analysis. Within the analysis and totaling from both QIMR cohorts, we included 7,499 samples, with 7,209 passing QC filters.<\/p>\n Data from the PUMAS project are being deposited at the NIMH Data Archive (NDA). At the time of manuscript submission, genomic data from the first 10,000 samples have been submitted. All remaining samples from the PUMAS grant will be submitted at the end of the grant period. Data that are designated with NDA-GRU data use will be deposited into DNA collection #3805. Data designated with Disease-Specific (Mental Health) \u2013 DS (Mental Health) will be deposited into DNA collection #4538. Finally, data designated Health\/Medical\/Biomedical, NDA-HMB-MDS will be deposited into DNA collection #4539.<\/p>\n Exome quality filtering<\/p>\n We conducted QC filtering using Hail v0.28.128 to restrict to sites and variants with high confidence in exome data. We filtered out sites with more than six alleles that failed VQSR, were located in low-complexity regions, or fell outside the Twist target capture regions. Within an individual, genotype calls were filtered if read depth was <10\u00d7, genotype quality was <20, or allele balance was <0.2 or >0.8 in heterozygous calls, or <0.8 in homozygous alternate calls.<\/p>\n We filtered samples based on WES and WGS coverage, ancestry and WES sample quality metrics to restrict to high-quality samples for subsequent analysis. First, we removed samples with low WES or WGS coverage. Exome coverage per sample was determined by calculating the fraction of the exome target covered with at least 10\u00d7 read depth. Samples with exome fractions less than 90% or estimated WGS coverage of less than 1\u00d7 were removed. In addition, samples with chimeric or contamination read rates greater than 5% were removed, resulting in 853 removed samples.<\/p>\n Next, we determined genetic ancestry of the PUMAS samples using the quality-filtered WES portion of the BGE. We combined the PUMAS data with four reference panels with diverse ancestries: HGDP, 1kGP, AWI-GEN and the African Genome Variation Project (AGVP)18<\/a>,28<\/a>,29<\/a>. Throughout this Article, we use ancestry labels assigned by existing genomic reference panels, including EUR (European), AFR (African) and AMR (Admixed American\u2014an imprecise label introduced by the 1000 Genomes Project to describe individuals with recent admixture from multiple continents including Indigenous American ancestry). Sites overlapping between the PUMAS data and all reference datasets were filtered to biallelic variants with call rates >0.98 and MAFs >0.1%. We conducted linkage disequilibrium (LD) pruning to extract independent markers and calculated the top ten principal components across all samples using Hail\u2019s hwe_normalized_pca function. We fit a random forest algorithm using the top ten principal components to the reference datasets of known ancestry and then applied the algorithm to the PUMAS data. We required PUMAS samples to have a probability of at least 0.7 to assign ancestry. Samples with probabilities less than 0.7 were removed (N removed = 1,799).<\/p>\n We calculated sample quality metrics using Hail\u2019s sample_qc function. Within each ancestry and cohort combination, we removed samples with outlier values defined as values more than four median absolute deviations from the mean for the following metrics: N singletons, N insertions, N deletions, N transitions, N transversions, heterozygosity ratio, transition-to-transversion ratio and insertion-to-deletion ratio. Sample quality metrics filtering removed 2,317 samples.<\/p>\n Finally, we removed samples with discrepancies between imputed genetic sex and reported gender. Genetic sex was imputed separately for each genetic ancestry group by calculating the inbreeding coefficient on the X chromosome using common, independent markers after removing markers within the pseudoautosomal region. Samples with an inbreeding coefficient <0.6 were labeled \u2018female\u2019, and samples with an inbreeding coefficient >0.6 were called \u2018male\u2019. We removed samples whose imputed sex and reported gender did not match, resulting in 316 removed samples.<\/p>\n CNV calling and QCCNV calling<\/p>\n To call CNVs from the deep-coverage exome data, GATK-gCNV 4.1.0.0 was used on the exome intervals following a previously described pipeline10<\/a>. GATK-gCNV adjusts for known WES read-depth confounders such as GC content and mappability, while simultaneously adjusting for unspecified technical batch confounders such as sample extraction, sequencing, library preparation and mapping quality. We took the read-depth data as input from the BGE samples over a set of canonical transcript genomic intervals. In brief, the raw sequencing files were compressed into counts of the reads over the set of annotated exons and used as input into GATK-gCNV. A principal component analysis-based approach was used on the compressed observed counts to identify differences in the capture kit. Within each cohort, principal component analysis was used to define ancestrally similar batches of 1,000 samples; from each, a random subset of 200 samples were used to train a CNV-discovery model tailored to each batch. Subsequently, a distance-based and hybrid-density-based clustering approach was used to curate samples into batches to process in parallel. Once batching determination was completed, GATK-GCNV was used on each batch and metrics for filtering were produced by the Bayesian model underlying the data to balance recall and PPV.<\/p>\n CNV benchmarking<\/p>\n To benchmark BGE CNV exome data, we used data from 100 quads from the SSC30<\/a>,31<\/a>,32<\/a>, sequenced previously. Benchmarking was carried out on all rare CNVs (frequency <1%). The 100 quads have existing CNV calls31<\/a>,32<\/a> from deep WGS (~30\u00d7), and the CNV calls from WGS were used as a gold standard to benchmark the BGE CNVs. Recall was calculated by the proportion of WGS data sites that had a match in the BGE CNV callset. Moreover, for any given site, if at least 50% of samples that had that variant in the WGS data also had a GATK-gCNV call with consistent directionality (duplication or deletion) that overlapped at least 50% of captured intervals, this was considered a correct call. The optimal recall and PPV were observed for CNVs spanning more than four exons in the 100 quads. Samples were removed if there were more than ten rare CNVs, defined as <1% frequency across the cohort. Any samples with >more than ten high-quality CNVs were removed. We additionally removed samples that failed SNV QC. Subsequently, CNVs were filtered to have a quality score >200 and restricted to <1% frequency for each ancestry. We further benchmarked the performance of de novo CNVs derived from BGE data using GATK-gCNV as a function of the number of captured exons of canonical transcripts compared with validated WGS de novo CNVs from the quads. We additionally benchmarked the CNVs in saliva and blood against two independent constrained gene sets (top 1,000) based on LOEUF12<\/a> from gnomAD v2 and GISMO-mis13<\/a>.<\/p>\n Structural variant calling<\/p>\n The structural variants were identified using VISTA v114<\/a> for deletion and INSurVeyor v1.1.315<\/a> for insertions. VISTA is an ensemble tool that combines the best predictions from multiple individual callers, including Manta, Delly and others. We ran the individual callers and combined the results while applying the default parameters and quality filters. For InSurVeyor, we also applied the default parameters, which produced a high-quality callset that was used for analysis.<\/p>\n Imputation<\/p>\n To impute the low-coverage WGS data, we used the Genotype Likelihoods IMputation and Phasing (GLIMPSE2 v2.0.0) method16<\/a>. We used the HGDP+1kGP reference panel18<\/a> after filtering out singleton variants and indels, resulting in roughly about over 67 million variants available to be imputed.<\/p>\n To phase and impute this large-scale data in a cost-effective manner, we used Broad\u2019s Hail Batch Service to submit randomized batches of 200 individuals. Then, we merged imputed results from the batches using bcftools. In total, we imputed over 55,000 samples for just under US$20,000, averaging about US$0.36 per sample (Supplementary Table 5<\/a>).<\/p>\n AFs were corrected after merging by calculating a weighted mean of the batch allele frequencies (AFs) based on individual batch values (for i batches and N samples per batch):<\/p>\n $${\\mathrm{AF}}_{\\mathrm{cohort}}=\\frac{\\sum {\\mathrm{AF}}_{i}{N}_{i}}{{\\sum N}_{i}\\,}.$$<\/p>\n Similarly, INFO scores (I) were corrected post-merge via<\/p>\n $${I}_{\\mathrm{cohort}}=1\\,-\\,\\frac{\\sum (1-{I}_{i})\\times 2{N}_{i}\\,\\times \\,{\\mathrm{AF}}_{i}(1-{\\mathrm{AF}}_{i})}{2\\times \\sum {N}_{i}\\times \\frac{\\sum {\\mathrm{AF}}_{i}{N}_{i}}{\\sum {N}_{i}}\\left(1-\\frac{\\sum {\\mathrm{AF}}_{i}{N}_{i}}{\\sum {N}_{i}}\\right)}.$$<\/p>\n This process resulted in a final set of chromosome-specific BCF files, with all individuals of a cohort included and allele frequencies and INFO scores updated to reflect the full sample size.<\/p>\n QC of GSA data<\/p>\n Variants genotyped via the Illumina GSA platform were first filtered for a SNP call rate of at least 95%; then samples were filtered to have a call rate of at least 98%, heterogeneity F statistic or inbreeding coefficient |FHET| <0.20, and concordance between genetically inferred and reported sex. A second SNP filtering step with a call rate of at least 98% was performed after the sample QC, followed by filtering steps to remove SNPs with missingness differences >2% between cases and controls, and Hardy\u2013Weinberg equilibrium P value of 1\u2009\u00d7\u200910\u22126 for controls and 1\u2009\u00d7\u200910\u221210 for cases. We used King33<\/a> to calculate the genetic relationship matrix and removed one sample in any pairs of individuals with second-degree or closer relatives. Finally, we used the conform-gt tool provided by the Beagle V5.4 imputation software34<\/a> to align the strand and allele order within the GSA VCFs to the HGDP+1kGP reference panel18<\/a>.<\/p>\n Concordance between BGE and GSA GWAS array data<\/p>\n To evaluate BGE data accuracy post-imputation, we compared BGE with previously generated Illumina GSA data as ground truth on the subset of individuals with both data types (Supplementary Table 4<\/a>). Imputed variants were filtered to match those in the GSA. We computed nonreference concordance and aggregated R2 as a function of MAF for each cohort. To evaluate nonreference concordance3<\/a>, we define it as the number of true positives divided by the sum of the true positives, false negatives and false positives (excluding missing sites and SNPs with only reference alleles within the GSA subset), as shown in the formula below. We computed aggregated R2 per MAF bin by stacking all SNP dosages (from the imputed data) and genotype calls (from the GSA data) within a MAF bin, and calculating the squared Pearson\u2019s correlation coefficient. Note that, to avoid discrepancies in MAF for the NeuroGAP-Psychosis subsets (each with around 160 individuals), MAF values were based on the HGDP+1kGP AFR subset rather than in-sample allele frequencies. For the larger datasets, allele frequency values were taken from in-sample estimates.<\/p>\n $$\\begin{array}{l}{\\text{non-ref-concordance}}({{\\text{SNP}}}_{j})\\\\=\\displaystyle\\frac{{\\rm{\\#}}\\,{\\rm{true\\; positives}}}{{\\rm{\\#}}\\,{\\rm{true\\; positives}}+{\\rm{\\#}}\\,{\\rm{false\\; positives}}+{\\rm{\\#}}\\,{\\rm{false\\; negatives}}}.\\end{array}$$<\/p>\n Let XGSA be the matrix of genotypes for a given set of individuals and SNPs based on GSA array data, while XBGE, imputed is the matrix of imputed dosages for the same set of SNPs and individuals. We can calculate aggregated R2 per a given MAF bin for a SNP as follows35<\/a>:<\/p>\n $$\\begin{array}{l}\\mathrm{Aggregate}\\,{R}^{2}\\,\\mathrm{per}\\,\\mathrm{MAF}\\,\\mathrm{bin}\\\\={\\left({\\mathrm{Pearson}}^{^{\\prime} }{\\rm{s}}\\;\\mathrm{correlation}\\left(\\mathrm{vec}\\left({X}_{\\mathrm{GSA}}\\right),\\mathrm{vec}\\left({X}_{\\mathrm{BGE},\\mathrm{imputed}}\\right)\\right)\\right)}^{2}.\\end{array}$$<\/p>\n Concordance between BGE and GSA by local ancestry inference<\/p>\n We also evaluated post-imputation accuracy in BGE data stratified by each genotype\u2019s local ancestry background. We focus specifically on the Paisa and GPC populations, which have seen recent admixture between European, Native American Indigenous American and African continental ancestry. We first phased array genotypes via Beagle 5.436<\/a> using the full HGDP+1kGP reference panel (n\u2009=\u20094,099)18<\/a>. Next, we followed the Tractor37<\/a> tutorial and computed local ancestry with RFMix2 software22<\/a>. For the Paisa cohort, we modeled three-way local ancestry to capture admixture represented by EUR, AFR and AMR (as a proxy for Indigenous American ancestry) populations in the HGDP+1kGP reference panel. Specifically, we filtered the reference panel to only samples with AFR\/EUR\/AMR population labels. For local ancestry inference to be most useful, reference samples themselves should be representative of ancestral populations and non-admixed. Because AMR samples tend to themselves be recently admixed, we excluded 358\/549 reference AMR samples that have <90% Indigenous American ancestry as measured by RFMix2. As a result, the total number of reference samples Ntotal\u2009=\u20091,935 are split among 191 AMR, 752 EUR and 992 AFR. For the GPC population, we inferred two-way local ancestry to capture European and African admixture. When running RFMix2, we used the optional EM flag set to one iteration to account for residual admixture within the subsetted reference panel. Finally, we also included the \u2018-n 5\u2019 flag to account for reference panel sample size imbalance.<\/p>\n After local ancestry inference, two metrics were used to evaluate imputation accuracy for different ancestry backgrounds, including aggregate R2 and aggregate nonreference concordance, computed by stacking all SNPs within a MAF bin and evaluating a single nonreference concordance value. This is analogous to the aggregate R2 formula except that we compute concordance rather than squared correlation after stacking SNPs column-wise. Note that we favor the aggregate nonreference concordance metric as opposed to averaging SNP-by-SNP concordances due to the presence of rare SNPs. This is because a rare SNP may have only several copies of the minor allele, and by chance, they may all reside on one specific ancestry background. This leaves no copies of minor alleles for other ancestry backgrounds. Averaging SNP-by-SNP nonreference concordances will therefore include excess zeros for the remaining ancestries, artificially deflating their accuracy. Stacking rare SNPs rescues concordance for rare SNPs among all ancestry backgrounds by virtue of including more copies of the minor allele, ultimately producing more accurate and interpretable results.<\/p>\n Number of variants per data generation type<\/p>\n We evaluated a subset of samples from the NeuroGAP-Psychosis data (N\u2009=\u200979) that were assayed across three different data generation strategies, including BGE, Illumina GSA arrays and 30\u00d7 WGS. To count the number of variants per MAF category in the noncoding BGE data, we used the imputed data imputed with GLIMPSE2 after filtering for an INFO score of at least 0.8 and removing monomorphic variants. For the coding regions, we utilized the high-coverage exome data rather than the imputed data, filtering for VQSR PASS variants as well as excluding locus control regions and filtering for monomorphic variants. We required genotypes to have genotype quality >50 and depth >13, and normalized Phred-scaled likelihoods to be >50. These same filters were applied for both coding and noncoding regions in the WGS dataset. For the Illumina GSA GWAS array data, we used the imputed data generated previously3<\/a> using the phase 3 1000 Genomes Project reference panel with BEAGLE v5.1 for all categories, filtering for dosage R2 of at least 0.8 and removing monomorphic variants. Coding regions were defined by Twist target capture region intervals.<\/p>\n Reporting summary<\/p>\n