Roughly speaking, an AccuracyGoal of $g$ means absolute errors less than $10^{-g}$ are accepted/acceptable. When the magnitude of the solution gets less than this threshold, then errors greater in magnitude may cause the solution to change sign. I think of AccuracyGoal as telling NDSolve the threshold below which a solution may safely be treated as zero plus or minus an error. In this case, you don't want that. Theoretically raising AccuracyGoal so that $10^{-g}$ is much less than the magnitude of the solution is a fix. However, the corrector step in LSODA sometimes/often will exhibit a convergence problem, which it did for me. I could get that to go away by keeping PrecisionGoal modestly low.

So here's the first solution:

sol12a = NDSolve[{y1'[x] + 4 y1[x]/x == 0, 
     y2'[x] + (3/x) (((\[Sqrt](y1[x] + y2[x]))/((\[Sqrt](y1[x] + y2[x]) - 
               P1 ((y2[x] - y3)/y1[x])) - 
             P2 ((y2[x]^2 - y3^2)/y1[x])))) y2[x] == (-4 P2 (y2[x]^2 - y3^2) - 
         4 P1 (y2[x] - y3))/((\[Sqrt](y1[x] + y2[x]) - 
           P1 ((y2[x] - y3)/y1[x])) - P2 ((y2[x]^2 - y3^2)/y1[x]) x), 
     y1[10] == 10^-7, y2[10] == 0}, {y1, y2}, {x, 10, 10^12}, 
    Method -> "Automatic", PrecisionGoal -> 6, 
    WorkingPrecision -> 100, MaxSteps -> \[Infinity]];

Now, the first equation is self-contained subsystem with an exact solution. So here's a second method that avoid the numerical issue entirely:

y1sol = DSolve[{y1'[x] + 4 y1[x]/x == 0, y1[10] == 10^-7}, y1, x];
y2sol = 
  NDSolve[{y2'[x] + (3/x) (((\[Sqrt](y1[x] + 
               y2[x]))/((\[Sqrt](y1[x] + y2[x]) - 
               P1 ((y2[x] - y3)/y1[x])) - 
             P2 ((y2[x]^2 - y3^2)/y1[x])))) y2[x] == (-4 P2 (y2[x]^2 - y3^2) - 
         4 P1 (y2[x] - y3))/((\[Sqrt](y1[x] + y2[x]) - 
           P1 ((y2[x] - y3)/y1[x])) - P2 ((y2[x]^2 - y3^2)/y1[x]) x), 
     y2[10] == 0} /. First@y1sol,
   {y2}, {x, 10, 10^12}, Method -> "Automatic", 
   WorkingPrecision -> 30, MaxSteps -> \[Infinity]];
 sol12b = Join[y1sol, y2sol, 2];

I usually do something like Daniel's approach, but in this case I'd scale the variable y1, which is the misbehaving one. Transforming $x=e^t$ might also be convenient — or not. It often is when a variable ranges from $10$ to $10^{12}$. Here's a third solution:

eqns0 = {y1'[x] + 4 y1[x]/x == 0, 
   y2'[x] + (3/
        x) (((\[Sqrt](y1[x] + y2[x]))/((\[Sqrt](y1[x] + y2[x]) - 
             P1 ((y2[x] - y3)/y1[x])) - 
           P2 ((y2[x]^2 - y3^2)/y1[x])))) y2[
       x] == (-4 P2 (y2[x]^2 - y3^2) - 
       4 P1 (y2[x] - y3))/((\[Sqrt](y1[x] + y2[x]) - 
         P1 ((y2[x] - y3)/y1[x])) - P2 ((y2[x]^2 - y3^2)/y1[x]) x)};
eqns1 = eqns0 /. y1 -> Exp@*u1 // Simplify[#, u1[x] > 0] &;

sol12c = NDSolve[Join[eqns1, {u1[10] == Log[10^-7], y2[10] == 0}],  
    {u1, y2}, {x, 10, 1*^12}] /. (u1 -> usol_) :>  y1 -> Exp@*usol;

Unfortunately, at the second step you have replaced x by t instead of Exp[t] and y'[x] -> y2'[t]*Exp[-t], after correcting these minor typos the solution remains the same as previous, i.e. negative valued at some places and independent of P1 and P2

Theoretically the coupled equations are devised such that the y1, y2 should always remain positive and, with increasing x at first y2 should increase and dominate over y1 but at higher x, y1 should become higher than y2. I also have tried to keep the P1, P2 and initial condition practical.

   P1 = 10^7;
   P2 = 10^12;
   y3 = 10^-6*BesselK[2, x]*x^(-3/2);

   sol12 = N[DSolve[Rationalize[SetPrecision[{y1'[x] + 4 y1[x]/x == 0,
         y2'[
            x] + (3/x) (((\[Sqrt](y1[x] + 
                   y2[x]))/((\[Sqrt](y1[x] + y2[x]) - 
                   P1 ((y2[x] - y3)/y1[x])) - 
                 P2 ((y2[x]^2 - y3^2)/y1[x])))) y2[
             x] == (-4 P2 (y2[x]^2 - y3^2) - 
             4 P1 (y2[x] - y3))/((\[Sqrt](y1[x] + y2[x]) - 
               P1 ((y2[x] - y3)/y1[x])) - P2 ((y2[x]^2 - y3^2)/y1[x]) x),
         y1[0.1] == 10^7, y2[0.1] == 0.10}, 1000]], {y1, y2}, {x, 0.1,
        10^10}, 
      Method -> {"StiffnessSwitching", "ExtraPrecision" -> 1000, 
        InterpolationOrder -> \[Infinity]}], {25, 
      500000}](*,WorkingPrecision\[Rule]50,MaxSteps\[Rule]\[Infinity],\
   InterpolationOrder\[Rule]All]*)

by evaluating this I am getting at higher x the issue

Can you explain why you find the numerical solution "wrong"?

Your numerical solution y1 has small fluctuations around 0. In the intervals where y1 is negative, the LogLogPlot will leave a gap. P1 and P2 do not affect y1. Luckily, the differential equation for y1 can be solved symbolically:

DSolve[{y1'[x] + 4 y1[x]/x == 0, y1[10] == 10^-7}, y1, x]