We use the Negele-Vautherin density matrix expansion to derive a quasi-local density functional for the description of systems of fermions interacting with short-ranged interactions composed of arbitrary finite-range central, spin-orbit, and tensor components. Terms that are absent in the original Negele-Vautherin approach owing to the angle averaging of the density matrix are fixed by employing a gauge invariance condition. We obtain the Kohn-Sham interaction energies in all spin-isospin channels, including the exchange terms, expressed as functions of the local densities and their derivatives up to second (next to leading) order. We illustrate the method by determining the coupling constants of the Skyrme functional or Skyrme force that correspond to the finite-range Gogny central force. The resulting self-consistent solutions reproduce the Gogny-force binding energies and radii within the precision of 1-2%.