Symbolic integration and numerical integration yield very different results

This is a topic that has been discussed before, but after reading other posts I could not find a solution that applies for my specific case.

I have the following integral, which I need to evalutate at $n=1,2,3,4...$ and $a=0.1$:

$$ I(n)=\int_{-\pi}^{\pi} {\rm{d}}x \frac{\Gamma(n,0,2ix)}{(1+iax)^n}.$$

The problem is that Integrate and NIntegrate seem to yield different results for sufficiently large $n$. Numerical integration with $a=0.1$ leads to

a = 1/10; Table[NIntegrate[Gamma[n, 0, 2 I x]/(1 + I a x)^n, {x, -\[Pi], \[Pi]}], {n, 1,6}] // Chop {5.82505, 9.23019, 42.4493, 117.67, -303.323, -4233.64} 

However, integrating first and then substituting $a=0.1$ returns the following values:

T = Table[Integrate[Gamma[n, 0, 2 I x]/(1 + I a x])^n, {x, -\[Pi], \[Pi]}, Assumptions->a>0], {n, 1, 6}]; a=1/10; T//N//Chop {5.82505, 9.22998, 41.7814, 0.951343, -26457.8, -3.06652*10^6} 

which only agree with numerical integration for $n=1$ and then they become very different. Interestingly, using Integrate with $a=0.1$ yields another different set of values:

a = 1/10; TT = Table[Integrate[Gamma[n, 0, 2 I x]/(1 + I a x)^n, {x, -\[Pi], \[Pi]}, Assumptions -> a > 0], {n, 1, 6}]; TT // N // Chop {5.82505, 9.23058, 42.4145, 164.712 + 1.19844*10^-9 I, 0, 701973.} 

All these evaluations occur without mathematica returning any errors. What is happening here? Which values are the correct ones?

Things that I know/ tried:

1. The integral is real, as the imaginary part of the integrand is an odd function of x.
2. Assumptions->a>0 is necessary, since for arbitrary $a$ the default solution mathematica returns is valid in the complex plane of $a$ minus the positive real axis.
3. From external arguments, the values that I am most inclined to trust are those from exact integration for arbitrary a>0 and then numerical evaluation at a=0.1 (the second set of numbers).
4. The integrand oscillates more and more the higher $n$ is, and I believe NIntegrate may have problems getting the right answer for highly oscillatory integrands. However, since it does not return any error I do not know if this is the issue.
5. The problem persists if I use GammaRegularized(n,0,2ix) to work with smaller numbers.
6. The problem persists if I increase numerical precision (I am not very familiar with how numerical precision should be managed, maybe I am doing it the wrong way: I increase $MaxExtraPrecision and WorkingPrecision in NIntegrate)

I've read similar posts in which integration methods for NIntegrate such as ExtrapolatingOscillatory and DoubleExponentialOscillatory where suggested, but as far as I know those only work for integrals with infinite range.