Developing a computer-aided diagnostic system for detecting various skin malignancies from images has attracted many researchers. Unlike many machine learning approaches such as Artificial Neural Networks, Genetic Programming (GP) automatically evolves models with flexible representation. GP successfully provides effective solutions using its intrinsic ability to select prominent features (i.e. feature selection) and build new features (i.e. feature construction). Existing approaches have utilized GP to construct new features from the complete set of original features and the set of operators. However, the complete set of features may contain redundant or irrelevant features that do not provide useful information for classification. This study aims to develop a two-stage GP method where the first stage selects prominent features, and the second stage constructs new features from these selected features and operators such as multiplication in a wrapper approach to improve the classification performance. To include local, global, texture, color, and multi-scale image properties of skin images, GP selects and constructs features extracted from local binary patterns and pyramid-structured wavelet decomposition. The accuracy of this GP method is assessed using two real-world skin image datasets captured from the standard camera and specialized instruments, and compared with commonly used classification algorithms, three state-of-the-art, and an existing embedded GP method. The results reveal that this new approach of feature selection and feature construction effectively helps improve the performance of the machine learning classification algorithms. Unlike other black-box models, the evolved models by GP are interpretable, therefore the proposed method can assist dermatologists to identify prominent features, which has been shown by further analysis on the evolved models.