Table of Contents
The goal of AdmixPoly is to provide functions to perform global and local admixture inference from bi- and multi-allelic marker dosages (discrete or continuous) in polyploid species.
You can install AdmixPoly directly from the CRAN (not available at the moment):
install.packages(pkg='AdmixPoly',repos='https://cran.r-project.org/')
or from the source repository (GitLab) using remotes:
install.packages(pkg='remotes', repos='https://cran.r-project.org/')
remotes::install_git('https://gitlab.cirad.fr/agap/seg/admixpoly.git')
When installing from the GitLab repository, first make sure to have a compilation environment installed on your system.
Download and install the appropriate 64-bits version of Rtools and git.
Debian and derivative (e.g. Ubuntu, Mint):
apt-get install g++ git libicu-dev
Redhat and derivative (e.g. Fedora, Centos, RockyLinux):
dnf install gcc-c++ git libicu-devel
Check if git and the clang compiler are installed with these commands in a terminal window:
clang --version
git --version
If not, run this command:
xcode-select --install
This is a basic example illustrating how to perform admixture inference using AdmixPoly.
First, an admixed dataset is simulated using the SimulatePop function,
considering N=100 individuals of ploidy P=8, M=1000 markers with
L=10 alleles distributed on C=5 chromosomes, and K=4 ancestral
groups.
DataSim <- AdmixPoly::SimulatePop(K=3L, N=100L, P=8L, M=1000L, C=2L, L=10L)
The resulting DataSim object includes:
1) Geno: list of genotyping matrices over markers
head(DataSim$Geno[[1]])
#> Allele1 Allele2 Allele3 Allele4 Allele5 Allele6 Allele7 Allele8 Allele9
#> Ind1 0 1 0 0 1 1 1 0 1
#> Ind2 0 2 3 0 1 0 1 0 0
#> Ind3 0 5 2 0 0 0 1 0 0
#> Ind4 0 1 2 2 0 0 1 2 0
#> Ind5 0 3 1 1 0 0 2 0 0
#> Ind6 0 1 1 1 1 2 1 1 0
#> Allele10
#> Ind1 3
#> Ind2 1
#> Ind3 0
#> Ind4 0
#> Ind5 1
#> Ind6 0
2) Prop: matrix of admixture proportions
head(DataSim$Prop)
#> K1 K2 K3
#> Ind1 0.106375 0.431750 0.461875
#> Ind2 0.587625 0.394000 0.018375
#> Ind3 0.680500 0.217500 0.102000
#> Ind4 0.509500 0.490500 0.000000
#> Ind5 0.859500 0.033125 0.107375
#> Ind6 0.181125 0.267750 0.551125
3) Freq: list of matrices of allele frequencies in ancestral groups
DataSim$Freq[[1]]
#> Allele1 Allele2 Allele3 Allele4 Allele5 Allele6 Allele7
#> K1 0.02999937 0.3216599409 0.2827303 0.06252135 0.04569980 0.01259175 0.1560145
#> K2 0.06062462 0.0007481235 0.3483231 0.01478613 0.04326888 0.04442491 0.1591765
#> K3 0.00337067 0.0736675715 0.3034438 0.06786697 0.17255898 0.11904414 0.1678521
#> Allele8 Allele9 Allele10
#> K1 0.01215938 0.019408693 0.05721497
#> K2 0.19739411 0.001450772 0.12980285
#> K3 0.03792193 0.036377344 0.01789647
4) GeneticMap: genetic map dataframe
head(DataSim$GeneticMap)
#> Marker Chromosome Distance
#> 1 Marker1 Chr1 0.0000000
#> 2 Marker2 Chr1 0.2004008
#> 3 Marker3 Chr1 0.4008016
#> 4 Marker4 Chr1 0.6012024
#> 5 Marker5 Chr1 0.8016032
#> 6 Marker6 Chr1 1.0020040
From simulated data, admixture proportions can be estimated using the
AdmixGlobal function:
ResGlobalAdmix <- AdmixPoly::AdmixGlobal(Geno=DataSim$Geno, K=3L,MaxIter=1000,Verbose=F)
The estimated admixture proportions are available from:
head(ResGlobalAdmix$Prop)
#> K1 K2 K3
#> Ind1 0.42751014 0.09196272 0.48052714
#> Ind2 0.39607077 0.60392823 0.00000100
#> Ind3 0.18982802 0.71856785 0.09160413
#> Ind4 0.47932675 0.52067225 0.00000100
#> Ind5 0.01352721 0.92827717 0.05819562
#> Ind6 0.24174073 0.19640609 0.56185317
They can be represented using the
ggplot2-based GlobalPlot
function:
AdmixPoly::GlobalPlot(ResGlobalAdmix$Prop)

