The lowest-lying glueballs are investigated in lattice QCD using $N_f=2$ clover Wilson fermion on anisotropic lattices. We simulate at two different and relatively heavy quark masses, corresponding to physical pion mass of $m_pisim 938$ MeV and $650$ MeV. The quark mass dependence of the glueball masses have not been investigated in the present study. Only the gluonic operators built from Wilson loops are utilized in calculating the corresponding correlation functions. In the tensor channel, we obtain the ground state mass to be 2.363(39) GeV and 2.384(67) GeV at $m_pisim 938$ MeV and $650$ MeV, respectively. In the pseudoscalar channel, when using the gluonic operator whose continuum limit has the form of $epsilon_{ijk}TrB_iD_jB_k$, we obtain the ground state mass to be 2.573(55) GeV and 2.585(65) GeV at the two pion masses. These results are compatible with the corresponding results in the quenched approximation. In contrast, if we use the topological charge density as field operators for the pseudoscalar, the masses of the lowest state are much lighter (around 1GeV) and compatible with the expected masses of the flavor singlet $qbar{q}$ meson. This indicates that the operator $epsilon_{ijk}TrB_iD_jB_k$ and the topological charge density couple rather differently to the glueball states and $qbar{q}$ mesons. The observation of the light flavor singlet pseudoscalar meson can be viewed as the manifestation of effects of dynamical quarks. In the scalar channel, the ground state masses extracted from the correlation functions of gluonic operators are determined to be around 1.4-1.5 GeV, which is close to the ground state masses from the correlation functions of the quark bilinear operators. In all cases, the mixing between glueballs and conventional mesons remains to be further clarified in the future.