Admixtools qpWave: Who do you form a clade with?

[Jovialis]

This is basically how you model yourself with a single population.

# Full F2 extraction + qpWave clade test script

# Load necessary libraries
library(admixtools)
library(tidyverse)
library(crayon)  # for colored console output

# ------------------------------------------------------------------
# 1) Define your data directories
# ------------------------------------------------------------------
prefix <- "C:/Users/Jovialis/Documents/Bioinformatics/Jovialis_HO_merge_PLINK/merged_HO"
f2_dir <- "C:/Users/Jovialis/Documents/Bioinformatics/Jovialis_HO_merge_PLINK/f2_blocks"

# ------------------------------------------------------------------
# 2) Specify populations
# ------------------------------------------------------------------
target    <- "Jovialis"
left      <- c(target, "Italian_South.HO")
outgroups <- c(
  "Ethiopia_4500BP.AG",     "Russia_UstIshim_IUP.DG",
  "Italy_Epigravettian.AG.BY.AA", "Russia_YuzhniyOleniyOstrov_Mesolithic.AG",
  "Russia_MA1_UP.SG",        "Georgia_Satsurblia_LateUP.SG",
  "Jordan_PPNB.AG"
)
all_pops <- c(left, outgroups)

# ------------------------------------------------------------------
# 3) Extract F2 blocks with relaxed filters
#    - auto_only = TRUE for autosomes only
#    - maxmiss   = 1    to allow all SNPs regardless of missingness
#    - overwrite = TRUE to refresh existing files
# ------------------------------------------------------------------
extract_f2(
  prefix,
  f2_dir,
  pops      = all_pops,
  overwrite = TRUE,
  auto_only = TRUE,
  maxmiss   = 1,
  blgsize   = 0.05
)

# ------------------------------------------------------------------
# 4) Load the precomputed F2 blocks
# ------------------------------------------------------------------
f2_blocks <- f2_from_precomp(
  f2_dir,
  pops   = all_pops,
  afprod = TRUE
)

# ------------------------------------------------------------------
# 5) Run qpWave for rank = 0 (clade test)
# ------------------------------------------------------------------
qpwave_res <- qpwave(
  f2_blocks,
  left    = left,
  right   = outgroups,
  verbose = TRUE
)

# ------------------------------------------------------------------
# 6) Display the rank=0 summary
# ------------------------------------------------------------------
print(qpwave_res$rankdrop)

# ------------------------------------------------------------------
# 7) Print populations tested for record-keeping in blue text
# ------------------------------------------------------------------
cat(blue("\n*** Populations tested ***\n"))
cat(blue(paste0("Target: ", target, "\n")))
cat(blue(paste0("Left group(s): ", paste(left, collapse = ", "), "\n")))
cat(blue(paste0("Outgroups: ", paste(outgroups, collapse = ", "), "\n")))

 

[baeticvs]

you can do this using one way models in qpadm, no? how it differs?

 

[Jovialis]

you can do this using one way models in qpadm, no? how it differs?

In most aDNA studies, qpWave is used to test if a target and a single reference form the same clade, and qpAdm is only used when fitting mixtures of two or more sources.

 

[Jovialis]

I form a very tight fit for a clade with C6_Central_Italy_Medieval_Central_Mediterranean.

 

[Jalisciense]

I made it in the Classic ADMIXTOOLS, the process is actually different from the ADMIXTOOLS 2 but I was able to figured out how to do it and run it in the classical one.

This is a test one tho but it was able to be statistically robust (p-value/taildiff ≥0.05):

qpWave: parameter file: parqpWave
### THE INPUT PARAMETERS
##PARAMETER NAME: VALUE
S1: Jalisciense_merged_with_v62HO
indivname: Jalisciense_merged_with_v62HO.ind
snpname: Jalisciense_merged_with_v62HO.snp
genotypename: Jalisciense_merged_with_v62HO.geno
popleft: leftwavepops.txt
popright: rightwavepops.txt
details: YES
inbreed: NO
allsnps: YES
##  qpWave -p parqpWave
## qpWave version: 1570
inbreed set NO
inbreed set NO
tmp: /tmp/fsx.4679

left pops:
Mexico_Colonial_European.AG
          Spanish.DG

right pops:
            Mbuti.DG
  Ethiopia_4500BP.SG
Morocco_Iberomaurusian.AG
Serbia_EBA_Yamnaya.AG
Italy_Epigravettian.AG.BY.AA
  Israel_Natufian.AG
Iran_TepeAbdulHosein_N.SG
    Russia_MA1_UP.SG

