We report computer simulations of a three-dimensional system of semiflexible polymers consisting of hard spherocylinders connected by joints of variable flexibility. In these simulations, we have studied the influence of molecular flexibility on the location of the isotropic-nematic phase transition. A comparison of our numerical results with available theoretical predictions indicates that the existing theories systematically overestimate the density of the coexisting phases. We observe that our simulation data agree well with the available experimental data.