Merge pull request #67 from EvolEcolGroup/run_admixture

Run admixture

.Rbuildignore

@@ -5,3 +5,4 @@
 ^\.github$
 pkgdown/
 _pkgdown.yml
+admixture/

.github/workflows/R-CMD-check.yaml

@@ -2,7 +2,7 @@
 # Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help
 on:
   push:
-    branches: [main]
+    branches: [main, run_admixture]
   pull_request:
     branches: [main]
 
@@ -46,6 +46,17 @@ jobs:
           extra-packages: any::rcmdcheck
           needs: check
 
+      - name: Get ADMIXTURE for linux
+        if: matrix.config.os == 'ubuntu-latest'
+        run: |
+             mkdir admixture
+             cd admixture
+             wget https://dalexander.github.io/admixture/binaries/admixture_linux-1.3.0.tar.gz
+             tar -xvzf admixture_linux-1.3.0.tar.gz
+             cd ..
+             echo "$PWD/admixture/dist/admixture_linux-1.3.0" >> "$GITHUB_PATH"
+
+
       - uses: r-lib/actions/check-r-package@v2
         with:
           upload-snapshots: true

.github/workflows/test-coverage.yaml

@@ -29,6 +29,16 @@ jobs:
           extra-packages: any::covr, any::xml2
           needs: coverage
 
