The choice of which stars you want to plot can vary and when you start plotting this many objects, you need to optimize where possible. The code for querying for the data is about as optimized as it can get, but if you want to plot several thousand stars and maintain their approximate color (not taking brightness into account), the following would be an extension of your code. I’ve hard coded the PlotRange to about 2500 light years. You can make out a nearly spherical distribution of stars with the more distant ones flattening out along the galactic plane and being biased towards the blue end. This is because the stars we can see that are that far away are very luminous stars and tend to be blue. Closer to us, there tend to be more low mass red stars. These fainter red stars are harder to see that far away and so gives an observation bias favoring blue stars at further distances.
stars = {AstronomicalData[#, "Name"],
AstronomicalData[#, "RightAscension"],
AstronomicalData[#, "Declination"],
AstronomicalData[#, "DistanceLightYears"],
AstronomicalData[#, "EffectiveTemperature"]} & /@
AstronomicalData["NakedEyeStar"];

goodstars =
Cases[stars, {nam_, ra_?NumberQ, dec_?NumberQ, dist_?NumberQ,
temp_?NumberQ}];

Plotting points is more efficient than spheres and using a single Point primitive with the VertexColors option allows for more efficient rendering.

Graphics3D[{PointSize[0],
Point[equatorialToCartesian[#[[2 ;; 4]]] & /@ goodstars,
VertexColors -> (ColorData["BlackBodySpectrum"][#[[5]]] & /@
goodstars)]}, PlotRange -> 2500, SphericalRegion -> True,
Background -> Black]