Reconstructing Past Admixture Processes from Local Genomic Ancestry Using Wavelet ...

36 Pages·2015·1.89 MB·English

HIGHLIGHTED ARTICLE
GENETICS | INVESTIGATION
Reconstructing Past Admixture Processes from Local
Genomic Ancestry Using Wavelet Transformation
JeanSanderson,*,1Herawati Sudoyo,†TatianaM.Karafet,‡Michael F.Hammer,‡,§andMurrayP.Cox*,2
*StatisticsandBioinformaticsGroup,InstituteofFundamentalSciences,MasseyUniversity,PalmerstonNorth4442,NewZealand,
†EijkmanInstituteforMolecularBiology,Jakarta,Indonesia,‡DivisionofBiotechnology,ArizonaResearchLaboratoriesand
§DepartmentofAnthropology,UniversityofArizona,Tucson,Arizona85721
ABSTRACTAdmixturebetweenlong-separatedpopulationsisadefiningfeatureofthegenomesofmanyspecies.Themosaicblock
structureofadmixedgenomescanprovideinformationaboutpastcontactevents,includingthetimeandextentofadmixture.Here,
wedescribeanimprovedwavelet-basedtechniquethatbettercharacterizesancestryblockstructurefromobservedgenomicpatterns.
principal components analysis is first applied to genomic data to identify the primary population structure, followed by wavelet
decomposition to develop a new characterization of local ancestry information along the chromosomes. For testing purposes, this
method is applied to human genome-wide genotype data from Indonesia, as well as virtual genetic data generated using genome-
scalesequentialcoalescentsimulationsunderawiderangeofadmixturescenarios.Timeofadmixtureisinferredusinganapproximate
Bayesian computation framework, providing robust estimates of both admixture times and their associated levels of uncertainty.
Crucially,wedemonstratethatthisrevisedwaveletapproach,whichwehavereleasedastheRpackageadwave,providesimproved
statistical poweroverexisting wavelet-basedtechniques andcanbeusedtoaddressabroadrangeofadmixturequestions.
KEYWORDSwavelets;principalcomponentanalysis(PCA);admixture;localancestry;dating
ADMIXTURE occurs when previously separated popula- onrawgenomicdata(Falushetal.2003;Priceetal.2009).
tionsinteractandmerge.Thisprocesshasbeeninstru- Somemethodsassignancestrydirectly.Forinstance,HAPMIX
mentalinhumanhistory,withmostglobalgroupsshowingat uses a hidden Markov model to estimate the break points
least some signals of population merger (Hellenthal et al. ofancestryblocks,whileotherapproachesdefineancestry
2014). The admixture process produces “mosaic” genomes blocks using simple empirical criteria, such as strings of
with alternatingblocksofDNAfromeachancestralpopula- shared vs. nonshared polymorphisms (Pool and Nielsen
tion.Overtime,recombinationdecreasesthelengthofthese 2009)orthedifferentialpresenceofpopulation-specificvar-
ancestryblocks,andthereforethedistributionofblocksizesis iants(BrownandPasaniuc2014).Anothersetofmethods
informativeaboutthetimeofadmixture.However,theextent is more indirect. ROLLOFF (Moorjani et al. 2011), LAMP
to which these patterns can provide additional information (Baranetal. 2012),andALDER (Loh etal. 2013)allsearch
about historic admixture processes is still a young area of for rapid changes in linkage disequilibrium to define the
exploration. bordersofancestryblocks,whileotherapproachesassign
Arangeofmethodshavebeendevelopedtopartitionthe ancestryforpredefinedgenomicwindowsusingconditional
genomeofanadmixedindividualintoancestryblocksbased random fields (Maples et al. 2013) or principal component
analysis (PCA) (Gravel 2012).
Copyright©2015bytheGeneticsSocietyofAmerica Thesemethodsvaryintheireffectiveness.Simpleempir-
doi:10.1534/genetics.115.176842 icalcriteriaperformsurprisinglywellforadmixturebetween
Manuscript received October 29, 2014; accepted for publication April 3, 2015;
publishedEarlyOnlineApril7,2015. species (as for the mouse admixture zone studied by Pool
Availablefreelyonlinethroughtheauthor-supportedopenaccessoption. andNielsen2009).Similarly,mostofthesemethodstendto
Supportinginformationisavailableonlineathttp://www.genetics.org/lookup/suppl/
behighlyaccurateforrecentadmixturebetweenwell-separated
doi:10.1534/genetics.115.176842/-/DC1.
1Presentaddress:SchoolforHealthandRelatedResearch,UniversityofSheffield, human groups (such as African Americans or American La-
Sheffield,UnitedKingdom.
tinos). Indeed, in these settings, subtleties such as multiple
2Correspondingauthor:InstituteofFundamentalSciences,MasseyUniversity,Private
Bag11222,PalmerstonNorth4442,NewZealand.E-mail:[email protected]. waves of admixture have even be detected (Gravel 2012).
Genetics,Vol.200,469–481 June2015 469
However, reconstructing complex demographic features for directly to the SNP-level data, without the need to compute
much older admixture events (i.e., thousands rather than averagesoverlocalizedgenomicwindowsasimplementedin
hundredsofyearsinthepast)remainsextremelychallenging related procedures (Pugach et al. 2011; Brisbin et al. 2012).
(Moorjanietal.2011).Whilemethodshaveinprinciplebeen Instead,windowingisperformednaturallyandobjectivelyas
proposed to detect multiple waves of ancient admixture, in part of the wavelet decomposition procedure. We show that
manyrealisticsettingstheyarestillrestrictedtosingleadmix- thisnewmethodprovidesrobustestimatesofadmixturetime
ture events (Loh et al. 2013), although some evidence for (includingimprovedcontrolofuncertaintyestimates),aswell
multiple ancient admixture events has been presented for as recognizing other aspects ofadmixture processes that pre-
several Indianpopulations (Moorjanietal. 2011). vious wavelet-based methods have not been able to identify
Other indirect methods look increasingly promising in with anyaccuracy.
this “old admixture” space. Approaches based on principal
componentsanalysisandwaveletshavebeenemployedwith Methods
some success. PCA is a nonparametric data-reduction tech-
Generalframework
nique, which has been used widely to identify patterns of
population structure in genetic data (Patterson et al. 2006; Initially,weconsiderasimpleadmixturescenariowheretwo
Novembre and Stephens 2008; McVean 2009; Bryc et al. ancestral populations P and P merged T generations ago
A B
2010; Ma and Amos 2012). Dispersion of admixed individ- to form the admixed population P . The ancestral popula-
C
ualsalongthefirstprincipalcomponentconnectingancestral tions contribute to the admixed population with probabili-
populations can be used as a diagnostic for two-way admix- tiespand12p.Thesizesofthepopulations,theadmixture
ture (Patterson et al. 2006; Mcvean 2009). For instance, time, and the admixture proportions are free to vary.
PCAdmix employs PCA to assign ancestry to localized win- To quantify patterns of genomic block size variation,
dows along the genome for each individual (Brisbin et al. athree-stepanalysisprocedurewasused:(i)PCAwasapplied
2012).Pugachetal.(2011)alsousePCA,butdonotdirectly tothegenomicdatatodescribepopulationstructure;(ii)the
assignancestrytogenomicregions,insteadapplyingawavelet wavelet variance was computed to provide a scale-by-scale
transform to obtain an indirect measure of the average ad- decomposition of the variance for each population; and (iii)
mixtureblocklength.Whilethisapproachhasbeenshownto theportionofthismeasurethatisinformativeforadmixture
be powerful for dating old admixture events, there remains processes was extracted relative to background levels ob-
considerable scope for (i) the development of more sophisti- served in the ancestral populations.
cated wavelet constructions, (ii) examining the resulting
Datasimulation
waveletdecompositionsingreaterdetail(particularlytoiden-
tifyaspectsofnon-time-relatedinformationinthetransformed Genome-wideSNPdataweresimulatedusingthesequential
data),and(iii)toprovideamoreuser-friendlysoftwaresolu- coalescent simulator MaCS (Chen et al. 2009). Because our
tionforwaveletanalysis. primaryinterestisintheadmixturehistoryofIslandSouth-
Wavelettechniquesthemselvesareanactiveandevolving east Asia (see Real genomic data section below), we chose
area,withmuchpotentialfornovelapplicationinpopulation parameter settings that produce genomic data that broadly
genetics, as highlighted in the review article by Liò (2003). fitobservedpatternsofgeneticdiversityinthisstudyregion
Waveletscanbethoughtofaslocalized,oscillatoryfunctions (Coxetal.2008).Thedemographicmodel,parameters,and
and are particularly useful for representing data that has informationsourcesaredescribedinmoredetailintheSup-
localfeaturessuchassharpchangesanddiscontinuities.In porting Information (Figure S1). We emphasize, however,
thecontextofgenome-widesinglenucleotidepolymorphism thatthemethodwedescribeisgeneralandcanbeappliedto
(SNP)data,waveletscanbeusedtorepresentthemosaic most admixed genomic systems.
pattern of ancestry blocks. A wavelet decomposition of
Datasetup
the data provides information on the size of the ancestry
blocks and, importantly, how they are distributed along GivenanadmixedpopulationP derivedfromtwoancestral
C
thechromosomes. Summarymeasures ofthe wavelet de- populations P and P , the number of individuals in the
A B
composition allow aspects of the admixture process to be analysis (i.e., present day samples) is n = n = n + n .
A B C
reconstructed,suchasthetimeofadmixtureandadmixture For each individual i, we observe a collection of T SNPs
proportions. along a chromosome. Thus the raw data matrix X is
Here, we present a substantially revised wavelet-based a T3n matrix with T genotype counts in columns and n
approachtodescribepopulationadmixturethatbuildsonthe individualsinrows.TheSNPssareorderedbytheirphysical
work of Pugach et al. (2011). This new method has signifi- positions along the chromosome, with the cells of the data
cantly fewer model assumptions and allows us to identify matrixXs;itakingthevalue0ifheterozygous,andarbitrarily
more complex demographic processes, such as multiple 21or1forthealternativehomozygousstates.Priortoprin-
admixture events. As with previous methods, PCA is first cipal componentsanalysis, the datamatrixis centered such
employedtodescribethepopulationstructure.Themaximal that the column mean with respect to the ancestral refer-
overlapdiscretewavelettransform(MODWT)isthenapplied ence populations is zero, giving
470 J.Sandersonetal.
X
X9s;i ¼Xs;i2n þ1 n Xs;i:
A B i2PE;PB
Principal componentsanalysis
PCA is performed using only individuals from the ancestral
populations. Rather than performing PCA on all samples
combined,thisapproachhastheadvantagethatotherfeatures
of the admixed sample (such as admixture from additional
ancestral populations) will not influence the projection
(McVean2009).Thefirsteigenvectorv reflectstheprimary
1
populationstructure.Projectionofindividualsontothisaxis
of variation is given by
X
y1i ¼ Ts¼1X9s;iv1;s: (1)
TheproportionofancestryinheritedfrompopulationP can
A
be estimated for each individual (or population) using the
distancefromthecentroidsoftheancestralpopPulations;that
ics,¼pið¼1=ðncBÞ2Py1iÞ=ðyciBa2recAthÞe; wcehnetrroeidcsAo¼ftðh1e=nanAcÞestir2aPlEpyo1ipuanlad-
tiBonsalongBthei2fiPrBst1principalaxis(Brycetal.2010).Notethat
variation between individuals within a population is repre-
sentedbythesmallereigenvaluesandcorrespondingeigenvec-
tors(MaandAmos2010).
This representation of admixed individuals in PCA space, Figure1 Simulatedexamplewith13,000SNPs,15diploidindividualsin
as shown in Figure 1A, provides a genome-wide estimate of ancestralpopulations(P ,P ),and20diploidindividualsintheadmixed
A B
averageancestry,butdoesnotindicatehowadmixturetracts population(PC).Populationsareshowningreen(PA),blue(PB),andred
(P ).(A)PCAisusedtodescribetheprimarypopulationstructure;(B)raw
are distributed along the chromosomes. To obtain localized C
waveletvarianceforeachpopulationillustrateshighfrequencynoise;(C)
estimates,theprojectionisperformedattheSNPlevelrather
informativevariationintheadmixedpopulationafterstandardcorrection
thansummingoverthelengthofthegenomeasinEquation fornoiseestimatedfromtheancestralpopulations.Notethatthisexam-
1.The raw SNP-leveladmixture signalsare given by pleusesthedefaultthresholdm=1.
8 (cid:2) (cid:2)
Ysi ¼>>>><>>>>:2Xs9;ivð1Y;ssB22ðYYsBsAÞþYsA0Þ;; (cid:2)(cid:2)(cid:2)(cid:2)(cid:2)(cid:2)(cid:2)YYsBsB22YYsAsA(cid:2)(cid:2)(cid:2)(cid:2)(cid:2)(cid:2)(cid:2)$,ee; (2) cwSsiajtapifovtoenurlreeostcsfavttlhheasereiajhvn¼aicgre1hia;en(.scP.te.ef;rroJceifqvputarhoelenv1icdsy9ieg9spn5aaa)t.lts.ecTTranhhlesee,-brfiwyer-apssvtcreaessllceeeatndlteevinca(gorjim¼avnpe1croye)-
local information. Increasing the scale index provides suc-
P
wtehrmerseiYnGsE¼quða1t=ionnGÞ2 ein2PsGuXrse9;ivt1h;satfotrhGe s2igAn;alsB.aTrheenaodrmdiatiloizneadl ctoes“szivoeolmyicnogarosuert,”oornlotwheersfirgenqauleunnctyilintfhoermleavteiolno,fetqhueiveanlteinret
such that the mean of the ancestral populations are arbi- chromosome is reached. A plot of Sj vs. j indicates which
trarily1 and21.Thisnormalizationstepmakesthemea- scales are important contributors to the process variance
sure robust to uneven sample sizes, which can affect the and indirectly provides information about the distribution
structure of the PCA (Novembre and Stephens 2008; ofadmixturetracts.Forexample,recentadmixtureproduces
Mcvean2009).Stabilityofthesignalsismaintainedbyspec- a peak in the wavelet variance at a large wavelet scale,
ifyingatoleranceE forseparatingtheancestralpopulations reflectinglongadmixturetracts,whilemoreancientadmix-
at a given SNP. This ensures that SNPs with poor discrimi- tureeventsproducepeaksatlowerwaveletscales,reflecting
nation are treated as uninformative in the next step of the shorter admixture tracts.
analysis. The wavelet variance for an individual i is given by
Wavelet transform Si ¼1 XT (cid:2)(cid:2)di (cid:2)(cid:2)2; (3)
The resulting SNP-level admixture signals indicate how j T j;k
k¼1
ancestry varies along the genome, but they invariably P
exhibitahighnoise-to-informationratio.Tointerpretthe where dij;k ¼ sYsicj;s2k are the wavelet coefficients for
signal, its frequency content can be described using the the signal Yi constructed using the wavelet system c. To
AdmixtureProcessesfromLocalAncestry 471
appreciate the methodology, it is sufficient to understand Indonesian Institute of Sciences. Blood samples or buccal
that the wavelet variance reflects the frequency content of swabs were collected from consenting, closely unrelated,
the signal, but more detailed background material is pro- and seemingly healthy individuals by J. Stephen Lansing
vided in the Supporting Information (Figure S5, Figure S6, (University of Arizona) and Herawati Sudoyo (Eijkman In-
Figure S7, Figure S8, Figure S9, and File S1). Our imple- stituteforMolecularBiology,Indonesia),withtheassistance
mentation employs Daubechies’ least asymmetric wavelet ofIndonesianPublicHealthclinicstaff.Allsamplecollection
number 8 in the waveslim (Whitcher 2013) package of the followed protocols for the protection of human subjects
statisticalsoftwareR(RDevelopmentCoreTeam2014).We establishedbyboththeEijkmanInstituteandtheUniversity
emphasize, however, that the methods proposed here are of Arizona institutional review boards. Participant inter-
robust to other choices of analyzing wavelet (see Supporting viewsconfirmedlocalresidenceforatleasttwogenerations
Information,FigureS5,FigureS6,FigureS7,FigureS8,Figure into the past. Samples were genotyped with the Affymetrix
S9,andFileS1). Axiom chip, yielding 548,994 SNPs across the autosomes.
Population averages are computed as (Sex-linked markers were excluded from the analysis.) The
X SNPdatawerecleanedusingstandardprotocolsinPLINKv.
1
SC ¼ Si 1.07 (Purcell et al. 2007; Purcell 2009) and the wavelet
j n j
C i2PC transform performed as described above.
TheapproximateBayesiancomputationanalysisemployed
(andsimilarlyforpopulationsP andP ).Anexampleofthe
A B 1000 data sets with sample sizes and SNP numbers set to
average wavelet variance for each population is shown in
those of the real data. These data sets were simulated by
Figure1B.Thewaveletvarianceishighestatfinescales,but
drawing from a uniform prior of admixture times between
astheancestralpopulationsalsoshowthispattern,itshould
10 and 300 generations. The admixture proportion for
beconsideredbackgroundnoise.Itisintuitivethatthevery
the Bena population (used as our primary test case) was
finest wavelet scales are uninformative because small num-
setto0.6,asestimatedpreviouslyfromtherealdata.The
bers of SNPs should beinsufficient to differentiate between
ABS metric was calculated for each simulation, and the
populations. The raw wavelet variance is therefore consid-
multiplechromosomestructureofthedataweremimicked
ered as a combination of informative variation and back-
bysamplingeachindividualrepeatedlywithdifferentdata
ground noise
densities.
SC ¼ICþN : (4)
j j j
Results
To extract the informative variance IC, wesubtract the pro-
J Asproofofconcept,wefirstappliedourwaveletmethodto
portion that can be attributed to noise. This is estimated
simulated data. A range of admixture scenarios was ex-
from the variation observed in the ancestral populations;
Nb ¼m(cid:2)maxðSA;SBÞ;wheremisamultiplicativefactorthat plored by varying parameters of the demographic model,
j j j particularly the time of admixture, admixture proportion,
allows the degree of thresholding to be controlled. Under
almost all conditions, a default value of m¼1 may be as- and single vs. multiple admixture events. Fifty simulations
were performed for each scenario with modest (but there-
sumed, and this threshold should be raised only if the ad-
forerealistic)ancestralsamplesizesofn =n =15andan
mixture signals exhibit high levels of noise (see Supporting A B
admixed sample size n = 20.
Information, Figure S5, Figure S6, Figure S7, Figure S8, C
FigureS9,andFileS1fordetails).Populationcharacteristics Admixturetime
that influence noise levels in the admixture signals are ex-
Becausetheabilityofwaveletmethodstocalculatethetime
plored in the next section. The final measure of the infor-
ofadmixtureiswellknownfromearlierwork(Pugachetal.
mative variance is given by
2011),weexploredthisfeaturefirst.Simulationswereper-
(cid:3) (cid:4)
b
IC ¼max SC2N;0 ; (5) formed for admixture times ranging from 10 to 320 gener-
J J j ations(i.e.,fromtherecentpastto(cid:3)10,000yearsago,using
which describes the frequency content that is unique to the agenerationintervalof30years;Fenner2005).Admixture
admixedpopulation(incontrasttotheancestralpopulations). at 10 generations shows the highest informative wavelet
varianceatscale13,reflectingrelativelyfew,longadmixture
Real genomicdata
blocks (Figure 2). As the time of admixture occurs further
To illustrate that our method performs well in real-world backinthepast,thepeakinwaveletvarianceshiftstoward
situations, itwasappliedtoaSNPgenotypingchipdataset successivelylowerwaveletscales,reflectingever-smallerad-
of 394 individuals from 16 communities spread across the mixture blocks driven by cumulative recombination along
Indonesianarchipelago(Table1).EquivalentSNPdatafrom the chromosome. The average frequency content can be
SouthernHanChineseandPapuaNewGuineaHighlanders characterized by the average block size metric ABS,termed
were used as proxies for the ancestral populations. Permis- the “wavelet center” by Pugach et al. (2011), which as
sion to conduct research in Indonesia was granted by the shown later, can be used to date the admixture event
472 J.Sandersonetal.
Table1 Summaryofcasestudypopulationsdescribingsamplesize areneartheaverageofallchromosome-levelrecombination
(n), proportionof Asian ancestryas inferredby PCA(p), andthe rates(FigureS2).Theeffectoflowersamplesize(condition
averageblocksizemetric(ABS,foradmixedpopulationsonly)
3) was investigated by reducing the number of individuals
Averageblock sampled from each population by 5, thus yielding sample
Population n p sizemetric(ABS) sizesthatwouldbesmallerthanalmostanypublishedpop-
SouthernHanChinese 13 1.00 – ulationgeneticsstudy(n ,n =10,n =15).Theeffectof
A B C
Nias 28 0.87 4.26 more recent divergence between the ancestral populations
Mentawai 29 0.87 4.30
(condition 4) was investigated by decreasing T from
Ancestral
Java 21 0.84 4.41 2000 to 1200 generations ago (50,000–30,000 years ago).
Sumatra 30 0.83 4.49
Bali 19 0.83 4.90 The effect of using a misrepresentative modern population
Sulawesi 21 0.80 6.57 as a proxy for an ancestral population (condition 5) was
Sumba,Wunga 30 0.67 7.90 investigated by studying ancestral populations with mixed
Sumba,Anakalang 30 0.66 7.83 (rather than “pure”) ancestry. Rather than using samples
Flores,Rampasasa 12 0.66 8.05
fromthetrueancestralpopulationP ,anadmixedancestral
Flores,Bena 30 0.57 8.14 A
Flores,Bama 30 0.55 8.34 population P*A was employed instead (p¼ 0.1). A wide
Timor,UmanenLawalu 17 0.55 8.44 range of parameters was applied for sensitivity testing, but
Timor,Kamanasa 19 0.53 8.42 forclarity,onlyresultsforsingleparametervaluesareshown
Lembata 28 0.53 8.39
on Figure 4. These examples are representative of all the
Pantar 27 0.45 8.47
tests that were run.
Alor 23 0.42 8.46
PapuaNewGuineaHighlands 13 0.00 – Variationinsummarymeasuresbetweensimulationswas
compared by computing the relative standard deviation
(RSD)ateachadmixturetime.Foralloftheerrorconditions
P above,thecomputedABSmetricsareconsistentwiththe
j(cid:2)IAB referencecase(condition1),butwithslightlylargerrelative
ABS¼ Pj j : (6)
IAB standarddeviations.Onlyforonecase(condition4;reduced
j j
divergencebetweentheancestralpopulationsandadmixture
at 320 generations) are the ABS metrics biased, with the
Admixture proportion
meanfallingoutsidetherangeofvaluesobservedforthe
Admixture proportions were varied between 0.5 (equal reference example. We emphasize that this is expected:
ancestryfromP andP )and0.025(ancestrypredominately admixtureshouldbemoredifficulttodetectwhenitoccurs
A B
fromP ).Forthisanalysis,thetimeofadmixturewasfixedat betweentwoancestralpopulationsthatdivergedonlyrecently.
A
160 generations. As the proportion of admixture decreases, StabilityoftheABSmetricsinthisparticularscenariocouldbe
the raw wavelet variance exhibits increasing levels of noise improvedbyapplyingahigherlevelofthresholding.However,
relative to informative variation. This is shown by the thedefaultvalueofm¼1wasretainedheretoprovideconsis-
reduced magnitude of the informative wavelet variance tencyacrossscenarios,todemonstratethedeteriorationinres-
(Figure3)andemphasizesthat,asexpected,itisincreasingly olution,andtoillustratethatthethresholdingparametercanbe
difficult to extract informative variation at low admixture ignoredforallbutthemostextremeadmixturecases.
proportions (small p) even where the signal is technically Inalloftheseexamples,includingthestandardreference
present.Inthisexample,informativeestimateswereobtained case, the localized admixture signals provide a noisy in-
for admixture proportions as low as 2.5%, although in dication of how ancestry varies along the chromosome.
general,therangeofpforwhichthismethodisapplicable Indeed, the inherent stochasticity of the block structure is
willalsodependonothercharacteristicsofthedata,such the primary reason why other sources of variance, such as
as the SNP density and sample size, as considered in the the cases discussed above, have relatively little additional
next section. effect on the overall results. This noise is addressed using
wavelets to capture the distribution of block sizes, coupled
Sensitivity analyses
with a correction based on the ancestral populations to
The sensitivity of the method to a wide range of data distinguish informative signals from background variation.
characteristicswasconsideredbyrepeatingtheresultsofthe Thecasesconsideredaboveallslightlyincreasenoiselevels
admixture time example with a large number of simulated relativetoinformativevariation,which,asdemonstratedby
data sets. Results are summarized in Table 2 and Figure 4. the admixture proportion example in Figure 3, reduces the
Condition 1 shows the original results, exactly as de- magnitudeoftheextractedinformativewaveletvariance.As
scribed above. New simulations were then performed to noiseincreases,itnaturallybecomesmoredifficulttoextract
mimicrealisticlinkagedisequilibrium(LD)(condition2).To informativevariation. However,this increaseinnoiselevels
do so as accurately as possible, we applied the real re- is minimal for all but the most extreme confounds, thus
combination rates observed along the first 100 Mb of allowingthetechniquetobeappliedrobustlytoaverywide
chromosome 1, as recombination rates for chromosome 1 range of scenarios.
AdmixtureProcessesfromLocalAncestry 473
Figure2 Informativewaveletvarianceforeachtimeofadmixture(10–320generationsusingdefaultthresholdingm¼1).Shadedbarsrepresentthe
averageover50simulationsateachadmixturetime;blackbarsrepresenttherangeacrossindividualsimulations.Theaverageblocksizemetricforeach
scenarioisindicatedbyadottedblueline.
The effect of SNP density (which is always a known only minimal overlap.This illustrates thatadwave offers in-
variable) is demonstrated by down sampling the data creasedpowertodifferentiatebetweenolderadmixturesce-
(conditions6–8).Theoriginaldensityof4306SNPs(condi- narios, with substantially reduced uncertainty in dating.
tion 6) was chosen to correspond to the size of our real We also emphasize that adwave requires far fewer user
chromosome 22 data set. Reducing the SNP density of this specifications with regard to runtime options. The only vari-
data set means that the resulting wavelet decomposition is ableforadwaveisthethresholdingparameter,andasshown
givenover11waveletscalesratherthantheearlier13,and above, the defaultvalue of m¼1 should beusedfor almost
so as expected, the computed mean ABS metrics are corre- all admixture scenarios. In contrast, the StepPCO results re-
spondinglymuchsmaller.However,thishasnoeffectonthe quiredasignallengthparameter(K¼1024),awindowsize
inference,asthedatasizeisalwaysknownandsimulations parameter(l¼5),andtwothresholdingparameters(thresh-
are simply run to match the size of the observed data. Fur- old = 0.1, maxlevel = 6) (all notations from Pugach et al.
ther reductions in SNP density to 3250 SNPs (condition 7) 2011).Adetaileddemonstrationofthismethodcomparison,
and 1625 SNPs (condition 8) are also shown. Note that with explanation of the settings chosen for StepPCO, is pro-
althoughtheabsolutevaluesoftheABSmetricsareshifted, vided in FigureS3.
the trend with admixture time remains consistent.
AdmixtureinIndonesianpopulations
Method comparison
Populations across Indonesia show genomic admixture
The original StepPCO method (Pugach et al. 2011) has al- between Asian and Melanesian ancestral sources (Cox
ready been tested extensively against other admixture de- et al. 2010), which has been dated using other methods to
tection methods, particularly HAPMIX (Price et al. 2009). an admixture event (cid:3)4000 years ago ((cid:3)130 generations)
We therefore focus here on comparing our improved wave- (Xuetal.2012).Wecalculatedwaveletsummarymeasures
let method against the StepPCO procedure. Figure 5 shows for16communitiesacrosstheIndonesianarchipelagousing
that the summary measure (wavelet center) used in 548,994autosomalSNPsscreenedin394individuals(Table
StepPCO is comparable to the adwave ABS metrics, as both 1). Equivalent data from Southern Han Chinese and Papua
exhibitastrongtrendwithtimeofadmixture.However,the New Guinea Highlanders was used as proxies for ancestral
dispersion is consistently smaller for the adwave ABS met- populations, as described in Cox et al. (2010).
rics. For example, the wavelet centers (StepPCO) computed The PCA for all individuals, where only the ancestral
for T ¼320 and T ¼160 show substantial overlap, while populationswereusedtodefinetheaxes,isshowninFigure
the ABS metrics (adwave) for the same populations show 6. Admixed individuals dispersed along the first principal
474 J.Sandersonetal.
The ABC inference procedure allows us to capture un-
certainty in admixture time estimates more robustly than
earlier wavelet dating approaches (Pugach et al. 2011; Xu
etal.2012).Toillustratethisprocess,datingwasperformed
on the Bena population of Flores in eastern Indonesia,
resulting in an estimated median admixture time of 147
generations (95% credible region: 122–178 generations),
or 4410 years before present (95% CR: 3660–5340 years
BP). This almost exactly matches earlier point estimates of
the admixture time (Xu et al. 2012) and is consistent with
ourcurrentunderstandingofIslandSoutheastAsianprehis-
tory (Bellwood 2007).
TherelationshipbetweentimeofadmixtureandtheABS
metricacrossallsimulationsisillustratedinFigure7A.ABC
Figure3 Relationship between proportion of admixture and informa-
was implemented using the R package abc (Csilléry et al.
tivewaveletvariance.Forthisexampleonly,anondefaultvalueforthe
threshold m¼1:1was used to account for increased noise in the ad- 2012), and the posterior distribution of admixture time
mixturesignalsduetolowproportionsofadmixture,asdescribedinthe was computed using a local linear regression (Beaumont
text.Themagnitudeofthewaveletvariancedecreaseswiththeadmix- et al. 2002) with a tolerance rate of 0.2. Cross validation
tureproportion,shownascoloredbarsfromblack(P=0.50)toyellow
wasusedtoevaluatetheaccuracyofthisestimate:thepre-
(P=0.025).
diction error was low (0.038) and insensitive to the exact
tolerance value. For future research focusing on parameter
component illustrate the primary genomic signal, a strong inference, this procedure could be modified to use a larger
gradient in Asian-Melanesian ancestry that has previously number of simulated data sets and a lower tolerance rate.
been observed across the region (Cox et al. 2010). The in- However, this simple example clearly illustrates that the
formative wavelet variance was computed separately for adwavemethodhasgoodstatisticalpowertodateadmixture
each chromosome and individual and subsequently com- using a relatively small number of simulations.
binedtoprovideasinglemeasureforeachpopulation(Fig-
Multipleadmixtureevents
ure S4). To combine information across chromosomes,
which vary considerably in size, the raw admixture signals Another aim of this work is to show that our improved
were windowed: all signals were reduced to the size of the waveletapproachcanbeusedtostudyotherfeaturesofthe
smallest chromosome (importantly without discarding any admixture process beyond the well-explored question of
data) by computing averages over a window of SNPs admixture time. In the examples covered thus far, it has
(details of the windowing procedure are provided in Sup- been assumed that admixture occurred as a single event.
portingInformation,FigureS5,FigureS6,FigureS7,Figure However, additional waves of admixture will result in the
S8, Figure S9, and File S1). The SNP density and window introduction of new ancestry tracts, replacing a proportion
sizeforeachchromosomeareshowninTableS1.Thiswin- of older, shorter ancestry blocks with newer, longer ones.
dowingprocedureisusedonlytostandardizechromosomes Pugachetal.(2011)brieflyconsideredtheeffectsofcontin-
tothesamelengthandutilizesveryshortwindowsofSNPs uous admixture within a wavelet setting, showing that this
(unlike the approach of Pugach et al. 2011). leads to underestimated admixture times in their original
The average block size metrics calculated for each pop- methodologicalframework.Incontrast,weinsteadconsider
ulation are shown in Table 1. The first six Indonesian popu- scenarioswithtwodistinctadmixtureevents.Weshowthat
lations (Nias, Mentawai, Java, Sumatra, Bali, and Sulawesi) this process creates distinctive patterns in the observed in-
exhibit predominantly Asian ancestry, with high-frequency formative variation, which can be used to reconstruct more
noiseinthesignalscausingsomebiasinthecomputedABS complex demographic processes (as opposed to being trea-
metrics (Figure S4). The remaining Indonesian populations ted solely as a potential source of bias).
exhibit less extreme Asian ancestry proportions (42–67%), Inthefollowingdual-admixturescenarios,thefirstadmix-
with the resulting ABS metrics appearing broadly similar ture event always occurs at 160 generations. To investigate
between populations. theeffectofseparationbetweenadmixtureevents,thesecond
Undertheassumptionofasingleadmixturetime(relaxed admixtureeventvariesbetween10and80generations.Inthe
in later sections), the average blocksizemetriccanbe used extremecaseofadmixtureat160and10generationsago,the
to date the time of admixture using approximate Bayesian localizedadmixturesignalscontaintwodominantfrequencies.
computation (ABC). A general introduction to ABC can be Single admixture events at 160 and 10 generations lead to
found in Csilléry et al. (2010) and Sunnåker et al. (2013), peaks in the informative wavelet variance at wavelet scales
while ABC in the context of parameter estimation forpopu- of 9 and 13, respectively. When two admixture events occur,
lationadmixturehasbeenconsideredbySousaetal.(2009) the informative wavelet variance is instead spread between
andRobinsonet al. (2014). thesescales(Figure8A).Astheadmixtureeventsoccurcloser
AdmixtureProcessesfromLocalAncestry 475
Table2 Sensitivityoftheadwavemethodtoarangeofdatalimitations
Datalimitations Admixturetime(generations)
Condition Description 10 20 40 80 160 320
1 Reference 11.55(0.69) 11.14(0.58) 10.65(0.7) 10.13(0.66) 9.54(1.12) 8.84(1.72)
2 RealisticLD 11.71(0.84) 11.31(0.88) 10.85(1.08) 10.32(1.04) 9.82(1.23) 9.09(2.28)
3 Reducedsamplesize 11.66(0.90) 11.25(0.63) 10.74(0.88) 10.20(0.85) 9.58(1.19) 8.85(1.78)
4 Reduceddivergencebetween 11.64(0.83) 11.24(0.69) 10.76(0.82) 10.21(1.02) 9.49(1.35) 8.32(3.59)
ancestralpopulations
5 Non-representativeancestral 11.73(1.34) 11.23(1.42) 10.77(1.51) 10.21(1.45) 9.58(1.94) 8.60(3.06)
populations
6 SNPdensityT=4036 10.47(0.91) 10.04(0.73) 9.55(0.92) 8.98(1.22) 8.37(1.62) 7.60(3.13)
7 SNPdensityT=3250 9.96(0.76) 9.60(0.60) 9.19(0.76) 8.68(1.08) 8.12(1.41) 7.37(2.6)
8 SNPdensityT=1625 8.99(1.23) 8.63(1.09) 8.20(1.50) 7.68(2.29) 7.10(3.15) 6.33(5.86)
Meanaverageblocksizevalues(relativestandarddeviationinparentheses)areshownforeachadmixturetime.ReferencedataweresimulatedwithT=13,000SNPs,
populationssizesofnA,nB=15,nC=15,anddivergencebetweentheancestralpopulationsatTAncestral=2000generationsago.
together,thisspreadintheobservedinformativewaveletvar- usingtheinformativewaveletvariance,while theABSmet-
iancedecreases (Figures8,B–D). ric correctly classifies only 60% of cases. For real multiple
For one admixture event, a single dominant peak is admixture situations, this classification framework could be
observed in the informative wavelet variance, and the ABS extended to a more complex inferential setting (such as
metric therefore provides a convenient summary measure. ABC), but this simple example demonstrates the potential
Formultipleadmixtureevents,theABSmetricdescribesthe forreconstructingcomplexadmixturescenariosfromthefull
averageadmixturetime, butprovides noinformationabout wavelet variance profile.
thedurationoverwhichadmixtureoccurred.Incontrast,the
informative wavelet variance should provide additional in-
Discussion
formationaboutthepeakdispersion.Toexplorethepotential
for identifying more complex admixture scenarios, a simple Wavelettechniquesprovideinformationontheancestryblock
classification rule was implemented. An admixed population structure of admixed genomes and hence can be used to
P isassignedtooneoftwogroupsG ; G ,whicharechar- reconstructtheprocessesinvolvedinpastadmixtureevents.
C 1 2
acterizedbythesummarymeasuresM ; M .Thisschemeis Ancestryblocksarestrictlyunobservableandcanbeinferred
1 2
described with abstract choice of summary measure, but be- only from the data. Wavelets provide indirect information
low, we consider how different summary measures (taking on the block structure, thus providing an alternative over
M ; M to be either the ABS metric or wavelet informative methods that assign ancestry directly (Sankararaman et al.
1 2
variance)affect thesuccess ofclassification. 2008; Price et al. 2009). A growing body of methods now
The classification rule is implemented as follows: assignancestryindirectlyusingvariousunrelatedapproaches
(Moorjani et al. 2011; Baran et al. 2012; Gravel 2012; Loh
1. The“true” summarymeasuresM , M arecomputedfor
1 2 et al. 2013; Mapleset al. 2013;BrownandPasaniuc 2014),
each group using values obtained from the first 25
but here we extend the use of wavelet techniques as intro-
simulations.
duced by Pugach et al. (2011). Importantly, our implemen-
2. Foreachoftheremaining25trialdatasets(s=1,...,25),
tation differs markedly from the original StepPCO program,
c
estimated summary statistics M are calculated. The di-
s withthemaindifferencesateachstageoftheanalysishigh-
vergence measures are defined as
lighted below:
XS (cid:2)(cid:2) (cid:2)(cid:2) (cid:129) Localized admixture signal formation: StepPCO (Pugach
D ¼ (cid:2)Mc2M(cid:2); (7)
i s i et al. 2011) uses large windows of SNPs to produce an
s¼1
averaged admixture signal in localized windows along
the genome. Our work demonstrates that wavelet meth-
for i¼1;2.
odsareequallyapplicabletotherawunwindowedsignals,
with the windowing procedure performed intrinsically as
3. If D ,D , classify to G ; otherwise classify to G .
1 2 1 2 part of the wavelet analysis, and therefore not requiring
TheclassificationratesareshowninTable3forscenario arbitrary a prioridecisions on window size.
1 (a single admixture event at 60 generations; mean ABS (cid:129) Wavelet analysis: The wavelet methods we describe are
10.47, range 10.30–10.64) and scenario 2 (two admixture basedontheMODWT,whichoffersmoreflexibilityinits
events at 160 and 10 generations; mean ABS 10.57, range applicationsincethereisnorestrictiononthelengthofthe
10.35–10.90). With a sample size of just 10 individuals for signals.Conversely,StepPCOemploysthediscretewavelet
the admixed population, perfect classification is achieved transform (DWT), which has the strict requirement that
476 J.Sandersonetal.
Figure4 Sensitivitytoarangeofrealisticdatalimi-
tations. Comparison to reference data (condition 1)
simulated with T = 13,000 SNPs, populations sizes
n ;n ¼15,n ¼15,andancestralpopulationdiver-
A B C
genceatT ¼2000generations.Thegrayarea
Ancestral
showstherangeofABSmetricsobservedunderthe
standardreferencecondition.(A)Potentialsourcesof
error(conditions2–5);(B)varyingSNPdensities(con-
ditions6–8).Notethatthedeclineinabsolutevalues
oftheABSmetricsinBisexpected;theseareeasily
accountedforinaninferencesettingbecausetheSNP
densityisalwaysaknownvariable.Conditiondescrip-
tionsandnumericvaluesarepresentedinTable2.
signalsbeoflength2n.Datamustthereforebewindowed, The work presented here also makes several other
ordiscarded,tomeettherestrictivelengthrequirementsof advances. Theaverageblocksizemetric has previouslybeen
the DWT framework. Another advantage of the MODWT shown to capture the time of admixture. Here, we have
is that the resulting wavelet coefficients are translation implemented a more formal dating procedure using ABC
equivariant, meaning that circularly shifting the data undertheassumptionofasingleadmixtureevent.Inreality,
results in the same shifting of the coefficients. Said differ- populationsmayhaveexperiencedmultipleadmixtureevents
ently, changing the starting point—for instance, to avoid leading to complex patterns of genetic variation. We have
a poor quality SNP—does not affect the resulting wavelet shown that the wavelet variance contains additional in-
coefficients,whereasthisisnottrueundertheDWTframe- formationtoidentifythesemorecomplexadmixturescenar-
work. This property is particularly important if the results ios.Thishighlightsthepotentialofwavelet-basedtechniques
are to be used for specific localized genomic regions (as tobecoupledwithformalstatisticalinferenceproceduresto
discussedbrieflybelow)andthusprovidesasolidstatistical robustly distinguish between the range of scenarios that
foundation forfuturework. could have resultedin any observedgenetic pattern.
(cid:129) Extraction of relevant information: The portion of the Method performance for the StepPCO procedure has al-
resulting wavelet decomposition that is informative ready been tested against other admixture detection meth-
about the admixture process is extracted in a simple ods,mostextensivelywithHAPMIX(Priceetal.2009),with
procedure with reference to the ancestral populations, favorableresults.Thisisespeciallytrueforolderadmixture
offeringgreatersimplicityandobjectivitythanthemul- events (Pugach et al. 2011). While an in-depth comparison
tistage thresholding procedure described by Pugach with other local ancestry detection methods would be of
et al. (2011). great interest (Moorjani et al. 2011; Baran et al. 2012;
(cid:129) Software: The adwave software, which implements the Gravel 2012; Loh et al. 2013; Maples et al. 2013; Brown
method described in this article, is an official package in and Pasaniuc 2014), such an analysis is beyond the scope
the R project (http://cran.r-project.org/web/packages/ of this manuscript. We have therefore focused instead on
adwave/index.html).Thisallowsextremelyeasyinstalla- showing how adwave markedly improves on the original
tion and use, as well as providing a series of simple waveletmethodimplementedinStepPCO.Asshownabove,
worked examples as a learning exercise. The adwave adwaveoffersimprovedstatisticalpowertodifferentiatebe-
package is also faster than the existing StepPCO code tweenadmixturescenarios,offersmuchreduceduncertainty
and offers more flexibility in the choice of analyzing inmodelparameterestimates,andimportantly,isfareasier
wavelet (unlike StepPCO, which employs only the sim- tousethan StepPCO, especially byrequiringfarfewer user-
plest “square-shaped” Haar wavelet). specified runtime parameters.
Figure5 ComparingStepPCOandadwaveshowingthe
relationship between wavelet transform summaries and
timeofadmixture.(A)Adwaveusingm¼1;(B)StepPCO
usingK¼1024,l¼5,threshold=0.1,andmaxlevel=
6.Numbersindicatetherelativestandarddeviation(RSD,
%) foreach admixturetime. Note thedifferenceindis-
crimination power between the two methods for older
admixture events (95% confidence intervals as dashed
blueandgreenhorizontallines).
AdmixtureProcessesfromLocalAncestry 477
Figure6 PCAofautosomalSNPdatafromIndone-
sian populations, with Southern Han Chinese (blue
circles) and Papua New Guinea Highlanders (green
circles) employed as proxy ancestral populations.
Numbersgivecalculatedadmixtureproportions.
In the future, considering the full wavelet periodogram, surrogate populations become more different from the true
rather than the genome-wide summary measures used in ancestralpopulations,thelocalizedadmixturesignalsbecome
bothadwaveandStepPCO,mayyieldpromisingresultswher- increasingly noisy. Although this ultimately leads to a loss
ever the distribution of ancestry tracts along the genome is ofidentifiabilityinextremecases,themethodisremarkably
substantially nonstationary. Bryc et al. (2010) use their for- robust to moderate deviations from these assumptions.
mulation of localized admixture signals to address whether As shown above for low admixture proportions, through
regions of the genome show predominant ancestry from judicial choice of the thresholding parameter even ex-
agivenpopulation.Waveletsarewellsuitedtodistinguishing tremely noisy data can still provide meaningful estimates
local features in data and could be helpful in this regard, (the only situation in which we encourage deviation from
identifying features that may not be easily detected by con- the default setting).
sideringthe localized admixturesignalsintheir raw form. Samplesize(bothintermsofSNPsandindividuals)isalso
Other prospective areas for further work include the important and affects the PCA step of the procedure. The
extension of these methods to the more general case of purposeofthePCAstepistosummarizetheoverallvariability
multipopulationadmixture.MaandAmos(2012)describethe among individuals, which includes both between-population
use of PCA as a diagnostic in this setting, and PCA has been and within-population variability. In reconstructing popula-
usedtoassignmultipopulationancestryinthesoftwarePCAd- tion ancestry, we aim to describe between-population vari-
mix(Brisbinetal.2012).Thewaveletmethodsdescribedhere ation, while ignoring within-population variation. This is
could be extended in a similar way by considering pairwise achieved by selecting the first principal component, as long
combinations ofany numberofancestralpopulations. asthesamplesizesaresufficientlylarge.Within-population
In contrast, key restrictions that determine our ability to fluctuationsofindividualcoordinatesonthePCAscatterplot
reconstruct admixture events include the degree of differ- can be caused by subtle population substructure. Assuming
entiation between the ancestral populations and the repre- that no such substructure is present, these fluctuations
sentativenessofsamplesusedassurrogateancestralgroups. decrease as the total sample size increases, and an asymp-
As the ancestral populations become more similar or the toticallystablepatternoftheeigenvectorplotresults(Maand
Figure 7 Dating time of admixture for Bena (Flores,
eastern Indonesia) using approximate Bayesian com-
putation. (A) Relationship between admixture time
and average block size metric for all simulations; (B)
weightedposteriordistributionofadmixturetime.Me-
dian estimated time of admixture, indicated by the
blue line, is 147 generations (95% credible region:
122–178generations).
478 J.Sandersonetal.