Clustering by genetic ancestry using genome-wide SNP data

16 Pages·2011·2.48 MB·English

Solovieffetal.BMCGenetics2010,11:108
http://www.biomedcentral.com/1471-2156/11/108
METHODOLOGY ARTICLE Open Access
Clustering by genetic ancestry using
genome-wide SNP data
Nadia Solovieff1*, Stephen W Hartley1, Clinton T Baldwin2, Thomas T Perls3, Martin H Steinberg4, Paola Sebastiani1
Abstract
Background: Population stratification can cause spurious associations in a genome-wide association study (GWAS),
and occurs when differences in allele frequencies of single nucleotide polymorphisms (SNPs) are due to ancestral
differences between cases and controls rather than the trait of interest. Principal components analysis (PCA) is the
established approach to detect population substructure using genome-wide data and to adjust the genetic
association for stratification by including the top principal components in the analysis. An alternative solution is
genetic matching of cases and controls that requires, however, well defined population strata for appropriate
selection of cases and controls.
Results: We developed a novel algorithm to cluster individuals into groups with similar ancestral backgrounds
based on the principal components computed by PCA. We demonstrate the effectiveness of our algorithm in real
and simulated data, and show that matching cases and controls using the clusters assigned by the algorithm
substantially reduces population stratification bias. Through simulation we show that the power of our method is
higher than adjustment for PCs in certain situations.
Conclusions: In addition to reducing population stratification bias and improving power, matching creates a clean
dataset free of population stratification which can then be used to build prediction models without including
variables to adjust for ancestry. The cluster assignments also allow for the estimation of genetic heterogeneity by
examining cluster specific effects.
Background test statistic for each SNP is divided by the same value,
In GWAS, varying ancestral backgrounds lead to popu- irrespective of whether the particular SNP is structured,
lation stratification inflating the type I error rate [1,2]. and thus results in a loss of power [4]. A wide range of
Even European Americans are affected by population alternative methods have been proposed to adjust for
stratification bias [3]. population stratification by applying an adjustment to
PCA [4], spectral graph theory [5], structured associa- the test statistic, using permutation tests or by perform-
tion analysis [6,7] and genomic control [8] can detect ing stratified analyses [9-11].
underlying population substructure with SNP data. Alternatively,populationstratificationcanbeobviated
Investigators typically use PCA to detect population by matching cases and control with respect to genetic
structure and then adjust the estimates of genetic effects ancestry[12].Matchingdoesnotrequiretheinvestigator
for the top number of principal components (PCs) to toaccountforpopulation stratificationbiasbyadjusting
control for population stratification bias [4]. Genomic forPCsordividingteststatisticsbythegenomiccontrol
control divides the test statistic for each SNP by the inflationfactor.Guanetal[13]matchcasesandcontrols
genomic control inflation factor (l), defined as the using a genetic similarity score computed directly from
observed median test statistic across all tests genome- genotypedataasaweightedidentitybystateestimate.
wide divided by the expected median test statistic. The We propose a different algorithm in which we identify
clusters of varying genetic ancestry using the results of
PCA from genome-wide data. The algorithm includes a
*Correspondence:[email protected]
1DepartmentofBiostatistics,BostonUniversitySchoolofPublicHealth, novel approach to choosing the appropriate number of
Boston,MA,02118,USA informative PCs, a clustering step to group subjects into
Fulllistofauthorinformationisavailableattheendofthearticle
©2010Solovieffetal;licenseeBioMedCentralLtd.ThisisanOpenAccessarticledistributedunderthetermsoftheCreative
CommonsAttributionLicense(http://creativecommons.org/licenses/by/2.0),whichpermitsunrestricteduse,distribution,and
reproductioninanymedium,providedtheoriginalworkisproperlycited.
Solovieffetal.BMCGenetics2010,11:108 Page2of16
http://www.biomedcentral.com/1471-2156/11/108
clusters of genetic diversity and a novel “scoring index” their top PCs. Since k-means does not provide a metric
(SI) to choose the appropriate number of clusters. Luca for choosing the appropriate number of clusters, we
et al [14] implemented a similar algorithm for matching developed an algorithm and scoring index to identify
cases and controls with the main difference being that the optimal number of clusters. The algorithm performs
parametric tests that may not be robust to departures M executions of k-means clustering for each cluster
from model assumptions are used while our algorithm size, k = 2,3,...K, and computes a measure of accuracy,
is non-parametric and does not require strict stability and between cluster distance for each execution
assumptions. and cluster size. We use these 3 measures because an
We test the sensitivity and specificity of our algorithm optimal cluster assignment should accurately assign
on simulated data and provide applications on African individuals based on the PCs, should be consistent from
populations from the Human Genome Diversity Project execution to execution for a particular number of clus-
[15] and a cohort of centenarians with European ances- ters and should also maximize the distance between
try [16]. We demonstrate that by matching cases and individuals in different clusters. We compare these three
controls within each cluster the inflation in test statistics measures computed in the data to those expected under
caused by population stratification bias substantially random cluster allocation using permutation analysis,
decreases as compared to unmatched cases and controls and summarize the gain of accuracy, stability and dis-
and in certain situations matching is more powerful tance into a score. The derivation of the components is
than adjusting for the top number of PCs. More impor- described in detail in the methods. In the following sec-
tantly, our method allows for the creation of a dataset tions we describe the results of the experiments we con-
free of population stratification which can then be used ducted to evaluate the properties of the algorithm.
for prediction modeling [17]. Prediction models contain-
ing principal components as covariates are study specific Simulations
and cannot easily be generalized to a different study in We first simulated genotype data containing no popula-
which the values of the PCs are unknown. The clusters tion substructure (see Methods). The heatmap of the
produced by the algorithm also allow for the exploration top PCs shows no distinct pattern and the scree plot is
of locus heterogeneity and we provide an example of a flat (Figure 1A) suggesting that the data do not contain
SNP in APOE in which the odds ratio varies widely substructure and clustering is unnecessary. To assess
across ethnic groups in Europe. the sensitivity of the method, we simulated genotype
data with the number of clusters ranging from 2 to 10
Results clusters under two scenarios, equal and unequal sample
Algorithm sizes across clusters (see Methods). The power of the
The algorithm that we developed consists of the follow- scoring index to identify the correct cluster size and
ing three components: 1) selection of informative PCs allocate subjects to the correct cluster was greater than
for the cluster analysis; 2) clustering to discover popula- 95% in most cases and in some cases 100% (Table 1).
tion substructure; 3) a novel score to automatically When the maximum scoring index was used to identify
select the best number of clusters that satisfy a variety the cluster size, the power was low for a cluster size of
of criteria. We select the informative PCs by visually 2 with equal sample sizes (78%) and a cluster size of
inspecting a heatmap and scree plot of the top PCs (Fig- 2 and 3 with unequal sample sizes (56%, 64%). In these
ure 1). A heatmap displays the top PCs in the columns situations the scoring index remained relatively constant
and the subjects in the rows. The intensity of the color across several cluster sizes and the maximum SI tended
corresponds to the value of the PC and the subjects are to overestimate the true number of clusters. However,
reordered so that individuals with similar PC values are using the optimal SI increased the power to 100% for a
displayed next to each other in the heatmap. The heat- cluster size of 2 with equal sample sizes and 99% and
map allows one to visually inspect the joint pattern pro- 97% for a cluster size of 2 and 3 with unequal sample
duced by a set of PCs and to choose the informative sizes. The optimal SI is defined as the smallest cluster
PCs. The scree plot graphs the natural logarithm of the size for which the estimated SI falls within the 95% con-
eigenvalue versus the component number. The informa- fidence interval of the maximum SI.
tive number of PC is identified by a “kink” in the plot
after which a straight line is observed. Based on the top Real Data with African Ancestry
PCs, we cluster individuals using k-means clustering. We then tested the sensitivity of our algorithm to group
K-means clustering assigns subjects to a pre-specified real data from individuals with known ancestry from
number of clusters by maximizing the distance between 7 African populations in the Human Genome Diversity
subjects in different clusters. Distance is defined as the Panel. The heatmap identified the top 7 PCs as most
Euclidean distance between subjects with respect to informative while the scree plot identified the top 5 PCs
Solovieffetal.BMCGenetics2010,11:108 Page3of16
http://www.biomedcentral.com/1471-2156/11/108
(cid:36) (cid:37) (cid:38)
(cid:19)(cid:20)(cid:21)(cid:22)(cid:23)(cid:24)(cid:25)(cid:26)(cid:27)(cid:28)(cid:19) (cid:19)(cid:20)(cid:21)(cid:22)(cid:23)(cid:24)(cid:25)(cid:26)(cid:27)(cid:28)(cid:19) (cid:19)(cid:20)(cid:21)(cid:22)(cid:23)(cid:24)(cid:25)(cid:26)(cid:27)(cid:28)(cid:19)
(cid:20)(cid:20)(cid:21)(cid:21)(cid:22)(cid:22)(cid:23)(cid:23)(cid:24)(cid:24)(cid:25)(cid:25)(cid:26)(cid:26)(cid:27)(cid:27)(cid:28)(cid:28)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:21)(cid:21) (cid:20)(cid:20)(cid:21)(cid:21)(cid:22)(cid:22)(cid:23)(cid:23)(cid:24)(cid:24)(cid:25)(cid:25)(cid:26)(cid:26)(cid:27)(cid:27)(cid:28)(cid:28)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:21)(cid:21) (cid:20)(cid:20)(cid:21)(cid:21)(cid:22)(cid:22)(cid:23)(cid:23)(cid:24)(cid:24)(cid:25)(cid:25)(cid:26)(cid:26)(cid:27)(cid:27)(cid:28)(cid:28)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:20)(cid:21)(cid:21)
(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38) (cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38) (cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)(cid:38)
(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51) (cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51) (cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)(cid:51)
Figure1ScreePlotandHeatmapofTop20PCs:Thescreeplotsandheatmapsarepresentedforthesimulateddatawithnostructure(left),
HGDPAfricans(middle),NECS(right).Thescreeplotgraphsthenaturallogarithmoftheeigenvalue(y-axis)versusthecomponentnumber(x-
axis).Akinkisnotobservedinthescreeplotforthesimulateddatawithnopopulationsubstructure(a).FortheHGDPAfricans(b)andthe
NECS(c)thescreeplotidentifiesakinkatthe5thand4thPCs,respectively.Intheheatmaps,eachcolumnisaPCandeachrowisanindividual.
OriginalPCsarestandardizedbyrowandthecolorintensitiescorrespondtothestandardizedvalueofthePCforeachindividual(green:higher
thanaverage,red:loweraverage)andaresortedbythecorrespondingeigenvalues,androwsaresortedbyhierarchicalclustering.Whileno
patternisfoundintheheatmapofthesimulatedexamplewithoutpopulationstructure,adistinctpatternisobservedfortheHGDPAfricansand
theNECS.ThepatternintheHGDPAfricansisthemostdistinctsincethesepopulationsarewelldefinedandareverydifferentfromeachother.
ThepatternfortheNECSismoresubtlebecausethevariabilityinsubjectsofEuropeanancestryismuchlowerthaninAfricans.However,one
canstillobservevariabilitybeyondthefirst2principalcomponents.
Table 1PowerofScoring Index (Figure 1). Since the pattern for the top 7 PCs is quite
Scenario1:EqualN Scenario2:UnequalN strong in the heatmap, we used the top 7 PCs instead of
K MaxSI OptimalSI MaxSI OptimalSI only the top 5 for clustering. The SI identified 8 clusters
2 0.78 1 0.56 0.99 as the optimal cluster size (Figure 2). The Biaka, Man-
3 0.97 1 0.64 0.97 denka, Mbuti Pygmy and San populations are assigned
4 0.99 1 0.86 0.99 to their own cluster. The Yorubans are clustered with
5 1 1 0.99 1 the Bantu and 2 of the Mozabite subjects in the 8th clus-
6 0.99 1 0.99 0.99 ter and show a stronger similarity with each other than
7 1 1 0.99 0.99 with the other Africans in this analysis. The plots of the
8 1 1 0.99 0.99 PCs support the clustering of these subjects (Figure 3)
9 1 1 0.99 0.99 and are consistent with the origin of the Bantu near the
10 1 1 1 1 southern boundary of modern Nigeria and Cameroon
and the known genetic similarity between Bantu and
Wereporttheproportionofdatasetsforwhichthecorrectnumberofclusters
isidentifiedandthesubjectsarecorrectlyallocatedtheirrespectivecluster. Yorubans [18,19]. In the plot of PC6 and PC7 (Figure 3)
ThemaximumSIcorrespondstotheclustersizeatwhichtheSIismaximized
we see a distinct separation of the Mandenka from the
andtheoptimalSIcorrespondstothesmallestclustersizeforwhichtheSI
fallswithinthe95%CIofthemaximumSI. Yorubans and Bantu which indicates that these PCs are
Solovieffetal.BMCGenetics2010,11:108 Page4of16
http://www.biomedcentral.com/1471-2156/11/108
Figure2Accuracy,Stability,BetweenClusterDistanceandScoringIndex(SI)forHGDPAfricans.Accuracy,stabilityandbetweencluster
distanceonthek-meansclusterassignmentsaredisplayedinredandthemeasuresonthepermutedclusterassignmentsaredisplayedinblue.
TheSIaveragestherelativegainintheaccuracy,stabilityandbetweenclusterdistancesandmaximizesat8clusters(0.931,CI:0.923,0.938)).
Notethatthegraphicaldisplayofaccuracy,stabilityanddistanceshowsthatnoneofthescorecomponentswouldbesufficienttoidentifythe
correctnumberofclusters.Forexample,highaccuracyisnotsufficienttoconcludethattheclusteringisoptimal.Inthisexample,theaccuracy
isnearlyperfectwhenthenumberofclustersislessthan7suggestingthatanynumberofclusterslessthan7isequallyoptimal,howeverthese
numbersarenotallequallyoptimalasdemonstratedbythemeasuresofstabilityandbetweenclusterdistancesthatcontinuallyincrease.
important to increase the sensitivity of clustering. The confounding in GWAS and confirms the utility of the
remaining Mozabite subjects were split among 3 clus- heatmap plot of PCs.
ters. This genetic diversity is consistent with the history
of the Mozabites and recent studies that reveal genetic Real Data with European Ancestry
ancestry from sub-Saharan Africa, the Middle East and In a cohort of 1051 centenarians and 290 controls from
Europe [18]. Previously reported analyses of these data the New England Centenarian Study, the algorithm
that were based on inspection of the first two PCs failed detected many specific European ethnic groups. The
to detect this level of diversity based on purely genetic heatmap and the scree plot identified the top 4 PCs as
data. Patterson et al [19] could only detect differences informative (Figure 1) and the algorithm identified
between San and Bantu/Yorubans combined. Our analy- 9 clusters (Figure 4). Based on partial information of
sis shows that a larger number of PCs is necessary to ancestry (birth places and native language of grandpar-
detect finer population structure that could introduce ents) for the subjects in the NECS, we found that the 9
Solovieffetal.BMCGenetics2010,11:108 Page5of16
http://www.biomedcentral.com/1471-2156/11/108
Figure3PairwisePCPlotsofHGDPAfricansColoredbyClusterAssignment.Thelegendreportstheclusterassignment,populationand
samplesizeforeachcluster.TheBiaka,Mandenka,MbutipygmyandSanareseparatedintothefirst4clusters.Theeighthcluster(black)
containstheYorubans,Bantuand2subjectsfromtheMozabitepopulationandthesesubjectsclustertogetherineachoftheaboveplots
indicatingthatthesecannotbedifferentiatedbasedonthefirst7PCs.TheremainingMozabitearesplitbetweenclusters5-7andshow
variabilityinPC5,PC6andPC7suggestingthattheseindividualsaremoreheterogeneousthantheotherAfricanpopulations.
clusters correlated with distinct ethnic groups in Europe formed by PC1 and PC2 in Figure 5 that identifies
(Figure 5) demonstrating that the clustering algorithm AshkenaziJewsandaNorthwesttoSoutheastclineacross
can detect fine distinctions among subjects of European Europe [3,20]. The pattern formed by PC3 and PC4,
descent. The population substructure in subjects of generallynotreported,alsocorrelatesclearlywithspecific
European descent is typically identified as the pattern ethnic groups in Europe such as the Anglo-Saxons,
Solovieffetal.BMCGenetics2010,11:108 Page6of16
http://www.biomedcentral.com/1471-2156/11/108
Figure4Accuracy,Stability,BetweenClusterDistanceandScoringIndexforNECS.Accuracy,stabilityandbetweenclusterdistanceonthe
k-meansclusterassignmentsaredisplayedinredandthemeasuresonthepermutedclusterassignmentsaredisplayedinblue.Thescoring
indexmaximizesat11clusters(SI=0.886,95%CI:(0.876,0.896)).Themeanscorefor9clusters(0.879)fallswithinthe95%confidenceintervalof
11clustersandthusisasgoodas11clustersbutismoreparsimonious.andisalsotheoptimalsolution.Ingeneral,weexpectthepermuted
stabilitytoincreasewithanincreasingnumberofclustersbecausealargernumberofpairswillbeconsistentlyassignedtodifferentclusters
sincetherearemoreavailableclusters.However,notethepeculiartrendofclusterstabilityinthisexample:thepermutedstabilitydecreases
whenmovingfromk=2to4clustersandthengraduallyincreases.Thisoccursbecauseoneclustercontainsthemajorityofthesubjectsfor
theseclusterssizesandthusthenumberofpairsofsubjectsassignedtothesamecluster,bychance,ismuchhigherthaninthecasewhere
thereareequalsamplesizesineachcluster.Thusitisimportanttoaccountforthepermutedmeasurementswhenselectingtheoptimal
numberofclusters.
Celto-Germanic, Scandinavians, Saxons, Irish-Celtic, matched subjects from a sample of approximately 3,600
Italians, Ashkenazi Jewish, Slavic, Slavo-Germanic and subjects from the Illumina database. We conducted a
Celtic. combined PCA of the NECS centenarians and controls,
and ~3,600 Illumina controls and the heatmap of the
Effect of Matching on the Power of a GWAS top 20 PCs suggested that the top 4 PCs are the most
Understanding the genetic diversity of the NECS allows informative (Figure 6). The clustering algorithm identi-
one to select genetically matched controls and reduce fied 8 clusters as the optimal solution (SI = 0.919, 95%
population stratification bias. For example, in our CI: (0.907, 0.931)). We matched cases to controls within
GWAS of exceptional longevity (EL) [17] we set out to each of the 8 clusters to balance the overall proportion
expand the NECS control sample by drawing cluster- of cases and controls across the clusters, resulting in the
Solovieffetal.BMCGenetics2010,11:108 Page7of16
http://www.biomedcentral.com/1471-2156/11/108
Figure5PlotofFirst4PCsofNECSColoredbyClusterAllocation.PC1andPC2clearlyseparatethesubjectsofJewishdescent,shownin
purple,fromtherestofthesubjectsandcreateanorthwesterntosoutheasternEuropeancline.PC3andPC4capturefurthervariabilitywithin
Europeanethnicgroups.TheclusterassignmentsfoundbythealgorithmcloselycorrespondtoestablishedethnicgroupsacrossEurope.Based
onpartialinformation(birthplacesandnativelanguageofgrandparents)ofthesubjects’ancestry,theclusterswerelabeledaccordingtothe
mostcommonethnicgroupintheregion[37].
addition of 745 Illumina controls to the existing NECS PC adjusted analysis (l = 1.04) is slightly lower than for
controls set. We then performed two sets of GWAS for the matched analysis, however the difference is not
EL, one with the matched data and a second analysis substantial.
with an equal number of randomly selected Illumina
controls, using additive models. The QQ plot in Figure 6 Effect of Matching on Power
displays the results for SNPs genome-wide and demon- To investigate the effect of matching on the power of a
strates the striking reduction in the inflation of test sta- GWAS, we simulated genotype data to compare the
tistics for the matched controls as compared with the power and false positive rate of a GWAS in which we
unmatched controls. The genomic control factor l match cases and controls and a GWAS in which the
decreases from 1.44 for the unmatched set to 1.07 for subjects are unmatched but the analysis is adjusted for
the matched set. Matching does not require the investi- the top 10 PCs (see Methods). The false positive rate
gator to account for population stratification bias by across all scenarios was close to the nominal level
adjusting for PCs or by dividing test statistics with the (Table 2). The GWAS using the matched subjects was
genomic control inflation factor which can decrease the more powerful than the PC adjusted unmatched analysis
power of detecting true associations [4]. For example, an with the same sample size (Figure 8). The gain in power
unmatched PC-adjusted analysis with the same number for matched analysis ranged from 1.0% to 23.4%. When
of controls randomly chosen from the Illumina database matching is performed, a number of controls are
would lose almost all of the top associations suggesting excluded due to the inability to match them. To exam-
a loss in power. (Figure 7) Furthermore, in the QQ plot ine the effect of adding in all controls versus using a
(Figure 6) the tail of the distribution of p-values for the smaller pool of matched controls, we added a varying
adjusted unmatched GWAS is lower than the tail for number of controls to the unmatched analysis (see
the matched GWAS which again suggests a loss in Methods). As expected the power of the PC adjusted
power. The genomic control factor for the unmatched unmatched analysis increases as more controls are
Solovieffetal.BMCGenetics2010,11:108 Page8of16
http://www.biomedcentral.com/1471-2156/11/108
Table 2False Positive Rate ofMatched andPCadjusted
(cid:50)(cid:69)(cid:86)(cid:72)(cid:85)(cid:89)(cid:72)(cid:71)(cid:3)(cid:177)(cid:79)(cid:82)(cid:74)(cid:20)(cid:19)(cid:11)(cid:83)(cid:12) (cid:104)(cid:87)(cid:68)(cid:104)(cid:18)(cid:374)(cid:374)(cid:258)(cid:3)(cid:373)(cid:373)(cid:4)(cid:410)(cid:272)(cid:282)(cid:258)(cid:258)(cid:346)(cid:361)(cid:410)(cid:410)(cid:856)(cid:286)(cid:272)(cid:272)(cid:282)(cid:346)(cid:346)(cid:286)(cid:286)(cid:282)(cid:282)(cid:3) (cid:644)(cid:644)(cid:644)(cid:1089)(cid:1089)(cid:1089)(cid:1005)(cid:1005)(cid:1005)(cid:856)(cid:856)(cid:856)(cid:1008)(cid:1004)(cid:1004)(cid:1008)(cid:1011)(cid:1008) (cid:51)(cid:38)(cid:20)(cid:51)(cid:38)(cid:21)(cid:51)(cid:38)(cid:22)(cid:51)(cid:38)(cid:23)(cid:51)(cid:38)(cid:24)(cid:51)(cid:38)(cid:25)(cid:51)(cid:38)(cid:26)(cid:51)(cid:38)(cid:27)(cid:51)(cid:38)(cid:28)(cid:51)(cid:38)(cid:20)(cid:19)(cid:51)(cid:38)(cid:20)(cid:20)(cid:51)(cid:38)(cid:20)(cid:21)(cid:51)(cid:38)(cid:20)(cid:22)(cid:51)(cid:38)(cid:20)(cid:23)(cid:51)(cid:38)(cid:20)(cid:24)(cid:51)(cid:38)(cid:20)(cid:25)(cid:51)(cid:38)(cid:20)(cid:26)(cid:51)(cid:38)(cid:20)(cid:27)(cid:51)(cid:38)(cid:20)(cid:28)(cid:51)(cid:38)(cid:21)(cid:19) UnUmnmaSMatctcacehthnceeahddreido+GPWCAsS Cas22e00s00/00C//o2200n00tr00ols FalseP85o..00s00iEEti--v00e77Rate*
Unmatched+PCs 2000/2250 9.00E-07
(cid:40)(cid:91)(cid:83)(cid:72)(cid:70)(cid:87)(cid:72)(cid:71)(cid:3)(cid:177)(cid:79)(cid:82)(cid:74)(cid:20)(cid:19)(cid:11)(cid:83)(cid:12) Unmatched+PCs 2000/2500 1.30E-06
Unmatched+PCs 2000/2750 1.20E-06
Unmatched+PCs 2000/3000 9.00E-07
(cid:54)(cid:70)(cid:82)(cid:85)(cid:76)(cid:81)(cid:74)(cid:3)(cid:44)(cid:81)(cid:71)(cid:72)(cid:91) UUUnnnmmmaaatttccchhheeeddd+++PPPCCCsss 222000000000///333257505000 189...000000EEE---000677
Unmatched+PCs 2000/4000 8.00E-07
(cid:49)(cid:88)(cid:80)(cid:69)(cid:72)(cid:85)(cid:3)(cid:82)(cid:73)(cid:3)(cid:38)(cid:79)(cid:88)(cid:86)(cid:87)(cid:72)(cid:85)(cid:86) *Falsepositiveratewithap-valuecut-offof1E-06
Figure6Q-QPlotofGWAresults,ScoringIndexandHeatmap WereportthefalsepositiverateofamatchedGWASandcompareittoPC-
forNECSwithMatchedandUnmatchedControls.TheQQplot adjustedunmatchedGWASswithvaryingnumberofcontrols.Thefalse
positiverateiscomputedusingap-valuecut-offof1E-06.
(topleft)showsthesubstantialdecreaseoftheinflationintheGWA
resultswhencontrolsarematchedtocaseswithineachclusterwith
thegenomiccontrolfactor,l,decreasingfrom1.44to1.07.When
added and surpasses the power of the matched analysis
adjustingforthetop4PCsintheunmatchedanalysis,ldecreases
when 1000 to 1500 subjects are added (Figure 8). This
to1.04,however,manyofthetopassociationsaremissed(Figure
7).Thescoringindex(bottomleft)forthecombinedanalysisofthe result suggests that using all controls and adjusting for
NECSandcontrolsfromtheIlluminadatabasemaximizedand PCs as opposed to matching will become more powerful
optimizedat8clusters(SI=0.919,95%CI:(0.907,0.931)).The if the number of unmatched controls is substantially lar-
heatmapdisplaysapatternforthetop4PCs.
ger than the number of matched controls.
Locus Heterogeneity
Breaking the subjects into clusters also has the advan-
tage of examining cluster specific allele frequencies and
genetic effects. As an example, the allele frequencies for
rs405509, a SNP near APOE known to be associated
with exceptional longevity and other age-related diseases
[21,23], varies greatly among European ethnic groups
(Figure 9). Interestingly the association between EL and
this SNP also varies between clusters with odds ratios
ranging from 0.38 in the Italians to 1.09 in the Irish Cel-
tics. (Figure 9) Small differences in allele frequencies can
affect the power to detect a true main effect when the
SNP is part of a SNP- SNP interaction [24] and thus
close examination of cluster specific allele frequencies
between discovery and replication sets can be useful for
replication.
Discussion
Our analyses suggest that the algorithm works well in
Figure 7 GWAS Results of EL for Matched and Unmatched
Subjects.Thex-axisreportsthe-log10p-valuesfortheGWASin simulated datasets and in real data of Africans and Cau-
whichsubjectswerematchedbasedon8clusters.They-axis casians of European ancestry and can reduce or elimi-
reportsthe-log10p-valuesfortheGWAScontrollingforthetop4 nate population stratification. Choosing the appropriate
PCsinwhichanequalnumberofsubjectswererandomlychosen
number of PCs for clustering is a critical choice. We
fromtheIlluminacontrols.Thedottedlineisdrawnatavalueof7.
chose the number of PCs by evaluating patterns
TheSNPsinthebottomrighthandrectangleareSNPswhichwould
beselectedasthetopSNPsinthematchedandwouldbemissed observed in a heatmap and in a scree plot; although
intheunmatchedPCadjustedanalysis.Thereare19SNPswithap- subjective, in many situations a clear choice is evident.
valuelessthan10E-07forthematchedanalysisforwhichthep- We emphasize that investigators should carefully exam-
valueisgreaterthan10E-07fortheunmatchedPCadjusted
ine the results of the PCA prior to clustering to ensure
analysis.
that population substructure in fact exists in the data. If
Solovieffetal.BMCGenetics2010,11:108 Page9of16
http://www.biomedcentral.com/1471-2156/11/108
Figure 8 Power of Matched versus Unmatched PC Adjusted GWAS. We compared the power of a matched GWAS (2000 cases/2000
controls)toanunmatchedPCadjustedGWASwithavaryingnumberofcontrols(2000to4000).Thex-axisreportsthenumberofcontrols
addedtotheunmatchedcontrolgroup(n=2000)andtheletterMcorrespondstothematchedGWAS.PowerwasassessedforaMAFranging
between0.1to0.5andarelativerisk(RR)rangingfrom1.2to1.5.Inallscenariosthematchedanalysisismorepowerfulthantheunmatched
analysiswhenthesampleisthesameforbothanalyses.Asweaddmorecontrolstotheunmatchedanalysis,thepowerincreasesandsurpasses
thepowerofthematchedanalysiswhen1000to1500additionalcontrolsareadded.Notethatthepowerisadjustedforthevaryingfalse
positiverateacrossscenarios.
population substructure is absent, clustering is unneces- informative clusters are created and will increase by a
sary and will lead to overfitting. Although we use smaller amount when the clustering method begins to
k-means clustering in our algorithm [25], alternative create uninformative clusters [25]. However, one does
techniques including k-medoids, hierarchical clustering, not always observe a distinct “kink” (see Figure 4). An
and model based clustering could be used. Model based alternative statistic called the gap statistic [26] monitors
clustering requires parametric assumptions about the the within cluster distance, compares it to a null distri-
distribution of the data and it can outperform k-means bution and chooses the number of clusters for which
clustering when the assumptions are correct but per- the deviation of the within cluster distance from the
form poorly when they are not. null distribution is maximized. In practice, a maximum
The SI provides a novel statistic for choosing the best may never be attained for a reasonable range of cluster
number of clusters by combining measures of accuracy, sizes as demonstrated in the examples of the HGDP
stability and between cluster distances. Maximizing only Africans and NECS (Figure 10).
one or two of these measures is insufficient for choosing We have shown that matching is an effective method
an optimal number of clusters. For example, the cluster- to reduce the inflation in test statistics due to popula-
ing may be accurate and stable for a given number of tion stratification both in simulated cases and real data.
clusters, however there may be a different number of Matching in some situations is more powerful than
clusters for which the accuracy and stability are compar- using an unmatched dataset and adjusting for PCs. Our
able but the between cluster distance is better. Some method also allows one to create a dataset free of popu-
advocate choosing the number of clusters at which a lation stratification which can be used to build predic-
“kink” occurs in the between cluster distance since the tion models. We created a clean dataset for the NECS
measure will increase by a larger amount when which was used to build a model of exceptional
Solovieffetal.BMCGenetics2010,11:108 Page10of16
http://www.biomedcentral.com/1471-2156/11/108
Furthermore as the availability of external controls
increases through organizations like dbSNP, investiga-
tors can use this method to match the controls to cases
with respect to ancestry and reduce the issue of popula-
tion stratification. Zhuang et al [27] show that expand-
ing the control group can improve power. Of course,
investigators must take precautions to avoid introducing
other forms of bias due to the use of external controls.
As next generation sequencing and the discovery and
analysis of rare variants continue to emerge, population
stratification will remain an obstacle. Our algorithm can
easily be implemented to match cases and controls
when selecting subjects for sequencing and insure that
the selection adequately represents all ethnicities in the
sample thus increasing the chance of finding true novel
polymorphisms.
Conclusions
Figure9Frequencyofrs405509inCentenariansin8Clusters. We developed an algorithm to cluster individuals into
Theallelefrequenciesofrs4005509,aSNPnearAPOE,varygreatly ethnic groups based on PCs computed from SNP data.
betweenthe8clustersofcentenariansidentifiedintheanalysisof We showed that the algorithm works well in real data
theNECSandIlluminasubjects.Theoddsratioforeachclusteris
of Africans and Caucasians with European ancestry and
calculatedusingaseparatelogisticregressionmodelofexceptional
also in simulated data. The cluster assignments can be
longevityandtheSNPgenotypecodedasthenumberofAalleles.
Theoddsratiosarequitedifferentamongvariousethnicgroups used to effectively match cases and controls in GWAS
indicatingthatAPOEmayhaveadifferentimpactindifferentethnic to reduce and even eliminate population stratification
groups. bias. Matching can also aid in genetic risk prediction
models by creating a dataset free of population stratifi-
cation. Furthermore, the cluster assignments can be
longevity and we replicated the model in an indepen-
used to study locus heterogeneity and identify SNPs and
dent study with high accuracy [17]. The advantage of
genes that have a different effect on the phenotype in
this approach is that we did not need to account for
different ethnic groups.
population substructure in the model and were able to
generalize it to an independent study. Including PCs in
Methods
a prediction model causes the model to be study specific
Simulated Datasets
and may be difficult to reproduce.
Dataset with No Population Substructure
To test the algorithm on a dataset containing no popu-
lation substructure, we simulated 100,000 SNPs in link-
age equilibrium for 2000 individuals with no underlying
pattern. The allele frequencies for each SNP were ran-
domly simulated with values ranging between 0.05 and
0.95. Genotype frequencies were computed from the
allele frequencies assuming Hardy Weinberg equilibrium
and individuals were randomly assigned genotypes
according to the genotype frequencies.
Dataset with Known Structure
We simulated genotype data in which the number of
clusters ranged from k = 2 to 10. Each dataset contained
Figure10GapStatisticAcrossVaryingClusterSizesforHGDP 1000 individuals and 10,000 SNP in linkage equilibrium.
AfricansandNECS.Thegraphplotsthemeangapstatistic(see
Each cluster was simulated with Fst = 0.01 correspond-
Methods)ateachclustersize.ThegapstatisticfortheHGDP
ing to differences observed among divergent European
Africanssteadilyincreasesto21clustersatwhichpointthestatistic
cannotbecomputedduetosingletonclusters.FortheNECS,the populations [28]. For each SNP the founder allele fre-
statisticagainincreaseswithanincreasingclustersizeanddoesnot quency, p, was selected from a uniform(0.05,0.50) distri-
reachapointafterwhichthestatisticcontinuouslydecreases.In bution, the allele frequency for each cluster, p , was
k
bothsituationschoosinganoptimalclustersizeisnotobvious.
selected from a beta distribution with shape parameter