Identification of condition-specific protein interaction subnetworks has emerged as an attractive research field to reveal molecular mechanisms of diseases and provide reliable network biomarkers for disease diagnosis. Several methods have been proposed, which integrate gene expression and protein-protein interaction (PPI) data to identify subnetworks. However, existing methods treat differential expression of genes and network topology independently, which is an oversimplified assumption to model real biological systems. In this paper, we propose a sampling-based subnetwork identification approach to take into account the dependency between gene expression and network topology. Specifically, we apply Markov random field (MRF) theory to model the dependency of genes in PPI network using a Bayesian framework, followed by a Markov Chain Monte Carlo (MCMC) approach to identify significant subnetworks. The MCMC approach estimates the posterior distribution of genes' significant scores and network structure iteratively. Experimental results on both synthetic data and real breast cancer data demonstrated the effectiveness of the proposed method in identifying subnetworks, especially several functionally important, aberrant subnetworks associated with pathways involved in the development and recurrence of breast cancer.