A computer-implemented method for identifying a ligand binding site of a protein, comprising: a. obtaining sequence information for a plurality of sequence variants of the protein; b. obtaining a measure of the biological activity of each of the plurality of the sequence variants; c. using a three-dimensional model of the protein's structure to define patches on the protein structure containing a plurality of residues within close proximity of each other, wherein the residue at the center of the patch is determined based on solvent accessibility, and wherein the residues are within 20 angstroms of the residue at the center; d. correlating the sequence variability within each of the patches with the biological activity data for each sequence variant, wherein the information on the sequence variability within each patch and the biological activity data for each sequence variant is represented in a non-binary or binary form, compiled as a vector file, and processed by two or more computational methods testing each patch for an ability of the sequence variability within the patch to predict changes in the biological activity; wherein the two or more computational methods comprise i) at least one method involving multivariate analysis, at least one method involving univariate analysis, and at least one machine learning algorithm, or ii) two or more of multiple regression analysis, logistic regression analysis, support vector machine, Fischer's Exact test, Hidden Markov Models, Neural Networks, Random Forest, Decision Trees, or Bayesian Networks; and wherein a residue is identified as a potential binding site residue if variation at the residue exhibits statistically significant correlation with the biological activity; and e. identifying the residues forming the binding site by comparing the potential binding site residues identified with the two or more computational methods and identifying residues based on the agreement between the two or more computational methods. a. obtaining sequence information for a plurality of sequence variants of the protein; b. obtaining a measure of the biological activity of each of the plurality of the sequence variants; c. using a three-dimensional model of the protein's structure to define patches on the protein structure containing a plurality of residues within close proximity of each other, wherein the residue at the center of the patch is determined based on solvent accessibility, and wherein the residues are within 20 angstroms of the residue at the center; d. correlating the sequence variability within each of the patches with the biological activity data for each sequence variant, wherein the information on the sequence variability within each patch and the biological activity data for each sequence variant is represented in a non-binary or binary form, compiled as a vector file, and processed by two or more computational methods testing each patch for an ability of the sequence variability within the patch to predict changes in the biological activity; wherein the two or more computational methods comprise i) at least one method involving multivariate analysis, at least one method involving univariate analysis, and at least one machine learning algorithm, or ii) two or more of multiple regression analysis, logistic regression analysis, support vector machine, Fischer's Exact test, Hidden Markov Models, Neural Networks, Random Forest, Decision Trees, or Bayesian Networks; and wherein a residue is identified as a potential binding site residue if variation at the residue exhibits statistically significant correlation with the biological activity; and e. identifying the residues forming the binding site by comparing the potential binding site residues identified with the two or more computational methods and identifying residues based on the agreement between the two or more computational methods.