mTADA is a framework for identifying risk genes from de novo mutations in multiple traits

Joint analysis of multiple traits can result in the identification of associations not found through the analysis of each trait in isolation. Studies of neuropsychiatric disorders and congenital heart disease (CHD) which use de novo mutations (DNMs) from parent-offspring trios have reported multiple putatively causal genes. However, a joint analysis method designed to integrate DNMs from multiple studies has yet to be implemented. We here introduce multiple-trait TADA (mTADA) which jointly analyzes two traits using DNMs from non-overlapping family samples. We first demonstrate that mTADA is able to leverage genetic overlaps to increase the statistical power of risk-gene identification. We then apply mTADA to large datasets of >13,000 trios for five neuropsychiatric disorders and CHD. We report additional risk genes for schizophrenia, epileptic encephalopathies and CHD. We outline some shared and specific biological information of intellectual disability and CHD by conducting systems biology analyses of genes prioritized by mTADA.

T he analysis of multiple traits can help characterize the genetic architectures of complex disorders 1 . One approach is to meta-analyze results derived from separate single-trait studies 2 . However, joint analysis of multiple traits can better accommodate heterogeneity of genetic effects of the same variants or genes across traits 3,4 . Numerous studies have jointly analyzed two or more traits and successfully identified shared commonvariant associations [5][6][7][8] . In addition, additional risk loci have been identified using these approaches 7,9 . However, none of these studies has examined rare variation from case-control data, or de novo variants for which mutation rates should be taken into account. For these rare variants, gene based tests have identified several genes associated with different disorders [10][11][12][13] . Some recent studies have shown that there are multiple risk genes that are shared between neurodevelopmental disorders 10,14,15 , and also with congenital heart disease (CHD) 16,17 . These results are based on the intersection among the top prioritized genes from each disorder; therefore, reported numbers of genes shared by multiple disorders remain low 10,17 . Development of multi-trait rare-variant methods for neuropsychiatric disorders (NPDs) and related disorders will facilitate the understanding of this important aspect of genetic architecture for these phenotypes.
Currently, there is still a limitation in the risk gene identification for a single trait of NPDs and relevant disorders from parentoffspring trio studies. One reason is that risk gene discovery is underpowered when sample sizes are limited, as well as when relative risks are not large 10,11 . Multiple risk genes have been reported for undiagnosed developmental disorders (DD), intellectual disability (ID) and autism spectrum disorder (ASD) 12,18,19 thanks to large sample sizes and/or relative risks 10 . However, there are a few risk genes identified for schizophrenia (SCZ), epileptic encephalopathies (EE) and other disorders because of small genelevel relative risks or small sample sizes 10,20,21 . Increasing sample sizes will increase power to identify additional risk genes, but this is an expensive solution and may not be feasible for some studies. If there are genetic overlaps, methods that can leverage the information from one trait to increase power for risk-gene identification of another trait could help in obtaining additional genes for these disorders.
Here, we have developed a new statistical model, mTADA (multi-trait transmission and de novo association test), that jointly analyzes de novo mutations (DNMs) of two traits in order to estimate the gene-level genetic overlap of the two traits, and to identify additional risk genes for each analyzed trait as well as shared and specific risk genes. First, we utilize simulation data and demonstrate that, compared with a single-trait method, mTADA substantially improves the power of risk-gene identification when genetic overlaps increase, especially for traits with smaller sample sizes or smaller relative risks. For example, mTADA is able to statistically increase evidence for multiple genes in a tested trait which shows 1) marginally statistical evidence in that trait, and 2) strong evidence in the other trait if the two traits have a high genetic overlap. To illustrate the advantage of the new pipeline over its previous single-trait version, we apply the method to large data sets of different NPDs and CHD (>13,000 parent-offspring trios) and identify shared genes between each pair of these disorders. mTADA identifies additional risk genes for each disorder by borrowing the information of other traits. We validate these results in an independent cohort of 1,241 trios with CHD, 197 trios with EE, and 4,877 SCZ cases and 6,203 controls. In addition, we demonstrate that mTADA's results could be used to better understand the shared and specific biological information for two tested disorders by using multiple systems biology approaches to test the top prioritized risk genes of the CHD-ID pair. CHD-specific genes are specific to certain biological pathways.

