We present the halo number counts and its two-point statistics, the observable angular power spectrum, extracted for the first time from relativistic N-body simulations. The halo catalogues used in this work are built from the relativistic N-body code gevolution, and the observed redshift and angular positions of the sources are computed using a non-perturbative ray-tracing method, which includes all relativistic scalar contributions to the number counts. We investigate the validity and limitations of the linear bias prescription to describe our simulated power spectra. In particular, we assess the consistency of different bias measurements on large scales, and we estimate up to which scales a linear bias is accurate in modelling the data, within the statistical errors. We then test a second-order perturbative bias expansion for the angular statistics, on a range of redshifts and scales previously unexplored in this context, that is 0.4 \le \bar{z} \le 2 up to scales \ell_\mathrm{max} \sim 1000. We find that the angular power spectra at equal redshift can be modelled with high accuracy with a minimal extension of the number of bias parameters, that is using a two-parameter model comprising linear bias and tidal bias. We show that this model performs significantly better than a model without tidal bias but with quadratic bias as extra degree of freedom, and that the latter is inaccurate at \bar{z} \ge 0.7. Finally, we extract from our simulations the cross-correlation of halo number counts and lensing convergence. We show that the estimate of the linear bias from this cross-correlation is consistent with the measurements based on the clustering statistics alone, and that it is crucial to take into account the effect of magnification in the halo number counts to avoid systematic shifts in the computed bias.
