Maybe replace the algebraic equation by a function definition. Reduce the AccuracyGoal. I see no need for such a high setting, and it reduces the step size to zero. You could raise the PrecisionGoal if a higher quality solution is needed.

Following integrates over the whole interval:

Block[{P1},
 q[x_?NumericQ] := (830397. - 
     0.375 x^(87/296))/(4.51833*10^8 x^(34/37) - 1.02416*10^7 x - 
     1. x^(383/296));
 P1[x_] := 2.372463129699238`*^-19 Sqrt[y2[x] + y4[x] + y1[x] + y3[x]];
 sol1 = NDSolve[
   SetPrecision[{y1'[x] == y1[x] (-3 P1[x] - q[x] + G'[x]/G[x]), 
     3 P1[x] y2[x] + y2'[x] == 0.95 y1[x] q[x], 
     4 P1[x] y3[x] + (y1[x] G'[x])/G[x] + y3'[x] == 0, 
     4 P1[x] y4[x] + y4'[x] == 0.05 y1[x] q[x],
     (*P1[x]==2.372463129699238`*^-19 Sqrt[y2[x]+y4[x]+y1[x]+y3[x]],*)
     G'[x] == (-8.722718060807123`*^73 + 
         1.80559635200202`*^-76 G[x]^4 y3[x])/G[x]^2, 
     3 y1[10^-30] == 10^51 \[Pi]^2, y2[10^-30] == 0, 
     y3[10^-30] == 3.3*10^56 \[Pi]^2, y4[10^-30] == 0, 
     G[10^-30] == 5.468912167778173`*^27}, 100], {y1, y2, y3, y4, 
    G}, {x, 10^-30, 10^5}, Method -> "Automatic", AccuracyGoal -> 20, 
   PrecisionGoal -> 5, WorkingPrecision -> 80]
 ]

If you didn't suppress the error message with Quiet, you would see what it thinks the problem is: 6 equations (DEs), 7 unknowns (dependent variables).

Just as a simple quick experiment, can you

scale all your 10^-30 and 10^-76 and 10^-19 up nearer to 10^-2

scale all your 10^51 and 10^73 and 10^56 and 10^27 and 10^5 down nearer to 10^2

include your definition of q[x]

possibly eliminate your use of P[x] because that doesn't seem to be the usual style

eliminate all your /G[x] denominators by multiplying

remove the Quiet and DeleteCases

and then see whether it is able to produce plausible interpolating functions?

If so then gradually and gently begin slowly backing out each of those changes a little bit at a time and confirming that it still solves?