Statistical and Ontological Analysis of Adverse Events Associated with Monovalent and Combination Vaccines against Hepatitis A and B Diseases

Vaccinations often induce various adverse events (AEs), and sometimes serious AEs (SAEs). While many vaccines are used in combination, the effects of vaccine-vaccine interactions (VVIs) on vaccine AEs are rarely studied. In this study, AE profiles induced by hepatitis A vaccine (Havrix), hepatitis B vaccine (Engerix-B), and hepatitis A and B combination vaccine (Twinrix) were studied using the VAERS data. From May 2001 to January 2015, VAERS recorded 941, 3,885, and 1,624 AE case reports where patients aged at least 18 years old were vaccinated with only Havrix, Engerix-B, and Twinrix, respectively. Using these data, our statistical analysis identified 46, 69, and 82 AEs significantly associated with Havrix, Engerix-B, and Twinrix, respectively. Based on the Ontology of Adverse Events (OAE) hierarchical classification, these AEs were enriched in the AEs related to behavioral and neurological conditions, immune system, and investigation results. Twenty-nine AEs were classified as SAEs and mainly related to immune conditions. Using a logistic regression model accompanied with MCMC sampling, 13 AEs (e.g., hepatosplenomegaly) were identified to result from VVI synergistic effects. Classifications of these 13 AEs using OAE and MedDRA hierarchies confirmed the advantages of the OAE-based method over MedDRA in AE term hierarchical analysis.

