In this paper, compressible, stationary Navier-Stokes equations are considered. The model is well-posed, and there exist weak solutions in bounded domains, subject to inhomogeneous boundary conditions. The shape sensitivity analysis is performed for Navier-Stokes boundary value problems in the framework of small perturbations of the so-called approximate solutions. The approximate solutions are determined from the Stokes problem, and the small perturbations are given by the unique solutions to the full nonlinear model. The differentiability of small perturbations of the approximate solutions with respect to the coefficients of differential operators implies the shape differentiability of the drag functional. The shape gradient of the drag functional is derived in a form convenient for computations, and an appropriate adjoint state is introduced to this end. The shape derivatives of solutions to the Navier-Stokes equations are given by smooth functions; however, the shape differentiability of the solutions is shown in a weak norm. The proposed method of shape sensitivity analysis is general. The differentiability of solutions to the Navier-Stokes equations with respect to the data leads to the first order necessary conditions for a broad class of optimization problems. The boundary shape gradient as well as the boundary value problems for the shape derivatives of solutions to state equations and the adjoint state equations are obtained in the form of singular limits of volume integrals. This method of shape sensitivity analysis seems to be new and is appropriate for nonlinear problems. It is an important contribution in the field of numerical methods of shape optimization in fluid mechanics.