Adhesively bonding has been increasingly used for numerous industrial applications to meet the high demand for lightweight and safer structures. Debonding of adhesively bonded joints is a typical mixed mode failure process. It is highly depended on the interactional effects of material properties and geometric definitions of the constituents, which is very complicated. The existing studies in identifying fracture modes of joints based on either experiments or finite element analysis are often prohibitively time and computational expensive. This paper proposed an innovate method by combining Finite Element Analysis (FEA), Latin Hypercube Sampling (LHS) and Genetic Programming (GP) to understand the effect of the physical attributes on the fracture modes of adhesively single lap joints. A dataset of 150 adhesive joint samples has been generated using LHS, including different combinations of adherend and adhesive’s material properties and thicknesses. The mixed mode ratios of the 150 samples are calculated using Strain Energy Release Rate (SERR) outputs embedded in Linear Elastic Fracture Mechanics (LEFM), which has been validated by experimental tests. Finally, a GP model is developed and trained to provide an extracted explicit expression used for evaluating the early-state failure modes of the adhesively bonded joints against the design variables.