We examine the lepton-specific 2HDM as a solution of muon $g-2$ anomaly under various theoretical and experimental constraints, especially the direct search limits from the LHC and the requirement of a strong first-order phase transition in the early universe. We find that the muon $g-2$ anomaly can be explained in the region of $32<\tan \beta <80$, $10\text{GeV}<m_A<65\text{GeV}$, $260\text{GeV}<m_H<620\text{GeV}$ and $180\text{GeV}<m_{H^{\pm }}<620\text{GeV}$ after imposing the joint constraints from the theory, the precision electroweak data, the 125 GeV Higgs data, the leptonic/semi-hadronic τ decays, the leptonic Z decays and $\text{Br}(B_s\to {\mu }^{+}{\mu }^{-})$. The direct searches from the $h\to AA$ channels can impose stringent upper limits on $\text{Br}(h\to AA)$ and the multi-lepton event searches can sizably reduce the allowed region of $m_A$ and $\tan \beta$ ($10\text{GeV}<m_A<44\text{GeV}$ and $32<\tan \beta <60$). Finally, we find that the model can produce a strong first-order phase transition in the region of $14\text{GeV}<m_A<25\text{GeV}$, $310\text{GeV}<m_H<355\text{GeV}$ and $250\text{GeV}<m_{H^{\pm }}<295\text{GeV}$, allowed by the explanation of the muon $g-2$ anomaly.