Congenital defects

Users can prepare a prioritization in four simple stages (see Figure 1). In the first stage, the species with which to work is selected among the six available species (H. sapiens, M. musculus, R. norvegicus, D. melanogaster, C. elegans and D. rerio). Then, the set of seed genes (or training genes) is prepared in the second stage. These are the genes that are already associated with the process of interest (e.g. genes already associated with congenital heart defects). Typically, users would like to have at least 5 but no more than 40 seed genes, but larger seed sets can still provide reliable results. The third stage consists in selecting the data sources to use for the prioritization. Users are advised to select the resources that best suit their specific problem and to avoid using conjointly redundant resources. The list of usable data sources depends on the selected species. There is no limit in the number of data sources to select, besides its impact on computing time. Finally, the set of candidate genes is defined in the last stage. The hypothesis is that at least one of these candidate genes is associated with the biological process of interest, but this association is yet unknown. This set can for example be derived from previous experiment (e.g. genes from a deletion frequently observed in patients with congenital heart defects). However, there is no size limit regarding this set, so users who have no a priori regarding relevant candidate genes can instead prioritize the whole genome. The single output is an ordered list of the candidate genes with on top the most promising candidate genes. In an ideal situation, the gene yet unknown associated with the process of interest should be ranked first or at least among the top ranking genes. Each candidate gene is also given a P-value that represents the significance of this combination of rankings. In addition, rankings for each individual data source are also available as to better understand the global ranking (e.g. to identify the sources that contributed the most to prioritize a given gene).
The algorithm behind Endeavour prioritizes genes in three simple steps (see Figure 1). In the first step, it trains a model of the biological process of interest, using the seed genes provided by the user. The model contains one sub-model per user-selected data source and is trained using simple statistics. For instance, for annotation-based data sources (e.g. Gene Ontology) only the features over-represented within the set of seed genes are kept in the sub-model. In the second step, the candidate genes are scored using the model built in the first step. More precisely, for a given data source, the scores are computed using the associated sub-model and represent how well the candidate genes fit this sub-model. For instance, for annotation-based data sources, the score of a candidate gene is the Fisher's omnibus combination of the P-values of its annotations that are present in the sub-model (i.e. annotations not present in the sub-model are ignored). At the end of this step, each data source is associated with a ranking of the candidate genes, with the most promising candidate genes on top. In the third step, these rankings, which correspond to prioritizations made using different data sources, are then integrated using order statistics to obtain a single global ranking (13 (link)). This method allows candidate genes with little data (e.g. poorly annotated genes) to be fairly compared with candidate genes that have a lot more data. The algorithm then outputs a ranked list of candidate genes, with P-values that represent the significance of this ranking. In addition, the prioritization results are now displayed graphically using a parallel coordinate representation so that users can easily check which data sources contributed the most to the global ranking.
Endeavour integrates 75 data sources for six species that can be classified into broad categories that describe what we know about genes. They are briefly described below, and a detailed list of all available databases is available in Supplementary Material 1. The ‘Gene and protein function’ category includes resources such as Gene Ontology (14 (link)) and InterPro (15 (link)), which are usually organized as annotations between ontologies (that describe function in a broad sense) and biological entities (i.e. genes or gene products). In addition, the category ‘Biomolecular pathways’ includes pathway databases such as Reactome (16 (link)), which are complementary to the purely ontological annotations described above. Several resources that describe gene or protein interaction networks, such as BioGrid (17 (link)) and IntAct (18 (link)), are also integrated and classified into the ‘Interaction networks’ category. Chemical data sets are also integrated in Endeavour within the ‘Chemical information’ category. These contain annotations between gene or gene products and other entities such as drugs. One example is the DrugBank database (19 (link)). The ‘Phenotypic information’ category gathers all resources that collect associations between genes or gene products and phenotypes or diseases, possibly using model organism data. Examples are the OMIM database for human Mendelian disorders (20 (link)) and the Rat Disease Ontology from the Rat Genome Database (21 (link)). Expression data are split into two categories depending on whether complete expression profiles are available (category ‘Expression profiles’), such as large expression data sets, or whether the data have already been summarized into annotations, such as within PaGenBase (22 (link)) (category ‘Expression ontologies’). The last category ‘Sequence-based features’ contains all resources that are based on gene, transcript or protein sequences. This includes protein sequence similarities computed using BLAST or predicted miRNA regulation using transcript sequences (23 (link),24 ).
The web server is based on a three-tier architecture, developed using Microsoft .Net technology and Microsoft SQL database. The core Endeavour library is written in Perl and the data are stored in a MySQL database. We have introduced a parallel coordinate representation to visualize our multidimensional numerical results. This gives the users the opportunity to easily recognize patterns in the data and identify possible correlation among the different sources. Uploaded data and prioritization results are kept private and not viewable by other users, and are in any case deleted after 30 days. The web site is free and open to all and there is no login requirement. Two examples derived from the literature are available for users who simply want to try out Endeavour. In addition, a manual describes how to run a simple prioritization and a help page contains hints on how to solve the more frequent issues.

