tutorial.R
In admixr: An Interface for Running 'ADMIXTOOLS' Analyses

## ----setup, include = FALSE---------------------------------------------------
evaluate <- .Platform$OS.type == "unix" && system("which qpDstat", ignore.stdout = TRUE) == 0

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "",
  eval = evaluate
)

## ---- eval = FALSE------------------------------------------------------------
#  install.packages("devtools")
#  devtools::install_github("bodkan/admixr")

## ---- eval = FALSE------------------------------------------------------------
#  install.packages("tidyverse")

## ----libraries, message = FALSE, warning = FALSE------------------------------
library(admixr)
library(tidyverse)

## ----eigenstrat_path----------------------------------------------------------
(prefix <- download_data(dirname = tempdir()))

## ----eigenstrat_trio----------------------------------------------------------
list.files(path = dirname(prefix), pattern = basename(prefix), full.names = TRUE)

## ----ind_file, echo = FALSE---------------------------------------------------
cat(system(paste0("column -t ", prefix, ".ind"), intern = TRUE), sep = "\n")

## ----snp_file, echo = FALSE---------------------------------------------------
cat(system(paste0("head -n 3 ", prefix, ".snp"), intern = TRUE), sep = "\n")

## ----geno_file, echo = FALSE--------------------------------------------------
cat(system(paste0("head -n 3 ", prefix, ".geno"), intern = TRUE), sep = "\n")

## -----------------------------------------------------------------------------
snps <- eigenstrat(prefix)

## ---- comment = "#>"----------------------------------------------------------
snps

## ----pop_def1-----------------------------------------------------------------
pops <- c("French", "Sardinian", "Han", "Papuan", "Khomani_San", "Mbuti", "Dinka")

## ----d------------------------------------------------------------------------
result <- d(W = pops, X = "Yoruba", Y = "Vindija", Z = "Chimp", data = snps)

## ----eval = FALSE-------------------------------------------------------------
#  head(result)

## ----d_kable, echo = FALSE----------------------------------------------------
knitr::kable(head(result))

## ----d_plot, fig.width = 7, fig.height = 4------------------------------------
ggplot(result, aes(fct_reorder(W, D), D, color = abs(Zscore) > 2)) +
  geom_point() +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_errorbar(aes(ymin = D - 2 * stderr, ymax = D + 2 * stderr))

## ----f4-----------------------------------------------------------------------
result <- f4(W = pops, X = "Yoruba", Y = "Vindija", Z = "Chimp", data = snps)

## ----eval = FALSE-------------------------------------------------------------
#  head(result)

## ----f4_kable, echo = FALSE---------------------------------------------------
knitr::kable(head(result))

## ----f4ratio------------------------------------------------------------------
result <- f4ratio(X = pops, A = "Altai", B = "Vindija", C = "Yoruba", O = "Chimp", data = snps)

## ---- eval=FALSE--------------------------------------------------------------
#  head(result)

## ----f4ratio_kable, echo=FALSE------------------------------------------------
knitr::kable(head(result))

## ----f4ratio_plot, fig.width = 7, fig.height = 4------------------------------
ggplot(result, aes(fct_reorder(X, alpha), alpha, color = abs(Zscore) > 2)) +
  geom_point() +
  geom_errorbar(aes(ymin = alpha - 2 * stderr, ymax = alpha + 2 * stderr)) +
  geom_hline(yintercept = 0, linetype = 2) +
  labs(y = "Neandertal ancestry proportion", x = "present-day individual")

## ----pops2--------------------------------------------------------------------
pops <- c("French", "Sardinian", "Han", "Papuan", "Mbuti", "Dinka", "Yoruba")

result <- f3(A = pops, B = pops, C = "Khomani_San", data = snps)

## ---- eval=FALSE--------------------------------------------------------------
#  head(result)

## ----f3_kable, echo=FALSE-----------------------------------------------------
knitr::kable(head(result))

## ----f3_plot, fig.width = 8, fig.height = 6-----------------------------------
# sort the population labels according to an increasing f3 value relative to French
ordered <- filter(result, A == "Mbuti", B != "Mbuti") %>% arrange(f3) %>% .[["B"]] %>% c("Mbuti")

# plot heatmap of pairwise f3 values
result %>%
  filter(A != B) %>%
  mutate(A = factor(A, levels = ordered),
         B = factor(B, levels = ordered)) %>%
  ggplot(aes(A, B)) + geom_tile(aes(fill = f3))

## -----------------------------------------------------------------------------
result <- qpWave(
 left = c("French", "Sardinian", "Han"),
 right = c("Altai", "Yoruba", "Mbuti"),
 data = snps
)

## ---- eval = FALSE------------------------------------------------------------
#  result

## ---- echo = FALSE------------------------------------------------------------
knitr::kable(result)

