Probabilistic Boolean Networks have been proposed for estimating the behaviour of dynamical systems as they combine rule-based modelling with uncertainty principles. Inferring PBNs directly from gene data is challenging however, especially when data is costly to collect and/or noisy, e.g., in the case of gene expression profile data. In this paper, we present a reproducible method for inferring PBNs directly from real gene expression data measurements taken when the system was at a steady state. The steady-state dynamics of PBNs is of special interest in the analysis of biological machinery. The proposed approach does not rely on reconstructing the state evolution of the network, which is computationally intractable for larger networks. We demonstrate the method on samples of real gene expression profiling data from a well-known study on metastatic melanoma. The pipeline is implemented using Python and we make it publicly available.