Study population. The NBDPS recruits cases from population-based, active surveillance congenital anomaly registries in nine U.S. states and includes live births and stillbirths > 20 weeks gestation or at least 500 g, as well as elective terminations of prenatally diagnosed defects when available (Yoon et al. 2001 (link)). Arkansas, Iowa, and Massachusetts ascertain cases statewide, whereas California, Georgia, New York, North Carolina, Texas, and Utah ascertain cases in select counties. Cases are reviewed by clinical geneticists using standardized study protocols to determine study eligibility and classification, and cases with chromosomal/microdeletion disorders and disorders of known single-gene deletion causation are excluded. Controls are unaffected livebirths who are randomly selected from vital records or hospital records, depending upon study center. The NBDPS has been approved by the institutional review boards (IRBs) of all participating centers, and all participants provided written or oral informed consent before participation. These analyses were reviewed and approved by the University of North Carolina IRB.
For this analysis, the study population consisted of all controls and eligible cases with a simple, isolated CHD (i.e., a single CHD with no extra-cardiac birth defects present) and an estimated date of delivery (i.e., due date) from 1 October 1997 through 31 December 2006. During this time period, the participation response was 69% among all cases and 65% for controls. Within the NBDPS, a team of clinicians with expertise in pediatric cardiology reviewed information abstracted from the medical record and centrally assigned a single, detailed cardiac phenotype to each case whose diagnosis was confirmed by echocardiography, cardiac catheterization, surgery, or autopsy and documented in the medical record. Phenotypes were then aggregated into individual CHDs and defect groupings (Botto et al. 2007 (link)). The following additional groups were created because of limited sample size of individual defects: a) other conotruncal defects, which included common truncus, interrupted aortic arch–type B (IAA-type B), interrupted aortic arch–not otherwise specified (IAA-NOS), double outlet right ventricle associated with transposition of the great arteries (DORV-TGA) and not associated with TGA (DORV-other), and conoventricular septal defects (VSD-conoventricular); and b) atresias that included both pulmonary and tricuspid atresia. Simple, isolated CHD cases represented 64% (n = 12,383) of the total CHD cases. We restricted the analysis to offspring with a single CHD to create more etiologically homogeneous case groups, although this limits the generalizability of our findings. Women who reported having pregestational diabetes (types 1 and 2) during their pregnancy were excluded owing to the established association with CHD (Correa et al. 2008 ). Women living > 50 km from a pollutant-specific air monitor were excluded from that analysis.
Exposure assignment. Each woman reported the due date that was provided by her clinician during pregnancy to obtain the gestational age of the infant at birth. Using the gestational age to estimate the date of conception, we assigned calendar dates to each week of pregnancy. Women’s residential addresses during pregnancy were centrally geocoded to ensure consistency across study centers. Each geocoded address during weeks 2–8 of pregnancy was matched to the closest air monitor for each pollutant, with > 50% of the data available using ArcGISv10 (ESRI, Redlands, CA) and monitor locations obtained from U.S. EPA’s Air Quality System (U.S. EPA 2013 ). Participants from 1996–1998 were excluded from the analysis of PM_2.5 because monitoring began in 1999.
We used the daily maximum hourly measurement for CO, NO₂, and SO₂, the daily maximum 8-hr average for O₃, and 24-hr measurements of PM₁₀ and PM_2.5 to assign exposure. We averaged over the daily maximum or 24-hr measurements for weeks 2–8 of pregnancy to assign a 7-week and also 1-week averages of the daily values. We included week 2 in addition to the standard window of cardiac development, because of the potential for lag effects of air pollution (van den Hooven et al. 2012 (link)). If only a single measurement was taken during a given week, it was assigned as the weekly exposure. Ambient levels of each pollutant except O₃ were categorized into the following categories, using the distribution of pollutant concentration among controls: less than the 10th centile (referent), 10th centile to less than the median, the median to less than the 90th centile, and greater than or equal to the 90th centile. These categories captured the departure from linearity observed in initial, exploratory analyses (data not shown). For similar reasons, O₃ was categorized into quartiles. Centiles were calculated separately for the 7-week and 1-week measures of exposure.
Statistical analysis. The following variables obtained from the maternal interview were identified as potential confounders through directed acyclic graph analysis (Greenland et al. 1999 (link)) and included in the final adjustment set: maternal age, race/ethnicity, educational attainment, household income, tobacco smoking in the first month of pregnancy, alcohol consumption during the first trimester, and maternal nativity. Maternal age was represented as a single, continuous term, measured at the time of conception. Race/ethnicity was self-reported and categorized into the following groups: white non-Latino, black non-Latino, Latino, Asian or Pacific Islander, and other. Educational attainment was collapsed into six categories: 0–6 years of education, 7–11 years, high school graduate or equivalency, 1–3 years of college or trade school, 4 years of college or completion of a bachelor’s degree, and an advanced degree. Household income was self-reported as < $10,000 annually, > $50,000 annually, or in-between. We adjusted for any tobacco use in the first month of pregnancy and differentiated between some alcohol consumption (less than four drinks) and binge drinking (four or more drinks) during the first trimester. Maternal nativity was defined as self-report of being born outside the United States.
To account for potential differences in case ascertainment by study center, models were also adjusted for the center-specific ratio of septal defects to total CHDs. Identifying septal defects often depends on method of case ascertainment (Martin et al. 1989 (link)). All potential confounders, as well as distance to major roadway, prepregnancy body mass index (BMI), and maternal occupation status during pregnancy were assessed for effect measure modification by constructing logistic regression models with and without interaction terms and conducting likelihood ratio tests using an a priori alpha level of 0.1. Distance to the closest major road—defined as an interstate, U.S. highway, state, or larger county highway—was constructed using ArcGISv10 and then dichotomized at 50 m. Prepregnancy BMI was defined using self-reported maternal height and weight and categorized according to National Institutes of Health (1998) guidelines into underweight (BMI < 18.5), normal weight (18.5 ≤ BMI < 25), overweight (25 ≤ BMI < 30), and obese (BMI ≥ 30). Maternal occupation status was defined as ever working outside the home during any time during pregnancy.
For each pollutant, models were constructed to explore individual defects and defect-groupings. If a woman did not have at least one monitoring value for each week of exposure, she was excluded from the weekly analysis. We explored the relationships between all weeks and all defects because of uncertainty in pregnancy dating when using an estimated date of conception and the lack of clearly elucidated mechanisms by which cardiac development could be disrupted by exposure to air pollution. Animal research suggests that exposures outside the typical period of development for an individual heart structure could also be etiologically relevant (Morgan et al. 2008 (link)).
Because we simultaneously assessed multiple weeks of exposure and multiple defects/groupings, we constructed two-stage hierarchical regression models to account for the correlation between estimates and partially address multiple inference (Greenland 1992 (link); Witte et al. 1998 (link)). The first-stage, represented in Equation 1, was an unconditional, polytomous logistic regression model of individual CHDs on exposure (x) defined as either all 1-week averages of maximum or 24-hr pollutant values or the single 7-week average, and the full adjustment set (w) detailed above.
b_d is the vector of regression coefficients corresponding to pollutant exposure for an individual CHD (d), c_d is the vector of regression coefficients corresponding to the covariates for a given defect, and m is the total number of individual types of CHDs. The second-stage model, which defines how the first-stage betas are associated, is given in Equation 2:
β_i = Z_ir + δ_i, [2]
where Z_i is a row in the design matrix that includes an intercept term and then indicator variables for type of defect, broader defect grouping, and exposure week/level for the ith β, r is the vector of coefficients corresponding to the variables included in the design matrix, and δ_i are independent normal random variables with a mean of zero and a variance of τ² that describe the residual variation in β_i. The obtained second-stage coefficients, r, are used to estimate values toward which the first-stage coefficients will be shrunk, with the magnitude of the shrinkage depending on the precision of the maximum-likelihood estimate obtained in stage 1 and the value of the second-stage variance, τ² (Greenland 1992 (link); Witte et al. 1998 (link)). We fixed τ² at 0.5, corresponding to a prior belief with 95% certainty that the residual odds ratio (OR) will fall within a 16-fold span.
To assess whether our results were robust to changes in model specification, we conducted sensitivity analyses by setting the value of τ² to 0.25, corresponding to a 7-fold OR span, as well as to a value of 1, corresponding to a 50-fold span. We also explored different specifications for the design matrix, in turn defining the prior value as a common mean for all defects, a common mean for each defect, or a common mean for each exposure week/level, across defects. Individual defects with > 10 but < 100 cases were excluded from hierarchical models and explored using Firth’s penalized maximum-likelihood method to address the quasi-complete separation that occurred due to small sample size (Heinze and Schemper 2002 (link)). These defects included the individual defects collapsed into the other conotruncals and atresia categories described above; Ebstein’s anomaly, which was part of the right ventricular outflow tract obstruction (RVOTO) defect grouping; and muscular ventricular septal defects (VSD_muscular), which was part of the septal defect-grouping. IAA-type A and partial anomalous pulmonary venous return had < 10 cases each and were excluded from all individual analyses, but were included in the left ventricular outflow tract obstruction (LVOTO) and anomalous pulmonary venous return (APVR) defect groupings, respectively. To assess whether pollutant–defect relationships conformed to a monotonic dose response, we reanalyzed the data using incremental coding which compares each category of exposure to its immediate predecessor. If the incremental ORs are all above (or below) 1, the relationship conforms to a monotonic dose response (Maclure and Greenland 1992 (link)).
To explore associations with CHDs within a multipollutant context, a principal component analysis (PCA) was conducted among participants who lived within 50 km of each type of monitor. PCA is used to reduce the number of correlated variables into a smaller number of artificial variables that capture most of the variance of the original variables while being uncorrelated with each other (Hatcher 1994 ). This allows the resulting factors to be included within the same model, reducing issues of multicollinearity. Applying PCA, we retained components that accounted for at least the same or more variance than one of the original pollutant variables. We then applied a varimax rotation and calculated factor scores for each participant. These factor scores were categorized using the 10th, 50th, and 90th centiles and used to assign exposure in hierarchical models.

