In the past two decades, the biologists are able to identify gene signatures associated with the phenotypes through monitoring gene expressions using high-throughput biotechnologies. The gene signatures have been successfully applied to drug development, disease prevention, crop improvement, etc. However, ignoring the interactions among genes has weaken the prediction power of gene signatures in real applications. The gene regulatory network, in which genes are present by nodes and the associations between genes are present by edges, are typically constructed to analyze and visualize the gene interactions. Particularly, we proposed to measure the strength of (direct or indirect) associations by the coefficient of intrinsic dependence (CID) to capture possible nonlinear gene relationships. While encountering pathways analysis in a larger scale, the stepwise gene (variable) selection may help to identify relevant genes with correct order from upstream to downstream in a pathway. In this study, we propose to perform the stepwise pathway analysis on microarray expression data using the CID along with the partial coefficient of intrinsic dependence (pCID). The proposed method aims to reduce the high false-positive rates using the CID along in stepwise variable selections. The method was examined using the simulated networks, and the well known CBF-COR pathway under cold stress using Arabidopsis microarray data. It was also practiced on construction of bHLH gene regulatory pathways under abiotic stresses using rice microarray data. The proposed method can efficiently decipher the gene regulatory pathways and achieve higher prediction power in real applications.