Inter-chromosomal variation in the pattern of human population genetic structure

Data mining, processing and description

We downloaded the HapMap SNP data (http://www.hapmap.org, release # 24, on NCBI B36 assembly, dbSNP b126). The HapMap project contains genotypes from 60 unrelated Caucasians from the USA with northern and western European ancestry (CEU), 60 unrelated Yoruba individuals from Ibadan, Nigeria (YRI), 45 Japanese individuals from Tokyo (JPT) and 45 Han Chinese individuals from Beijing (CHB) (http://www.hapmap.org). Two criteria were used to filter the SNPs included in the analysis: (i) locus call rate ≥95 per cent (ie we excluded all SNPs with more than 5 per cent missing data) and (ii) the SNP should be shared among populations so that the same sets of SNPs were used throughout in the population comparisons. A computer program using Python (http://www.python.org) was written to export and pre-process SNP genotype information from the databases. Genotypes were summarised for each population. For each dataset, the number of alleles per locus (SNP) was coded to a string of numbers to obtain a full design matrix of alleles (the cells give the number of copies of each major allele for each individual: zero, one or two). Figure 1 depicts our approach to SNP mining, multivariate chromosomal and population diversity and network analysis strategies. Of the total 3.7 million SNPs in the HapMap data release,[23, 25] 809,000 SNPs fulfilled the criteria and were used in this analysis.

Statistical analysis

Multivariate statistical techniques (namely, principal component [PC], cluster, discriminant, network analyses and F_ST statistics) were used to examine chromosomal structure within and between populations and associated functional networks by estimating chromosomal overall differentiation values. The analysis was carried out either using all SNPs together or separately for each chromosome. Because PC analysis (PCA) does not take into account group differences in reducing the dataset to a few representative variables, and it can be difficult to make appropriate inference about population relationships from the PC scatter plot, we further analysed the data using cluster analysis (CA) to classify individuals into mutually exclusive groups with high homogeneity within clusters and with low homogeneity between clusters. In other words, CA provides a visual assessment and identifies individuals who are similar (or dissimilar) to one another. To further confirm the grouping obtained in CA, discriminate analysis (DA) was performed. DA consists of the separation of a priori given classes for each individual. The variance-covariance between classes is maximised and the variance-covariance within classes is minimised under simultaneous consideration of all analysed data.

PCA was done using the EIGENSOFT software package (http://genepath.med.harvard.edu/~reich/Software.htm) either on all SNPs simultaneously (all loci together) or separately per each chromosome. The analysis follows singular value decomposition, a procedure that produces eigenvectors, corresponding eigenvalues and proportions of eigenvalues, as well as the scores of the PCs [26]. Using PCA, we estimated axes of variation corresponding to ancestry. The first eigenvector separates the samples in a way that explains the largest amount of variability, while the second and subsequent ones explain lesser amounts of variability. The spatial relationships of populations in each chromosome and all chromosomes were presented by plotting the scores of the first and second PCs. The numbers of significant PCs (at the level of p < 0.05) were tested using Tracy-Widom statistics. Pairwise population genetic diversity was determined by calculation of Wright's F_ST using EIGENSOFT. F_ST values indicate how much of the genetic variability between individuals from different populations is due to population affiliation.

Hierarchical clustering of molecular variance was followed using the similarity for qualitative data (SIMQUAL) module with the first 10 PCs that account for most of the variation. Average taxonomic distance matrices (DIST) were computed as a measure of genetic distance. This matrix was subjected to unweighted pair-group method analysis (UPGMA) to generate a dendrogram using the Sequential, Agglomerative, Hierarchical and Nested (SAHN) module. Both numerical taxonomic analyses were performed using the Numerical Taxonomy and Multivariate Analysis System Program, version 2.11f (NTSYS-pc) [27]. The cophenetic correlation coefficient was calculated, and Mantel's test [28] was performed, to check the goodness of fit of a CA.

In DA, a linear combination of features that best separates two or more groups of objects is sought. The discriminant functions are determined based on the maximisation of the ratio of the external (between populations) to the internal (between individuals within the same population) variability [29]. The values of Wilks' lambda (λ) and their X² statistics are used to evaluate the number of significant discriminant functions. In turn, to determine the most important features of the objects, partial Wilks' λ and its Fisher statistics were utilised [29]. Discriminant function analysis [30] was done following the SAS system [31] DISCRIM, CANDISC and STEPDISC procedures, and significance was tested using Wilks' λ [32]. In order to avoid the limitation of a large number of alleles compared with the number of observations and the correlation occur in allele frequencies, we ran discriminant analysis using the uncorrelated SNPs in the top significant PCs. This ensures that variables submitted to DA are perfectly uncorrelated and that their number is lower than that of analysed individuals. Linear discriminant analysis is similar to logistic regression and is useful for building a predictive model of group membership based on observed characteristics. The procedure yields a set of discriminant functions based on the linear combinations of variables that provide the best discrimination between groups.

In the final set of analyses, a dataset containing a total of 126 genomic regions linked to SNPs that differed between populations (F_ST ≥ 0.5) was uploaded into the Ingenuity Pathways Analysis (IPA) 8.7 network analysis (Ingenuity Systems, Redwood City CA, USA). The network generated from the 126 input genes (called focus genes) uses both direct and indirect relationships/connectivity. These networks were ranked by scores that measured the probability that the genes were included in the network by chance alone. Networks with scores of three or more were classified as not being generated by random chance [33]. The significance threshold for Fisher's exact test to determine the probability that each biological function and/or disease assigned to that network is due to chance alone was 0.05 or less. Canonical pathways associated with input genes were elucidated with a statistical significance value. The gene ontology (GO) analysis was used to identify functional commonalities between the genes based on the number of shared ancestors in gene products (http://gostat.wehi.edu.au).

Estimates of F_STdiffer between chromosomes and populations

The empirical genome-wide distribution of F_ST showed heterogeneity in chromosomal ancestry across the genome (Figure 2). The average F_ST values for autosomes and sex chromosomes were significantly different (0.120 and 0.210, respectively; t-test, t = 16.1, p < 10^-15). The higher average F_ST for chromosome X compared with autosomes might indicate differences in inheritance mechanisms that potentially affect sex chromosomes and autosomes differently [34–36]. Similarly, statistically significant differences in F_ST estimates (F_ST analysis per chromosome- F_ST_SNP_CHROM) among autosomal chromosomes were detected. The variability that drives F_ST distribution among autosomes could be due to variations in natural selection and/or recombination rates during meiosis [37]. Based on Wright's qualitative guidelines, F_ST statistics range from 0 (no differentiation) to 1 (fixed difference between populations for different alleles). Values of F_ST less than 0.05 represent low or little population genetics differentiation, values between 0.05 and 0.15 represent moderate population divergence, values between 0.15 and 0.25 indicate large population differentiation and F_ST values greater than 0.25 represent very large population divergence [38]. Usually, an F_ST > 0.5 is considered sufficient for ancestry differentiation.

The global pairwise F_ST value estimated for the 210 worldwide samples using all loci together (F_ST_SNP_ALL) was 0.130 (p < 10^-6). As F_ST increases, populations become more distant and/or unrelated to each other [39]. We observed an average genetic differentiation between CEU and YRI (F_ST = 0.153), CEU and CHB (F_ST = 0.110), CEU and JPT (F_ST = 0.111), YRI and CHB (F_ST = 0.190) and YRI and JPT (F_ST = 0.192) (Figure 2). It is evident that more divergence has occurred between YRI and each of the three other populations than between the other pairs of populations. Genetic distances between CHB and JPT populations were low (mean F_ST = 0.007), indicating that substantial gene flow compensates for the effects of genetic drift. This low F_ST value as a result of high similarities in allele frequencies between CHB and JPT samples motivates researchers to analyse CHB and JPT populations jointly, as a single panel [25].

Regardless of the populations compared, most of the variation was observed within populations (average, 87 per cent versus 13 per cent variation observed between populations). Within-population diversity reflects the number of different types in the population, taking into account their frequencies. By contrast, between-population differentiation measures variation based on the relative frequencies of types within these subpopulations and, ideally, measures the average distance of subpopulations from their respective lumped remainders. The fact that only 13 per cent of the total genetic variation results from differences between populations indicates that alleles present in one population are also present in other populations [40–43]. The remaining 87 per cent represents the average difference between members of the same population. One way to interpret this number is to say that the expected genetic difference between unrelated individuals from distant continents exceeds by 13 per cent the expected difference between members of the same community [44]. An interesting common feature in population genetics studies of humans, animals, plants and other types of species is that within-population diversity is greater than between-population diversity [45–48]. This estimate is highly consistent for protein polymorphisms, blood groups, microsatellites, SNPs and morphological/phenotypic markers [49, 50]. Therefore, it is necessary to quantify and control population structure, not only for major population differences, but also for subtle variation/structure arising within populations.

Significant numbers of PCs vary by chromosome

As shown in Figure 3, the first component -- which accounted for 50.2 per cent of variation -- separated YRI populations from CEU and CHB/JPT, while the second (accounting for 24.2 per cent of the total variance) could be associated with an Africa/Europe gradient. The overlap between the CHB and JPT populations suggests a low level of genetic differentiation, as shown by the pairwise F_ST divergence (Figure 2). Overall, these results demonstrate a clear partition of the West African populations considered. On a genome-wide average scale, about 74 per cent of the diversity in human samples was explained by the first two PCs. The eigenvalues for PC3-PC10 showed a plateau, suggesting that the first two PCs account for most of the populations' average substructure in this analysis. Such a genome-wide level of population structure may lead to an erroneous conclusion that the samples are genetically homogeneous. Thus, correction for population structure is only as good as the level of structure (at a finer or coarser level) that one wishes to correct.

In order to obtain fine-scale resolution of population membership, PCA was performed on each chromosome. Our analysis showed that the contribution of the first two PCs in classifying geographical regions varied among chromosomes, ranging from 65 per cent (Chr X) to 76 per cent (Chr 15). The contribution of PC1 ranged from 47 per cent (Chr X) to 51 per cent (Chr 3, Chr 8). The contribution of PC2 to the total variation ranged from 18 per cent (Chr X) to 27 per cent for Chr 15 (Supplementary Figure S1 (Figure 6)). As shown in Figure 4, on a finer scale, the number of significant PCs accounting for population differentiation varies from 2 (Chr 2) to 31 (Chr X) among chromosomes. The higher number of significant PCs on Chr X explains why it has the lowest chromosome-wise contribution to the first two PCs (Supplementary Figure S1 (Figure 6)).

We next characterised the genetic relationships existing among the four different populations. The diagrammatic output of CA (constructed from PCs) for the mean of 210 individuals indicated that these individuals could be clustered into groups that basically coincided with their geographical distribution (data not shown). The analysis confirmed the distinctiveness of the CEU and YRI populations and the close average genetic distance between the CHB and JPT populations (Supplementary Figure S2 (Figure 7)). The results of the chromosomal-based CA (data not shown) were comparable to those of the PCA, and both methods classified racial populations into separate groups.

DA predicts population membership for 70 per cent of individuals

The eigenvalue is the ratio of the between-groups sum of squares to the within-groups sum of squares [52]. This value measures the spread of the group centroids in the dimension of multivariate space (eigenvalue 10587.53; p < 0.0001). The canonical correlation measures the association between discriminant scores and groups. This association appeared to be statistically significant (stepwise elimination, R² = 0.99, 0.96; p < 0.0001), and the data were subjected to the stepwise procedure. The first canonical discriminant function had a high eigenvalue, accounting for more than 81 per cent of the total variance. The first and second functions together accounted for 100 per cent of the variance.

Functional networks and pathways in highly stratified genomic regions

To characterise the main functional networks/pathways underlying genes with substantial population differentiation, we carried out network analysis (see Materials and methods) for between-population comparisons (between CEU-YRI, CEU-CHB/JPT and YRI-CHB/JPT). A summary of networks, molecular and cellular functions, diseases and disorders and canonical pathways associated with the genomic regions are presented in Table 3. A total of 126 genes were significantly differentiated among populations and eligible for network analysis, which led to the identification of five significant networks (Figure 5). Network 1 was centred on the nuclear factor (NF)-κB complex and had 19 focus genes; network 2 was centred on tumour necrosis factor (TNF) and had 14 focus genes; network 3 was centred on v-myc myelocytomatosis viral oncogene homologue (MYC) and had 13 focus genes; network 4 was centred on interferon-γ (IFNG) and had 12 focus genes; and network 5 was centred on hepatocyte nuclear factor-4α and had 12 focus genes. Interestingly, although no genes were shared among all these five different networks, two networks (networks 1 and 2) contained chloride intracellular channel 4 (CLIC4), which plays a role in apoptosis, differentiation and diabetes. The overlap of this gene suggests that similar biological pathways were targeted by selection in these populations. In addition, the gene involved in skin, hair and eye pigmentation -- including oculocutaneous albinism II (OCA2), hect domain and RLD 2 (HERC2), ectodysplasin A receptor (EDAR) and solute carrier family 45, member 2 (SLC45A2) were over-represented in our GO analysis (http://gostat.wehi.edu.au) (Table 4). Enriched GO biological function terms include cytoskeletal protein binding (p = 1.1 × 10^-6), actin binding (p = 1.0 × 10^-6) and fibroblast growth factor receptor antagonist activity (p = 6.2 × 10^-5).

	IPA categories	Statistical measures	Associated gene(s)
No	Top networks	Network score*	No. of candidate genes
1	Cancer, cell death, dermatological diseases and conditions	33	19
2	Carbohydrate metabolism, dermatological diseases and conditions, lipid metabolism	24	14
3	Post-translational modification, embryonic development, tissue development	24	13
4	Inflammatory response, immunological disease, carbohydrate metabolism	20	12
5	Cell cycle, hair and skin development and function, nervous system development and function	19	12
	Molecular and cellular functions	p value +	No. of candidate genes
1	Cell death	2.55E-04 - 3.60E-02	19
2	Cell-to-cell signalling and interaction	6.33E-04 - 3.25E-02	6
3	Cellular assembly and organisation	6.33E-04 - 3.36E-02	18
4	Cellular compromise	6.33E-04 - 3.25E-02	8
5	Gene expression	1.50E-03 - 3.25E-02	7
	Associated disease and disorders	p value	No. of candidate genes
1	Inflammatory disease	1.43E-07 - 3.25E-02	47
2	Gastrointestinal disease	1.22E-06 - 3.39E-02	29
3	Genetic disorder	1.22E-06 - 3.25E-02	69
4	Endocrine system disorders	1.72E-05 - 3.25E-02	35
5	Metabolic disease	1.72E-05 - 1.31E-02	39
	Top canonical pathways	p value	No. of candidate genes
1	Androgen and oestrogen metabolism	2.15E-02	3
2	Neuroprotective role of THOP1 in Alzheimer's disease	2.85E-02	2
3	Alanine and aspartate metabolism	3.39E-02	2
4	Retinol metabolism	3.39E-02	2
5	Pentose and glucuronate interconversions	4.76E-02	2
	Physiological system development and function	p value	No. of candidate genes
1	Hair and skin development and function	6.26E-05 - 3.25E-02	9
2	Nervous system development and function	7.22E-04 - 3.36E-02	15
3	Connective tissue development and function	1.17E-03 - 3.25E-02	6
4	Skeletal and muscular system development and function	1.17E-03 - 3.25E-02	11
5	Tissue development	1.17E-03 - 3.25E-02	14

*Networks with scores ≥3 have a 99.9 per cent confidence of not being generated randomly.

⁺The IPA computes p values of statistically significant findings by comparing the number of molecules of interest relative to the total number of occurrences of these molecules in all functional/pathway annotations stored in the Ingenuity Pathways Knowledge Base (Fisher's exact test with p value adjusted using the Benjamini-Hochberg multiple testing correction).

In conclusion, our approach offers a complementary statistical strategy for summarising overall variability and global versus chromosomal structure, assessing population structure and identifying genomic regions driving genetic divergence among populations. We first used PCA (to reduce data dimensionality); however, because PCA does not take into account group differences in reducing the dataset to a few representative variables, we further analysed the data using CA to classify individuals into mutually exclusive groups with high homogeneity within clusters and low homogeneity between clusters. To further confirm and predict group membership, DA was performed using the top significant PCs. PCs were used to ensure that variables submitted to DA were perfectly uncorrelated, and that their number was lower than that of analysed individuals. Finally, using F_ST (to study population differentiation) analysis, we described the importance of chromosome-based population genetic structure to identify differing genomic regions driven by natural selection. We followed the target genomic regions using network/pathway analysis to elucidate their roles and functional implications in human genetic variations and diseases.

Principal component analysis (PCA)

where T is the score matrix summarising the X-variables and P' is the loading matrix showing the influence of the variables on the projection model. E is the residual matrix expressing the deviations between the original values and the projections. In general, PCA transforms a number of correlated allele frequencies into a smaller number of uncorrelated synthetic variables, or PCs.

Cluster analysis (CA)

where l is the number of loci for which both individuals have been tested.

Discriminant function analysis

Discriminant functions based on population grouping were obtained by the stepwise inclusion of SNPs to minimise Wilks' lambda (λ) between groups, as described by Rechner [29] and as follows: L = b₁x₁ + b₂x₂ + b₃x₃ ... + b_zx_z; where x₁ through x_z represent the various predictor variables (SNPs); b₁ through b_z represent the weight associated with each of the predictor variables; and L is the object's resultant qualitative discrimination score, with a cut-off score to assign objects to one group or another. Objects with L > × are assigned to one group, and those with L < × are assigned to another group,[51] based on allele frequency differences. L represents classifying variables.

Fixation index (F_ST) estimates between Populations

Network and gene ontology analysis of genes showing differentiation between populations

Ingenuity Pathways Analysis (IPA) was used to organise genes showing evidence of selection into networks of interacting genes and to identify pathways containing functionally related genes [33]. More precisely, network analysis consists of searching for direct and indirect interactions between candidate genes and all other molecules (genes, gene products or small molecules) contained in the Ingenuity Pathways Knowledge Base (IPKB). The complete list of gene identifiers was uploaded into IPA, and each was mapped to its corresponding IPKB gene object. Candidate genes are eligible for network generation if there is at least one wild-type IPKB interacting molecule. Based on the information available for eligible candidate genes (focus genes), IPA further constructs networks by maximising the number of focus genes and their inter-connectivity within the limit of 35 molecules per network. Note that additional highly connected non-focus molecules are also included. Finally, for each network, a right-tailed Fisher exact test is implemented to evaluate how likely it is that the focus genes it contains might be found together by chance. Only those networks with a score (-log[p value]) greater than three were considered as significant [75]. In addition, networks might be inter-connected (ie sharing at least one molecule), which strengthens the importance for the underlying biological functions. Networks are graphically represented by nodes with various shapes (according to the molecule type) and edges (according to their biological relationships).

The likelihood for a gene pair to be regulated in the same manner increases with the similarity of their gene ontology (GO) description. The GO similarity score between two gene products is based on the number of shared ancestors. As a gene product might be assigned with multiple GO terms, we seek the maximum similarity score between all possible combinations. As we seek to discover gene-gene interactions, we reformulate the GO approach as follows. Let gene i and gene j be assigned h_i and h_j GO terms, respectively. Then, the GO similarity for the gene (i, j) pair is taken to be the maximum number of shared ancestors for all combinations of the h_i and h_j [76]. IPA essentially evaluates the enrichment of particular biological processes and molecular functions of gene sets by examining information collected by databases such as GO, Kyoto Encyclopedia of Genes and Genomes or the IPKB.