Estimated proportions are very close to simulated ones (up to an arbitrary labeling of groups), as indicated by the root mean square error (RMSE):
group_corresp <- apply(cor(DataSim$Prop,ResGlobalAdmix$Prop),2,which.max)
sqrt(mean((DataSim$Prop-ResGlobalAdmix$Prop[,group_corresp])^2))
#> [1] 0.02168584
Estimated allele frequencies are available from:
ResGlobalAdmix$Freq[[1]]
#> Allele1 Allele2 Allele3 Allele4 Allele5 Allele6 Allele7
#> K1 0.06565551 0.00000100 0.3477740 0.00000100 0.03057859 0.06962224 0.1208821
#> K2 0.03109042 0.29182439 0.2713478 0.06719518 0.07038353 0.01340646 0.1872236
#> K3 0.00000100 0.05668937 0.2049396 0.07104907 0.13114363 0.14552308 0.2449462
#> Allele8 Allele9 Allele10
#> K1 0.19902928 0.000001000 0.16645533
#> K2 0.00000100 0.009649829 0.05787779
#> K3 0.06641198 0.034075409 0.04522066
The convergence can be checked from the log-likelihood:
plot(ResGlobalAdmix$LogLik, xlab="Iterations", ylab="LogLik", type="o")

Let us consider a second admixed population from the same ancestral groups with a higher admixture level by specifying a higher ‘AlphaProp’ value.
DataSim2 <- AdmixPoly::SimulatePop(K=3L, N=100L, P=8L, M=1000L, C=2L, L=10L,
AlphaProp = 10, Freq = DataSim$Freq)
The higher admixture level can be illustrated using a barplot of simulated amdixture proportions.
AdmixPoly::GlobalPlot(DataSim2$Prop)

In this case, it can be beneficial to initialize the algorithm with ancestral allele frequencies estimated using the first population ‘FreqInit=ResGlobalAdmix\$Freq’ and only update admixture proportions ‘ParamToUpdate=“Prop”’.
ResGlobalAdmix2 <- AdmixPoly::AdmixGlobal(Geno=DataSim2$Geno, K=3L, MaxIter=1000, Verbose=F,
FreqInit=ResGlobalAdmix$Freq, ParamToUpdate="Prop")
Again, estimated proportions are very close to simulated ones, as indicated by the RMSE:
sqrt(mean((DataSim2$Prop-ResGlobalAdmix2$Prop[,group_corresp])^2))
#> [1] 0.02526292
Based on global admixture results, local admixture can be estimated for
a given individual using the AdmixLocal function:
Ind <- "Ind3"
ResLocalAdmix <- AdmixPoly::AdmixLocal(Geno=DataSim$Geno,ResAdmixGlobal = ResGlobalAdmix,
GeneticMap = DataSim$GeneticMap,Ind = Ind,P = 8L,
SmoothParam = 1,Verbose=F)
Local ancestry dosages based on posterior probabilities are available from:
head(ResLocalAdmix$Posterior$Chr1)
#> K1 K2 K3
#> Marker1 0.02838614 6.966770 1.004844
#> Marker2 0.02683719 6.966129 1.007034
#> Marker3 0.02464581 6.965626 1.009728
#> Marker4 0.02337629 6.966975 1.009649
#> Marker5 0.02250143 6.968826 1.008673
#> Marker6 0.02199903 6.970135 1.007866
Alternatively, local ancestry dosages based on the Viterbi algorithm are available from:
head(ResLocalAdmix$Viterbi$Chr1)
#> K1 K2 K3
#> Marker1 0 7 1
#> Marker2 0 7 1
#> Marker3 0 7 1
#> Marker4 0 7 1
#> Marker5 0 7 1
#> Marker6 0 7 1
The RMSE of estimated local ancestry proportions (ancestry dosages divided by the ploidy) is:
sqrt(mean((DataSim$Ancestry[[Ind]]$Chr1/8-ResLocalAdmix$Posterior$Chr1[,group_corresp]/8)^2))
#> [1] 0.01208612
sqrt(mean((DataSim$Ancestry[[Ind]]$Chr1/8-ResLocalAdmix$Viterbi$Chr1[,group_corresp]/8)^2))
#> [1] 0.01118034
Local ancestry dosages can be represented using the
ggplot2-based LocalPlot
function:
AdmixPoly::LocalPlot(Ancestry = ResLocalAdmix$Posterior, GeneticMap = DataSim$GeneticMap)

