Usually, a mesh with "values" at the vertices is the result of some Galerkin scheme or interpolation.
I say "values", because those are really degrees of freedom in a basis of functions. It just so happens that a very commonly used basis of functions is the Lagrange basis with the property that the coefficients in the basis are evaluations of the function. This is not true of all bases. Take the Bernstein (Bézier) basis, for instance. In this case, you'd have to compute linear combinations of the DoFs to obtain values at the vertices. Another case is with discontinuous representations, where solutions are described as polynomial over each element but only $L^2$ over the entire domain (e.g. Discontinuous Galerkin schemes); this means values at the vertices are multiple (one per adjacent element). Yet another example is Hermite FE where some DoFs relate to gradients directly.
What this means is that the solution field over the mesh is partial information: the other half is knowing which basis $\{\phi_i\}_{i=1}^{N_P}$, where $N_P$ is the number of DoFs, is used. Then your solution field writes
$f_h = \sum f_i \phi_i$
where $f_i$ are the values supplied by the discrete solution field. Now, to compute partial derivatives, you simply differentiate the basis functions:
$\partial_j f_h = \sum_i f_i \partial_j \phi_i$
By convention, the basis is the Lagrange basis. This is by far the most common one. Only if another basis is used would a person note it explicitely. The answer by @ConvexHull gives you the expression in that case.
Gradients at vertices
These bases are defined element-by-element and, if you only have DoFs at the vertices, are bases of affine functions. As such, derivatives are discontinuous accross elements. This means you cannot compute the derivative at a vertex.
To solve this, there are several approaches. One is to enrich the FE space and resolve the new degrees of freedom by $L^2$ projection. The end result, for the gradient at a point $P$ belonging to elements $\mathcal{B}(P)$, is
$\nabla f_{h^+}(P) = \frac{\sum_{K\in\mathcal{B}(P)} |K| \nabla (f_h)_{|K}}{\sum_{K\in\mathcal{B}(P)} |K|}$
with $|K|$ the element measure. As often, an intimidating theoretical setting leads to the mundane: a weighted average. One thing of note is that the gradient at $P$ then involves values at all points belonging to adjacent elements (degree 1 neighbourhood). You can also reconstruct the Hessian, in which case you'll involve the degree 2 neighbourhood (up to neighbours of adjacent elements).
For a more Finite-Differences-friendly approach, you could also approximate the gradient by least-squares (or other norm) regression. What is the gradient if not the linear component of the Taylor expansion of $f$ at $P$? Then seek $f$ of the form
$f(x) = f(P) + b^T (x-P)$
and minimize the energy
$\sum_{i\ \text{s.t.} \ P_i \ \text{neighbour of} \ P}||f(x_i) - f_i||^2$
You could also fit the more lax $f(x) = a + b^T (x-P)$ to possibly improve the gradient, then add the term $||f_{i_0} - a||^2$ with $f_{i_0}$ the DoF at $P$.