The goal of this work is to develop a general theory for non-local singular operators of the type $$ L^{\mathcal{B}}_αf(x)=\lim_{ε\to 0} \int_{D,\, |y-x|>ε}\big(f(y)-f(x)\big) \mathcal{B}(x,y)|x-y|^{-d-α}\,dy, $$ and $$ L f(x)=L^{\mathcal{B}}_αf(x) - κ(x) f(x), $$ in case $D$ is a $C^{1,1}$ open set in $\mathbb{R}^d$, $d\ge 2$. The function $\mathcal{B}(x,y)$ above may vanish at the boundary of $D$, and the killing potential $κ$ may be subcritical or critical. From a probabilistic point of view we study the reflected process on the closure $\overline{D}$ with infinitesimal generator $L^{\mathcal{B}}_α$, and its part process on $D$ obtained by either killing at the boundary $\partial D$, or by killing via the killing potential $κ(x)$. The general theory developed in this work (i) contains subordinate killed stable processes in $C^{1,1}$ open sets as a special case, (ii) covers the case when $\mathcal{B}(x,y)$ is bounded between two positive constants and is well approximated by certain Hölder continuous functions, and (iii) extends the main results known for the half-space in $\mathbb{R}^d$. The main results of the work are the boundary Harnack principle and its possible failure, and sharp two-sided Green function estimates. Our results on the boundary Harnack principle completely cover the corresponding earlier results in the case of half-space. Our Green function estimates extend the corresponding earlier estimates in the case of half-space to bounded $C^{1, 1}$ open sets.