Change of variables in a linear differential operator

Question

Suppose I start with a differential operator in two variables x1 and x2 as follows:

dop[x1_, x2_] = a[x1,x2]*D[#, {x1, 2}] + b[x1,x2]*D[#, x1, x2] + c[x1,x2]*D[#, {x2, 2}] +
                d[x1,x2]*D[#, x1] + e[x1,x2]*D[#, x2] + f[x1,x2] &;  

Now, I want to transform this operator expressed as a set of two transformed variables: y1, y2 with the following definitions:

pdx1[y1_, y2_] = 
  1/D[x1[y1, y2], y1]*D[#, y1] + 1/D[x1[y1, y2], y2]*D[#, y2] &;
pdx2[y1_, y2_] = 
  1/D[x2[y1, y2], y1]*D[#, y1] + 1/D[x2[y1, y2], y2]*D[#, y2] &;

Now I just want to substitute these definitions for partial derivatives wrt x1 and x2 in terms of y1 and y2 into my linear differential operator. I try the following:

dop1[y1_, y2_] = 
  Replace[Replace[dop[x1, x2], 
     Derivative[order__][#][x1, x2] & -> 
      Nest[pdx2[y1, y2], 
       Nest[pdx1[y1, y2], f[y1, y2], {order}[[1]]], {order}[[2]]], 
     2] // ExpandAll, 
   Derivative[order__][f][y1, y2] -> Derivative[{order}][#][y1, y2] &,
    2];

The level specification up to 2 is to allow for Plus(level 0), Times(level 1) before hitting the derivative terms(level<=2) where I want to replace by those in terms of new coordinates.

However, this code is not working at all. What am I doing wrong?

Edit:
Following Kuba's suggestion, I am revising my code as follows:

ClearAll["Global`*"];
dop[x1_, x2_] = 
  a[x1[y1, y2], x2[y1, y2]]*D[#, {x1, 2}] + 
    b[x1[y1, y2], x2[y1, y2]]*D[#, x1, x2] + 
    c[x1[y1, y2], x2[y1, y2]]*D[#, {x2, 2}] + 
    d[x1[y1, y2], x2[y1, y2]]*D[#, x1] + e[x1[y1, y2], x2]*D[#, x2] + 
    f[x1[y1, y2], x2[y1, y2]] &;
dop1[y1_, y2_] = 
  Block[{f}, 
    Function[Evaluate[dop[x1, x2]@f[y1[x1, x2], y2[x1, x2]]]] /. 
     Derivative[n_, m_][f][_y1, _y2] :> D[#, {y1, n}, {y2, m}]] /. 
   Derivative[n_, m_][y : y1 | y2][x1, x2] ->  Nest[(1/\!\(
\*SubscriptBox[\(\[PartialD]\), \(y1\)]\(x2[y1, y2]\)\) \!\(
\*SubscriptBox[\(\[PartialD]\), \(y1\)]\((#)\)\) + 1/\!\(
\*SubscriptBox[\(\[PartialD]\), \(y2\)]\(x2[y1, y2]\)\) \!\(
\*SubscriptBox[\(\[PartialD]\), \(y2\)]\((#)\)\)) &, Nest[(1/\!\(
\*SubscriptBox[\(\[PartialD]\), \(y1\)]\(x1[y1, y2]\)\) \!\(
\*SubscriptBox[\(\[PartialD]\), \(y1\)]\((#)\)\) + 1/\!\(
\*SubscriptBox[\(\[PartialD]\), \(y2\)]\(x1[y1, y2]\)\) \!\(
\*SubscriptBox[\(\[PartialD]\), \(y2\)]\((#)\)\)) &, 1/\!\(
\*SubscriptBox[\(\[PartialD]\), \(y\)]\(x1[y1, y2]\)\), n], m];

Why is still the Nest function not working as desired? It is just supposed to determine the $\frac{\partial y}{\partial x}$'s in terms of $\frac{\partial x}{\partial y}$'s by repeated differentiation.