Cancer remains a leading cause of death worldwide, with its heterogeneous and patient-specific nature complicating treatment. Immunotherapy, especially immune checkpoint blockades (ICB), shows promising potential by reactivating the immune system to detect and destroy cancer cells. However, their effectiveness largely depends on the patient’s clinical and molecular characteristics, making them unsuitable for everyone. This highlights the importance of developing predictive and interpretable models that could distinguish potential responders from non-responders. This study aims to utilize gene expression data to predict patients’ responses to ICB and identify key genes whose expression is associated with treatment outcome using probabilistic Bayesian networks (BNs). To address the high-dimensional nature of gene expression data, we propose and compare the use of machine learning-based feature selection (FS) and clustering by weighted gene co-expression network analysis (WGCNA). To our knowledge, this is the first study to apply BNs along with FS and WGCNA dimensionality reduction techniques to tackle ICB response using transcriptomic data. We achieved a maximal classification accuracy of 75% and identified eight genes and two gene clusters directly influencing ICB response. These results may offer more biological insights and mechanisms underlying ICB response. The proposed pipeline will contribute to guiding treatment decisions, minimizing harmful side effects, and reducing the financial costs associated with ineffective treatments.