Context. Multi-messenger observations of binary neutron star mergers can provide information on the neutron star’s equation of state (EOS) above the nuclear saturation density by directly constraining the mass-radius diagram. Aims. We present a Bayesian framework for joint and coherent analyses of multi-messenger binary neutron star signals. As a first application, we analyze the gravitational-wave GW170817 and the kilonova (kN) AT2017gfo data. These results are then combined with the most recent X-ray pulsar analyses of PSR J0030+0451 and PSR J0740+6620 to obtain new EOS constraints. Methods. We extend the bajes infrastructure with a joint likelihood for multiple datasets, support for various semi-analytical kN models, and numerical-relativity (NR)-informed relations for the mass ejecta, as well as a technique to include and marginalize over modeling uncertainties. The analysis of GW170817 used the TEOBResumS effective-one-body waveform template to model the gravitational-wave signal. The analysis of AT2017gfo used a baseline multicomponent spherically symmetric model for the kN light curves. Various constraints on the mass-radius diagram and neutron star properties were then obtained by resampling over a set of ten million parameterized EOSs, which was built under minimal assumptions (general relativity and causality). Results. We find that a joint and coherent approach improves the inference of the extrinsic parameters (distance) and, among the intrinsic parameters, the mass ratio. The inclusion of NR-informed relations marks a strong improvement over the case in which an agnostic prior is used on the intrinsic parameters. Comparing Bayes factors, we find that the two observations are better explained by the common source hypothesis only by assuming NR-informed relations. These relations break some of the degeneracies in the employed kN models. The EOS inference folding-in PSR J0952-0607 minimum-maximum mass, PSR J0030+0451 and PSR J0740+6620 data constrains, among other quantities, the neutron star radius to R1.4TOV = 12.30− 0.56+ 0.81 km(R1.4TOV = 13.20− 0.90+ 0.91 km) and the maximum mass to MmaxTOV = 2.28− 0.17+ 0.25M⊙(MmaxTOV = 2.32− 0.19+ 0.30M⊙), where the ST+PDT (PDT-U) analysis of Vinciguerra et al. (2024, ApJ, 961, 62) for PSR J0030+0451 was employed. Hence, the systematics on the PSR J0030+0451 data reduction currently dominate the mass-radius diagram constraints. Conclusions. We conclude that bajes delivers robust analyses in line with other state-of-the-art results in the literature. Strong EOS constraints are provided by pulsars observations, albeit with large systematics in some cases. Current gravitational-wave constraints are compatible with pulsar constraints and can further improve the latter.