Skip to content

Commit

Permalink
small fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
mijangos81 committed Jan 24, 2022
1 parent 7ad64c3 commit 9162d13
Show file tree
Hide file tree
Showing 15 changed files with 155 additions and 225 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,7 @@ biocViews:
Imports: stats, methods, utils, plyr, tidyr, MASS, stringr, ape, vegan,
SNPRelate, StAMPP, sp, PopGenReport, robustbase,
gridExtra, foreach, dplyr, crayon, devtools,reshape2,
patchwork, data.table, hierfstat
patchwork, hierfstat, data.table
Suggests:
tibble,
HardyWeinberg,
Expand Down Expand Up @@ -62,7 +62,6 @@ Suggests:
plotly,
ggthemes,
ggrepel,
raster,
strataG,
ggtern,
markdown,
Expand All @@ -71,6 +70,7 @@ Suggests:
SIBER,
fields,
snpStats,
raster,
boot
License: GPL-2
LazyData: true
Expand Down
4 changes: 2 additions & 2 deletions NAMESPACE
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
# Generated by roxygen2: do not edit by hand

export(AIC)
export(AIC.comp)
export(CLR)
export(AIC_comp)
export(fsc.bootstraps)
export(fsc.estimate)
export(fsc.multiple.estimate)
export(gi2gl)
Expand Down
99 changes: 59 additions & 40 deletions R/fastsimcoal.R
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@