f4info:
f4rank: 0 dof:      7 chisq:    12.830 tail:         0.0763716673 dofdiff:      0 chisqdiff:     0.000 taildiff:                    1

f4info:
f4rank: 1 dof:      0 chisq:     0.000 tail:                    1 dofdiff:      7 chisqdiff:    12.830 taildiff:         0.0763716673
B:
          scale   170.255
Ethiopia_4500BP.SG     0.162
Morocco_Iberomaurusian.AG     0.063
Serbia_EBA_Yamnaya.AG     2.507
Italy_Epigravettian.AG.BY.AA     0.727
Israel_Natufian.AG    -0.236
Iran_TepeAbdulHosein_N.SG    -0.064
Russia_MA1_UP.SG     0.311
A:
          scale     1.000
     Spanish.DG     1.000

removing /tmp/fsx.4679
##end of qpWave:        0.006 seconds cpu        0.131 Mbytes in use

So in summary: Mexico_Colonial_European.AG and Spanish.DG form a clade.

 

[Jalisciense]

In most aDNA studies, qpWave is used to test if a target and a single reference form the same clade, and qpAdm is only used when fitting mixtures of two or more sources.

I didn't know that, thanks for let us know bro, from now on I'll use qpWave to model pops in 1-way.

The 2-way, 3-way....and so on will be modeled just in qpAdm.

 

[baeticvs]

In most aDNA studies, qpWave is used to test if a target and a single reference form the same clade, and qpAdm is only used when fitting mixtures of two or more sources.

i have seen qpadm 1 way models in studies

 

[Jovialis]

You're right,

I've composed a single-left pop qpAdm, that produces the same exact result as the qpwave:

qpwave result:

One‐way F₂‐block qpAdm script:

# One‐way F₂‐block qpAdm
#
# Matches the extraction/filtering settings of the full qpWave script
# (auto_only = TRUE, maxmiss = 1, overwrite = TRUE, blgsize = 0.05)

# 1) Load libraries
library(admixtools)
library(tidyverse)
library(crayon)   # for colored console output

# 2) File paths
prefix <- "C:/Users/jovialis/Documents/Bioinformatics/Jovialis_HO_merge_PLINK/merged_HO"
f2_dir <- "C:/Users/jovialis/Documents/Bioinformatics/Jovialis_HO_merge_PLINK/f2_blocks"

# 3) Define populations
target     <- "Jovialis"
one_source <- "C6_Central_Italy_Medieval_Central_Mediterranean"
outgroups  <- c(
  "Ethiopia_4500BP.AG",
  "Russia_UstIshim_IUP.DG",
  "Italy_Epigravettian.AG.BY.AA",
  "Russia_YuzhniyOleniyOstrov_Mesolithic.AG",
  "Russia_MA1_UP.SG",
  "Georgia_Satsurblia_LateUP.SG",
  "Jordan_PPNB.AG"
)
all_pops <- c(target, one_source, outgroups)

# 4) Extract F₂ blocks with the same filters as the full script
extract_f2(
  prefix,
  f2_dir,
  pops      = all_pops,
  overwrite = TRUE,
  auto_only = TRUE,
  maxmiss   = 1,
  blgsize   = 0.05
)

# 5) Load the precomputed F₂ blocks
f2_blocks <- f2_from_precomp(
  f2_dir,
  pops   = all_pops,
  afprod = TRUE
)

# 6) Print populations tested
cat(blue("\n*** Populations tested for one‐way qpAdm ***\n"))
cat(blue(paste0("Target:      ", target, "\n")))
cat(blue(paste0("Left source: ", one_source, "\n")))
cat(blue(paste0("Outgroups:   ", paste(outgroups, collapse = ", "), "\n\n")))

# 7) Run one‐way, F₂‐block–based qpAdm
res_f2_oneway <- qpadm(
  f2_blocks,
  left    = one_source,
  right   = outgroups,
  target  = target,
  verbose = TRUE
)

# 8) Print results
cat("\n— One-way F₂-block qpAdm Results —\n")
cat("Admixture weight (≈1 if clade holds):\n")
print(res_f2_oneway$weights)
cat("\nClade-test p-value (popdrop):\n")
print(res_f2_oneway$popdrop)

 

[Jovialis]

I'm assuming a single-qpAdm model would be most appropate if you desire to employ the ALLSNPS=TRUE flag, which is not possible with qpWave, as it requires pre-computed F2 stats

ALLSNPS=TRUE Single-Left pop qpAdm:

prefix <- "C:/Users/jovialis/Documents/Bioinformatics/Jovialis_HO_merge_PLINK/merged_HO"

