Let $W^D$ be a killed Brownian motion in a domain $D\subset {\mathbb R}^d$ and $S$ an independent subordinator with Laplace exponent $φ$. The process $Y^D$ defined by $Y^D_t=W^D_{S_t}$ is called a subordinate killed Brownian motion. It is a Hunt process with infinitesimal generator $-φ(-Δ|_D)$, where $Δ|_D$ is the Dirichlet Laplacian. In this paper we study the potential theory of $Y^D$ under a weak scaling condition on the derivative of $φ$. We first show that non-negative harmonic functions of $Y^D$ satisfy the scale invariant Harnack inequality. Subsequently we prove two types of scale invariant boundary Harnack principles with explicit decay rates for non-negative harmonic functions of $Y^D$. The first boundary Harnack principle deals with a $C^{1,1}$ domain $D$ and non-negative functions which are harmonic near the boundary of $D$, while the second one is for a more general domain $D$ and non-negative functions which are harmonic near the boundary of an interior open subset of $D$. The obtained decay rates are not the same, reflecting different boundary and interior behaviors of $Y^D$.