The data were derived from the Population-based Birth Defects Surveillance System in Hunan Province, China, which is run by the Hunan Provincial Health Commission. In 2008, the Hunan Provincial Health Commission selected Liuyang County and Shifeng District as population-based birth defects surveillance sites, which had undergone a comprehensive evaluation process by experts before the decision. These two places are in the central Hunan Province, with a resident population of about 1.5-2 million and approximately 20,000 live births per year, and the location, demographics, economic conditions, and healthcare facilities genuinely mirror those of the entire province.
The surveillance population included all live births, stillbirths, infant deaths, and legal termination of pregnancy between 28 weeks gestation and 42 days postpartum, whose mothers lived in Liuyang County and Shifeng District between 2014 and 2020. Surveillance data included demographic characteristics such as sex, residence, parents’ age, and other key information. In this study, almost all available demographic characteristics that may be associated with birth defects were chosen for analysis, including sex, residence, number of births, paternal age, maternal age, number of pregnancies, parity, and maternal household registration.
The maternal and child healthcare workers at the community health service centers in urban areas and village doctors in rural areas are responsible for collecting surveillance data. They follow up with the live infants until 42 days after birth. According to the “Maternal and Child Health Monitoring Manual in Hunan Province (2013 Edition)” formulated by the Hunan Provincial Health Commission, diagnostic methods for birth defects included clinical examination, ultrasonography, biochemical examination, chromosomal analysis, genetic testing, autopsy, and other appropriate examination. All birth defects should be diagnosed by medical institutions above the district or county level as soon as possible. Each quarter, county-level surveillance centers will collect the surveillance data and submit it to municipal surveillance centers for review, which then submit it to the provincial surveillance centers (the Hunan Provincial Maternal and Child Health Care Hospital) for review.
Birth defects are coded according to the WHO International Classification of Diseases (Tenth Revision, ICD-10, codes Q00–Q99). The ICD codes of common specific defects are as follows: congenital heart defects (Q20-Q26), congenital metabolic disorders, congenital metabolic disorders (E03, E25.0, E70-E90, D55.0), congenital limb defects (Q69-Q74), congenital ear defects (Q16-Q17), congenital kidney and urinary defects (Q60-Q64), cleft lip and/or palate (Q35-Q37), chromosomal abnormalities (Q90-Q99), congenital digestive system defects (Q38-Q45), congenital nervous system defects (Q00-Q07), and other unclassified defects (Q00–Q99, excluding the codes mentioned above).

