Get["RISC`HolonomicFunctions`"];
(
vars = Array[x, 6];
expr = Total[Times @@@ Subsets[vars, {3}]]^2 / (Times @@ vars);
(* We want the constant coefficient of expr^n, which is obtained by integrating
   expr^n / (x1 .. x6) over some contour. *)
anns = {Annihilator[expr^n / (Times @@ vars), Append[Der /@ vars, S[n]]]};
Do[
  res = Timing[FindCreativeTelescoping[anns[[i]], Der[x[i]]][[1]]];
  AppendTo[anns, res[[2]]];
  Print[]; Print[StringForm["Integration `1` took `2`s.", i, res[[1]]]];
, {i, 6}];
InputForm[ApplyOreOperator[Factor[anns[[-1,1]]], a[n]]]
)