Are angles ($\theta$ and $\phi$) in spherical coordinates treated as operators in quantum mechanics? Position is specifically considered as an operator in quantum mechanics. I want to know if $\theta$ and $\phi$ are explicitly considered as operators in quantum mechanics for solutions to 3D Schrodinger equation. Also, how do they commute with each other and with operators like $r$ or $J$?
Also, if they are considered as operators then how do we calculate expectation value for, say, solution to electron wave function in Hydrogen atom potential. For as far as I understand, the potential is spherically symmetric so there should be no preference for any particular angle.
(This seems to be the case for s-orbitals but for p and higher orbitals, the wave function is not spherically symmetric. So what is the origin of this asymmetry?)
 A: For simplicity, let's go down by one dimension, so that wavefunctions are of the form $\psi(x, y)$. You can of course change the wavefunction to polar coordinates, getting $\psi(r, \theta)$, and define
$$\hat{\theta} |\psi \rangle = |\psi' \rangle, \quad \psi'(r, \theta) = \theta \psi(r, \theta).$$
There's been a lot of discussion of the phase operator in the liteature, but the fundamental problem is simply that $\theta$ is only defined up to multiples of $2\pi$, so $\psi(r, \theta) = \psi(r, \theta + 2 \pi)$. In order for the action of the phase operator to be well-defined, you need to try something like
$$\psi'(r, \theta) = (\theta - 2\pi \lfloor \theta/2\pi \rfloor) \psi(r, \theta)$$
which adjusts the phase to lie in $[0, 2 \pi)$. But such an ugly definition screws up the commutation relations. Of course the phase operator still commutes with the radius operator $\hat{r}$, which just multiplies the wavefunction by $r$, but it no longer has canonical commutation relations with the angular momentum operator $\hat{L} = - i \partial / \partial \theta$. Instead, you get nasty additional terms from differentiating the floor function.

Also, if they are considered as operators then how do we calculate expectation value for, say, solution to electron wave function in Hydrogen atom potential.

Under the "ugly" definition I gave above, the expectation value of $\hat{\theta}$ for any rotationally symmetric state would be $\pi$, since it's in the middle of the interval $[0, 2\pi)$. But this depends on the interval you choose, which seems artificial. That's again just showing how the idea of a phase operator becomes ugly, when the phase can wrap all the way around the circle, and thus why we don't often use it.

For as far as I understand, the potential is spherically symmetric so there should be no preference for any particular angle. (This seems to be the case for s-orbitals but for p and higher orbitals, the wave function is not spherically symmetric. So what is the origin of this asymmetry?)

That's actually a completely separate question. It's like how space has no preferred direction, but if we describe space in Cartesian coordinates, then there suddenly are preferred directions, namely those of the coordinate axes. But if we're consistent, then at the end of a calculation all the dependence on the choice of axes will drop out, leaving a rotationally symmetric result.
Similarly, while an individual $p$-orbital $|p_i \rangle$ is not rotationally symmetric, the subspace spanned by the $p$-orbitals $\{|p_x \rangle, |p_y \rangle, |p_z \rangle \}$ is, just like how space in Cartesian coordinates is.
A: You are probably really asking about  spherical harmonics,  ideally specified in a fine text, such as Modern Quantum Mechanics by Sakurai and Napolitano, for instance.  You don't treat these angles as operators, as you normally work in their coordinate representation, where they are pure numerical angle variables (periodic ones).
To avoid confusion with the $\hat \bullet $ notation, let us keep using that for unit vectors and use Capitals for operators, instead, and lower case for numerical variables, such as $r,\theta,\phi$.
In the spherical coordinate representation,
\begin{align}
  \mathbf L &= i \hbar \left(\frac{\hat{\boldsymbol{\theta}}}{\sin(\theta)} \frac{\partial}{\partial\phi} - \hat{\boldsymbol{\phi}} \frac{\partial}{\partial\theta}\right), \\      
        L^2 &= -\hbar^2 \left(\frac{1}{\sin(\theta)} \frac{\partial}{\partial\theta} \left(\sin(\theta) \frac{\partial}{\partial\theta}\right) + \frac{1}{\sin^2(\theta)}\frac{\partial^2}{\partial\phi^2}\right), \\
        L_z &= -i \hbar \frac{\partial}{\partial\phi}~~~.
