Let X n = {x j } n j=1 be a set of n points in the d-cube, and Φ n = {ϕ j } n j=1 a family of n functions on I d . We consider the approximate recovery of functions f onerror of sampling recovery is measured in the norm of the space L q (I d )-norm or the energy quasi-norm of the isotropic Sobolev space W γ q (I d ) for 1 < q < ∞ and γ > 0. Functions f to be recovered are from the unit ball in Besov type spaces of an anisotropic smoothness, in particular, spaces B α,β p,θ of a "hybrid" of mixed smoothness α > 0 and isotropic smoothness β ∈ R, and spaces B a p,θ of a nonuniform mixed smoothness a ∈ R d + . We constructed asymptotically optimal linear sampling algorithms L n (X * n , Φ * n , ·) on special sparse grids X * n and a family Φ * n of linear combinations of integer or half integer translated dilations of tensor products of B-splines. We computed the asymptotic order of the error of the optimal recovery. This construction is based on B-spline quasi-interpolation representations of functions in B α,β p,θ and B a p,θ . As consequences we obtained the asymptotic order of optimal cubature formulas for numerical integration of functions from the unit ball of these Besov type spaces.