How do I find co-occurring labels?

I have a matrix that is of the shape n_samples by n_features, where n_samples ~ 7000 and n_features ~ 800. To clarify, n_features are categorical labels, each representing an amino acid position.

I'd like to know which columns have amino acid positions which are highly correlated with each other. For example, one might have a protein where amino acid position 255 = H when 288 = Y ("HY" combination), and 255 = N when 288 = P ("NP" combination). This is a completely hypothetical scenario, btw. In the data set, I might find, for example, 3500 instances of HY, 100 instances of HP, 3400 instances of NP. In this case, we would say that position 255 and 288 are correlated with each other.

The problem is, I'm not sure how to compute this. Using scikit-learn, I have transformed the data into numerical categorical labels, where an amino acid position is now represented by numbers (i.e. from 0 to 20). So now, it's rows=each sequence, and columns=each position, and the numbers in each cell=0 to 20, each letter represented by a number.

What would your advice be for going from this data matrix to finding columns that "covary" together?

I'm familiar with one branch of the literature. The early literature on coevolution formed positional covariance using $$\chi^2$$ measures. Someone got the idea to compare positional covariance to structural contact maps derived from pdbs. When compared in this way normalized mutual information handedly beats $$\chi^2$$ metrics. Several other metrics, such as SCA(which means a number of things, but probably means inner product), sum-of-squares-methods et al. have fell by the wayside. Normalized mutual information was improved upon by z-scoring away phylogenetic background for positions. I'm not sure where the literature is currently.