We discuss massive quark effects in the endpoint region $x \to 1$ of inclusive deep inelastic scattering, where the hadronic final state is collimated and thus represents a jet. In this regime heavy quark pairs are generated via secondary radiation, i.e. due to a gluon splitting in light quark initiated contributions starting at $\mathcal{O}(\alpha_s^2)$ in the fixed-order expansion. Based on the factorization framework for massless quarks in Soft Collinear Effective Theory (SCET), we construct a variable flavor number scheme that deals with arbitrary hierarchies between the mass scale and the kinematic scales exhibiting a continuous behavior between the massless limit for very light quarks and the decoupling limit for very heavy quarks. We show that the threshold matching corrections for all gauge invariant components at the mass scale are related to each other via consistency conditions. This is explicitly demonstrated by recalculating the known threshold correction for the parton distribution function at $\mathcal{O}(\alpha_s^2 C_F T_F)$ within SCET. The latter contains large rapidity logarithms $\sim \ln(1-x)$ that can be summed by exponentiation. Their coefficients are universal which can be used to obtain potentially relevant higher order results for generic threshold corrections at colliders from computations in deep inelastic scattering. In particular, we extract the $\mathcal{O}(\alpha_s^3)$ threshold correction multiplied by a single rapidity logarithm from results obtained earlier.