+      - name: Get ADMIXTURE for linux
+        run: |
+             mkdir admixture
+             cd admixture
+             wget https://dalexander.github.io/admixture/binaries/admixture_linux-1.3.0.tar.gz
+             tar -xvzf admixture_linux-1.3.0.tar.gz
+             cd ..
+             echo "$PWD/admixture/dist/admixture_linux-1.3.0" >> "$GITHUB_PATH"
+
+
       - name: Test coverage
         run: |
           cov <- covr::package_coverage(

DESCRIPTION

@@ -50,6 +50,7 @@ Suggests:
     LEA,
     rmarkdown,
     readr,
+    reticulate,
     testthat (>= 3.0.0),
     vcfR
 Remotes:

NAMESPACE

@@ -8,13 +8,15 @@ S3method(augment,gt_dapc)
 S3method(augment,gt_pca)
 S3method(augment,q_matrix)
 S3method(augment_loci,gt_pca)
+S3method(autoplot,gt_admix)
 S3method(autoplot,gt_cluster_pca)
 S3method(autoplot,gt_dapc)
 S3method(autoplot,gt_pca)
 S3method(autoplot,gt_pcadapt)
 S3method(autoplot,q_matrix)
 S3method(autoplot,qc_report_indiv)
 S3method(autoplot,qc_report_loci)
+S3method(c,gt_admix)
 S3method(count_loci,tbl_df)
 S3method(count_loci,vctrs_bigSNP)
 S3method(dplyr_col_modify,grouped_gen_tbl)
@@ -74,7 +76,7 @@ S3method(show_loci,tbl_df)
 S3method(show_loci,vctrs_bigSNP)
 S3method(show_ploidy,tbl_df)
 S3method(show_ploidy,vctrs_bigSNP)
-S3method(summary,q_matrix_list)
+S3method(summary,gt_admix)
 S3method(summary,rbind_report)
 S3method(summary,vctrs_bigSNP)
 S3method(tbl_sum,gen_tbl)
@@ -92,6 +94,7 @@ export(filter_high_relatedness)
 export(gen_tibble)
 export(get_p_matrix)
 export(get_q_matrix)
+export(gt_admixture)
 export(gt_as_genind)
 export(gt_as_genlight)
 export(gt_as_geno_lea)
@@ -145,6 +148,7 @@ export(q_matrix)
 export(qc_report_indiv)
 export(qc_report_loci)
 export(rbind_dry_run)
+export(read_q_files)
 export(scale_fill_distruct)
 export(select_loci)
 export(select_loci_if)

R/autoplot_gt_admix.R

@@ -0,0 +1,48 @@
+#' Autoplots for `gt_admix` objects
+#'
+#' For `gt_admix`, the following types of plots are available:
+#' - `cv`: the cross-validation error for each value of `k`
+#' - `barplot` a standard barplot of the admixture proportions
+#'
+#' `autoplot` produces simple plots to quickly inspect an object. They are
+#' not customisable; we recommend that you use `ggplot2` to produce publication
+#' ready plots.
+#'
+#' @param object an object of class `gt_admixture`
+#' @param type the type of plot (one of "cv", and "boxplot")
+#' @param k the value of `k` to plot (for `barplot` type only)
+#' param repeat the repeat to plot (for `barplot` type only)
+#' @param run the run to plot (for `barplot` type only)
+#' @param ... not used.
+#' @returns a `ggplot2` object
+#' @name autoplot_gt_admix
+#' @export
+autoplot.gt_admix <- function(object,
+                            type=c("cv", "barplot"),
+                            k = NULL,
+                            run = NULL,
+                            ...)
+{
+  type <- match.arg(type)
+  if (type== "cv") {
+    if (is.null(object$cv)){
+      stop("No cross validation error available")
+    }
+    ggplot2::ggplot(data.frame(k=object$k, cv=object$cv), ggplot2::aes(x=.data$k, y=.data$cv)) +
+      ggplot2::geom_point() +
+      ggplot2::geom_line() +
+      ggplot2::labs(x="k", y="Cross validation error")
+  } else if (type == "barplot") {
+    # check that k is specified
+    if (is.null(k)){
+      stop("You must specify a value for k")
+    }
+    # check that run is specified
+    if (is.null(run)){
+      stop("You must specify a value for repeat")
+    }
+    # get the Q matrix for the specified k and repeat
+    Q <- get_q_matrix(object, k = k, run = run)
+    autoplot(Q)
+  }
+}

R/autoplot_gt_pcadapt.R

@@ -1,9 +1,9 @@
 #' Autoplots for `gt_pcadapt` objects
 #'
 #' For `gt_pcadapt`, the following types of plots are available:
-#' - `qq`: a qunatile-quantile plot of the p-values from pcadapt
+#' - `qq`: a qunatile-quantile plot of the p-values from `pcadapt`
 #' (wrapping [bigsnpr::snp_qq()])
-#' - `manhattan` a manhattan plot of the p-values from pcadapt
+#' - `manhattan` a manhattan plot of the p-values from `pcadapt`
 #' (wrapping [bigsnpr::snp_manhattan()])
 #'
 #' `autoplot` produces simple plots to quickly inspect an object. They are

R/gt_admix_methods.R

@@ -0,0 +1,80 @@
+#' Combine method for gt_admix objects
+#'
+#' @param ... A list of `gt_admix` objects
+#' @return A `gt_admix` object with the combined data
+#' @export
+
+c.gt_admix <- function(...) {
+  # check that all the objects are of class gt_admix
+  if (!all(sapply(list(...), function(x) inherits(x, "gt_admix")))) {
+    stop("All the objects must be of class gt_admix")
+  }
+
+  combined_obj <- list()
+  # combine all the elements from each list
+  combined_obj$k <- sapply(list(...), function(x) x$k)
+  combined_obj$Q <- sapply(list(...), function(x) x$Q)
+  # if we have a P element in any of the objects, combine it
+  if (all(sapply(list(...), function(x) !is.null(x$P)))) {
+    combined_obj$P <- sapply(list(...), function(x) x$P)
+  }
+  # if we have a log_lik element in any of the objects, combine it
+  if (all(sapply(list(...), function(x) !is.null(x$loglik)))) {
+    combined_obj$loglik <- unlist(sapply(list(...), function(x) x$loglik))
+  }
+
+  # if we have a log element in any of the objects, combine it
+  if (all(sapply(list(...), function(x) !is.null(x$log)))) {
+    combined_obj$log <- sapply(list(...), function(x) x$log)
+  }
+  # if we have a cv element in any of the objects, combine it
+  if (all(sapply(list(...), function(x) !is.null(x$cv)))) {
+    combined_obj$cv <- unlist(sapply(list(...), function(x) x$cv))
+  }
+  # if the first object has an id element, use it in the combined object
+  if (!is.null(list(...)[[1]]$id)) {
+    combined_obj$id <- list(...)[[1]]$id
+  }
+  # if the first object has a group element, use it in the combined object
+  if (!is.null(list(...)[[1]]$group)) {
+    combined_obj$group <- list(...)[[1]]$group
+  }
+  # set the class of the object
+  class(combined_obj) <- c("gt_admix", "list")
+  return(combined_obj)
+}
+
+#' Summary method for gt_admix objects
+#'
+#' @param object a `gt_admix` object
+#' @param ... unused (necessary for compatibility with generic function)
+#' @return A summary of the `gt_admix` object
+#' @export
+summary.gt_admix <- function(object, ...) {
+  cat("Admixture results")
+  # if we only have one element, give the k
+  if (length(object$k) == 1) {
+    cat(" for k = ", object$k, "\n")
+  } else {
+    tab_sum <- table(object$k)
+    tab_sum <- rbind(as.numeric(names(tab_sum)),tab_sum)
+    rownames(tab_sum) <- c("k", "n")
+    colnames(tab_sum) <- rep("", ncol(tab_sum))
+    cat(" for multiple runs:")
+    print(tab_sum)
+  }
+  cat("with slots:\n")
+  cat("$Q for Q matrices\n")
+  # if there is a lot P in the object, print it
+  if ("P" %in% names(object)){
+    cat("$P for  matrices\n")
+  }
+  # if there is a slot log in the object, print it
+  if ("log" %in% names(object)){
+    cat("$log for logs from the algorithm\n")
+  }
+  # if there is a slot cv in the object, print it
+  if ("cv" %in% names(object)){
+    cat("$cv for cross validation error\n")
+  }
+}