A total of 10,542 (response rates: 73.42%) birth cohort study was conducted (2010–2012) at the Gansu Provincial Maternity and Child Care Hospital in Lanzhou, northwest of China. Portions of this text were previously published as part of a preprint (Li et al., 2022 (link)). The study excluded individuals with multiple gestation, stillbirths, non-CHDs birth defects and the case with incomplete information, resulting in 10,090 live-born singleton participants. Data were collected as previously described in Li et al. (2021) (link), Sun et al. (2022) (link) (as showed in Fig. 1). Ethical approval was obtained from the Gansu Provincial Maternity and Child Care Hospital’s the Human Investigation Committees.
All the CHDs were classified into “isolated,” “multiple,” or “syndrome” based on the defect’s complexity of heart. “Isolated CHDs” refers to the case with either an isolated CHD; “Multiple CHDs” means the case with one more cardiac malformation, (e.g., atrial septal defect with patent ductus arteriosus) and “Syndrome defects” refer the infants with distinct CHDs associated with any non-CHD congenital anomalies. The cases were then further classified into three main subtypes based on anatomical structures: septal defects, atrial septal defect (ASD) and patent ductus arteriosus (PDA). Subtypes with insufficient cases were not discussed in our study. CHDs cases and non-CHDs controls were matched by maternal age (≤2 years) and living area that living radius is within 2 km. Among 10,090 individuals, a total of 97 mothers gave birth to live, singleton infants with CHDs (cases), and 194 selected matched mothers with healthy children (controls) were recruited. The cases were grouped into isolated defects (43 cases), multiple defects (46 cases), and syndrome defects (eight cases). More detailed information on the study population selection can be found in our previous studies (Huang et al., 2022 (link); Mao et al., 2017 (link); Sun et al., 2022 (link)).