The easiest way to see how you can recover the powers of $\hbar$ is by looking at the path integral, which (in natural units) is given by
\begin{align}
\mathcal{Z}=\int \left[\mathcal{D}\phi\right]e^{iS}\end{align}
with the action $S$. Going back to dimensionful units, $S$ will have dimension of an action, so the same as $\hbar$. To make the argument of the exponential dimensionless, we therefore need the path integral to be
\begin{align}
\mathcal{Z} = \int \left[\mathcal{D}\phi\right] e^{\frac{i}{\hbar}S},\end{align}
so the action is rescaled by $\frac{1}{\hbar}$. This means that every internal propagator gets a factor of $\hbar$ (because it is the inverse of the kinetic term) and every vertex gets a factor of $\frac{1}{\hbar}$. You can probably continue from here on your own.
One thing to note is, that $D_F$ and $S_F$ don't need to have the same dimensions, as one of those has the dimensions of a bosonic field squared and the other those of a fermion field squared (remember that they are given as expectation values of time-ordered field products).
You've got 5 base dimensions here, so you can associate any physical dimension with a $5$-tuple of real numbers. For example, we could write
$$M^\alpha L^\beta T^\gamma I^\delta \Theta^\epsilon\rightarrow \pmatrix{\alpha\\\beta\\\gamma\\\delta\\\epsilon}\in \mathbb R^5$$
When two quantities are multiplied or divided, the $5$-tuples representing their dimensions are added or subtracted. Similarly when a quantity is raised to a power, the corresponding $5$-tuple is multiplied by that power. For instance,
$$L/T \rightarrow\pmatrix{0\\1\\-1\\0\\0}, \qquad ML/T^2 \rightarrow \pmatrix{1\\1\\-2\\0\\0}$$
In this sense, physical dimensions can be understood as elements of the vector space $\mathbb R^5$. In order to come up with a system of "natural units," you just need five quantities whose dimensions are linearly independent of one another.
Your example fails because you only provide four - their dimensions obviously don't span $\mathbb R^5$, so there will be some dimensions which cannot be expressed as a combination of them. You can make it work by choosing one additional constant whose dimension is linearly independent of the four you have already provided.
If we use the four physical constants I give in my question, plus the Coulomb constant K (dimension M¹ L³ T⁻⁴ I⁻² ϴ⁰) to get a full set of five, the resulting system of equations is unsolvable, but I don't see how any of the five constants are linearly dependent with another.
Let's make the following identifications
$$[c]=\pmatrix{0\\1\\-1\\0\\0}\equiv \vec v_1, \quad [\hbar]=\pmatrix{1\\2\\-1\\0\\0}\equiv \vec v_2, \quad [e]=\pmatrix{0\\0\\1\\1\\0}\equiv \vec v_3$$
$$[k_B]=\pmatrix{1\\2\\2\\0\\-1}\equiv \vec v_4, \quad [k_C]=\pmatrix{1\\3\\-4\\-2\\0}\equiv \vec v_5$$
The series of equations obtained by writing $\vec v_5 = A\vec v_1+B\vec v_2 + C\vec v_3 + D\vec v_4$ is
$$\matrix{B+D=1\\A+2B+2D=3\\-A-B+C+2D = -4\\C = -2 \\-D = 0}$$
This set of equations can be solved easily to yield $(A,B,C,D)=(1,1,-2,0)$, implying that $\vec v_5 = \vec v_1+\vec v_2-2\vec v_3$. In other words, the dimensions of $k_C$ are equal to the dimensions of $c\hbar/e^2$ - namely
$$\left[\frac{c\hbar}{e^2}\right] = \frac{L}{T}\frac{ML^2}{T} \frac{1}{I^2T^2} = \frac{ML^3}{T^4I^2} = \left[k_C\right]$$
Best Answer
For reference, here is table 2 from page 4 in reference [$1$]:
Yes. The source of confusion here seems to be a misunderstanding of what the first line in table 2 means. It means that in a system of units in which the speed of light ($c$) and Newton's constant ($G$) are both equal to $1$, the SI unit "$1$ meter" can also be expressed as a number of kilograms, specifically $1.3466\times 10^{27}$ kg. The table is not starting with $1$ kg and then ending up with $1.3466\times 10^{27}$ kg. Instead, it is showing how to express "$1$ meter" in kilograms. Here's the explicit calculation, using $c=G=1$, keeping only two significant digits for simplicity: \begin{align} 1\text{ meter} &= 1\text{ m }\times\frac{c^2}{G} \\ &\approx 1\text{ m }\times\frac{ (3.0\times 10^8\text{ m/s})^2 }{ 6.7\times 10^{-11}\text{ m}^3/(\text{kg}\cdot\text{s}^2) } \\ &\approx 1\text{ m }\times\frac{1.3\times 10^{27}\text{ kg}}{1\text{ m}} \\ &\approx 1.3\times 10^{27}\text{ kg}. \end{align} The first step simply multiplies by $1$, expressed as $c^2/G$. Since we're using units in which $c=G=1$, we could also multiply by $1$ expressed as, say, $c^{42}G^{7/3}$, if we wanted to, because that's also equal to $1$ in these units, and the result would still be legitimate. However, that would give an awkward combination of the SI units "meter" and "kilogram" on the right-hand side. The reason for multiplying by $c^2/G$ is that all of the meters cancel, leaving only kilograms, so we can use this to express a given number of meters as some number of kilograms, or conversely. For example, the mass $M$ of the sun is $2.0\times 10^{30}$ kg, which can be expressed in meters like this: \begin{align} M\approx 2.0\times 10^{30}\text{ kg} &\approx 2.0\times 10^{30}\text{ kg}\times \frac{1\text{ m }}{1.3\times 10^{27}\text{ kg}} \\ &\approx 1.5\times 10^3\text{ m}. \end{align} The Schwarzschild radius of the sun is $R=2GM/c^2$, which can be written simply as $R=2M$ in units where $c=G=1$. Either way, it comes out to be $R\approx 3$ km.
Of course, $1$ Joule and $1$ GeV are two entirely different amounts of energy. They are not equivalent. On the contrary, $1$ GeV is equivalent to $\approx 1.6\times 10^{-10}$ Joules, according to page 126 in reference [2] and also acknowledged explicitly on page 2 in reference [$1$]. Page 2 in reference [$1$] is saying that if we use units in which $c$ and Planck's constant $\hbar$ are both equal to $1$ (that is, $c=\hbar=1$), then we can express kilograms, meters, and seconds all in units of energy. Once we have a quantity in units of energy, we can express it either in GeV or in Joules — with different numeric values, of course, because $1$ GeV $\approx 1.6\times 10^{-10}$ Joules. That's all reference [$1$] means by "where $E$ is an arbitrarily chosen energy unit" below equation (2); nothing novel, just the usual freedom to express a given amount of energy using either GeV, or Joules, or ergs ($1$ erg $=10^{-7}$ J), or kilowatt-hours, or whatever happens to be convenient.
References:
[$1$] Myers, "NATURAL SYSTEM OF UNITS IN GENERAL RELATIVITY," https://www.seas.upenn.edu/~amyers/NaturalUnits.pdf
[2] "The International System of Units (SI), 8th edition," International Bureau of Weights and Measures (BIPM), http://www.bipm.org/utils/common/pdf/si_brochure_8_en.pdf