library(admixtools)
library(tidyverse)

target <- c('Jovialis')
left <- c('C6_Central_Italy_Medieval_Central_Mediterranean')
right <- c("Ethiopia_4500BP.AG", "Russia_UstIshim_IUP.DG", "Italy_Epigravettian.AG.BY.AA", "Russia_YuzhniyOleniyOstrov_Mesolithic.AG", "Russia_MA1_UP.SG", "Georgia_Satsurblia_LateUP.SG", "Jordan_PPNB.AG")

results <- qpadm(prefix, left, right, target, allsnps = TRUE)

print(results$weights)
print(results$popdrop)

 

[Jovialis]

Essentially both approaches complement one another to determine I indeed am not-rejected as being modeled from a single source (C6 Medieval Central Italy). The difference is the examination of SNPs. ALLSNPS=TRUE maximizing the availability; while Pre-comp F2 curates the SNPs to a higher standard reducing the possibility of bias.

 

[Jalisciense]

I form a very tight fit for a clade with C6_Central_Italy_Medieval_Central_Mediterranean.

View attachment 18493

If I put any other basal population other than Mbuti, I can't run qpAdm, so how did you do to put Ethiopia_4500BP, bro?

 

[thunorr]

I form a clade with British Islanders, and other Northwestern Europeans.

This is basically how you model yourself with a single population.


View attachment 18456

View attachment 18457

View attachment 18458

# Full F2 extraction + qpWave clade test script

# Load necessary libraries
library(admixtools)
library(tidyverse)
library(crayon)  # for colored console output

# ------------------------------------------------------------------
# 1) Define your data directories
# ------------------------------------------------------------------
prefix <- "C:/Users/Jovialis/Documents/Bioinformatics/Jovialis_HO_merge_PLINK/merged_HO"
f2_dir <- "C:/Users/Jovialis/Documents/Bioinformatics/Jovialis_HO_merge_PLINK/f2_blocks"

# ------------------------------------------------------------------
# 2) Specify populations
# ------------------------------------------------------------------
target    <- "Jovialis"
left      <- c(target, "Italian_South.HO")
outgroups <- c(
  "Ethiopia_4500BP.AG",     "Russia_UstIshim_IUP.DG",
  "Italy_Epigravettian.AG.BY.AA", "Russia_YuzhniyOleniyOstrov_Mesolithic.AG",
  "Russia_MA1_UP.SG",        "Georgia_Satsurblia_LateUP.SG",
  "Jordan_PPNB.AG"
)
all_pops <- c(left, outgroups)

# ------------------------------------------------------------------
# 3) Extract F2 blocks with relaxed filters
#    - auto_only = TRUE for autosomes only
#    - maxmiss   = 1    to allow all SNPs regardless of missingness
#    - overwrite = TRUE to refresh existing files
# ------------------------------------------------------------------
extract_f2(
  prefix,
  f2_dir,
  pops      = all_pops,
  overwrite = TRUE,
  auto_only = TRUE,
  maxmiss   = 1,
  blgsize   = 0.05
)

# ------------------------------------------------------------------
# 4) Load the precomputed F2 blocks
# ------------------------------------------------------------------
f2_blocks <- f2_from_precomp(
  f2_dir,
  pops   = all_pops,
  afprod = TRUE
)

# ------------------------------------------------------------------
# 5) Run qpWave for rank = 0 (clade test)
# ------------------------------------------------------------------
qpwave_res <- qpwave(
  f2_blocks,
  left    = left,
  right   = outgroups,
  verbose = TRUE
)

# ------------------------------------------------------------------
# 6) Display the rank=0 summary
# ------------------------------------------------------------------
print(qpwave_res$rankdrop)

# ------------------------------------------------------------------
# 7) Print populations tested for record-keeping in blue text
# ------------------------------------------------------------------
cat(blue("\n*** Populations tested ***\n"))
cat(blue(paste0("Target: ", target, "\n")))
cat(blue(paste0("Left group(s): ", paste(left, collapse = ", "), "\n")))
cat(blue(paste0("Outgroups: ", paste(outgroups, collapse = ", "), "\n")))

 

[thunorr]

I form a clade with British Islanders, and other Northwestern Europeans.View attachment 18590View attachment 18591View attachment 18592View attachment 18593
View attachment 18594
View attachment 18595
View attachment 18596

However, the problem with the British Isles is that they're all genetically identical when it comes to Neolithic results, which this model is evidently set up to prepare. Iron Age results can vary a whole lot inbetween various British Isles (from now on abbreviated as BI) groups.