Inflated p-values in quantitative trait analysis I am performing a quantitative trait association between the expression of one gene and ~400,000 methylation values. First, both variables are rank inverse normal transformed, adjusted for confounders (technical and for family structure) in a mixed-effects model and we use the residuals for the association. Then we do a linear regression between the residuals of expression and each methylation site, adding other four confounders in the model (age, BMI, etc.). The sample size is of ~500. For the association we are using MatrixEQTL.

When we check the p-values, they are clearly inflated. Is it something that I have to be concerned about? Did I probably miss something in the analysis? Thanks in advance.
Permutation as suggested by @StupidWolf's comment is essential to understand what's going on. If permutation makes this pattern go away, then you have a problem with your model specification, there's something uncorrected.

If your data are weird, well, that's just how they are. But this argues to me that something else is going on confounding your associations. Cryptic population structure in the data would lead to this kind of pattern, e.g. here. That is the first thing that most reviewers would say, I'd guess.

Some random thoughts:


What does the distribution of your methylation data look like after normalization? 
Your observed distribution is basically never on the expected line, which argues strongly for confounding. 
How are you adjusting for non-independence of methylation types? Family-wise correction is likely not specific enough. You may have to estimate an identity-by-state SNP similarity matrix for your population and fit that as part of your model, or at least fit 10-20 principal components to account for population structure.
How are you representing your methylation status? If it's binary you might look at this paper.