Admixtools 2-way Italy_BA/Crete_EMBA-cline

[Jovialis]

The outgroups and left groups work for all Italians and create a cline between Italy_BA.SG and Greece_Crete_HgCharalambos_EMBA.AG. Although this model is mostly historically plausible only for southern regions (particularly Apulia), I think it hints at something broader.

library(admixtools)
library(tidyverse)

prefix   <- "C:/Users/jovialis/Documents/Bioinformatics/Jovialis_HO_merge/merged_HO"
f2_dir   <- "C:/Users/jovialis/Documents/Bioinformatics/Jovialis_HO_merge/f2_blocks"

target   <- "Italian_South.HO"
left     <- c(
    "Italy_BA.SG",
    "Greece_Crete_HgCharalambos_EMBA.AG"
)

outgroups <- c(
    "Ethiopia_4500BP.AG", "Russia_UstIshim_IUP.DG", "http://Italy_Epigravettian.AG.BY.AA", "Russia_YuzhniyOleniyOstrov_Mesolithic.AG", "Russia_MA1_UP.SG", "Georgia_Satsurblia_LateUP.SG", "Israel_Natufian.AG")

mypops <- c(target, left, outgroups)

extract_f2(
    prefix, f2_dir,
    pops      = mypops,
    overwrite = TRUE,
    auto_only = TRUE,
    blgsize   = 0.05
)

f2_blocks <- f2_from_precomp(
    f2_dir,
    pops   = mypops,
    afprod = TRUE
)

qpwave_results <- qpwave(
    f2_blocks,
    left    = left,
    right   = outgroups,
    verbose = TRUE
)
print(qpwave_results$rankdrop)

min_p <- min(qpwave_results$rankdrop$p, na.rm = TRUE)

if (!http://is.na(min_p) && min_p > 0.05) {
    res_geno <- qpadm(
        prefix,
        left      = left,
        right     = outgroups,
        target    = target,
        allsnps   = TRUE,
        auto_only = TRUE,
        verbose   = TRUE,
        return_f4 = TRUE
    )
    print(res_geno$weights)
    print(res_geno$popdrop)
} else {
    res_geno <- qpadm(
        prefix,
        left      = left,
        right     = outgroups,
        target    = target,
        allsnps   = TRUE,
        auto_only = TRUE,
        verbose   = TRUE,
        return_f4 = TRUE
    )
    res_f2 <- qpadm(
        f2_blocks,
        left      = left,
        right     = outgroups,
        target    = target,
        verbose   = TRUE
    )
    print(res_geno$weights)
    print(res_geno$popdrop)
    print(res_f2$weights)
    print(res_f2$popdrop)
}