Local admixture inference over the five first individuals can be run using:
ResLocalAdmix_list <- lapply(paste0("Ind",1:5),function(i){
res_i <- AdmixPoly::AdmixLocal(Geno=DataSim$Geno,ResAdmixGlobal = ResGlobalAdmix,
GeneticMap = DataSim$GeneticMap,Ind = i,P = 8L,
SmoothParam = 1,Verbose=F)
return(res_i)
})
When computing time is large, each individual can be run in parallel on a high-performance computing cluster.
Genotyping data formatted as a variant call format (VCF) can be imported as:
vcf_path <- system.file("extdata", "Test.vcf", package = "AdmixPoly")
DataVCF <- AdmixPoly::ReadVCF(File = vcf_path)
The resulting DataVCF object includes:
1) MarkerInfo: first eight columns of the VCF
head(DataVCF$MarkerInfo)
#> # A tibble: 6 × 9
#> MARKER CHROM POS ID REF ALT QUAL FILTER INFO
#> <chr> <chr> <int> <chr> <chr> <chr> <chr> <chr> <chr>
#> 1 Chr1_15892318 Chr1 15892318 . A T . . .
#> 2 Chr1_31209423 Chr1 31209423 . T C . . .
#> 3 Chr1_39962954 Chr1 39962954 . C G . . .
#> 4 Chr1_40515463 Chr1 40515463 . G A . . .
#> 5 Chr1_44904209 Chr1 44904209 . A T . . .
#> 6 Chr1_47542399 Chr1 47542399 . T C . . .
2) Geno: list of genotyping matrices
head(DataVCF$Geno[[1]])
#> A T
#> IND1 4 4
#> IND2 5 3
#> IND3 1 7
#> IND4 4 4
#> IND5 0 8
#> IND6 4 4
3) GeneticMap: genetic map dataframe in which physical between first
and last marker of each chromosome are scaled to 100 cM
head(DataVCF$GeneticMap)
#> # A tibble: 6 × 3
#> Chromosome Marker Distance
#> <chr> <chr> <dbl>
#> 1 Chr1 Chr1_15892318 0
#> 2 Chr1 Chr1_31209423 20.4
#> 3 Chr1 Chr1_39962954 32.1
#> 4 Chr1 Chr1_40515463 32.9
#> 5 Chr1 Chr1_44904209 38.7
#> 6 Chr1 Chr1_47542399 42.2
Instead of estimated dosages, read depth ratios can be imported by specifying the allele read depths field of the FORMAT column:
DataVCF2 <- AdmixPoly::ReadVCF(File = vcf_path,AlleleDepthField = "AD")
head(DataVCF2$Geno[[1]])
#> A T
#> IND1 0.41935484 0.5806452
#> IND2 0.50000000 0.5000000
#> IND3 0.08108108 0.9189189
#> IND4 0.68181818 0.3181818
#> IND5 0.03125000 0.9687500
#> IND6 0.53125000 0.4687500
Please note that read depth ratios are not scaled to the ploidy level, which must then be informed in the inference functions.
Another function can be used to read the haplotype presence-absence (HPA) format obtained from HaploCharmer:
hpa_path <- system.file("extdata", "Test.hpa", package = "AdmixPoly")
DataHPA <- AdmixPoly::ReadHPA(File = hpa_path)
head(DataHPA$Geno[[1]])
#> hap1 hap2 hap3
#> IND1 0.54545455 0.1818182 0.2727273
#> IND2 0.48888889 0.2444444 0.2666667
#> IND3 0.08823529 0.6764706 0.2352941
#> IND4 0.08108108 0.4054054 0.5135135
#> IND5 0.00000000 0.3684211 0.6315789
#> IND6 0.45454545 0.1818182 0.3636364
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.