\end{align}
Atomic wavefunctions are of the form $\langle  x,y,z|n,l,m\rangle=R_{nl}(r) Y^m_l(\theta, \phi) $, and we normally use the normalization $\langle \hat r|l,m\rangle=\langle \theta,\phi|l,m\rangle =Y^m_l(\theta,\phi)$, periodic functions. Consequently, using the standard "abuse" of notation for operator action on coordinate representations,
$$
\langle \hat r|L^2|l,m\rangle= \hbar^2 l(l+1) Y^m_l(\theta,\phi) \\
\\=  -\hbar^2 \left(\frac{1}{\sin(\theta)} \frac{\partial}{\partial\theta} \left(\sin(\theta) \frac{\partial}{\partial\theta}\right) + \frac{1}{\sin^2(\theta)}\frac{\partial^2}{\partial\phi^2}\right)  Y^m_l(\theta,\phi),\\
\langle \hat r|L_z|l,m\rangle=\hbar m Y^m_l(\theta,\phi)=  -i \hbar \frac{\partial}{\partial\phi} Y^m_l(\theta,\phi).
$$
To address your questions then, the notional spherical operators you are talking about, $R,\Theta,\Phi$, have continuous c-number eigenvalues $r,\theta,\phi$, in principle,  but they need not be used explicitly. When acting on the spherical coordinate representation wavefunctions we considered, they commute with each other and produce $\Theta Y^m_l(\theta,\phi)= \theta Y^m_l(\theta,\phi)$, $\Phi Y^m_l(\theta,\phi)= \phi Y^m_l(\theta,\phi)$, so they are out of the picture---taking away troubled formal underpinnings with them.

*

*They most emphatically do not commute with $L^2$ and (for $\Phi$) $L_z$, in general (except for $Y^0_0$), as evident above!

I also wouldn't spend any time on the gradients $-i\hbar(\partial_r,\partial_\theta, \partial_\phi)$ which  underlie the respective canonical momentum operators, as you again never need use them, in practice.
The spherical harmonics $Y^m_l$ obey well-known orthogonality relations, where you must recall the angular measure involves a factor of $\sin\theta$. What expectation value do you have in mind? Typically, $\langle l,m| \Theta|l,m\rangle= 2\pi \int_0^\pi \!\!d\theta ~(\sin\theta) ~~\theta ~~|Y^m_l|^2$.

NB Response to comment question.
The $Y^m_l$ are eigenfunctions of $L^2$ and $\Theta $, but $\theta Y^m_l$ is an eigenfunction of $\Theta$ but not $L^2$, in general, so $[L^2,\Theta]\neq 0$. This is because θ  is both an eigenvalue and a variable, so on the $Y^m_l$ basis where $L^2$ is represented by the above double gradient,
$$
[L^2,\theta] Y^m_l= L^2 (\theta Y^m_l) - \theta \hbar^2 l(l+1) Y^m_l\neq 0.
$$
To bypass unfamiliar features of this formalism, you may revert to the harmonic oscillator, with rationalized unit hamiltonian H, quadratic, $2H= P^2+X^2$, so that $[H,X]\neq 0$. Now, the Hermite functions $\psi_n(x)=\langle x|\psi_n\rangle$ are eigenfunctions of the  Hamiltonian, $H\psi_n(x)=  E_n \psi_n(x)$, and, of course, X. But the functions $X\psi_n(x)=x\psi_n(x)$ are not eigenfunctions of the Hamiltonian, anymore, and so
$$
\tfrac{1}{2} [-\hbar^2\partial_x^2 + x^2,x] \psi_n(x) =-\hbar^2\partial_x \psi_n(x) \neq 0~.
$$
You see then that the set of all functions of x, all eigenfunctions of X, is much larger than the set of eigenfunctions of H, and the two operators fail to commute.
A: Yo can define such an operator. Let $\{|\mathbf{r}(r,\theta,\varphi)\rangle\}$ be the eigenstates of the position operator $\mathbf{r}$:
$$
\mathbf{r}|\mathbf{r}(r,\theta,\varphi)\rangle = \mathbf{r}(r,\theta,\varphi) |\mathbf{r}(r,\theta,\varphi)\rangle.
$$
Let us try the following definition of an "angular-operator" $\hat{\theta}$:
$$
\hat{\theta}|\mathbf{r}(r,\theta,\varphi)\rangle  = \theta |\mathbf{r}(r,\theta,\varphi)\rangle .
$$
But you will definitely agree with $|\mathbf{r}(r,\theta,\varphi)\rangle = |\mathbf{r}(r,\theta+2 \pi,\varphi)\rangle$. But then we have
$$
\hat{\theta}|\mathbf{r}(r,\theta,\varphi)\rangle  = (\theta+ n 2 \pi) |\mathbf{r}(r,\theta,\varphi)\rangle , \quad n \in \mathbb N_0.
$$
Hence, this operator, or the action of this operator, is not well-defined. An alternative would be a "modulo" definition like this:
$$
\hat{\theta}|\mathbf{r}(r,\theta,\varphi)\rangle  = (\theta + \mathrm{mod}\, 2\pi) |\mathbf{r}(r,\theta,\varphi)\rangle.
$$
:-)
