#' Interfaces with the COLONY pedigree assignment program.
#'
#' Interfaces with the command-line version of the COLONY pedigree program.
#' Requires a snpRdata object containing offspring genotypes and can optionally
#' take snpRdata objects containing maternal or paternal genotypes, or both. No
#' facet support.
#'
#' Requires that the COLONY program is installed locally if running within snpR.
#' Text files exported from write_colony_input command can also be imported to
#' command line versions of COLONY (eg. on a compute cluster). Input and output
#' files will be stored in a colony folder created in the current working
#' directory. The functions documented here can write input files, call them
#' using COLONY, and parse some parts of the results given the original
#' snpRdata object. Note that no facet support is currently available here due
#' to the complexity and number of possible input parameters that are difficult
#' to handle across multiple facets and facet levels. Facet support for the
#' basic, default operation with few additional options may be added in the
#' future. For now, facets can be handled during parentage and pedigree creation
#' in snpR using the \code{\link{run_sequoia}} function, which runs a notably
#' simpler (with respect to implementation) pedigree toolkit.
#'
#' These functions include many commonly used options but not all possible
#' parameters for COLONY inputs. See COLONY User Guide for using extra features.
#'
#' This is still in development. The defaults and the most commonly used options
#' have been tested and work well, but some of the more esoteric options haven't
#' been fully tested yet.
#' @param x snpRdata object containing offspring genotypes.
#' @param outfile character, default "colony_input". Output file name. A file
#' path may be provided (e.g. "colony/colony_run_1.txt").
#' @param method character, default "FPLS". Pedigree reconstruction method. For
#' more details see the Colony User Guide. Options: \itemize{\item{"FPLS":
#' }{Pure pairwise likelihood method, combines the full likelihood and
#' pairwise likelihood methods. A good compromise between speed and
#' accuracy.}{\item{"FL":}{Full Likelihood. More accurate than PLS but more
#' computationally intensive and slow to run, especially with large complex
#' datasets.}}\item{"PLS": }{Pairwise likelihood score. Less accurate but less
#' computationally intensive than FL.}}
#' @param run_length numeric in c(0,1,2,3), default 2. Length of run:
#' short/medium/long/very long.
#' @param sampleIDs character, default NULL. Name of a column in the sample
#' metadata that designates sample identifications/"names". Each name must be
#' unique!
#' @param sibship_prior numeric in c(0, 1, 2, 3, 4), default 0. Strength the
#' sibship size prior (no prior, weak, medium, strong, OR determine from known
#' prior values). Values other than 0 require additional parameters. Option
#' for supplying value of 4, currently not implemented in snpR. See Colony
#' User Guide for more details.
#' @param paternal_sib_size numeric, default NULL. Minimum value is 0. The
#' number of offspring in the candidate pool that are known to share the same
#' father. If this value is not zero, then you must include a file with the
#' known paternal sibship/paternity.
#' @param maternal_sib_size numeric, default NULL. Minimum value is 0. The
#' number of offspring in the candidate pool that are known to share the same
#' mother. If this value is not zero, then you must include a file with the
#' known paternal sibship/maternity.
#' @param nruns integer, default 1. A number of replicate runs for the dataset.
#' @param seed integer, default NULL. Supply a four digit integer (eg: 1234,
#' 9876) as a starting point for the algorithm.
#' @param maternal_genotypes snpRdata object containing maternal genotypes.
#' @param paternal_genotypes snpRdata object containing paternal genotypes.
#' @param maternal_inclusion_prob numeric in 0:1, default 0. Probability the
#' mother is in the dataset ranging from 0 to 1.
#' @param paternal_inclusion_prob numeric in 0:1, default 0. Probability the
#' father is in the dataset ranging from 0 to 1.
#' @param update_af character, default TRUE. Should Colony update the allele
#' frequencies used in the calculations?
#' @param dioecious character, default TRUE. Is this species diploid/dioecious?
#' FALSE = haploid/monoecious. Colony does not work with more ploidy.
#' @param inbreeding character, default TRUE. Should Colony assume inbreeding in
#' the calculations?
#' @param male_monogamous character, default FALSE. Should Colony assume males
#' are monogamous?
#' @param female_monogamous character, default FALSE. Should Colony assume
#' females are monogamous?
#' @param clone_inference character, default FALSE. Should Colony infer clones
#' in the sample set?
#' @param sibship_scaling character, default TRUE. Should Colony scale sibling
#' groups?
#' @param known_af character, default FALSE. If TRUE snpR will calculate and
#' supply mafs, else, supply a numeric vector containing the known maf for
#' each locus.
#' @param precision integer in c(0,1,2,3), default 2. Low/Medium/High/Very High
#' for calculating the maximum likelihood.
#' @param dropout numeric vector where each value is in 0:1, default 0.01.
#' Supply a flatrate value for all markers, or a vector corresponding to the
#' allelic droput rate for each marker.
#' @param genotyping_error numeric vector where each value is in 0:1, default
#' 0.01. Supply a flatrate value for all markers, or a vector corresponding to
#' the genotyping error rate for each marker.
#' @param known_maternal_dyads character, default NULL. Supply matrix or
#' dataframe with known maternal-offspring dyads. Offspring ID in column 1,
#' Maternal ID in column 2.
#' @param known_paternal_dyads character, default NULL. Supply matrix or
#' dataframe with known paternal-offspring dyads. Offspring ID in column 1,
#' Paternal ID in column 2.
#' @param known_maternal_max_mismatches integer in c(0,1,2:nsample), default 0.
#' @param known_paternal_max_mismatches integer in c(0,1,2:nsample), default 0.
#' @param known_maternal_sibships character, default NULL. Data frame or matrix
#' with sibship size followed by single column containing all of the sibling
#' IDs separated by spaces.
#' @param known_paternal_sibships character, default NULL. Data frame or matrix
#' with sibship size followed by single column containing all of the sibling
#' IDs separated by spaces.
#' @param maternal_exclusions character, default NULL. Data.frame or matrix with
#' column 1 the offspring ID, column 2 the number of excluded females, column
#' 3 the IDs of excluded females separated by spaces.
#' @param paternal_exclusions character, default NULL. Data.frame or matrix with
#' column 1 the offspring ID, column 2 the number of excluded males, column 3
#' the IDs of excluded males separated by spaces.
#' @param excluded_maternal_siblings character, default NULL. Data.frame or
#' matrix with column 1 the offspring ID, column 2 the number of excluded
#' siblings, column 3 the IDs of excluded siblings separated by spaces.
#' @param excluded_paternal_siblings character, default NULL. Data.frame or
#' matrix with column 1 the offspring ID, column 2 the number of excluded
#' siblings, column 3 the IDs of excluded siblings separated by spaces.
#' @param infile character. Path to the pre-written colony input file to be run.
#' @param colony_path character. Path to the colony executable.
#' @param path character. Path to the directory containing colony results.
#' @param x snpRdata object from which metadata for colony results can be found.
#' @param prefix character. The prefix for the colony files to be parsed.
#' @param update_bib character or FALSE, default FALSE. If a file path to an
#' existing .bib library or to a valid path for a new one, will update or
#' create a .bib file including any new citations for methods used. Useful
#' given that this function does not return a snpRdata object, so a
#' \code{\link{citations}} cannot be used to fetch references.
#' @param cleanup logical, default FALSE. If TRUE, colony files will be removed
#' after parsing.
#' @param verbose Logical, default TRUE. If TRUE, colony progress will be
#' reported to console.
#'
#' @author William Hemstrom
#' @author Melissa Jones
#'
#' @aliases write_colony_input call_colony parse_colony run_colony
#'
#' @name colony_interface
#'
#' @references Jones,O.R. and Wang,J. (2010) COLONY: a program for parentage and
#' sibship inference from multilocus genotype data. Mol. Ecol. Resour., 10,
#' 551–555.
#'
#' @examples
#'
#' # A simple example for running all individuals in the snpR object as siblings
#' # in colony. Not run to avoid clutter.
#' \dontrun{
#' write_colony_input(x = stickSNPs, outfile = "stk.col")
#' }
#'
#' # A more complex example requires 1) creating and adding variables to the
#' # stickSNPs example dataset and 2) creating subset snpR objects.
#' a <- 2013:2015 #create a vector of possible birthyears
#' b <- c("M", "F", "U") #create a vector of possible sexes
#' stk <- stickSNPs
#' sample.meta(stk)$BirthYear <-
#' sample(x = a, size = nsamps(stk), replace = TRUE) #create birthyears
#' sample.meta(stk)$ID <- 1:nsamps(stk) #create unique sampleID
#' sample.meta(stk)$Sex <-
#' sample(x= b, size = nsamps(stk), replace = TRUE) # create sexes
#' #generating snpR objects for male and female potential parents and offspring
#' # (no U sexes in the potential parents in for this example)
#' # get list of samples which are now "M" for subsetting
#' lsir <- which(stk@sample.meta$Sex =="M" &
#' stk@sample.meta$BirthYear == "2013") #list sires
#' ldam <- which(stk@sample.meta$Sex =="F" &
#' stk@sample.meta$BirthYear == "2013") #list dams
#' loff <- which(stk@sample.meta$BirthYear %in% c("2014","2015")) #list offspr
#'
#'
#' #creating new snpR objects for Colony formatting
#' sir <- subset_snpR_data(x = stk, .samps = lsir)
#' dam <- subset_snpR_data(x = stk, .samps = ldam)
#' off <- subset_snpR_data(x = stk, .samps = loff)
#' # not run to avoid clutter
#' \dontrun{
#' write_colony_input(x = off, outfile = "parents_example.col",
#' maternal_genotypes = dam, paternal_genotypes = sir)
#' }
#'
#' # running a simple model
#' \dontrun{
#' ## intentionally shorter run, with a small subset of the samples
#' asp <- which(stickSNPs@sample.meta$pop == "ASP")
#' test_dat <- subset_snpR_data(stickSNPs, .samps = asp)
#' run_colony(x = test_dat, colony_path = "/usr/bin/colony2s.exe",
#' method = "PLS", run_length = 1)
#' }
NULL
#' @describeIn colony_interface Create a colony input file
#' @export
write_colony_input <- function(x, outfile = "colony_input", method = "FPLS", run_length = 2, sampleIDs = NULL,
sibship_prior = 0, paternal_sib_size = NULL, maternal_sib_size = NULL,
nruns = 1, seed = NULL, maternal_genotypes = NULL, paternal_genotypes = NULL,
maternal_inclusion_prob = 0, paternal_inclusion_prob = 0,
update_af = TRUE, dioecious = TRUE, inbreeding = TRUE, male_monogamous = FALSE,
female_monogamous = FALSE, clone_inference = FALSE, sibship_scaling = TRUE, known_af = FALSE,
precision = 2, dropout = 0.01, genotyping_error = 0.01,
known_maternal_dyads = NULL, known_paternal_dyads = NULL,
known_maternal_max_mismatches = 0, known_paternal_max_mismatches = 0,
known_maternal_sibships = NULL, known_paternal_sibships = NULL,
maternal_exclusions = NULL, paternal_exclusions = NULL,
excluded_maternal_siblings = NULL, excluded_paternal_siblings = NULL,
update_bib = FALSE){
#=====================initialize===============
# initialize storage directory
if(basename(getwd()) != "colony"){
if(!dir.exists("colony")){
dir.create("colony")
setwd("colony")
}
else{
setwd("colony")
}
}
# write anything preliminary to output
proj <- outfile
if(nchar(proj) > 40){
proj <- substr(proj, 1, 40)
}
outfile <- paste0(proj, ".dat")
#=====================sanity checks============
# check if file exists
if(file.exists(outfile)){
file.remove(outfile)
}
#=====================write first stuff=============
# project name
write(proj, outfile)
# outfile
write(proj, outfile, append = T)
# other info:
write(ncol(x), outfile, append = T) # sample size
write(nrow(x), outfile, append = T) # nloci
if(is.null(seed)){
seed <- sample(1000000, 1)
}
write(seed, outfile, append = T) # seed
write(as.numeric(update_af), outfile, append = T) # update allele frequency in bayesian model
if(dioecious){
dioecious <- 2
}
else{
dioecious <- 1
}
write(dioecious, outfile, append = T) # dio or monoecious
write(as.numeric(inbreeding), outfile, append = T) # inbreeding
write(0, outfile, append = T) # diploid only for now
write(paste(as.numeric(male_monogamous), as.numeric(female_monogamous), collapse = " "), outfile, append = T) # polygamy
write(as.numeric(clone_inference), outfile, append = T) # infer clones
write(as.numeric(sibship_scaling), outfile, append = T) # scale sibships
# sibship scaling priors
sib.prior <- sibship_prior
if(!is.null(paternal_sib_size)){
sib.prior <- paste(sib.prior, paternal_sib_size, maternal_sib_size)
}
write(sib.prior, outfile, append = T)
#=====================reformat genotypes, write info on them==================
colony_genotypes <- format_snps(x, "colony")
if(!is.null(sampleIDs)){
colony_genotypes[,1] <- x@sample.meta[,sampleIDs]
}
# allele frequencies
## if provided with minor allele frequencies
if(length(known_af) == nrow(x) & is.numeric(known_af)){
write(1, outfile, append = T)
write(rep(2, nrow(x)), outfile, append = T, sep = " ", ncolumns = nrow(x)) # number of alleles per locus, should all be two.
maj <- 1 - known_af
afs <- cbind(maj, known_af)
afs_tab <- matrix(0, nrow = 2*nrow(afs), ncol = 2)
afs_tab[seq(2, nrow(afs_tab), by = 2),] <- afs
afs_tab[seq(1, nrow(afs_tab), by = 2),] <- c(rep(1, nrow(afs)), rep(2, nrow(afs)))
utils::write.table(afs_tab, outfile, sep = " ", append = T, row.names = F, col.names = F)
}
# otherwise, check the logical and calculate minor allele frequencies if true
else if(is.logical(known_af[1]) & length(known_af) == 1){
write(as.numeric(known_af), outfile, append = T)
if(known_af){
write(rep(2, nrow(x)), outfile, append = T, sep = " ", ncolumns = nrow(x)) # number of alleles per locus, should all be two.
# afs
x <- calc_maf(x)
afs <- x@stats[x@stats$facet == ".base",]$maf
afs <- cbind(1-afs, afs)
afs_tab <- matrix(0, nrow = 2*nrow(afs), ncol = 2)
afs_tab[seq(2, nrow(afs_tab), by = 2),] <- afs
afs_tab[seq(1, nrow(afs_tab), by = 2),] <- c(rep(1, nrow(afs)), rep(2, nrow(afs)))
utils::write.table(afs_tab, outfile, sep = " ", append = T, row.names = F, col.names = F)
}
}
else{
file.remove(outfile)
stop("known_af must either be TRUE, FALSE, or a vector of minor allele frequencies.\n")
}
write(nruns, outfile, append = T) # nruns
write(run_length, outfile, append = T) # run length
write(0, outfile, append = T) # monitor method
write(10000, outfile, append = T) # monitor interval
write(0, outfile, append = T) # no GUI
# analysis method
analysis_methods <- c("PLS", "FL", "FPLS")
method <- match(method, analysis_methods) - 1
write(method, outfile, append = T)
write(precision, outfile, append = T) # precision
write("MK@", outfile, append = T) # marker names
write("0@", outfile, append = T) # marker types
# genotyping error and dropout
## dropout
if(length(dropout) == 1){
write(rep(dropout, nrow(x)), outfile, append = T, ncolumns = nrow(x), sep = " ")
}
else{
write(dropout, outfile, append = T, ncolumns = nrow(x), sep = " ")
}
## error rates
if(length(genotyping_error) == 1){
write(rep(genotyping_error, nrow(x)), outfile, append = T, ncolumns = nrow(x), sep = " ")
}
else{
write(genotyping_error, outfile, append = T, ncolumns = nrow(x), sep = " ")
}
# offspring genotypes
utils::write.table(colony_genotypes, outfile, T, quote = F, sep = " ", row.names = F, col.names = F)
#===================maternal and paternal genotype things==============
# number of male and female candidates, maternal and paternal inclusion probabilities in OMS and PMS
#there needs to be a newline between the offspring genotypes and the parent data #eg
# offspring genotype end
# "\n"
# 0 0 # M P
# 0 0 # M P
# "\n"
## then a newline before first parental genotypes
write("\n", outfile, append = T) # adding a newline between offspring genotypes and parents
par_nums <- numeric(2)
if(is.snpRdata(paternal_genotypes)){
par_nums[1] <- ncol(paternal_genotypes) #ncol not nrow!
}
else{
paternal_inclusion_prob <- 0
}
if(is.snpRdata(maternal_genotypes)){
par_nums[2] <- ncol(maternal_genotypes) #ncol not nrow!
}
else{
maternal_inclusion_prob <- 0
}
write(c(paternal_inclusion_prob, maternal_inclusion_prob), outfile, append = T, sep = " ") # write inclusion probs
write(par_nums, outfile, append = T, sep = " ") # OPS and OMS sample sizes
write("\n", outfile, append = T) #added newline after number parents and probs parents
# paternal and maternal genotypes
if(is.snpRdata(paternal_genotypes)){
male_colony <- format_snps(paternal_genotypes, "colony")
utils::write.table(male_colony, outfile, T, quote = F, sep = " ", row.names = F, col.names = F)
write("\n", outfile, append = T) #added newline for separating genotypes
}
write("\n", outfile, append = T) #added newline
if(is.snpRdata(maternal_genotypes)){
female_colony <- format_snps(maternal_genotypes, "colony")
utils::write.table(female_colony, outfile, T, quote = F, sep = " ", row.names = F, col.names = F)
write("\n", outfile, append = T) #added newline for separating stuff
}
# paternal dyads:
if(!is.null(known_paternal_dyads[1,1])){
n_p_d <- nrow(known_paternal_dyads)
write(c(n_p_d, known_paternal_max_mismatches), outfile, append = T, sep = " ") # write the number known + mismatch max
utils::write.table(known_paternal_dyads, outfile, T, quote = F, sep = " ", row.names = F, col.names = F) # write dyads. They should be a matrix or data frame, with the offspring ID in the first column and the paternal ID in the second
}
else{
write(c(0, 0), outfile, append = T, sep = " ") # no dyads
}
# maternal dyads:
if(!is.null(known_maternal_dyads[1,1])){
n_m_d <- nrow(known_maternal_dyads)
write(c(n_m_d, known_maternal_max_mismatches), outfile, append = T, sep = " ") # write the number known + mismatch max
utils::write.table(known_maternal_dyads, outfile, T, quote = F, sep = " ", row.names = F, col.names = F) # write dyads. They should be a matrix or data frame, with the offspring ID in the first column and the maternal ID in the second
}
else{
write(c(0, 0), outfile, append = T, sep = " ") # no dyads
}
# known paternal sibships
if(!is.null(known_paternal_sibships[1,1])){
n_p_s <- nrow(known_paternal_sibships)
write(n_p_s, outfile, append = T) # write the number known
utils::write.table(known_paternal_sibships, outfile, T, quote = F, sep = " ", row.names = F, col.names = F) # write sibsips. Data frame or matrix with sibship size followed by single column containing all of the sibling IDs.
}
else{
write(0, outfile, append = T) # no sibships
}
# known maternal sibships
if(!is.null(known_maternal_sibships[1,1])){
n_m_s <- nrow(known_maternal_sibships)
write(n_m_s, outfile, append = T) # write the number known
utils::write.table(known_maternal_sibships, outfile, T, quote = F, sep = " ", row.names = F, col.names = F) # write sibsips. Data frame or matrix with sibship size followed by single column containing all of the sibling IDs.
}
else{
write(0, outfile, append = T) # no sibships
}
# paternal exclusions
if(!is.null(paternal_exclusions[1,1])){
npe <- nrow(paternal_exclusions)
write(npe, outfile, append = T) # write number of exclusions
utils::write.table(paternal_exclusions, outfile, T, quote = F, sep = " ", row.names = F, col.names = F) # write the paternal exclusions. Data.frame or matrix with col one the offspring ID, col 2 the number of excluded males, col 3 the IDs of excluded males separated by spaces.
}
else{
write(0, outfile, append = T) # no exclusions
}
# maternal exclusions
if(!is.null(maternal_exclusions[1,1])){
npe <- nrow(maternal_exclusions)
write(npe, outfile, append = T) # write number of exclusions
utils::write.table(maternal_exclusions, outfile, T, quote = F, sep = " ", row.names = F, col.names = F) # write the maternal exclusions. Data.frame or matrix with col one the offspring ID, col 2 the number of excluded males, col 3 the IDs of excluded females separated by spaces.
}
else{
write(0, outfile, append = T) # no exclusions
}
# excluded paternal sibships
if(!is.null(excluded_paternal_siblings[1,1])){
n_eps <- nrow(excluded_paternal_siblings)
write(n_eps, outfile, append = T) # number of exclusions
utils::write.table(excluded_paternal_siblings, outfile, T, quote = F, sep = " ", row.names = F, col.names = F) # write the excluded paternal sibships. Data.frame or matrix with col one the offspring ID, col 2 the number of excluded siblings, col 3 the IDs of excluded siblings separated by spaces.
}
else{
write(0, outfile, append = T) # no exclusions
}
# excluded maternal sibships
if(!is.null(excluded_maternal_siblings[1,1])){
n_eps <- nrow(excluded_maternal_siblings)
write(n_eps, outfile, append = T) # number of exclusions
utils::write.table(excluded_maternal_siblings, outfile, T, quote = F, sep = " ", row.names = F, col.names = F) # write the excluded maternal sibships. Data.frame or matrix with col one the offspring ID, col 2 the number of excluded siblings, col 3 the IDs of excluded siblings separated by spaces.
}
else{
write(0, outfile, append = T) # no exclusions
}
#=================return to original wd===========
setwd("..")
return(paste0("Colony input saved to: ./colony/", outfile, collapse = ""))
}
#' Call a prepared colony infile.
#' @describeIn colony_interface Call a colony executable to run a prepared colony input file.
#' @export
call_colony <- function(infile, colony_path, update_bib = FALSE, verbose = TRUE){
# check system type
sys.type <- Sys.info()["sysname"]
# save the call
call <- paste0(colony_path," IFN:", infile)
# call
system(call, show.output.on.console = verbose)
# move results
pattern <- gsub(".+/", "", infile)
pattern <- gsub("\\..+$", "", pattern)
files <- list.files(pattern = pattern)
if(!dir.exists("colony")){
dir.create("colony")
}
file.rename(files, paste0("colony/", files))
.yell_citation("Jones2010", "parentage", "parentage or sibship reconstruction with COLONY2", update_bib)
}
#' Parse colony data.
#' @describeIn colony_interface Parse a previously run colony analysis.
#' @export
parse_colony <- function(prefix, x, path = "colony/", sampleIDs = NULL, cleanup = FALSE){
#===============full and half sibs==============================
# read in half and full sibs
fsd <- data.table::fread(paste0(path, prefix, ".FullSibDyad"))
if(file.exists(paste0(path, prefix, ".HalfSibDyad"))){
hsd <- data.table::fread(paste0(path, prefix, ".HalfSibDyad"))
## combine
if(nrow(fsd) > 0 & nrow(hsd) > 0){
dyads <- rbind(cbind(fsd, type = "FullSib", stringsAsFactors = F),
cbind(hsd, type = "HalfSib", stringsAsFactors = F))
}
else if(nrow(fsd) > 0){
dyads <- cbind(fsd, type = "FullSib", stringsAsFactors = F)
}
else if(nrow(hsd) > 0){
dyads <- cbind(hsd, type = "HalfSib", stringsAsFactors = F)
}
}
else{
dyads <- cbind(fsd, type = "FullSib", stringsAsFactors = F)
}
colnames(dyads)[1:2] <- c("Sample1", "Sample2")
## add in non-sibs
all_pairs <- as.data.frame(t(utils::combn(colnames(x), 2)), stringsAsFactors = F)
colnames(all_pairs) <- c("Sample1", "Sample2")
cdyads <- rbind(dyads, dyads[,c(2,1,3,4)])
all_pairs <- merge(all_pairs, cdyads, by = c("Sample1", "Sample2"), all.x = T)
all_pairs$type[which(is.na(all_pairs$type))] <- "none"
#=============add male and female parent ID info to the sample meta===========
clusters <- data.table::fread(paste0(path, prefix, ".BestCluster"))
colnames(clusters)[2] <- "ClusterProbability"
if(is.null(sampleIDs)){
x@sample.meta$.offspringID <- colnames(x)
x@sample.meta <- merge(x@sample.meta, clusters, by.x = ".offspringID", by.y = "OffspringID", sort = F)
x@sample.meta$.offspringID <- NULL
}
else{
x@sample.meta <- merge(x@sample.meta, clusters, by.x = sampleIDs, by.y = "OffspringID", sort = F)
}
.sidcol <- which(colnames(x@sample.meta) == ".sample.id")
x@sample.meta <- x@sample.meta[,c((1:ncol(x@sample.meta))[-.sidcol], .sidcol)]
#=============cleanup and return================
if(cleanup){
if(Sys.info()["sysname"] == "Windows"){
shell(paste0("rm -r ./", path))
}
else{
system(paste0("rm -r ./", path))
}
}
return(list(x = x, dyads = dyads, all_pairs = all_pairs, clusters = clusters))
}
#' Create and infile, run, and gather some basic results using the COLONY parentage analysis program.
#'
#' Create a COLONY infile using snpRdata sets containing offspring and possibly paternal genotypes given
#' specified parameters. Requires that the COLONY program is installed locally. Input and output files
#' will be stored in a colony folder created in the current working directory.
#'
#' @export
#' @describeIn colony_interface run colony on a snpRdata object, start to finish
run_colony <- function(x, colony_path, outfile = "colony_input", method = "FPLS", run_length = 2, sampleIDs = NULL,
sibship_prior = 0, paternal_sib_size = NULL, maternal_sib_size = NULL,
nruns = 1, seed = NULL, maternal_genotypes = NULL, paternal_genotypes = NULL,
maternal_inclusion_prob = 0, paternal_inclusion_prob = 0,
update_af = TRUE, dioecious = TRUE, inbreeding = TRUE, male_monogamous = FALSE,
female_monogamous = FALSE, clone_inference = FALSE, sibship_scaling = TRUE, known_af = FALSE,
precision = 2, dropout = 0.01, genotyping_error = 0.01,
known_maternal_dyads = NULL, known_paternal_dyads = NULL,
known_maternal_max_mismatches = 0, known_paternal_max_mismatches = 0,
known_maternal_sibships = NULL, known_paternal_sibships = NULL,
maternal_exclusions = NULL, paternal_exclusions = NULL,
excluded_maternal_siblings = NULL, excluded_paternal_siblings = NULL, update_bib = FALSE,
cleanup = FALSE, verbose = TRUE){
msg <- character(0)
if(!file.exists(colony_path)){
msg <- c(msg, "Cannot find colony executable.\n")
}
if(!is.snpRdata(x)){
msg <- c(msg, "x is not a snpRdata object.\n")
}
.check.installed("readr")
if(length(msg) > 0){
stop(msg)
}
# prep
write_colony_input(x = x,
outfile = outfile,
method = method,
run_length = run_length,
sampleIDs = sampleIDs,
sibship_prior = sibship_prior,
paternal_sib_size = paternal_sib_size,
maternal_sib_size = maternal_sib_size, nruns = nruns, seed = seed,
maternal_genotypes = maternal_genotypes,
paternal_genotypes = paternal_genotypes,
maternal_inclusion_prob = maternal_inclusion_prob,
paternal_inclusion_prob = paternal_inclusion_prob,
update_af = update_af,
dioecious = dioecious,
inbreeding = inbreeding,
male_monogamous = male_monogamous,
female_monogamous = female_monogamous,
clone_inference = clone_inference,
sibship_scaling = sibship_scaling,
known_af = known_af,
precision = precision,
dropout = dropout,
genotyping_error = genotyping_error,
known_maternal_dyads = known_maternal_dyads,
known_paternal_dyads = known_paternal_dyads,
known_maternal_max_mismatches = known_maternal_max_mismatches,
known_paternal_max_mismatches = known_paternal_max_mismatches,
known_maternal_sibships = known_maternal_sibships,
known_paternal_sibships = known_paternal_sibships,
maternal_exclusions = maternal_exclusions,
paternal_exclusions = paternal_exclusions,
excluded_maternal_siblings = excluded_maternal_siblings,
excluded_paternal_siblings = excluded_paternal_siblings)
# run
call_colony(infile = paste0("colony/", outfile, ".dat"),
colony_path = colony_path,
update_bib = update_bib, verbose = verbose)
# parse some basics
return(parse_colony(prefix = outfile,
x = x,
path = "colony/",
sampleIDs = sampleIDs,
cleanup = cleanup))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.