generate accurate and early prediction of vaccine-vaccine interactions (VVIs) to prevent the VVI-induced newfound or serious AEs. Unfortunately, rare studies have paid attention to VVIs.
In order to monitor post-licensure vaccine AEs, the Centers for Disease Control and Prevention (CDC) and the FDA established the Vaccine Adverse Event Report System (VAERS) surveillance program in 1990 9 . VAERS typically receives over 30,000 case reports annually from various submitters including professional healthcare providers, vaccinees, vaccinees' relatives, and vaccine manufacturers. However, as a passive surveillance system, VAERS is subject to a number of well-described limitations such as high variability in report quality, biased reporting, underreporting and the inability in VAE causality determination 10 . Nevertheless, with statistical support, VAERS has been used widely to identify potential vaccine safety issues 4,11,12 . In VAERS, each case report is annotated by certified coders who manually assign individual AEs with the codes of the Medical Dictionary for Regulatory Activities (MedDRA), a standard vocabulary nomenclature 13 . Until now, there are 9,893 MedDRA terms used in VAERS. While MedDRA has played a central role in standardizing vocabulary in the scope of AE reporting, there exist several issues related to the MedDRA usage 14 . First, MedDRA does not provide any term definitions, which may cause confusion. Second, MedDRA has a poorly defined hierarchical structure, making it difficult to use for advanced AE clustering analysis. To improve the AE representation and organization, the community-based Ontology of Adverse Events (OAE) has recently been developed 15 . OAE organizes AE terms in a logical hierarchy based on pathological processes of AE symptoms. In OAE, an AE term denotes a pathological bodily process in a patient that occurs after a medical intervention (e.g., vaccination). Since MedDRA is the standard method for AE representation in VAERS, OAE terms have been mapped to corresponding MedDRA terms to support VAERS AE classification. After the MedDRA-OAE term mapping, the VAERS contents associated with the MedDRA terminology can be analyzed by an OAE-based AE classification method. Compared to MedDRA, OAE has been empirically showed to have superior performance in classifying different AEs associated with live attenuated or killed inactivated influenza vaccines 11 .
Hepatitis A and B are vaccine-preventable epidemic diseases that are caused by highly contagious hepatitis A and B viruses, respectively. The CDC Advisory Committee on Immunization Practices (ACIP) recommends hepatitis A and B vaccinations for all children at age 1 year 16 . For adults, especially healthcare workers, patients with chronic liver disease, injection drug users, prisoners, travelers to endemic areas, and men who have sex with men should also get the vaccines timely 17 . Meanwhile, for the use of the hepatitis A and B vaccines, it is recommended to pay attentions to the risk of post-vaccination AEs, particularly those SAEs such as autoimmune diseases (e.g., multiple sclerosis, arthritis, and myelitis) that have been postulated to be caused by hepatitis vaccines 12 .
GlaxoSmithKline Biologicals manufactures three vaccines, Havrix, Engerix-B, and Twinrix, against hepatitis A and B diseases. Havrix is a monovalent vaccine that contains 1,440 ELISA units (EL.U) of the inactivated hepatitis A virus. Engerix-B is a monovalent vaccine containing 20 μ g of recombinant protein hepatitis B surface antigen (HbsAg). Twinrix is a combination vaccine that is a mixture of the Havrix (half dose) and Engerix-B (full dose), containing 720 EL.U of inactivated hepatitis A virus and 20 μ g of recombinant HbsAg 18 . Given that the three licensed vaccines are all manufactured by the same company and the formulation of Twinrix is a combination vaccine out of the mixture of Havrix and Engerix-B, the investigation of the responses stimulated by the three vaccines provides an ideal use case of VVI studies. Clinical studies have indicated that Twinrix is well tolerated and displays no increased reactogenicity compared to monovalent vaccines Havrix and Engerix-B administered concurrently [18][19][20] . However, it was found that the titers of antibodies to hepatitis A and hepatitis B viruses were significantly higher after administration of the combined or mixed vaccines than after separate injections of the monovalent vaccines 19,20 , suggesting a potential interaction between the two monovalent vaccines. The up-regulated antibody titers were likely due to the increased local production of cytokines and consequent enhancement of macrophage activity that result from the interactions between the components of Havrix and Engerix-B 21 . Nowadays, these monovalent and combination hepatitis A and B vaccines have been used in market more than ten years, and there are many AE case reports associated with these vaccines recorded in VAERS. However, systematic comparison and analysis the AE profiles of these monovalent and combination hepatitis A and B vaccines with large scale post-licensure AE case report data has not been performed yet. The analysis of VAERS AE data related to Havrix, Engerix-B, and Twinrix provides us an ideal case to study the AEs associated with hepatitis A and B vaccines and their combination vaccine.
In this study, we investigated the comparative profiles of AEs associated with Havrix, Engerix-B, and Twinrix. The statistically significant AEs were then classified and analyzed using OAE-based methods. Given Twinrix being a combination vaccine composed of Havrix and Engerix-B, we hypothesized that some Twinrix-associated AEs would be induced by the interaction between Havrix and Engerix-B. To address this hypothesis, we designed and implemented a novel statistical method that combines a logistic regression modeling method with a Markov Chain Monte Carlo (MCMC) sampling 22 .