Results
The mTADA framework. The mTADA method is gene-based and requires input data of the number of DNMs and mutation rate per gene. If the DNMs are stratified on the basis of predicted effect (e.g., 'missense', 'nonsense', etc.), then each geneannotation category should have its own mutation rate that reflects the predicted effects of the mutations within. In summary, for each gene, we consider four models M j (j = 0..3) reflecting four alternative hypotheses: the gene is associated with neither trait (H 0 ), the first trait only (H 1 ), the second trait only (H 2 ), or both traits (H 3 ). We assume prior probabilities π j (j = 0..3) for the four models and these π j are estimated from data and single-trait studies. DNMs are modeled using Poisson distributions with mean relative risks, mutation rates and sample sizes as main parameters 10 (Methods). For each gene, four posterior probabilities (PP), which are abbreviated as PP0, PP1, PP2 and PP3 respectively, are used to infer the status of the gene for the four models. To summarize the evidence for association with a given trait, we use the sum of PPs of models including the risk gene hypotheses for that trait, i.e., PP1 + PP3 for trait one and PP2 + PP3 for trait two (Fig. 1, Table 1, Methods).
Results of mTADA on simulated data. To validate the new method, we conducted simulation studies by using genetic parameters from real-data analyses of previous studies (Methods).
Power for single-trait risk gene discovery. We compared gene numbers identified by mTADA and our previous single-trait method, extTADA, using the same threshold PP > 0.8. For π 3 = 0 (no overlapping information), mTADA and extTADA reported nearly the same positive gene numbers (Fig. 2). However, mTADA identified more genes than extTADA when π 3 increased. In addition, mTADA's gene counts were also higher than those of extTADA when higher mean relative risks were used.
Comparison of risk-gene classification for single traits. We designed a simulation experiment to assess the performance in the classification of risk and non-risk genes. We applied extTADA to single-trait data from our simulated data. We then calculated areas under the Receiver Operating Characteristic curves (AUCs) for mTADA and extTADA using classification results from single-trait data. AUCs of both were equal when π 3 = 0 ( Fig. 2). However, AUCs of mTADA were higher than those of extTADA when π 3 's values were larger. In addition, mTADA also performed better than extTADA with larger mean relative risks.
The proportion of false positive shared risk genes for two traits with non-genetic overlaps. We estimated this information for identifying shared risk genes (i.e. associated with both traits). We simulated data with π 3 = 0 and calculated the proportion of shared risk genes (per 19,358 tested genes) using different PP Statistical models for four hypotheses in mTADA for one category of variants in each trait at the i th gene. mTADA assumes that the gene can be in one of four models M0..M3. πj (j = 0..3) is the prior probability of the j th model. xk and Nk (k = 1, 2) are the data and the sample size of the k th trait. μi is the mutation rate of the gene. For each trait, the relative risks of shared and specific genes (γk) are from a Gamma distribution with two parameters: γ k (mean relative risk) and βk (to control the variance of relative risks).   Fig. 3a).
Correlations between posterior probabilities and observed false discovery rates (FDRs). Since mTADA makes inference on risk genes using PPs, instead of commonly used FDR, we compared these two metrics. We calculated the correlation between PPs and observed false discovery rates (oFDRs) for all situations. For PP3, we found that PP = 0.8 and 0.5 approximately correspond to oFDR = 0.1 and 0.25, respectively. Small mean relative risks could lead to higher FDRs, but this inflation was modest (Fig. 3b). These results were similar for other situations: when genes were associated with only the first trait, only the second trait, single traits (e.g., Trait 1 or Trait 2 genes) ( Supplementary Figs. 1, 2).
The correlation between simulated and estimated values of π 3 was also assessed. For large mean relative risks, high correlations were observed for all sample sizes. For smaller mean relative risks (≤24), π 3 's values were over-or underestimated ( Supplementary  Fig. 3). However, these small differences did not affect the results of risk-gene identification (Fig. 2, Supplementary Figs. 1, 2).
The effects of misdiagnosis and ascertainment bias on the results. When sample phenotypes are misdiagnosed (a patient of one trait is mis-assigned to another), the estimated parameters of mTADA may be biased and this may affect the results. In another scenario, samples from one trait may contain a larger number of patients of the second trait than expected based on the comorbidity in the population. This ascertainment bias may also have an effect on mTADA's estimates. We tested the impact of these scenarios. Overall, π 3 and downstream results were not strongly affected when there was ascertainment bias. Similar results were also observed for misdiagnosis rates of 5-10% if the The correlation between posterior probabilities (x-axis) and observed false discovery rates (FDRs, y-axis). These are for the combination of different sample sizes (ntrio) and mean relative risks (mRR).
mean relative risks of the tested traits were not highly imbalanced. If the mean relative risks of one trait were substantially higher than those of the other trait, overestimation of π 3 might arise for misdiagnosis rates of ≥5%. Detailed results are in the Supplementary Note 1.
Application of mTADA to neuropsychiatric disorders and CHD data. mTADA was applied to DNM datasets of 15 pairs of six disorders: five neuropsychiatric disorders including DD, ID, ASD, SCZ, EE; and CHD. These DNMs were classified into different categories using annotation tools. Based on previous results 10, 22 , we used loss of function (LoF), missense damaging (MiD) DNMs for all disorders and also added synonymous DNMs within DNase I hypersensitive sites for SCZ (Methods). We defined the gene-level genetic overlap (gO) of two disorders as: gO ¼ 100% π 3 =ðπ 1 þ π 2 þ π 3 Þ. DD based results showed strong convergence with smaller credible intervals (CIs) because of its large sample size as well as high relative risks of DNMs ( Fig. 4a/b, Supplementary Table 1). As expected, high gOs were observed for pairs of DD, ID, and ASD (gO > 32%; π 3 > 0:018). CHD and EE had the lowest gO (gO = 2%, π 3 = 0.001) followed by SCZ-EE (gO = 4.6%, π 3 = 0.0023). Supplementary Fig. 4 shows sampling results of the proportions of overlapping risk genes for pairs of these traits, and Supplementary Fig. 5 shows the percentage of genetic overlaps for traits. The gO of ASD and SCZ which was approximately 16% (CI = 5.6-31.4%) was similar to previous studies (Supplementary Table 2).
We also compared mTADA and extTADA in the identification of risk genes for single traits using a threshold of PP > 0.8. For DD and ID, mTADA always performed better than extTADA (Fig. 4c). Similar results were observed for ASD; except for the pair ASD-SCZ in which mTADA was slightly better than extTADA for SCZ but extTADA was better than mTADA for ASD. For CHD, EE and SCZ, mTADA was better than extTADA when CHD was combined with DD.
Insights into top genes prioritized by mTADA. To better understand the top genes prioritized by mTADA, we extracted genes with PP > 0.8 for further analyses.
Overlapping genes between two traits. The highest number of overlapping genes was observed for DD and ID (89 genes) followed by ASD-DD (65 genes) and ASD-ID (47 genes). Eight genes (ARID1B, GABRB3, KCNQ2, STXBP1, CHD2, TLK2, POGZ, SCN2A) were observed for at least six pairs of disorders (Fig. 4d). POGZ and SCN2A were present in eight pairs of disorders. POGZ was significant for pairs relating to ASD, DD, ID, CHD and SZ while SCN2A was significant for pairs relating to ASD, DD, EE, ID and SCZ. We checked DNMs of these two genes. As expected, POGZ had no DNMs for EE, and SCN2A had no DNMs for CHD. Interestingly, in the latest CHD study 23 , POGZ was one of the top CHD genes while no DNMs were observed for SCN2A. In addition, in a recent study of 6,753 parent-offspring trios with neurodevelopmental disorders and epilepsy 24 , 16 DNMs were in POGZ, but only one DNM was from a patient who has both ID and epilepsy.
Significant genes of single traits. To demonstrate the application of mTADA in the identification of additional risk genes, we tested three disorders (CHD, EE, and SCZ) whose DN-based genes have not been reported as often as the three other disorders. We used DD-based results because the number of risk genes for the three disorders highly increased when their datasets were jointly analyzed with the DD dataset in Fig. 4c.
CHD. 33 genes were prioritized. 20/31 were not in the list of known CHD genes and in the meta-analysis results of a recent CHD study of Jin et al. 23 (Table 2). We validated these results by using different approaches. First, we tested the protein-protein interactions (PPIs) of these 33 genes by using the STRING database 25 . The number of edges was higher than expected between 33 protein nodes (PPI p = 6e-12, Fig. 5). Multiple protein products of novel and known genes interacted with each other. The number of interactions decreased when tested with only PPIs from experiments but was still significant (PPI p = 0.0174). Second, we tested these CHD genes from an independent data set which includes 1,241 trios and 226 cases 23 . From the 1,241 trios, three genes (CTNNB1, CUL3, LZTR1) of the 20 novel genes had LoF or MiD DNMs (Poisson-test p < 2.0e-4, Table 2). Each of these three genes had only one DNM in the primary analysis. In addition, these genes were not called significant genes by extTADA. Finally, we compared these 33 genes with the top 25 genes meta-analyzed by Jin et al. 23 . 8/33 were in the 25-gene list (permutation p < 9.99e-05; Table 2).
EE. There were 16 genes. Similar to top CHD genes, their protein products also had more interactions than expected by using the STRING database (PPI p = 3e-11, Supplementary  Fig. 6). Three genes HECW2, MLL, WDR19 were not in the list of EE genes on the Online Mendelian Inheritance in Man 26 . These three genes only had PP < 0.3 in extTADA. Interestingly, HECW2 had a DNM in a whole-genome-sequencing study recently 27 .
Biological insights into shared and specific genes from mTADA's analysis. To demonstrate the application of mTADA in helping to understand the shared and specific biological mechanism of two analyzed disorders. We extracted three gene lists (shared and specific genes) for each pair of disorders using a threshold of PP > 0.5. To increase the sample size for CHD, we combined both tested and independent datasets (Methods). We then focused on ID and CHD in this analysis because this pair of disorders had high numbers of risk genes for the three lists (30 shared genes, 40 ID-specific genes and 30 CHD-specific genes, Supplementary Data 2). Different systems biology approaches were used to test these three gene lists. First, we conducted gene-set enrichment analyses 30 using gene-ontology (GO) gene sets 31 . The majority of top enriched GO gene sets were related to heart/cardiocytedevelopment for CHD-specific genes, to chromatin modification or DNA binding for shared or ID-specific genes (Fig. 6a). Next, we used gene sets from a human single-cell RNA sequencing (scRNAseq) dataset of~4,000 cardiac cells from human embryos 32 . No overlaps were observed between shared or IDspecific genes with these gene sets, but interestingly the CHDspecific genes were enriched in multiple gene sets (Fig. 5b). We then tested the three gene lists by using mouse scRNAseq gene expression datasets from different brain regions 33 . The three gene lists were not significantly enriched in these cell types; however, for pyramidal cells, ID-specific genes were nominally significant while CHD-specific genes did not have the same direction (Fig. 5c). Finally, we used BrainSpan RNAseq gene expression data to cluster these three gene lists into spatiotemporal groups. Using eight time points and four regions as in recent studies 10,34 , shared and ID-specific genes were strongly expressed in the prenatal stages of the human brain while CHD-specific genes were expressed in both prenatal and postnatal stages for Region 3 including hippocampus, amygdala and striatum (Fig. 5d, Supplementary Fig. 7). The three other brain regions did not show strong differentiations between these three gene lists (Supplementary Fig. 7).

