In general, published papers don’t show all the steps of a tedious but straightforward calculation; they just give the final result.
You might as well calculate the more general result $[P_{;abcd}]$, which you can find here if it’s not in DeWitt’s book.
Expand $[(g^{mn}\sigma_{;m} P_{;n})_{;abcd}]=0$. There will be 16 terms.
One vanishes because $[\sigma_{;a}]=0$, six vanish because $[\sigma_{;abc}]=0$, and one vanishes because $[P_{;a}]=0$.
In four of the remaining eight terms, use known results for $[\sigma_{;abcd}]$ (in terms of $R_{abcd}$) and $[P_{;ab}]$ (in terms of $F_{ab}$). (These terms will vanish when you do your two contractions, either because you’ll contract a Riemann tensor on antisymmetric indices or you’ll get $R^{ab}F_{ab}$, which is a symmetric tensor times an antisymmetric one.)
Using the known result for $[\sigma_{;ab}]$ (in terms of $g_{ab}$), the other four reduce to a sum of permutations of $[P_{;abcd}]$. You’ll have to use the commutation rules for covariant derivatives to express $[P_{;bacd}]$, $[P_{;cabd}]$, and $[P_{;dabc}]$ in terms of $[P_{;abcd}]$ plus additional terms. (This is the hardest part of the calculation.) Many of the extra terms will vanish in the coincidence limit and others will vanish when you do your two contractions.
This hand-calculation should take less than six pages of paper, and you shouldn’t need to use a computer algebra system.
If you are mainly interested in the coincidence limits $[a_n]$ of the heat kernel expansion, then there is an alternate approach used by Gilkey. You can construct the most general possible form as a linear combination of a basis of invariants, and determine the coefficients by specializing to a particular geometry. He uses a torus to calculate the 46 terms of $[a_3]$.