## -----------------------------------------------------------------------------
result <- qpWave(
 left = c("Papuan", "French", "Sardinian", "Han"),
 right = c("Altai", "Yoruba", "Mbuti"),
 data = snps
)

## ---- eval = FALSE------------------------------------------------------------
#  result

## ---- echo = FALSE------------------------------------------------------------
knitr::kable(result)

## ----qpAdm--------------------------------------------------------------------
result <- qpAdm(
  target = c("Sardinian", "Han", "French"),
  sources = c("Vindija", "Yoruba"),
  outgroups = c("Chimp", "Denisova", "Altai"),
  data = snps
)

## ---- eval = FALSE------------------------------------------------------------
#  result$ranks

## ---- echo = FALSE------------------------------------------------------------
knitr::kable(result$ranks)

## ---- eval = FALSE------------------------------------------------------------
#  result$proportions

## ---- echo=FALSE--------------------------------------------------------------
knitr::kable(result$proportions)

## ---- eval = FALSE------------------------------------------------------------
#  result$subsets

## ---- echo=FALSE--------------------------------------------------------------
knitr::kable(result$subsets)

## ----orig_ind, echo = FALSE, comment = ""-------------------------------------
cat(system(paste0("column -t ", snps$ind), intern = TRUE), sep = "\n")

## ----relabel------------------------------------------------------------------
# paths to the original ind file and a new modified ind file, which will
# contain merged population labels
modif_snps <- relabel(
  snps,
  European = c("French", "Sardinian"),
  African = c("Dinka", "Yoruba", "Mbuti", "Khomani_San"),
  Archaic = c("Vindija", "Altai", "Denisova")
)

## ---- comment = "#>"----------------------------------------------------------
modif_snps

## ----modif_ind, echo = FALSE, comment = ""------------------------------------
cat(system(paste0("column -t ", modif_snps$group), intern = TRUE), sep = "\n")

## ----modif_d------------------------------------------------------------------
result <- d(W = "European", X = "African", Y = "Archaic", Z = "Chimp", data = modif_snps)

## ----eval = FALSE-------------------------------------------------------------
#  head(result)

## ----modif_d_kable, echo = FALSE----------------------------------------------
knitr::kable(head(result))

## ----present_snps, results = "hide"-------------------------------------------
count_snps(snps)

## ----present_snps_kable, echo = FALSE-----------------------------------------
knitr::kable(count_snps(snps))

## -----------------------------------------------------------------------------
bed <- file.path(dirname(prefix), "regions.bed")

## -----------------------------------------------------------------------------
# BED file contains regions to keep in an analysis
new_snps <- filter_bed(snps, bed)

# BED file contains regions to remove from an analysis
new_snps <- filter_bed(snps, bed, remove = TRUE)

## ---- comment = "#>"----------------------------------------------------------
new_snps

## ---- eval = FALSE------------------------------------------------------------
#  snps %>%
#    filter_bed("regions.bed") %>%
#    d(W = "French", X = "Mbuti", Y = "Vindija", Z = "Chimp")

## ---- eval = FALSE------------------------------------------------------------
#  new_snps <- transversions_only(snps)
#  
#  # perform the calculation only on transversions
#  d(W = "French", X = "Dinka", Y = "Altai", Z = "Chimp", data = new_snps)

## ---- eval = FALSE------------------------------------------------------------
#  snps %>%                                    # take the original data
#    filter_bed("regions.bed", remove = TRUE) %>%  # remove sites not in specified regions
#    transversions_only() %>%                      # remove potential false SNPs due to aDNA damage
#    d(W = "French", X = "Dinka", Y = "Altai", Z = "Chimp") # calculate D on the filtered dataset

## ----merge_eigenstrat, eval = FALSE-------------------------------------------
#  # this is just an example code - it will not run unless you specify the paths
#  merged <- merge_eigenstrat(
#      merged = <"prefix of the merged dataset">
#      a = first_EIGENSTRAT_object,
#      b = second_EIGENSTRAT_object
#  )

## ----d_log, comment = "#>"----------------------------------------------------
dres <- d(W = c("French", "Han", "Dinka"), X = "Yoruba", Y = "Vindija", Z = "Chimp", data = snps)

dres

## ----d_loginfo, comment = "#>"------------------------------------------------
loginfo(dres)

## ----qpAdm_log, comment = "#>"------------------------------------------------
qpadm_res <- qpAdm(
  target = c("Sardinian", "Han"),
  sources = c("Vindija", "Yoruba"),
  outgroups = c("Chimp", "Denisova", "Altai"),
  data = snps
)

qpadm_res

## ----qpadm_log_target, comment = "#>"-----------------------------------------
loginfo(qpadm_res, target = "Han")

## ----qpadm_log_save-----------------------------------------------------------
loginfo(qpadm_res, target = "Sardinian", save = TRUE, prefix = "qpAdm_Neandertal_ancestry")