Reasons for different parameters via MoM and MLE I got a study of 210 samples and I tried fitting gamma distribution to them. I used method of moments and maximum likelihood estimation to calculate the parameters, but parameters came out quite different (that much that the shape of curve is different).

This is the histogram with green graph of gamma PDF, calculated from MLE and red graph calculated from MoM. The data is this.
The derivation of MoM and MLE are certainly correct (I tried them both on samples of 1000 generated by gamma distribution and they are very precise), but if you want to check it out:

Method of momements
For gamma distribution we have $\text{E}(\bar{X}) = \frac{a}{\lambda}$ and $\text{E}(\bar{X}^2) = \frac{a(a+1)}{\lambda^2}$. On the other hand we have
formulas $\text{E}(\bar{X}) = \frac{1}{n} \sum_{i=1}^n x_i = \bar{x}$ and $\text{E}(\bar{X}^2) = \frac{1}{n} \sum_{i=1}^n x_i^2$. Combining these we get
\begin{equation}
\hat{a} = \hat{\lambda} \bar{x}~~~~\text{and}~~~~\hat{\lambda} = \frac{\bar{x}}{\frac{1}{n}\sum_{i=1}^n x_i^2 - \bar{x}^2}.
\end{equation}
From here I got
$$a \approx 0.799174075139 ~~~~\text{and}~~~~ \lambda \approx 1.31877117345. $$
Maximum likelihood method
For gamma distribution we have
$$f_{X_i}(x_i) = \frac{\lambda^a}{\Gamma(a)} x_i^{a-1} e^{-\lambda x_i}.$$
Taking the product from $1$ to $n$ we get the likelihood function
$$L = \prod_{i=1}^n \frac{\lambda^a}{\Gamma(a)} x_i^{a-1} e^{-\lambda x_i}$$
and log-likelihood function
$$\ell = \log L = \sum_{i=1}^n \left(a \log\lambda - \log\Gamma(a) + (a-1)\log x_i -\lambda x_i\right) =\\
            =na \log\lambda - n \log\Gamma(a) + (a-1)\sum_{i=1}^n \log(x_i) - \lambda \sum_{i=1}^n x_i.$$
Calculating partial derivatives with respect to $a$ and $\lambda$ and equating with $0$ we obtain
\begin{align}
\frac{\partial \ell}{\partial a} &= n \log \lambda - n \psi(a) + \sum_{i=1}^n \log(x_i) = 0,\\
\frac{\partial \ell}{\partial \lambda} &= \frac{na}{\lambda} - n \bar{x} = 0,
\end{align}
where $\psi(a)$ is a digamma function ($\psi(a) = \frac{\Gamma'(a)}{\Gamma(a)}$). Now I solved this system by plugging the second equation into first, threw everything on one side and using scipy.optimize.newton I found the root $a$. Plugging $a$ into second equation I got
$$a \approx 1.59540858049 ~~~~\text{in}~~~~ \lambda \approx 2.63269156405.$$

So what I would like to know is what are the reasons for this (e.g. small sample size, bad sample, etc.).

 A: One reason is, that the data does not seem to follow a gamma distribution. Comparing the ECDF with the one predicted by the gamma distribution using either parameter set shows this quite clearly (blue = MoM, red = MLE).

Therefore the different parameter estimation techniques (focusing on different aspects of the distribution) lead to quite different results.
A: The two methods generally look at different aspects of the data; if those aren't consistent with each other, they can result in very different estimates. 
In the case of the gamma the log-likelihood is a function of the mean and the geometric mean, while the MOM estimates are a function of the mean and the variance. 
Indeed the ML estimate of the shape parameter is a function of the ratio of the geometric mean to the arithmetic mean, while the MoM estimate is a function of the ratio of the arithmetic mean to the standard deviation, while for both cases the corresponding estimate of the rate parameter is simply the shape parameter estimate divided by the sample mean.
Clearly the variance (or the standard deviation) and the geometric mean look at different aspects of the sample. However, even when a moment statistic appears in the likelihood it may impact the two estimators in different ways because it enters the formula differently; for example, while the sample mean is present in both shape estimates it comes into the calculation differently.
There are multiple reasons the pieces of information you get may be inconsistent with each other. This includes that the gamma model may be a poor description (it looks like it may be the case here, and then you would definitely expect different fits) but it can also happen when samples are small or when the MoM estimates are not very efficient (or when they are badly biased - this may persist into fairly large samples).
You may not always be able to untangle contributions of the possible causes (at least not without something like a simulation study perhaps).
