This article deals with the study of the following singular quasilinear equation: $$\begin{aligned} (P) \left\{ \ -\Delta _{p}u -\Delta _{q}u = f(x) u^{-\delta },\; u>0 \text { in }\; \Omega ; \; u=0 \text { on } \partial \Omega , \right. \end{aligned}$$where \(\Omega \) is a bounded domain in \({\mathbb {R}}^N\) with \(C^2\) boundary \(\partial \Omega \), \(1< q< p<\infty \), \(\delta >0\) and \(f\in L^\infty _{loc}(\Omega )\) is a non-negative function which behaves like \(\text {dist}(x,\partial \Omega )^{-\beta },\) \(\beta \ge 0\) near the boundary of \(\Omega \). We prove the existence of a weak solution in \(W^{1,p}_{loc}(\Omega )\) and its behaviour near the boundary for \(\beta <p\). Consequently, we obtain optimal Sobolev regularity of weak solutions. By establishing the comparison principle, we prove the uniqueness of weak solution for the case \(\beta <2-\frac{1}{p}\). Subsequently, for the case \(\beta \ge p\), we prove the non-existence result. Moreover, we prove Hölder regularity of the gradient of weak solution to a more general class of quasilinear equations involving singular nonlinearity as well as lower order terms (see (1.6)). This result is completely new and of independent interest. In addition to this, we prove Hölder regularity of minimal weak solutions of (P) for the case \(\beta +\delta \ge 1\) that has not been fully answered in former contributions even for p-Laplace operators.