#' Run fastsimcoal to estimate parameters' values
#'
#' Run fastsimcoal assuming that in \code{dir.in} there are the tpl, est and obs
Expand All @@ -9,7 +8,6 @@
#' default value). If \code{ncpu}>0 and \code{nBatches=NULL}, \code{nBatches} is
#' set to twice ]code{nspu}.
#'
#'
#' @param dir.in The path where to run fsc
#' @param n The number of coalescent simulations to approximate the expected SFS
#' (-n option). This should be larger than 100,000.
Expand All @@ -24,13 +22,20 @@
#' @param fsc.path The path where fsc is installed or \code{"path"} if it is in
#' the PATH (that means that it can be called regardless of the working
#' directory)
#'@author Carlo Pacioni (Post to \url{https://groups.google.com/d/forum/dartr})
#' @author Carlo Pacioni (Post to \url{https://groups.google.com/d/forum/dartr})
#' @export
#' @references Excoffier L., Dupanloup I., Huerta-Sánchez E., Sousa V. C. and
#' Foll M. (2013) Robust demographic inference from genomic and SNP data. PLoS
#' genetics 9(10)
fsc.estimate <- function(dir.in, n=500000, L=100, maf=TRUE, ncpu=0, nBatches=NULL,
fsc.cmd="fsc2702", fsc.path="path") {
fsc.estimate <- function(dir.in,
n=500000,
L=100,
maf=TRUE,
ncpu=0,
nBatches=NULL,
fsc.cmd="fsc2702",
fsc.path="path") {

if(!is.integer(as.integer(ncpu))) stop("The number of CPUs has to be an integer")
if(ncpu == 0 & is.null(nBatches)) nBatches <- 12
if(ncpu > 0 & is.null(nBatches)) nBatches <- 2 * ncpu
Expand Down Expand Up @@ -62,27 +67,26 @@ fsc.multiple.estimate <- function(dir.in, n=500000, L=100, maf=TRUE, ncpu=0,

}


#' AIC
#'
#' compute the AIC given the likelihood in log10 scale (default fsc output) and
#' number of parameters, which is automatically read by the function if
#' \code{k=NULL).
#'
#' @name AIC
#' @title Compute the AIC given the likelihood in log10 scale
#' @description Compute the AIC given the likelihood in log10 scale (default fsc
#' output) and number of parameters, which is automatically read by the
#' function if \code{k=NULL}.
#' @param dir.in The directory where the analysis was conducted
#' @param k The number of parameters in the model. If \code{NULL} is read
#' automatically by the function
#' @details
#' It is important to note that this function will probably fail if there are
#' more than one .est within each directory as it will read these to evaluate
#' the number of parameters in the model. If there is more than one file, there
#' is no way to know which one is correct.
#'
#' @param dir.in The directory where the analysis was conducted
#' @param k The number of parameters in the model. If \code{NULL} is read
#' automatically by the function
#'@author Carlo Pacioni (Post to \url{https://groups.google.com/d/forum/dartr})
#' @author Carlo Pacioni (Post to \url{https://groups.google.com/d/forum/dartr})
#' @export
#' @references Excoffier L., Dupanloup I., Huerta-Sánchez E., Sousa V. C. and
#' Foll M. (2013) Robust demographic inference from genomic and SNP data. PLoS
#' genetics 9(10)
AIC <- function(dir.in, k=NULL) {
AIC <- function(dir.in,
k=NULL) {
log10toln<-function(l10) {
rlns=l10/log10(exp(1))
}
Expand All @@ -106,13 +110,14 @@ AIC <- function(dir.in, k=NULL) {
return(2*k-2*ln)
}

#' Compute and compare AIC for different fasimcoal models
#'
#' @name AIC_comp
#' @title Compute and compare AIC for different fastsimcoal models
#' @description
#' This functions expects that different fastsimcoal models were run in
#' different directories within the same parent directory (as for
#' \code{fsc.multiple.estimate}). Assuming that \code{dir.in} is the parent
#' directory, it will extract and compute AIC from all of them.
#'
#' @details
#' It is important to note that this function will probably fail if there are
#' more than one .est within each directory as it will read these to evaluate
#' the number of parameters in the model. If there is more than one file, there
Expand All @@ -121,15 +126,15 @@ AIC <- function(dir.in, k=NULL) {
#' cause problems in reporting the results even if the function completes.
#'
#' @param dir.in The parent directory within which the models have been run
#'@author Carlo Pacioni (Post to \url{https://groups.google.com/d/forum/dartr})
#' @author Carlo Pacioni (Post to \url{https://groups.google.com/d/forum/dartr})
#' @export
#' @references Excoffier L., Dupanloup I., Huerta-Sánchez E., Sousa V. C. and
#' Foll M. (2013) Robust demographic inference from genomic and SNP data. PLoS
#' genetics 9(10)
AIC.comp <- function(dir.in) {
AIC_comp <- function(dir.in) {
mod.nms <- list.dirs(dir.in, full.names=FALSE, recursive=FALSE)
ld <- list.dirs(dir.in,full.names=TRUE, recursive=FALSE)
lAICs <- sapply(ld, amplicR::AIC, USE.NAMES=FALSE)
lAICs <- sapply(ld, AIC, USE.NAMES=FALSE)
mod.rank <- rank(lAICs)
return(data.frame(Model=mod.nms, AIC=lAICs, Delta=lAICs - min(lAICs), Rank=mod.rank))
}
Expand Down Expand Up @@ -179,27 +184,41 @@ AIC.comp <- function(dir.in) {
#'@param boot.type The method to be used to compute the confidence intervals.
#' See \code{?boot.ci} for details. By default, percentile intervals are
#' computed.
#'@inheritParams fsc.estimate
#'@author Carlo Pacioni (Post to \url{https://groups.google.com/d/forum/dartr})
#'@return A list with the following elements \itemize{ \item Bootstr.stats:
#' Descriptive statistics from bootstraps (Median, lower and upper limit), an
#' the initial estimated parameters \item P.Rand.less.Obs: The probability that
#' a random value from the null Composite Likelihood Ratio distribution is less
#' than the observed CLR \item P.Rand.gt.Obs: The probability that a random
#' value from the null Composite Likelihood Ratio distribution is greater than
#' the observed CLR \item Sim: The estimates from the simulated data }
#' @inheritParams fsc.estimate
#' @author Carlo Pacioni (Post to \url{https://groups.google.com/d/forum/dartr})
#' @return A list with the following elements:
#' \itemize{
#' \item Bootstr.stats: Descriptive statistics from bootstraps (Median, lower and
#' upper limit), an the initial estimated parameters.
#' \item P.Rand.less.Obs: The probability that a random value from the null
#' Composite Likelihood Ratio distribution is less than the observed CLR.
#' \item P.Rand.gt.Obs: The probability that a random value from the null
#' Composite Likelihood Ratio distribution is greater than the observed CLR.
#' \item Sim: The estimates from the simulated data.
#' }
#' @export
#'@references Excoffier L., Dupanloup I., Huerta-Sánchez E., Sousa V. C. and
#' Foll M. (2013) Robust demographic inference from genomic and SNP data. PLoS
#' genetics 9(10)
#'@export

CLR <- MaxEstLhood <- MaxObsLhood <- Param <- NULL

fsc.bootstraps <- function(dir.in, nLoci=10000, nSim=100, maf=TRUE, ncpu=0,
nBatches=NULL, fsc.cmd="fsc2702", fsc.path="path",
par.indLoci="200000 0", par.nBlocks=1,
fsc.bootstraps <- function(dir.in,
nLoci=10000,
nSim=100,
maf=TRUE,
ncpu=0,
nBatches=NULL,
fsc.cmd="fsc2702",
fsc.path="path",
par.indLoci="200000 0",
par.nBlocks=1,
par.data="DNA 100 0 2.5e-8 0.33",
n=100000, L=50, nBoot=1000, conf=0.95, boot.type="perc") {
n=100000,
L=50,
nBoot=1000,
conf=0.95,
boot.type="perc") {

CLR <- MaxEstLhood <- MaxObsLhood <- Param <- NULL
#---------- Helper ---------------#
med.i <- function(x, i) median(x[i])

Expand Down
2 changes: 1 addition & 1 deletion R/gl.impute.r
Original file line number Diff line number Diff line change
Expand Up @@ -127,7 +127,7 @@ gl.impute <- function(x,
x3 <- NULL
# merge back populations
for (pop in pop_list) {
x3 <- rbind (x3, pop)
x3 <- rbind(x3, pop)
}
}
# if 1 population
Expand Down
12 changes: 8 additions & 4 deletions R/gl.report.hwe.r
Original file line number Diff line number Diff line change
Expand Up @@ -109,12 +109,16 @@
#'
#' Ternary plots can be used to visualise patterns of H-W proportions (plot.out
#' = TRUE). P-values and the statistical (non)significance of a large number of
#' bi- allelic markers can be inferred from their position in a ternary plot.
#' bi-allelic markers can be inferred from their position in a ternary plot.
#' See Graffelman & Morales-Camarena (2008) for further details. Ternary plots
#' are based on the function \code{\link[HardyWeinberg]{HWTernaryPlot}} from
#' the package HardyWeinberg. Loci that depart significantly from H-W
#' proportions are shown in red, and those not showing significant departure are
#' shown in green.
#' the package HardyWeinberg. Each vertex of the Ternary plot represents one of
#' the three possible genotypes for SNP data: homozygous for the reference
#' allele (AA), heterozygous (AB) and homozygous for the alternative allele
#' (BB). Loci deviating significantly from Hardy-Weinberg proportions after
#' correction for multiple tests are shown in pink. The blue parabola
#' represents Hardy-Weinberg equilibrium, and the area between green lines
#' represents the acceptance region.
#'
#' For these plots to work it is necessary to install the package ggtern.
#' @return A dataframe containing loci, counts of reference SNP homozygotes,
Expand Down
2 changes: 1 addition & 1 deletion R/gl.report.ld.map.r
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
#'
#' This function requires that SNPs to be mapped to a reference genome and the
#' information for SNP's position must be stored in the genlight accessor
#' @@position" and the SNP's chromosome name in the accessor @@chromosome.
#' position" and the SNP's chromosome name in the accessor chromosome.
#'
#' @param x Name of the genlight object containing the SNP data [required].
#' @param ld_max_pairwise Maximum distance in number of base pairs at which LD
Expand Down
11 changes: 7 additions & 4 deletions R/gl2sfs.R
Original file line number Diff line number Diff line change
Expand Up @@ -23,10 +23,13 @@
#'
#'It expects a dartR formatted genlight object, but it should also work with
#'other genlight objects.
#'
#'@param outfile_root The root of the name of the output file
#'@inheritParams gl.report.heterozygosity
#'@inheritParams gl2vcf
#' @param x Name of the genlight object containing the SNP data [required].
#' @param n.invariant.tags Number of invariant sites[default 0].
#' @param outfile_root The root of the name of the output file [default "gl2sfs"].
#' @param outpath Path where to save the output file [default tempdir()].
#' @param verbose Verbosity: 0, silent or fatal errors; 1, begin and end; 2,
#' progress log; 3, progress and results summary; 5, full report
#' [default 2, unless specified using gl.set.verbosity].
#'@return A list with two elements: the DAF and MAF.
#'@author Custodian: Carlo Pacioni (Post to
#' \url{https://groups.google.com/d/forum/dartr})
Expand Down
6 changes: 4 additions & 2 deletions man/AIC.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

8 changes: 4 additions & 4 deletions man/AIC.comp.Rd → man/AIC_comp.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

65 changes: 50 additions & 15 deletions man/CLR.Rd → man/fsc.bootstraps.Rd

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

Loading

0 comments on commit 9162d13

Please sign in to comment.