You can use the general local formula for the Laplace-Beltrami operator in terms of any local orthonormal frame:
$$\Delta = \sum_{i=1}^n X_i^2 +\mathrm{div}(X_i)X_i$$
where the $X_i$'s are seen as derivations on functions.
You can always find a local frame of vector fields $X_1,\ldots,X_n$ that are divergence-free at a given point $q$. In terms of this frame, the Laplacian at the point $q$ is just a "sum of squares".
Locally, the construction of a local divergence-free, orthonormal frame leads to a system of first order PDEs. The integrability conditions then give a local obstruction.
UPDATE
Unless your manifold is parallelizable you can't find a global orthonormal frame. Still, the above formula works also if the number of vector fields $W_1,\ldots,W_N$ is greater than the dimension of the manifold $n=\dim M \leq N$. To see this practically, pick a orthonormal frame $X_1,\ldots,X_n$ (local on $U \subset M$). We have
$$ W_I = \sum_{j=1}^n A_{Ij} X_i $$
for some smooth family of $N \times n$ matrix $A : U \to M_{N\times n}$. Assume that
$$A^T A = \mathbb{I}_n$$
on $U$. Then you can check that for any function $f$
$$ \sum_{I=1}^N W_I(f) W_I = \sum_{I=1}^N \sum_{i,j=1}^n A_{Ij}A_{Ii} X_j(f) X_i(f) = \sum_{i=1}^n X_i(f) X_i = \nabla f$$
where $\nabla$ is the Riemannian gradient. Then
$$\Delta f = \mathrm{div}(\nabla(f)) = \sum_{I=1}^N W_I^2(f) + \mathrm{div}(W_I)W_I(f) \tag{1}$$
That is the formula at the beginning of my answer. Observe that all of this (starting from the definition of the matrix $A$ is local (since the $X_i$'s) are Local, but clearly formula (1) is true globally.
More abstractly, the initial formula holds true for any set of vector fields $W_1,\ldots,W_N$ (local or global) such that the symbol (as a function on $T^*M$) is written
$$ \lambda \mapsto \frac{1}{2}\sum_{I=1}^N \lambda(W_I)^2, \qquad \lambda \in T^*M$$
THANKS to Jean Van Schaftingen for pointing out an imprecision in my previous answer.