Methods
Adverse event data extraction. Vaccinees are often administered with more than one vaccine during a short period of time. In this case, it is impossible to associate following AEs with individual vaccines. To improve the accuracy of predicting vaccine-specific AEs, especially AEs associated with the interaction between Havrix and Engerix-B, we only extracted the AE case reports from those VAERS cases where the vaccinees were administered with only one hepatitis vaccine (i.e., Havrix, Engerix-B, or Twinrix).
Our study considered the effects of different variables. Twinrix was approved for use in the USA in persons 18 years of age or older starting on May, 11, 2001 23 . In comparison, Havrix and Engerix-B had already been used for all ages before 2001. To ensure comparability among the VAEs of these three vaccines, we included only VAERS AE reports received from May 2001 to January 2015 for those vaccinees who were aged at least 18 years old.
The VAERS cases filtered based on the above three criteria (age, reporting time, and only one vaccine administered) were used for the following statistical analysis (e.g., PRR) to identify statistically significant AEs. Related Scientific RepoRts | 6:34318 | DOI: 10.1038/srep34318 AE case report details were extracted from the CDC Wonder VAERS data search website (http://wonder.cdc.gov/ controller/datarequest/D8).
Statistical data analysis with PRR, Chi-square (x 2 ) test, and base level filtration. The proportional reporting ratio (PRR) 24 is a main statistical method used in our data analysis. Basically, the VAERS database can be viewed as a contingency table with rows representing the MedDRA coding terms and columns representing the vaccine products. Each cell in the table contains a value that gives the number of reports for the AE of interest and the vaccine of interest. Using the contingency table data structure, PRR calculates the proportion of a specific AE for a vaccine of interest where the comparator is all other vaccines in the VAERS database 24 . As described above, the restriction of vaccinee age (18 years of age or older), case reporting period (May 2001 to January 2015), and single vaccine administration was applied in our data retrieval from the VAERS database. A large PRR score of a specific vaccine AE indicates that the AE has been disproportionately reported for that vaccine, compared with all the other vaccines. We applied the standard PRR score cutoff of 2, i.e., only those AEs with PRR ≥ 2 were further considered 24,25 .
In parallel with the PRR signal detection, we also applied a x 2 test to statistically analyze the likelihood of individual AE terms associated with specific vaccines 26 . The x 2 test calculation also relies on the contingency table described above. An AE is called significant when its x 2 score is greater than 4, which is equivalent to a p-value of approximately 0.05 or smaller 11 . More details about how to compute the scores of PRR and x 2 for each AE in each group based on a 2 × 2 contingency table is provided in Supplementary Table S1.
To filter out background noises effectively, when the total case report number is less than 1,500, we also applied a minimal sample size cutoff of 3 case reports as a base level filtration for each AE to be further considered. When the total case report number is greater than 1,500, the cutoff was set to be 0.2% of total reports for each AE. Such a base level filtration means that at least 2 out of 1,000 cases should report the AE of interest 11,27 . The combined use of PRR and a minimal sample size cutoff is also called screened PRR method (SPRR) 25 .
Note that the MedDRA controlled terminology system used by VAERS contains many terms indicating laboratory test result normal or negative (e.g., 'blood albumin normal' and 'hepatitis B test negative'). These normal or negative laboratory test results are not involved in any AEs. In addition, many MedDRA terms, such as 'lymphocyte percentage' , 'rheumatoid factor' , 'CSF cell count' and 'electroneurography' only represent plain variables or plain test methods and are actually not AEs. Such ambiguous and non-informative MedDRA terms were removed in our AE analysis.

OAE-based vaccine adverse event classification. The hierarchical structures of AE terms in OAE and
MedDRA are largely different. For better comparison and analysis, statistically significant MedDRA terms associated with Havrix, Engerix-B, and Twinrix were mapped to corresponding OAE terms. Through MedDRA-OAE term mapping, the AE terms were classified and analyzed using OAE and MedDRA hierarchical structures 11 . The OntoFox software program was used to automatically extract Havrix-, Engerix-, and Twinrix-specific AE terms, their parent terms, and associated hierarchies in OAE 28 . Note that the OAE version 1.1.316 was used in this study and the whole hierarchical structure of OAE can be viewed in Ontobee (http://www.ontobee.org/ontology/ OAE) 29 . The MedDRA hierarchies of selected AE terms were extracted from the MedDRA version obtained in the BioPortal website (http://bioportal.bioontology.org/ontologies/MEDDRA) on December 22, 2015. The hierarchical results were all visualized and compared using the Protégé-OWL editor 30 . SAE classification. According to FDA, an AE is any undesirable experience associated with the use of a medical product in a patient, and a SAE was defined as an AE at any dose that: (i) results in death, (ii) is life-threatening, (iii) results in hospitalization, (iv) results in disability or permanent damage, (v) is congenital anomaly or birth defect, (vi) requires intervention to prevent permanent impairment or damage, or other serious medical events 31 . In this paper, for AEs associated with the hepatitis vaccines, we classified an AE as a SAE based on the FDA's definition as well as the existing classification in previous publications 4,11,12,32,33 .

Identification of VVIs by developing and applying a new logistic regression model accompanied
with MCMC sampling. To analyze the vaccine-vaccine interactions (VVIs), we developed and applied a VVI-specific logistic regression model. In this method, we modeled the probability of an AE (present/absent) as a function of two vaccines: 1 2 In the above logistic regression, p denotes the probability of AE, and (x 1 , x 2 ) is the design vector for the two vaccines. For example, (1, 0) is for a case to have hepatitis A vaccine, (0, 1) is for hepatitis B vaccine, and (0, 0) is for hepatitis AB combination vaccine. We adopted Markov Chain Monte Carlo (MCMC) sampling to obtain posterior distributions of the model parameters (α , β , γ ) using a MCMClogit function in R 22,34 . Then these parameters (or their functions) were used to estimate probabilities of AEs. Specifically, the probability of AE for hepatitis A (p A ), hepatitis B (p B ), and hepatitis AB (p AB ) are shown in Equation (2~4).
Scientific RepoRts | 6:34318 | DOI: 10.1038/srep34318 For a particular AE, the synergistic effect of vaccine A and B is defined as the probability of an AE for the vaccine combination that is larger than the simple sum of AE probabilities of the two vaccines alone 35,36 . That is, the vaccine A and B has synergistic effect if the fold change (FC) = is large. Specifically, we calculated two posterior probabilities: p FC2 = P(FC > 2) and p FC1 = P(FC < 1). The strong synergistic effect is evidenced by a large p FC2 and a small p FC1 . In this study, we selected AE with a significant synergistic effect for hepatitis A and B vaccines if p FC2 > 0.80 and p FC1 < 0.05.

Results
The general project workflow shown in Fig. 1 outlines different steps in our study and the results out of each step. The details of these analysis processes are provided below.
Extracting AEs associated with Havrix, Engerix-B and Twinrix from VAERS. As 1). From these case reports, we identified 1,093, 2,118, and 1,851 AEs (coded with MedDRA terms) specifically associated with Havrix, Engerix-B, and Twinrix, respectively. Note that in these reported cases, the vaccinees were inoculated with only one hepatitis vaccine (i.e., Havrix, Engerix-B or Twinrix), and no other vaccines were co-administered.   (Table 2), and 82 Twinrix-specific AEs (Table 3) remained (Fig. 1). Based on a Venn diagram analysis ( Fig. 2A) 40 . To compare the AEs generated from the VAERS data and the FDA package insert knowledge, AEs related to the three vaccines represented in the OVAE were also extracted and compared (Fig. 2B). In OVAE, there are 7 AE types associated with Havrix, 4 AEs associated with Twinrix, and 2 AEs associated with Engerix-B (Fig. 2B). Among these AEs, two AE symptoms (i.e., headache and injection site redness) are shared between Havrix and Twinrix, one AE symptom (i.e., fatigue) shared between Engerix-B and Twinrix, and one AE symptom (i.e., injection site muscular soreness) shared among all three vaccines. Overall, compared to the VAERS results, there are less numbers of AEs recorded in FDA package inserts and represented in OVAE, and FDA-recorded AEs associated with these three vaccines are in general mild and self-limited.

AE hierarchical classification based on the OAE method.
After statistically significant AEs were detected, we analyzed these AEs using an OAE-based classification method. Tables 1-3 summarizes the Havrix-, Engerix-B-, and Twinrix-specific AEs after the OAE clustering analysis, respectively. Supplementary Figs S1-S3 records the OAE hierarchies of statistically significant AEs associated with these three vaccines. As shown in Tables 1-3 and Supplementary Figs S1-S3, the most frequently identified AE category is the behavioral and neurological AE category that includes 29 unique AEs for all three vaccines. The other commonly identified categories include immune system AEs and investigation result abnormal AEs. Specifically, for Havrix (Table 1), the AEs included in pregnancy, neonatal or perinatal disorder (e.g., abortion spontaneous and unintended pregnancy) were relatively frequent AEs. Engerix-B was associated with many nervous system AEs (e.g., demyelination, hyperreflexia, and optic neuritis), eye disorders (e.g., visual acuity reduced, double vision, and visual disturbance), and musculoskeletal or connective tissue AEs (e.g., fasciitis, fibrosis tendinous, and myofascitis) ( Table 2). For the combination vaccine Twinrix (Table 3), many commonly identified AEs occurred at the hepatobiliary system (e.g., jaundice and liver disorder), cardiovascular system (e.g., circulatory collapse and cardiovascular disorder), and nervous system (e.g., myelitis, optic neuritis, and polyneuropathy). It is remarkable that the number of AEs associated with Twinrix was more than that with Havrix or Engerix-B.

Thirteen AEs identified out of synergistic VVIs between Havrix and Engerix-B.
In this study, we developed a new statistical method by combining the established logistic regression model with MCMC sampling to identify the hidden interactions among the hepatitis vaccines. For 144 unique statistically significant AEs associated with Havrix, Engerix-B, and Twinrix, 13 AEs were satisfied with the defined threshold (p FC2 > 0.80 and p FC1 < 0.05) ( Table 4). As shown in Table 4, for Havrix, hepatic steatosis was reported only once, and the other 12 AEs were not reported. For Engerix-B, three AEs (i.e., hepatosplenomegaly, premature delivery, and sinus tachycardia) were not reported, and the other 10 AEs were reported but were not within the list of statistically significant AEs (Table 2). However, for Twinrix, each of the 13 AEs was reported for at least 4 times (Table 4), and these 13 AEs were all statistically significantly associated with Twinrix (Table 3). Our statistical analysis results further indicate that there is a high probability that these Twinrix-associated AEs are out of a synergistic interaction between the components of Havrix and Engerix-B. Among the 13 AEs (Table 4), hepatosplenomegaly, premature delivery, and hepatic steatosis belong to SAEs. Based on the OAE classification, the 13 AEs associated with VVIs are mainly involved in behavioral and neurological conditions (i.e., monoparesis, monoplegia, and hypoesthesia facial), immune system (i.e., allergic dermatitis and hepatosplenomegaly) and hepatobiliary condition (i.e., hepatic steatosis, and cholelithiasis) (Fig. 4A,B). These enriched AE categories are similar to the patterns found in the AE and SAE classification (Supplementary Figs S1-S3 and Fig. 3).

Confirmation of the OAE advantaged over MedDRA in AE classification.
A previous empirical study found that OAE provided better AE classification than MedDRA in AE classification 11 . To further confirm the results and illustrate the differences between OAE and MedDRA in their applications in AE classification, we applied both OAE and MedDRA to identify the hierarchical structures of the 13 VVI-associated AEs described above (Fig. 4). Many AE terms may be classified under two or more parent terms. For example, 'liver inflammation AE' is fit under 'inflammation AE' or 'liver AE' (Fig. 4A,B). The approach of asserting more than one parent terms in ontology is called multiple inheritance, which often makes an ontology difficult to manually maintain and update 41 . To avoid multiple inheritances, OAE asserts only one parent term, and allows the other parent term(s) to be obtained automatically by reasoning 15 . In the above example, the 'liver inflammation AE' was asserted as a subclass of 'inflammation AE' . After reasoning (based on internal logical axiom definitions), 'liver inflammation AE' was inferred to be a 'liver AE' as well (Fig. 4B). Such a feature does not exist in MedDRA.
Another difference in OAE and MedDRA exists in their strategies for basic hierarchical construction. As shown in Fig. 4C, MedDRA includes many terms ended with "NEC" (i.e., "not elsewhere classified"), for example, 'faecal abnormalities NEC' , 'neurological disorders NEC' , and 'respiratory disorders NEC' . Such an "NEC" term definition style is arbitrary and ambiguous, often leading to confusion and unclear classification results. For instance, the parent MedDRA term of 'abnormal feces' is 'faecal abnormalities NEC' , which is confusing and logically incorrect. In addition, MedDRA misses obvious parent-child term logic. For example, MedDRA   classifies 'feces pale' and 'abnormal feces' , or 'dermatitis allergic' and 'hypersensitivity' in the same hierarchical levels (Fig. 4C). These are logically incorrect since in reality, a 'feces pale' is a subclass of 'abnormal feces' , and 'dermatitis allergic' is a subclass of 'hypersensitivity' (Fig. 4A).
In summary, this comparative study further confirmed that OAE provides better classification outcomes than MedDRA. Since MedDRA is the default AE reporting terminology in VAERS, our approach of MedDRA-OAE term mapping followed by OAE hierarchy classification proved to be a valid method in VAERS AE studies.

Discussion
To the best of our knowledge, current study is the first to compare and analyze the AEs associated with hepatitis A vaccine Havrix, hepatitis B vaccine Engerix-B, and hepatitis A/B combination vaccine Twinrix. Our statistical and ontological analyses of the VAERS data found that the three hepatitis vaccines were associated with 144 unique AEs (including 29 SAEs), mainly occurring in behavioral and neurological, immune, and investigation result abnormal. We also analyzed the VVIs between the hepatitis A and B vaccines by developing and implementing a new method of logistic regression modeling accompanied with MCMC sampling. Our VVI analysis identified 13 AEs out of the synergistic interaction between hepatitis A and B vaccines. These VVI-associated AEs were mainly involved in behavioral and neurological, immune, and hepatobiliary conditions.
Since the VAERS passive surveillance system has many limitations including underreporting, incomplete information in many reports, and lack of a direct and unbiased comparison group 9,10,42 , direct and naïve usages of the VAERS data may result in wrong assertions of causal relations between vaccines and AEs. Nevertheless, the combinative usage of bioinformatics and statistical methods (e.g., PRR and Chi-square test) to retrieve and analyse the VAERS data can still generate many meaningful and interpretable results and draw sensible hypotheses between vaccines and AEs. For example, the research based on VAERS data by Sirarat et al. suggested that the live attenuated influenza vaccine (LAIV) had lower chance of inducing Guillain-Barre syndrome and paralysis than inactivated influenza vaccine (TIV) 11 . Additionally, VAERS reports of intussusception at 1-2 weeks after rotavirus vaccine administration helped to identify this potentially fatal adverse event 43 .
Our study identified 9 AEs associated with all three hepatitis vaccines, of which, 6 AEs (i.e., hepatitis, liver disorder, jaundice, gamma-glutamyltransferase level increased, transaminases level increased, and blood bilirubin level increased) are hepatitis-associated symptoms. It is often difficult to determine whether these AEs are caused by the hepatitis A and/or B vaccines or other factors. Since the three vaccines contain either inactivated hepatitis A virus and/or noninfectious hepatitis B virus surface antigen, it is impossible for the vaccinees to get infections using the vaccinations. However, before the vaccinations, the vaccinees might be exposed to virulent hepatitis A and B viruses, which cause the occurrence of hepatitis and associated symptoms (e.g., transaminases level increased and jaundice). Many vaccinees are the people who have potential risks for exposure to these two viruses. Since the infection of virulent hepatitis A and B viruses have relatively long incubation periods, such infections might not be detected before the vaccinations 32 .
In general, hepatitis A and/or B vaccines are highly safe, as indicated in FDA-approved vaccine package insert documents [37][38][39] and peer-reviewed studies [18][19][20] . However, our studies still found 29 SAEs following the vaccinations using the VAERS data (Tables 1-3 and Fig. 3). Compared to the spontaneous reported vaccine AE cases in VAERS, the AE records in FDA package inserts were generated from randomized and well-controlled studies and are thus more likely to be causal AEs associated with specific vaccines. The large number of vaccinees reported in VAERS had varied backgrounds (e.g., gender, race, age, and location) and pre-existing health conditions. In contrast, the randomized and well-controlled studies recorded in the package inserts were usually conducted in a small scale of healthy population. Therefore, there are usually less numbers of AEs recorded in FDA package inserts, and the analysis of VAERS data allows the identification of more AEs under special patient backgrounds and conditions. The 29 SAEs identified in our study are enriched in the area of immune system. How these SAEs are related to the hepatitis A and/or B vaccines under various conditions deserves further investigation.
With various data resources and data analysis methods, many groups have identified SAEs associated with hepatitis A vaccines 44,45 , hepatitis B vaccines [46][47][48][49] , and hepatitis A and B combination vaccines 32,50 in humans at different ages. Some of these studies also used VAERS data 32,44,[46][47][48][49] . In general, the findings from our study are consistent with previous results 12,32,46 . Below we focus our discussion on two important areas: autoimmune-related disorder and abortion AEs.
Epidemiological studies and retrospective reviews have shown an association between autoimmunity with hepatitis A and B vaccines 12,51-54 . By evaluating many clinical and laboratory findings of children following hepatitis A vaccination, Karali et al. found that none of the children developed autoimmune disorders although hepatitis A vaccine could induce the production of autoantibodies 52 . A case-control epidemiological study described by Geier et al. showed that hepatitis B vaccination to adults was associated with an increased risk of serious autoimmune adverse events (SAAEs) such as alopecia, thrombocytopenia, lupus erythematosus, and rheumatoid arthritis 12 . By conducting a nested case-control study within the General Practice Research Database (GPRD) in the United Kingdom, Miguel et al. discovered that vaccinees immunized with the hepatitis B vaccine would suffer an increased risk of multiple sclerosis, an autoimmune demyelinating disease 51 . A case report study by Csepregi A et al. suggested that Twinrix led to an acute exacerbation of an unrecognized autoimmune hepatitis 53 . In our VAERS study, no autoimmune-related AE associated with hepatitis A vaccine (Havrix) was identified. This result is consistent with the reviewed findings by Karali et al. 52 that people vaccinated with hepatitis A vaccine are unlikely to develop any autoimmune disorders. For hepatitis B and hepatitis AB vaccines, we identified 5 (i.e., multiple sclerosis, myelitis, polyarthritis, rheumatoid arthritis, and systemic lupus erythematosus) and 7 (i.e., autoimmune thyroiditis, ulcerative colitis, multiple sclerosis, myelitis, psoriasis, rheumatoid arthritis, and systemic lupus erythematosus) autoimmune-related AEs associated with Engerix-B and Twinrix, respectively. These findings echoed previous findings of the associations between hepatitis B or AB vaccines and autoimmune diseases 12,51,53 . However, to date, no clinical evidence exists regarding the causalities between these two types of vaccines and autoimmune diseases. Due to the frequent observations of autoimmune diseases developed after hepatitis B or AB vaccinations, it appears critical to further investigate possible causalities and underlying mechanisms.
Spontaneous abortion (SAB) is the most common pregnancy-specific vaccine AE [55][56][57][58] . From reviewing all VAERS reports for AEs during 1996-2013, Moro et al. did not identify any concerning pattern of AEs in pregnant women or their infants following maternal hepatitis A or hepatitis AB immunizations during pregnancy 44 . In our study, out of 1,624 Twinrix case reports, 5 cases included abortion AE (Table 3). Our further VAERS data investigation found that among the 5 cases of abortion following Twinrix vaccination, only 2 cases were spontaneous abortion, and the other 3 were elective termination (VAERS case IDs: 209240, 233067, and 245750). If we only consider the spontaneous abortion and exclude elective termination, the abortion AE would not be classified as significantly associated with Twinrix. In our study, abortion AE and two more specific abortion AEs (i.e., abortion missed and abortion spontaneous) were identified as statistically significant AEs associated with Havrix (Table 1). Specifically, out of 941 Havrix AE case reports, 4, 3, and 29 cases were reported to have abortion, missed abortion (i.e., silent miscarriage), and spontaneous abortion, respectively ( Table 1). The rate of abortion (36 over 941 or 3.8%) was considered high and calculated as statistically significant when it was compared with the abortion AE associated with all other vaccines. Such observations suggest that hepatitis A vaccine (Havrix) is more likely to induce abortion-related AE than hepatitis B vaccine (Engerix-B), possibly because Engerix-B is composed of a recombinant protein while Havrix uses the whole inactivated virus. To examine the details about these abortion cases, Supplementary Table S2 was generated to lay out the detailed information about these 36 cases. As shown in this supplementary Table, most (33/36) patients were labeled as coming from "foreign", and the gestational age at abortion varied from 4 weeks to 38 weeks. It is likely that the patients vaccinated in "foreign" territories were associated with other unexpected effects.
To our best knowledge, our study represents the first report of the analysis of vaccine-vaccine interactions (VVIs) using clinically reported AE case data. Our signal detection algorithm based on a logistic regression model accompanied with MCMC sampling is also the first to be used for detecting potential synergistic VVI or drug-drug interaction (DDI) effects. Many in silico models have been developed to predict potential DDIs. For example, a heterogeneous network-assisted inference (HNAI) framework was used to support the prediction of DDIs by integrating drug phenotypic, therapeutic, chemical, and genomic properties 59 . Logistic regression models were also used to predict DDIs using drug clinical AE case report data 60,61 . However, these DDI logistic regression models were not accompanied with MCMC sampling (which was used to model fitting) 60,61 . Although already reported in statistics 62,63 , the logistic regression model with MCMC sampling method is new in VVI and DDI studies. In this work, we first applied the logistic regression model accompanied with MCMC sampling method to estimate the synergistic effect of hepatitis A and B vaccines, which relies on a logistic regression model and makes inference based on the posterior distribution of the fold change (probability of AE for the combination over the summation of two vaccines alone). Our statistical analysis identified 13 significant AEs likely associated with VVIs (Table 4). These AEs were not present or weakly present in Havrix and Engerix-B case reports; however, each of them was strongly associated with the combination vaccine Twinrix. The results suggest that the vaccine contents in Havrix and Engerix-B have significant synergistic interactions which likely result in these 13 AEs in Twinrix-vaccinated patients. Further experimental verifications on these VVI-associated AEs would be important to evaluate the safety of the combinational usage of Havrix and Engerix-B.

Conclusions
In this paper, we investigated the AEs associated with two monovalent vaccines (i.e., Havrix and Engerix-B) and one combination vaccine (Twinrix) against hepatitis A and B diseases using the data from the VAERS database. The contributions of our study are multiple. First, by using three biometrical methods (PRR, Chi-square test, and base level filtration), we identified 144 unique statistically significant AEs related to the three vaccines. Among these 144 AEs, 29 were considered serious AEs (SAEs). Many SAEs, including autoimmune-related disorder AEs, deserve further investigation. Second, our study is the first time to statistically evaluate the impact of VVI synergistic effects on the AEs using clinically reported AE case data. Such VVI synergistic effects can increase the risk of some special AEs or exacerbate known adverse reactions. To support the VVI study, we developed a statistical method using logistic regression model with MCMC sampling. Thirteen VVI-associated AEs were identified in our study. Third, we compared the AE classification methods using OAE and MedDRA and confirmed the advantages of using OAE for AE classifications. Overall, our research methods and results facilitate vaccine safety surveillance and benefit rational design of more secure and effective vaccines.