Discussion
In this paper, we propose a method to jointly analyze two traits (mTADA) using de novo exome sequencing data. The method is an extension of our previous work for single traits 10,11 . mTADA estimates the proportion of overlapping risk genes (π 3 ) between two traits, and then uses this information to infer how many overlapping risk genes exist between two traits. The pipeline is also able to infer the number of risk genes for each trait by calculating posterior probabilities (PPs) of genes for each trait. On simulated data, mTADA performs better than a single-trait approach, extTADA, on the identification of risk genes (Fig. 2). We applied mTADA to more than 13,000 trios of five neuropsychiatric disorders and congenital heart disease, and reported overlapping genes between these disorders. We also saw that mTADA reported more risk genes for these disorders than extTADA (Fig. 4). This suggests that mTADA can help in the identification of additional risk genes, especially for disorders whose large sample sizes are challenging to obtain or whose mean relative risks are small. For such disorders, users can combine the data of the disorders with large public data sets (e.g., trio data of ASD, DD) to prioritize risk genes. Using one-trait information to leverage the information for other traits has been successful in fine-mapping 35 and common-variant 36 studies. Based on our best knowledge, mTADA is the first tool using this approach for de novo mutation data. We hope that mTADA (https://github.com/ hoangtn/mTADA) will be generally useful for analyzing de novo mutation data across complex traits.
By using mTADA for prioritizing top genes, multiple overlapping genes were observed for CHD, DD, ID and ASD. This replicates a recent study 37 in which high overlapping genes were observed for CHD and neurodevelopmental disorders. Interestingly, CHD did not show any overlapping information with Table 2 Information of genes prioritized for congenital heart disease (CHD). another neurodevelopmental disorder: EE. Two genes SCN2A and POGZ which have been reported as risk genes for some of these disorders 23,[38][39][40] are top overlapping genes from mTADA (Fig. 4d), but they show different trends. No SCN2A DNMs are in CHD data and no POGZ DNMs are in EE data. One possible reason is that the sample size of EE is small in this study (356 trios). Another hypothesis might be that they do not have strong overlapping biological pathways. We did not see any overlapping information between SCZ and CHD, or SCZ and EE. We analyzed in depth the top prioritized genes of CHD, EE and SCZ (Fig. 5, Supplementary Fig. 6). Some top risk CHD and EE genes from mTADA are also reported in recent studies 23,27 . Multiple top CHD genes have only one DNM, but have DNMs in independent data sets (Table 2). This suggests that they might be real risk-genes for this disorder. Interestingly, we identify 20 CHD genes (posterior probabilities >0.8) which are not in the list of 253 curated known human/mouse CHD genes. 3 of these 20 genes have DNMs in an independent data set. This shows the benefit of using mTADA in the prediction of risk genes for CHD by borrowing the information of DD (Fig. 5, Table 2). We used different systems biology approaches to understand the shared and specific risk gene lists of ID and CHD. Some specific information emerged from these analyses. CHD-specific genes were enriched for heart/cardiocyte pathways, cell types while shared and IDspecific genes were strongly expressed in the prenatal stages of the human brain and enriched in regulatory and binding pathways (Fig. 6). This suggests that a model-based approach as mTADA can help shed light on the shared and specific biological mechanism between disorders with larger sample sizes. Although mTADA performs better than the single-trait based extTADA, it does have some limitations. mTADA uses the parameters of single traits from extTADA to infer π 3 . Using parameters from extTADA makes mTADA much faster in its calculation, it means mTADA relies on the results of the singletrait pipeline extTADA that uses a full Bayesian approach. Also, mTADA as well as extTADA use de novo counts for each gene and divide these counts into different categories similar to other rare variant based studies 12,[41][42][43] . In this current pipeline, we estimated π 3 directly from data. However, common-variant based genetic co-heritabilities 44 and transcriptomic correlations 45 for multiple pairs of NPDs are available now. Other studies which are able to incorporate the annotation information of each mutation, prior information for π 3 from previous studies may increase the power of mTADA or similar tools. In the current version, users Result of protein-protein interaction analysis for genes associated with congenital heart disease (CHD). These genes were prioritized by using undiagnosed developmental disorders (DD) information. This is the top 33 genes, posterior probabilities > 0.8, identified by mTADA using the data set of Homsy et al. 16 . Novel genes have red background and known genes have green background. Additional information for these genes is in Table 2.
can set a prior or change a distribution for π 3 . Comorbid information might be used as prior information for π 3 in the analyses of mTADA. For example, our estimated gene-level genetic overlaps which are inferred from π 3 's estimations are very high for pairs of EE and ID (~31%, CI = 18.9-43.7%), ASD and ID (~38.1, CI = 30-46.6%). These three disorders are also highly comorbid 46,47 ; therefore, this information may be used as priors. However, genomic results and comorbid information might not always have the same trends. For example, the genetic overlapping information between ASD and SCZ are high in our study (~16%, CI = 5.1-31.4%), in previous common-variant based (r g = 0.16, se = 0.06, p = 0.0071) and transcriptomic studies (rho~0.5, p < 0.001, Supplementary Table 2), but the comorbidity of the two disorders might not be strong or almost zero in some recent studies 48,49 . In addition, ASD and SCZ can have overlapping copy number variable regions 50,51 ; however, duplications can be significantly seen for one disorder and deletions can highly present in the other disorder 51 . Finally, for all analyzed disorders, even though we observed multiple overlapping genes for pairs of disorders, the origin of these overlaps could be different. For example, for each pair, some overlapping genes could have more loss-of-function DNMs for one disorder and more missense damaging DNMs for the other disorder, and vice versa. Future rare-variant studies which are able to obtain comorbidity information from the overlapping samples and compare this with the genetic information will shed light on the genetic and clinical relationship of these disorders. Also, studies which are designed to understand in depth the information of variant categories for overlapping genes can elucidate the genomic mechanism of disorders.
Our analysis of de novo mutation data of neuropsychiatric disorders and CHD also has some limitations, in particular, overlapping phenotypes may lead to violation of mTADA assumptions. DD samples include people with different disorders 52 and some of CHD samples may have other neuropsychiatric disorders 23 . In a recent study, the DD dataset was combined with the ID dataset to create a larger ID dataset because of the high proportion of people with ID inside the DD cohort 53 . In this study, even though we analyzed DD and ID separately to better understand the gene-level genetic overlaps between ID and other disorders, overlapping phenotypes may still affect the current results. We tested possible scenarios of overlapping phenotypes (Supplementary Note 1). The proportion of overlapping risk genes was modestly affected by ascertainment bias or by low percentages of misdiagnosed cases (<20%). However, this metric might be overestimated if misclassification rates were substantially high and the gene-level mean relative risks of one disorder were greatly different from those of the other disorder (Supplementary Note 1). The inflation could have an impact on analysis results, especially shared and specific risk genes. Nevertheless, mTADA always performed better than extTADA in the identification of risk genes for single-trait analyses in tested scenarios. It is possible that mTADA would benefit by jointly modeling these biases and this will be a future extension of the method.
With further development, the mTADA approach can be generalized further to consider more than two traits simultaneously, and the increased information could increase the number of identified risk genes but at a cost of increased computational time. Currently, the number of hypotheses increases exponentially to 2 N with N being the number of traits. To reduce computational time, another approach which uses a small number of latent probability vectors 54 might be used for more than two-trait studies.
In conclusion, mTADA can be very useful for better understanding the genetic correlation across disorders (via the proportion of overlapping risk genes), and to prioritize additional risk genes for disorders. The approach of mTADA can be used to identify shared/specific risk genes for different categories of one trait (e.g., loss of function and missense de novo mutations). Genetic information of de novo mutations and rare case/control variants can be different 55 , mTADA might be adopted to pipelines which are able to apply to DNMs and rare case/control variants as two traits.
Methods mTADA: statistical models and parameter estimation. The mTADA is designed to jointly analyze two traits using DNMs. We use statistical models of extTADA, a single-trait method, to model DNM counts for each trait in mTADA as presented in Table 1. The likelihood of the data across all N genes can be computed as where x ki and ϕ kj are the i th gene data and the j th modelʼs parameters for trait k (k = 1, 2). In addition, if the data include multiple categories of variants then P k ij ¼ Q n C l¼1 P k l ij with n C being the number of categories. For gene i, the statistical support for the j th model is captured by its posterior probability (PP ij ¼ , abbreviated as PP0, PP1, PP2 or PP3 for a gene). We use our single-trait pipeline, extTADA, to estimate the proportions of risk genes (π S 1 and π S 2 ), mean relative risks ( γ S 1 and γ S 2 ) and dispersion parameters (β S 1 and β S 2 ) for each single trait (described as the superscript). We use these values inside mTADA: Therefore, we only estimate π 3 inside mTADA. Bayesian models are built using the rstan package 56 . We use Markov Chain Monte Carlo (MCMC) within rstan to estimate π 3 . Convergence is diagnosed by the estimated potential scale reduction statistic (R) and visualizing traces of results. The Locfit package 57 is used to obtain the mode, CI of π 3 . We use the mode as the estimated value of π 3 . We also tested a model with different mean relative risks for shared and specific risk genes. The model was more complex for the estimation process of parameters but did not improve the riskgene identification. Therefore, this complex model was not used in our analysis.
Generation and analyses of simulated data. We simulated DNMs for genes under the mTADA model in Table 1. All 19,538 genes and their mutation rates from our current real dataset were used. A gene was assigned to one of the four groups (four models) by using the probability (π 0 , π 1 , π 2 , π 3 ). We used π S 1 ¼ 0:05 and π S 2 ¼ 0:03 which are approximately equal to ASD, ID and DD results in our single-trait study 10 . π 3 was simulated with different values between 0 to min (π S 1 ; π S 2 ); and π 0 , π 1 and π 2 were calculated as described in the section above. A range of mean relative risks were simulated for each of the two traits. Two mutation categories were simulated for each trait; therefore, there were four mean relative risks for the two traits. We used results from our previous studies 10, 11 and other studies 58,59 for simulated values of mean relative risks. We simulated 100 values of each combination of π 3 and mean relative risks. We then calculated the mean of these 100 simulation results.
To calculate the proportion of false positive genes when there was not a genetic overlap between two tested disorders, we simulated different combinations of genetic parameters with π 3 = 0. For each PP threshold, we divided the number of identified overlapping genes by the total tested genes (n = 19,358 genes in our analysis).
We also used simulated data to assess the correlation between true and observed π 3 values and between PPs and oFDRs. An oFDR at a PP threshold was defined as the number of false positive genes divided by the number of identified genes. To use mTADA for single traits, for the i th gene, we calculated PP i1 + PP i3 and PP i2 + PP i3 for the first and second trait respectively.
To compare risk gene classification performance between mTADA and extTADA on single traits, we used AUCs. We calculated true and false positive rates for extTADA and mTADA across PP thresholds, and calculated the areas under these ROC curves.
Real datasets of de novo mutations and variants. For primary analyses, we used the DNM data collected by Nguyen et al. 10 and CHD data from Homsy et al. 16 . These data included 356 EE trios; 5,122 ASD trios; 4,293 DD trios; 1,012 ID trios; 1,077 SCZ trios; and 1,213 CHD trios. DNMs were annotated and classified into multiple categories as in our previous work 10 as follows. For EE, ASD, DD, ID, and CHD, we used two categories 10 : loss-of-function (LoF) and missense damaging (MiD) DNMs. The LoF category included nonsense, essential splice site, and frameshift DNMs defined by Plink/Seq 60 while the MiD category included DNMs annotated as missense by Plink/Seq and predicted damaging by each of seven methods 41 : SIFT, Polyphen2_HDIV, Polyphen2_HVAR, LRT, PROVEAN, Muta-tionTaster, and MutationAssessor. For SCZ, we used LoF, MiD and synonymous mutations within DNase I hypersensitive sites because this category showed significant DNM enrichment in SCZ probands 22 and non-null mean relative risk in extTADA 10 . Mutation rates were calculated as described by Fromer, et al. 60 and Nguyen, et al. 10 .
For the validation of mTADA's results and for better understanding the specific and shared risk genes between tested disorders, other datasets were used in the analysis. First, we used independent datasets to validate mTADA results. For CHD, we extracted variant data of 2,871 probands from Jin et al. 23 . These samples include 2,445 trios (1,204 trios are inside the data set of extTADA and used in the primary analysis of this study) and 226 singletons 23 . Only independent CHD samples were used in the validation process. For EE, we used the whole-genome-sequencing trio data of Hamdan et al. 27 . This dataset includes 197 trios not included in our mTADA analyses. For SCZ, a case/control independent SCZ dataset from Genovese et al. 41 was used. Disruptive and damaging ultra-rare variants from 4,877 cases and 6,203 controls were extracted from Table S3 of the study 41 .
Known risk-gene datasets. We extracted lists of known risk genes from two sources. 253 curated known human/mouse CHD genes were obtained from the supplementary data set 2 of Jin et al. 23 . A list of EE genes from the Online Mendelian Inheritance in Man 26 was downloaded on September 02, 2019 using keywords "epileptic encephalopathy" and "epileptic encephalopathies".
Gene expression datasets. Human scRNAseq expression datasets of 4,000 cardiac cells were from 18 human embryos which ranged from 5 weeks (5W) to 25W of gestation. These were classified into four major cell types (cardiomyocytes (CMs), cardiac fibroblasts, endothelial cells (ECs), and valvar interstitial cells (VICs)), and also filtered and clustered into nice clusters. Gene lists of the nine clusters were extracted from Table S2 of Cui et al. 32 . scRNAseq transcriptome datasets were obtained from Skene et al. 33 via the link: http://www.hjerling-leffler-lab.org/data/ scz_singlecell/ (Downloaded on August 01, 2018). These datasets included 9,970 single cells. These cells were clustered into 24 different cell types. Spatiotemporal transcriptome data were obtained from BrainSpan 61 , divided into eight developmental time points (four prenatal and four postnatal) 62 . The average expression at each spatiotemporal point was calculated across samples. For each gene, average expression values were standardized across spatiotemporal points to obtain z-scores 10,34 . Z-scores were used for visualizing gene lists.
Analysis of de novo mutations using mTADA. extTADA was used to obtain the proportions of risk genes and the mean relative risks of each category for each disorder. These values were then used as input for mTADA to estimate π 3 and then to calculate PP ij ði ¼ 1::N; j ¼ 0::3; N ¼ 19; 358 genesÞ for each pair of traits. The default algorithm, No-U-Turn Sampler (NUTS), in the rstan package was used to estimate π 3 . Two independent chains and 10,000 steps for each chain were used in the sampling process. Only 1,000 samples from each chain were chosen for further analyses.
For primary analysis, we applied mTADA to NPDs and 1,213 CHD trios. For understanding the specific and shared risk genes between tested disorders, we combined both tested and independent datasets of CHD (2,445 trios) in jointly analyzing with other disorders to increase power. Finally, we also applied mTADA to the two CHD datasets (tested and independent datasets) to test the performance of the method as described in Supplementary Note 2.
Other statistical methods for real data analyses. We used the EWCE package 63 to calculate the enrichment of our gene lists and the expression data from the 24 mouse cell types. To test the significance of the overlap of two gene sets, a permutation approach was used. We chose two random gene sets whose lengths are the same as the two tested gene sets from the background genes (19,358 genes from mTADA). This was carried out N times (N = 10,000 in this study) and the numbers of overlapping genes were recorded in a vector m. A p-value was calculated as ðlength m m>m 0 ½ ð Þþ1Þ=ðlength m ð Þ þ 1Þ) in which m 0 is the observed number of overlapping genes between the two tested gene sets. To conduct PPI analyses, we used the STRING database and the package STRINGdb 25 from the Bioconductor project 64 , and p-values of protein-protein interactions were extracted from these analyses. To examine expression information of identified genes, we used the package mclust 65 to cluster BrainSpan gene expression data (z-scores) in heatmap analyses. To test the significance for individual genes from DNMs, we used a Poisson test. The R function ppoisðy À 1; lambda ¼ 2 Ntrio μ; lower:tail ¼ FALSEÞ in which y and μ are the number of DNMs and the mutation rate of the tested gene; Ntrio is the number of trios. All analyses were carried out using the R software 66 .
Reporting summary. Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Data availability
All analyzed results are in Supplementary Data 1