Wow Thank you Gianluca so much. And equals zero!

Can I just check the extrema takes the range of integration from opK, has it made any assumption about t to avoid t+s==0? I had to put lots of "Assuming" to simplify explicit cases.

Mathematica does not do this kind of calculations automatically, but we can mimic what we would do by hand.

opK = Function[f, Function[t, Integrate[f[s]/(t + s), {s, 0, 1}]]];
opP = Function[u, 
   Function[t, D[((1 - t^2) t^2 u'[t]), t] - 2 t^2 u[t]]];
commutator = Simplify[opK[opP[f]][t] - opP[opK[f]][t]]

We can introduce some rules to bring the commutator into a desired form:

integralDistribute = 
  Integrate[func_, extrema_] /; Head[Expand[func]] == Plus :> 
   Map[Integrate[#, extrema] &, Expand[func]];
integralByParts[func_] = 
  Integrate[
    smthg_*Derivative[1][func][s], {s, a_, 
     b_}] :> (smthg func[s] /. s -> b) - (smthg func[s] /. s -> a) - 
    Integrate[func[s] Together[D[smthg, s]], {s, a, b}];
bringFactorsInside = 
  HoldPattern[fact_*Integrate[func_, extrema_]] :> 
   Integrate[fact*func, extrema];
integralCollect = 
  Integrate[func1_, extrema_] + Integrate[func2_, extrema_] :> 
   Integrate[func1 + func2, extrema];

Finally:

commutator /. integralDistribute /. integralByParts[f] /. 
      integralByParts[f'] /. integralByParts[f] /. 
    integralDistribute /. bringFactorsInside //. 
  integralCollect // Simplify