P-value mining on large number of combinations of variables
You need to investigate multiple hypothesis correction methods, like Bonferroni correction or Benjamini-Hochberg false discovery rate. The problem with this sort of analysis is that your associations are unexpected, so you don't have any a priori hypothesis. All you can do is test every combination, and then see what is statistically significant after accounting for all the tests.

As you do more tests, it becomes increasingly likely that you'll find an association just by chance alone. To combat this, multiple hypothesis correction methods make the p-value threshold for significance more conservative. If you do only 1 test, a p-value of 0.05 may indicate a significant result. If you do a million tests, however, many of them will have a p-value of 0.05 by chance, so you need to be more conservative. This practice is not generally frowned upon, since it appropriately takes into account the analysis methods. You can easily be dishonest with this approach, however, by doing many tests and only reporting the significant ones without correction.

Obligatory xkcd: https://xkcd.com/882/