In eukaryotes, hundreds of protein kinases (PKs) specifically and precisely modify thousands of substrates at specific amino acid residues to faithfully orchestrate numerous biological processes, and reversibly determine the cellular dynamics and plasticity. Although over 100,000 phosphorylation sites (p-sites) have been experimentally identified from phosphoproteomic studies, the regulatory PKs for most of these sites still remain to be characterized. Here, we present a novel software package of iGPS for the prediction of in vivo site-specific kinase-substrate relations mainly from the phosphoproteomic data. By critical evaluations and comparisons, the performance of iGPS is satisfying and better than other existed tools. Based on the prediction results, we modeled protein phosphorylation networks and observed that the eukaryotic phospho-regulation is poorly conserved at the site and substrate levels. With an integrative procedure, we conducted a large-scale phosphorylation analysis of human liver and experimentally identified 9719 p-sites in 2998 proteins. Using iGPS, we predicted a human liver protein phosphorylation networks containing 12,819 potential site-specific kinase-substrate relations among 350 PKs and 962 substrates for 2633 p-sites. Further statistical analysis and comparison revealed that 127 PKs significantly modify more or fewer p-sites in the liver protein phosphorylation networks against the whole human protein phosphorylation network. The largest data set of the human liver phosphoproteome together with computational analyses can be useful for further experimental consideration. This work contributes to the understanding of phosphorylation mechanisms at the systemic level, and provides a powerful methodology for the general analysis of in vivo post-translational modifications regulating sub-proteomes.