In this work we present an adaptive boundary element method for computing the electromagnetic response of wave interactions in hyperbolic metamaterials. One unique feature of hyperbolic metamaterial is the strongly directional wave in its propagating cone, which induces sharp transition for the solution of the integral equation across the cone boundary when wave starts to decay or grow exponentially. In order to avoid a global refined mesh over the whole boundary, we employ a two-level a posteriori error estimator and an adaptive mesh refinement procedure to resolve the singularity locally for the solution of the integral equation. Such an adaptive procedure allows for the reduction of the degree of freedom significantly for the integral equation solver while achieving desired accuracy for the solution. In addition, to resolve the fast transition of the fundamental solution and its derivatives accurately across the propagation cone boundary, adaptive numerical quadrature rules are applied to evaluate the integrals for the stiff matrices. Finally, in order to formulate the integral equations over the boundary, we also derive the limits of layer potentials and their derivatives in the hyperbolic media when the target points approach the boundary.