I found out that the problem was that isosurface does not accept irregular x,y and z coordinates. This is what I found to be the solution. I used scatteredInterpolant to interpolate the density and then created my own meshgrid.
input = importdata( 'denNa20.dat' );
input(1,:)=[];
[x,y,z] = sph2cart(input(:,3),pi/2 - input(:,2),input(:,1));
density = input(:,4);
F = scatteredInterpolant(x,y,z,density);
xRange = linspace(min(x),max(x),75);
yRange = linspace(min(y),max(y),75);
zRange = linspace(min(z),max(z),75);
[x,y,z] = meshgrid(xRange,yRange,zRange);
density = F(x,y,z);
axis([-5000,5000,-5000,5000,-5000,5000])
p = patch(isosurface(x,y,z,density,0.00001));
isonormals(x,y,z,density,p);
set(p,'FaceColor','red','EdgeColor','none','FaceAlpha',0.50);
daspect([1 1 1])
view(3)
grid on
camlight; lighting phong