Understanding how interactions among set of genes affect diverse phenotypes is having a greater impact on biomedical research, agriculture and evolutionary biology. Mapping and characterizing the isolated effects of single quantitative trait locus (QTL) is a first step, but we also need to assemble networks of QTLs and define non-additive interactions (epistasis) together with a host of potential environmental modulators. In this article, we present a full-QTL model with which to explore the genetic architecture of complex trait in multiple environments. Our model includes the effects of multiple QTLs, epistasis, QTL-by-environment interactions and epistasis-by-environment interactions. A new mapping strategy, including marker interval selection, detection of marker interval interactions and genome scans, is used to evaluate putative locations of multiple QTLs and their interactions. All the mapping procedures are performed in the framework of mixed linear model that are flexible to model environmental factors regardless of fix or random effects being assumed. An F-statistic based on Henderson method III is used for hypothesis tests. This method is less computationally greedy than corresponding likelihood ratio test. In each of the mapping procedures, permutation testing is exploited to control for genome-wide false positive rate, and model selection is used to reduce ghost peaks in F-statistic profile. Parameters of the full-QTL model are estimated using a Bayesian method via Gibbs sampling. Monte Carlo simulations help define the reliability and efficiency of the method. Two real-world phenotypes (BXD mouse olfactory bulb weight data and rice yield data) are used as exemplars to demonstrate our methods.