Package 'LocaTT'

Description

Performs adjustments to a taxonomy system according to a taxonomy edits file.

Usage

adjust_taxonomies(
  path_to_input_file,
  path_to_output_file,
  path_to_taxonomy_edits
)

Arguments

Value

No return value. Writes an output CSV file with adjusted taxonomies.

See Also

get_taxonomies.species_binomials for remotely fetching NCBI taxonomies from species binomials.

get_taxonomies.IUCN for formatting taxonomies from the IUCN Red List.

Examples

# Get path to input file.
path_to_input_file<-system.file("extdata",
                                "example_local_taxa_list.csv",
                                package="LocaTT",
                                mustWork=TRUE)

# Get path to taxonomy edits.
path_to_taxonomy_edits<-system.file("extdata",
                                    "example_taxonomy_edits.csv",
                                    package="LocaTT",
                                    mustWork=TRUE)

# Create temporary output file path.
path_to_output_file<-tempfile(fileext=".csv")

# Adjust taxonomies.
adjust_taxonomies(path_to_input_file=path_to_input_file,
                  path_to_output_file=path_to_output_file,
                  path_to_taxonomy_edits=path_to_taxonomy_edits)

Description

Performs binomial tests.

Usage

binomial_test(k, n, p, alternative = "greater")

Arguments

Details

Calls on the pbinom function in the stats package to perform vectorized binomial tests. Arguments are recycled as in pbinom. Only one-sided tests are supported, and only p-values are returned.

Value

A numeric vector of p-values from the binomial tests.

Examples

binomial_test(k=c(5,1,7,4),
              n=c(10,3,15,5),
              p=c(0.2,0.1,0.5,0.6),
              alternative="greater")

Description

Checks whether a BLAST program can be found.

Usage

blast_command_found(blast_command)

Arguments

Value

Logical. Returns TRUE if the BLAST program could be found.

Examples

blast_command_found(blast_command="blastn")

Description

Gets the version of a BLAST program.

Usage

blast_version(blast_command = "blastn")

Arguments

Value

Returns a string of the version of the BLAST program.

Examples

blast_version()

Description

Draws circle polygon.

Usage

circle(r, v = 1000, ...)

Arguments

Details

Draws a circle polygon of a given radius. The circle is drawn about the origin (i.e., x = 0, y = 0). Intended for use with template to generate detection and proportion plots.

Value

No return value.

See Also

sector for plotting sector polygons.

Examples

template(l=1)
circle(r=1)

Description

Checks whether DNA sequences contain wildcard characters.

Usage

contains_wildcards(sequences)

Arguments

Value

A logical vector indicating whether each DNA sequence contains wildcard characters.

Examples

contains_wildcards(sequences=c("TKCTAGGTGW","CATAATTAGG","ATYGGCTATG"))

Description

Generates coordinates along a circular path.

Usage

coordinates(a, r)

Arguments

Details

Calculates xy coordinates along a circular path given a vector of angles and a specified radius. The center of the circle is at the origin (i.e., x = 0, y = 0). This function is helpful for calculating the vertices of circle and sector polygons.

Value

A numeric matrix of xy coordinates.

See Also

circle for plotting circle polygons.

sector for plotting sector polygons.

Examples

coordinates(a=c(90,180,270,360),r=1)

Description

Derives covariance matrix from correlation matrix and standard deviation vector.

Usage

cor2cov(sd, R)

Arguments

Details

Given correlation matrix R and standard deviation vector sd, performs the operation diag(sd) %*% R %*% diag(sd) to derive the corresponding covariance matrix. This is a counterpart to stats::cov2cor, which scales a covariance matrix into the corresponding correlation matrix.

Value

Returns a numeric covariance matrix.

See Also

stats::cov2cor for scaling a covariance matrix into the corresponding correlation matrix.

Examples

# Define standard deviation vector.
sd<-c(9.655,1.157,1.128,2.925)

# Define correlation matrix.
R<-matrix(data=c(1.000,-0.80,0.64,-0.512,
                 -0.800,1.00,-0.80,0.640,
                 0.640,-0.80,1.00,-0.800,
                 -0.512,0.64,-0.80,1.000),
           ncol=4,byrow=TRUE)

# Derive covariance matrix.
cor2cov(sd=sd,R=R)

Description

Density function for the Gaussian copula.

Usage

dcopula(u, R, log = FALSE)

Arguments

Details

Computes the probability density of the Gaussian copula. Given uniformly-distributed margins u on the interval [0, 1], applies the inverse cumulative distribution function of the standard normal (i.e., stats::qnorm) to map uniform margins to normal scores. Then, uses equation 1 of Song (2000) with the normal scores to calculate the probability density of the Gaussian copula.

Value

Numeric vector of probability densities.

References

Song P. 2000. Multivariate dispersion models generated from Gaussian copula. Scandinavian Journal of Statistics, 27(2): 305-320. DOI: 10.1111/1467-9469.00191

See Also

dmvlogis for density of the multivariate logistic distribution.

Examples

# Define uniform margins.
u<-c(0.324,0.383,0.917,0.015)

# Define correlation matrix.
R<-matrix(data=c(1.000,-0.80,0.64,-0.512,
                 -0.800,1.00,-0.80,0.640,
                 0.640,-0.80,1.00,-0.800,
                 -0.512,0.64,-0.80,1.000),
           ncol=4,byrow=TRUE)

# Compute log probability density.
dcopula(u=u,R=R,log=TRUE)

Description

Density function for the Dirichlet-multinomial distribution.

Usage

ddirmult(x, p, theta, alpha, log = FALSE)

Arguments

Details

Computes the probability mass of the Dirichlet-multinomial distribution. Under the proportion parameterization, the alpha parameters of the conventional Dirichlet-multinomial distribution are derived as the product of a proportion vector (p) and an exponentiated precision parameter (exp(theta)). The precision parameter controls the degree of overdispersion relative to the multinomial distribution, where higher values of theta are associated with reduced overdispersion. When theta = 0, the alpha parameters of the conventional Dirichlet-multinomial distribution are equal to the proportion vector (p). To ensure a simplex, the values of p (vector or matrix records) are internally normalized to sum to one. If alpha is provided, then the conventional alpha parameterization of the Dirichlet-multinomial distribution is used.

Value

Numeric vector of probability densities.

References

Minka TP. 2000. Estimating a Dirichlet distribution.

See Also

waic for generic function to compute widely applicable information criterion.

dmWAIC for computing widely applicable information criteria for Dirichlet-multinomial regression models.

Examples

# Compute log probability density.
ddirmult(x=c(33,115,95,359),
         p=c(0.075,0.201,0.175,0.549),
         theta=4.027,log=TRUE)

Description

Decodes Phred quality scores in Sanger format from symbols to numeric values.

Usage

decode_quality_scores(symbols)

Arguments

Value

A numeric vector of Phred quality scores.

Examples

decode_quality_scores(symbols="989!.C;F@\"")

Description

Generates detection plots for multiple groups.

Usage

detection(
  x,
  r = 1,
  b = 0.025,
  v = 1000,
  w = 1,
  f = 0.5,
  c = "lightskyblue",
  m = 3,
  ...
)

Arguments

Details

Produces a pie-chart-like detection plot with grouping structure. Each circle represents a group. Each sector represents a sample, and each sub-sector represents a replicate. Filled replicates represent detections. Groups are sorted alphabetically (or inherit factor level ordering) and arranged from left to right and top to bottom. Samples are sorted alphabetically and arranged in a clockwise orientation (from angle zero). Samples are sorted independently for each group. This plot design is specialized for visualizing binary detection data.

Value

No return value.

References

A manuscript describing this plot design is in preparation.

See Also

singular.detection for singular detection plots.

proportion for grouped proportion plots.

Examples

set.seed(1234)
n.groups<-6
n.samples<-6
n.replicates<-3
data<-list(g=rep(x=LETTERS[1:n.groups],each=n.samples),
           s=rep(x=letters[1:n.samples],times=n.groups),
           r=rep(x=n.replicates,times=n.groups*n.samples),
           d=sample(x=0:n.replicates,size=n.groups*n.samples,
                    replace=TRUE))
detection(x=data)

Description

Compute Bray-Curtis dissimilarity from proportional abundances.

Usage

dissimilarity(p1, p2)

Arguments

Details

Calculates Bray-Curtis dissimilarity from proportional abundances using the formula sum(abs⁠(p1-p2))/2⁠. This is equivalent to the Bray-Curtis dissimilarity formula in the vegdist function of the vegan package when both communities have the same total counts (e.g., rarefied counts). This function is primarily intended to provide a method to compute Bray-Curtis dissimilarity from the posterior predictions of Dirichlet-multinomial regression models, which generate predictions of proportional abundances.

The dimensions of p1 and p2 must match. If p1 and p2 are matrices, then each record represents paired replicates for the communities (e.g., predictions from the same posterior draw). The elements (if p1 and p2 are vectors) or fields (if p1 and p2 are matrices) represent dimensions (e.g., taxa or species). If p1 and p2 are vectors, then a numeric scalar is returned. If p1 and p2 are matrices, then a numeric vector is returned whose elements correspond to the paired records of p1 and p2.

Value

Numeric scalar or vector of Bray-Curtis dissimilarity values.

References

Bray JR, and Curtis JT. 1957. An ordination of the upland forest communities of southern Wisconsin. Ecological Monographs, 27(4): 325-349. DOI: 10.2307/1942268

Legendre P, and Legendre L. 2012. Numerical Ecology: Third Edition. Elsevier.

Odum EP. 1950. Bird populations of the Highlands (North Carolina) Plateau in relation to plant succession and avian invasion. Ecology, 31(4): 587-605. DOI: 10.2307/1931577

See Also

dmreg for fitting Dirichlet-multinomial regression models.

dmpredict for generating predictions from Dirichlet-multinomial regression models.

diversity for computing Hill diversity from proportional abundances.

richness for computing species richness from occupancy probabilities.

Examples

# Compute Bray-Curtis dissimilarity.
dissimilarity(p1=c(0.15,0.25,0.4,0.2),
              p2=c(0.25,0.35,0.1,0.3))

Description

Compute Hill diversity from proportional abundances.

Usage

diversity(p, alpha = 2)

Arguments

Details

Calculates Hill diversity from proportional abundances as defined in Hill (1973), which provides a unifying theory for ecological diversity indices. When alpha = 0, Hill diversity is equal to species richness. When alpha = 1, Hill diversity is equal to the exponentiated Shannon's entropy. When alpha = 2 (the default), Hill diversity is equal to the inverse of Simpson's index. For any value of alpha, the Hill diversity of a community with uniform proportional abundances is equal to species richness. Hill diversity represents the effective number of species.

Value

Numeric scalar or vector of Hill diversity values.

References

Hill MO. 1973. Diversity and evenness: A unifying notation and its consequences. Ecology, 54(2): 427-432. DOI: 10.2307/1934352

See Also

dissimilarity for computing Bray-Curtis dissimilarity from proportional abundances.

richness for computing species richness from occupancy probabilities.

Examples

# Compute Hill diversity.
diversity(p=c(0.15,0.25,0.4,0.2))

Description

Density function for a joint species distribution model.

Usage

djsdm(x, psi, log = FALSE)

Arguments

Details

Computes the probability density of a joint species distribution model. The probability of observing a community is calculated as the product of the probabilities of observing each species. Observations for each species are Bernoulli-distributed, and species-specific probability densities are computed with stats::dbinom.

Value

Numeric vector of probability densities.

References

Wilkinson DP, Golding N, Guillera‐Arroita G, Tingley R, and McCarthy MA. 2021. Defining and evaluating predictions of joint species distribution models. Methods in Ecology and Evolution, 12(3): 394-404. DOI: 10.1111/2041-210X.13518

See Also

stats::dbinom for density of the binomial distribution.

mlWAIC for computing widely applicable information criteria for joint species distribution models.

Examples

# Define species occurrence.
x<-c(1,0,0,1)

# Define occupancy probabilities.
psi<-c(0.886,0.391,0.139,0.991)

# Compute log probability density.
djsdm(x=x,psi=psi,log=TRUE)

Description

Generate predictions for Dirichlet-multinomial regression models.

Usage

dmpredict(X, H, fit, names)

Arguments

Details

Generates posterior predictions for Dirichlet-multinomial regression models fit with the dmreg function. Predictions can either include or omit hierarchical effects, depending on whether argument H is provided. Returns a list where each element contains a matrix of posterior predictions for the respective record of X. Field names for the element matrices can optionally be provided with the names argument.

Value

A list whose elements contain numeric matrices of posterior predictions. Within the list, one element is returned for each record of X. Element names are taken from the row names of X.

See Also

dmreg for fitting Dirichlet-multinomial regression models.

dmWAIC for computing widely applicable information criteria for Dirichlet-multinomial regression models.

Examples

# Define example data file path.
path<-system.file("extdata",
                  "example_regression_data.rds",
                  package="LocaTT",
                  mustWork=TRUE)

# Read in example regression data.
data<-readRDS(file=path)

# Predict with fitted Dirichlet-multinomial regression.
out<-dmpredict(X=data$X,fit=data$fit,names=colnames(data$Y))

Description

Fit a Bayesian Dirichlet-multinomial regression model. Both fixed and hierarchical effects are supported. Installation of the rstan package is required to use this function.

Usage

dmreg(
  Y,
  X,
  H,
  ones = TRUE,
  priors = c(B.mu = 0, B.sd = 1, theta.mu = 0, theta.sd = 1, sigma2.alpha = 0.01,
    sigma2.beta = 0.01),
  control = list(adapt_delta = 0.95, max_treedepth = 20),
  ...
)

Arguments

Details

Fits the Bayesian Dirichlet-multinomial regression model of Goodwin et al. (2022) using the rstan interface to Stan (Carpenter et al. 2017). A stanfit object of the fitted model is returned, which can be used with standard rstan functions to evaluate model convergence (e.g., posterior trace plots, R-hat convergence diagnostics, and effective sample sizes). The model formulation is identical to that of Goodwin et al. (2022), except that the hard sum-to-zero constraint on hierarchical effects was removed to preserve the prior marginal variance of the final element. Up to four hierarchical variables are supported.

For each observation, counts are distributed according to the Dirichlet-multinomial distribution with alpha parameters defined as the product of an expected proportions vector and an exponentiated precision parameter. The precision parameter controls the degree of overdispersion relative to the multinomial distribution. The softmax function normalizes linear predictor combinations into expected proportions. For the model to be identifiable, the regression coefficients of the final dimension are set to zero. By default, weakly informative priors are used on the regression coefficients (B), precision parameter (theta), and hierarchical variances (sigma2). See the supplement of Goodwin et al. (2022) for details.

Value

Returns a stanfit object of the fitted Bayesian Dirichlet-multinomial regression model.

References

Carpenter B, Gelman A, Hoffman MD, Lee D, Goodrich B, Betancourt M, Brubaker M, Guo J, Li P, and Riddell A. 2017. Stan: A probabilistic programming language. Journal of Statistical Software, 76: 1-32. DOI: 10.18637/jss.v076.i01

Goodwin KB, Hutchinson JD, and Gompert Z. 2022. Spatiotemporal and ontogenetic variation, microbial selection, and predicted Bd-inhibitory function in the skin-associated microbiome of a Rocky Mountain amphibian. Frontiers in Microbiology, 13: 1020329. DOI: 10.3389/fmicb.2022.1020329

Harrison JG, Calder WJ, Shastry V, and Buerkle CA. Dirichlet-multinomial modelling outperforms alternatives for analysis of microbiome and other ecological count data. Molecular Ecology Resources, 20(2): 481-497. DOI: 10.1111/1755-0998.13128

See Also

dmpredict for generating predictions from Dirichlet-multinomial regression models.

dmWAIC for computing widely applicable information criteria for Dirichlet-multinomial regression models.

Examples

# Define example data file path.
path<-system.file("extdata",
                  "example_regression_data.rds",
                  package="LocaTT",
                  mustWork=TRUE)

# Read in example regression data.
data<-readRDS(file=path)

# Fit Dirichlet-multinomial regression.
out<-dmreg(Y=data$Y,X=data$X,H=data$H)

Description

Density function for the multivariate logistic distribution.

Usage

dmvlogis(x, location, scale, R, log = FALSE)

Arguments

Details

Computes the probability density of the multivariate logistic distribution. The multivariate logistic distribution is constructed using a Gaussian copula with logistic marginals. The probability density is the product of the densities of the logistic marginals, which is further multiplied by the density of a Gaussian copula of the transformed standard uniform margins (i.e., probability integral transformation of the logistic marginals with stats::plogis).

Value

Numeric vector of probability densities.

References

Decani JS, and Stine RA. 1986. A note on deriving the information matrix for a logistic distribution. The American Statistician, 40(3): 220-222. DOI: 10.2307/2684541

Song P. 2000. Multivariate dispersion models generated from Gaussian copula. Scandinavian Journal of Statistics, 27(2): 305-320. DOI: 10.1111/1467-9469.00191

See Also

stats::dlogis for density of the logistic distribution.

dcopula for density of the Gaussian copula.

Examples

# Define logistic margins.
x<-c(0.055,-1.625,0.329,-5.765)

# Define location parameters.
location<-c(0.477,-0.998,-0.776,0.064)

# Define scale parameters.
scale<-c(0.574,1.314,0.460,1.393)

# Define correlation matrix.
R<-matrix(data=c(1.000,-0.80,0.64,-0.512,
                 -0.800,1.00,-0.80,0.640,
                 0.640,-0.80,1.00,-0.800,
                 -0.512,0.64,-0.80,1.000),
           ncol=4,byrow=TRUE)

# Compute log probability density.
dmvlogis(x=x,location=location,
         scale=scale,R=R,
         log=TRUE)

Description

Computes the widely applicable information criterion (WAIC) for Dirichlet-multinomial regression models. Serves as a wrapper for dmreg, dmpredict, ddirmult, and waic for convenient WAIC calculations. Installation of the rstan package is required to use this function.

Usage

dmWAIC(
  Y,
  X,
  H,
  ones = TRUE,
  method = 2,
  priors = c(B.mu = 0, B.sd = 1, theta.mu = 0, theta.sd = 1, sigma2.alpha = 0.01,
    sigma2.beta = 0.01),
  control = list(adapt_delta = 0.95, max_treedepth = 20),
  ...
)

Arguments

Details

For convenience, wraps the steps involved in WAIC calculations for Bayesian Dirichlet-multinomial regression models. Begins by fitting a Bayesian Dirichlet-multinomial regression model with the dmreg function, then generates resubstitution posterior predictions using the dmpredict function. The pointwise log-likelihood is calculated with the ddirmult function given the response matrix, posterior predictions, and precision parameter. WAIC is calculated from the pointwise log-likelihood using the waic function.

Value

Returns numeric scalar of the widely applicable information criterion.

References

Carpenter B, Gelman A, Hoffman MD, Lee D, Goodrich B, Betancourt M, Brubaker M, Guo J, Li P, and Riddell A. 2017. Stan: A probabilistic programming language. Journal of Statistical Software, 76: 1-32. DOI: 10.18637/jss.v076.i01

Gelman A, Hwang J, and Vehtari A. 2014. Understanding predictive information criteria for Bayesian models. Statistics and Computing, 24(6): 997-1016. DOI: 10.1007/s11222-013-9416-2

Goodwin KB, Hutchinson JD, and Gompert Z. 2022. Spatiotemporal and ontogenetic variation, microbial selection, and predicted Bd-inhibitory function in the skin-associated microbiome of a Rocky Mountain amphibian. Frontiers in Microbiology, 13: 1020329. DOI: 10.3389/fmicb.2022.1020329

Harrison JG, Calder WJ, Shastry V, and Buerkle CA. Dirichlet-multinomial modelling outperforms alternatives for analysis of microbiome and other ecological count data. Molecular Ecology Resources, 20(2): 481-497. DOI: 10.1111/1755-0998.13128

Watanabe S. 2010. Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning Research, 11(116): 3571-3594.

See Also

dmreg for fitting Dirichlet-multinomial regression models.

dmpredict for generating predictions from Dirichlet-multinomial regression models.

ddirmult for probability mass function of the Dirichlet-multinomial distribution.

waic for generic function to compute widely applicable information criterion.

Examples

# Define example data file path.
path<-system.file("extdata",
                  "example_regression_data.rds",
                  package="LocaTT",
                  mustWork=TRUE)

# Read in example regression data.
data<-readRDS(file=path)

# Compute WAIC for Dirichlet-multinomial regression.
out<-dmWAIC(Y=data$Y,X=data$X,H=data$H)

Description

Extracts each taxonomic level from a vector of taxonomic strings.

Usage

expand_taxonomies(
  taxonomies,
  levels = c("Domain", "Phylum", "Class", "Order", "Family", "Genus", "Species"),
  full_names = TRUE,
  delimiter = ";",
  ignore
)

Arguments

Value

Returns a data frame of extracted taxonomic levels. One record for each element of taxonomies, and one field for each element of levels. Field names are inherited from levels. If a taxonomic level is not present in a taxonomic string, then the respective cell in the returned data frame will contain NA.

See Also

get_taxonomic_level for extracting a taxonomic level from taxonomic strings.

get_consensus_taxonomy for generating a consensus taxonomy from taxonomic strings.

Examples

expand_taxonomies(taxonomies=
   c("Eukaryota;Chordata;Amphibia;Caudata;Ambystomatidae;Ambystoma;Ambystoma mavortium",
     "Eukaryota;Chordata;Amphibia;Anura;Bufonidae;Anaxyrus;Anaxyrus boreas",
     "Eukaryota;Chordata;Amphibia;Anura;Ranidae;Rana;Rana luteiventris"),
                       full_names=FALSE,
                       delimiter=";")

Description

Filters DNA sequences by minimum read count within a PCR replicate, minimum proportion within a PCR replicate, and number of detections across PCR replicates.

Usage

filter_sequences(
  input_files,
  samples,
  PCR_replicates,
  output_file,
  minimum_reads.PCR_replicate = 1,
  minimum_reads.sequence = 1,
  minimum_proportion.sequence = 0.005,
  binomial_test.enabled = TRUE,
  binomial_test.p.adjust.method = "none",
  binomial_test.alpha_level = 0.05,
  minimum_PCR_replicates = 2,
  delimiter.read_counts = ": ",
  delimiter.PCR_replicates = ", "
)

Arguments

Details

For each set of input polymerase chain reaction (PCR) replicate FASTA files associated with a sample, writes out DNA sequences which are detected across a minimum number of PCR replicates (minimum_PCR_replicates argument). Detection within a PCR replicate is defined as a sequence having at least a minimum read count and exceeding a minimum proportion of reads (minimum_reads.sequence and minimum_proportion.sequence arguments, respectively). When binomial_test.enabled = TRUE, a sequence must significantly exceed the minimum proportion within a PCR replicate at the provided alpha level (binomial_test.alpha_level argument) based on a one-sided binomial test (i.e., binomial_test with alternative = "greater"). Within a PCR replicate, p-values can be adjusted for multiple hypothesis testing by setting the binomial_test.p.adjust.method argument (see stats::p.adjust.methods and p.adjust in the stats package). PCR replicates which contain fewer than a minimum number of reads are discarded (minimum_reads.PCR_replicate argument) and do not contribute detections to any sequence.

DNA sequences in the input FASTA files are assumed to be summarized by frequency of occurrence, with each FASTA header line beginning with "Frequency: " and followed by the sequence's read count. Output FASTA files from truncate_and_merge_pairs have this format and can be used directly with this function. Each input FASTA file is assumed to contain the DNA sequence reads for a single PCR replicate for a single sample.

For pipeline calibration purposes, a data frame containing unfiltered DNA sequences with their read counts, proportions, and p-values in each PCR replicate is invisibly returned (see return value section). While the primary output of this function is the written CSV file of filtered sequences (described below), the invisibly returned data frame of unfiltered sequences can be helpful when calibrating or troubleshooting filtering parameters. To aid in troubleshooting filtering parameters, the data frame is invisibly returned even if the error "Filtering removed all sequences" is received.

For the primary output, this function writes a CSV file of filtered DNA sequences with the following field definitions:

• Sample: The sample name.

• Sequence: The DNA sequence.

• Detections_across_PCR_replicates: The number of PCR replicates the sequence was detected in.

• Read_count_by_PCR_replicate: The sequence's read count in each PCR replicate the sequence was detected in.

• Sequence_read_count: The sequence's total read count across the PCR replicates the sequence was detected in. Calculated as the sum of the read counts in the Read_count_by_PCR_replicate field.

• Sample_read_count: The sample's total read count across all sequences detected in the PCR replicates. Calculated as the sum of the read counts in Sequence_read_count field associated with the sample.

• Proportion_of_sample: The proportion of sample reads comprised by the sequence. Calculated by dividing the Sequence_read_count field by the Sample_read_count field. Equivalent to the weighted average of the sequence's proportion in each PCR replicate, with weights given by the proportion of the sample's total reads contained in each PCR replicate.

Value

Invisibly returns a data frame containing unfiltered DNA sequences with their read counts, proportions, and p-values in each PCR replicate. While the primary output of this function is the written CSV file of filtered sequences described in the details section, the invisibly returned data frame of unfiltered sequences can be helpful when calibrating or troubleshooting filtering parameters. To aid in troubleshooting filtering parameters, the data frame is invisibly returned even if the error "Filtering removed all sequences" is received. Field definitions for the invisibly returned data frame of unfiltered sequences are:

• Sample: The sample name.

• PCR_replicate: The PCR replicate identifier.

• Sequence: The DNA sequence.

• Read_count.sequence: The sequence's read count within the PCR replicate.

• Read_count.PCR_replicate: The number of reads in the PCR replicate.

• Proportion_of_PCR_replicate.observed: The proportion of reads in the PCR replicate comprised by the sequence.

• Proportion_of_PCR_replicate.null (Field only present if binomial_test.enabled = TRUE): The null hypothesis for a one-sided binomial test (inherited from the minimum_proportion.sequence argument). See the p.value field below.

• p.value (Field only present if binomial_test.enabled = TRUE): The p-value from a one-sided binomial test of whether the proportion of reads in the PCR replicate comprised by the sequence exceeds the null hypothesis (i.e., binomial_test with alternative = "greater").

• p.value.adjusted (Field only present if binomial_test.enabled = TRUE): The p-value from the one-sided binomial test adjusted for multiple comparisons within each PCR replicate for each sample. See the p.value_adjustment_method field below.

• p.value_adjustment_method (Field only present if binomial_test.enabled = TRUE): The p-value adjustment method (inherited from the binomial_test.p.adjust.method argument).

References

A manuscript describing these methods is in preparation.

See Also

binomial_test for performing vectorized one-sided binomial tests.

truncate_and_merge_pairs for truncating and merging read pairs prior to sequence filtering.

local_taxa_tool for performing geographically-conscious taxonomic assignment of filtered sequences.

Examples

# Get example FASTA files.
input_files<-system.file("extdata",
                         paste0(rep(x=paste0("S0",1:3),
                                    each=3),
                                "P0",1:3,".fasta"),
                         package="LocaTT",
                         mustWork=TRUE)

# Create path for temporary output file.
output_file<-tempfile(fileext=".csv")

# Specify samples.
samples<-rep(x=paste0("S0",1:3),each=3)

# Specify replicates.
PCR_replicates<-rep(x=paste0("P0",1:3),times=3)

# Filter sequences.
filter_sequences(input_files=input_files,
                 samples=samples,
                 PCR_replicates=PCR_replicates,
                 output_file=output_file)

Description

Formats reference databases from MIDORI or UNITE for use with the local_taxa_tool function.

Usage

format_reference_database(
  path_to_input_database,
  path_to_output_database,
  input_database_source = "MIDORI",
  path_to_taxonomy_edits = NA,
  path_to_sequence_edits = NA,
  path_to_taxa_subset_list = NA,
  makeblastdb_command = "makeblastdb",
  ...
)

Arguments

Value

No return value. Writes formatted BLAST database files.

See Also

local_taxa_tool for performing geographically-conscious taxonomic assignment.

adjust_taxonomies for adjusting a taxonomy system.

Examples

# Get path to example reference sequences FASTA file.
path_to_input_file<-system.file("extdata",
                                "example_reference_sequences.fasta",
                                 package="LocaTT",
                                 mustWork=TRUE)

# Create a temporary file path for the output reference database FASTA file.
path_to_output_file<-tempfile(fileext=".fasta")

# Format reference database.
format_reference_database(path_to_input_database=path_to_input_file,
                          path_to_output_database=path_to_output_file)

Description

Gets the consensus taxonomy from a vector of taxonomic strings.

Usage

get_consensus_taxonomy(taxonomies, full_names = TRUE, delimiter = ";")

Arguments

Value

A character string containing the taxonomy agreed upon by all input taxonomies. If the input taxonomies are not the same at any taxonomic level, then NA is returned.

See Also

get_taxonomic_level for extracting a taxonomic level from taxonomic strings.

expand_taxonomies for extracting each taxonomic level from a vector of taxonomic strings.

Examples

get_consensus_taxonomy(taxonomies=
   c("Eukaryota;Chordata;Amphibia;Caudata;Ambystomatidae;Ambystoma;Ambystoma_mavortium",
     "Eukaryota;Chordata;Amphibia;Anura;Bufonidae;Anaxyrus;Anaxyrus_boreas",
     "Eukaryota;Chordata;Amphibia;Anura;Ranidae;Rana;Rana_luteiventris"),
                       full_names=TRUE,
                       delimiter=";")

Description

Gets the specified taxonomic level from a vector of taxonomic strings.

Usage

get_taxonomic_level(taxonomies, level, full_names = TRUE, delimiter = ";")

Arguments

Value

A character vector containing the requested taxonomic level for each element of the input taxonomies.

See Also

expand_taxonomies for extracting each taxonomic level from a vector of taxonomic strings.

get_consensus_taxonomy for generating a consensus taxonomy from taxonomic strings.

Examples

get_taxonomic_level(taxonomies=
   c("Eukaryota;Chordata;Amphibia;Caudata;Ambystomatidae;Ambystoma;Ambystoma_mavortium",
     "Eukaryota;Chordata;Amphibia;Anura;Bufonidae;Anaxyrus;Anaxyrus_boreas",
     "Eukaryota;Chordata;Amphibia;Anura;Ranidae;Rana;Rana_luteiventris"),
   level=5,
   full_names=TRUE,
   delimiter=";")

Description

Formats taxonomies from IUCN Red List taxonomy.csv and common_names.csv files for use with the local_taxa_tool function.

Usage

get_taxonomies.IUCN(
  path_to_taxonomies,
  path_to_common_names,
  path_to_output_file,
  domain = "Eukaryota",
  path_to_taxonomy_edits = NA,
  ...
)

Arguments

Value

No return value. Writes an output CSV file with formatted taxonomies.

See Also

get_taxonomies.species_binomials for remotely fetching NCBI taxonomies from species binomials.

adjust_taxonomies for adjusting a taxonomy system.

Examples

# Get path to example taxonomy CSV file.
path_to_taxonomies<-system.file("extdata",
                                "example_taxonomy.csv",
                                package="LocaTT",
                                mustWork=TRUE)

# Get path to example common names CSV file.
path_to_common_names<-system.file("extdata",
                                  "example_common_names.csv",
                                  package="LocaTT",
                                  mustWork=TRUE)

# Create a temporary file path for the output CSV file.
path_to_output_file<-tempfile(fileext=".csv")

# Format common names and taxonomies.
get_taxonomies.IUCN(path_to_taxonomies=path_to_taxonomies,
                    path_to_common_names=path_to_common_names,
                    path_to_output_file=path_to_output_file)

Description

Remotely fetches taxonomies from the NCBI taxonomy database for a list of species binomials. Installation of the taxize package is required to use this function.

Usage

get_taxonomies.species_binomials(
  path_to_species_binomials,
  path_to_output_file,
  path_to_taxonomy_edits = NA,
  print_queries = TRUE,
  ...
)

Arguments

Value

No return value. Writes an output CSV file with added taxonomies. Species which could not be found in the NCBI taxonomy database appear in the top records of the output file.

See Also

get_taxonomies.IUCN for formatting taxonomies from the IUCN Red List.

adjust_taxonomies for adjusting a taxonomy system.

Examples

# Get path to example input species binomials CSV file.
path_to_species_binomials<-system.file("extdata",
                                       "example_species_binomials.csv",
                                       package="LocaTT",
                                       mustWork=TRUE)

# Create a temporary file path for the output CSV file.
path_to_output_file<-tempfile(fileext=".csv")

# Fetch taxonomies from species binomials.
get_taxonomies.species_binomials(path_to_species_binomials=path_to_species_binomials,
                                 path_to_output_file=path_to_output_file,
                                 print_queries=FALSE)

Description

Trims DNA sequences to an amplicon region using forward and reverse primer sequences. Ambiguous nucleotides in forward and reverse primers are supported.

Usage

isolate_amplicon(sequences, forward_primer, reverse_primer)

Arguments

Details

For each DNA sequence, nucleotides matching and preceding the forward primer are removed, and nucleotides matching and following the reverse complement of the reverse primer are removed. The reverse complement of the reverse primer is internally derived from the reverse primer using the reverse_complement function. Ambiguous nucleotides in primers (i.e., the forward and reverse primer arguments) are supported through the internal use of the substitute_wildcards function on the forward primer and the reverse complement of the reverse primer, and primer regions in DNA sequences are located using regular expressions. Trimming will fail for DNA sequences which contain ambiguous nucleotides in their primer regions (e.g., Ns), resulting in NAs for those sequences.

Value

A character vector of DNA sequences trimmed to the amplicon region. NAs are returned for DNA sequences which could not be trimmed, which occurs when either primer region is missing from the DNA sequence or when the forward primer region occurs after a region matching the reverse complement of the reverse primer.

Examples

isolate_amplicon(sequences=c("ACACAATCGTGTTTATATTAACTTCAAGAGTGGGCATAGG",
                             "CGTGACAATCATGTTTGTGATTCGTACAAAAGTGCGTCCT"),
                 forward_primer="AATCRTGTTT",
                 reverse_primer="CSCACTHTTG")

Description

Performs taxonomic assignment of DNA metabarcoding sequences while considering geographic location.

Usage

local_taxa_tool(
  path_to_query_sequences,
  path_to_BLAST_database,
  path_to_output_file,
  path_to_local_taxa_list = NA,
  include_missing = FALSE,
  blast_e_value = 1e-05,
  blast_max_target_seqs = 2000,
  blast_task = "megablast",
  full_names = FALSE,
  underscores = FALSE,
  separator = ", ",
  blastn_command = "blastn",
  ...
)

Arguments

Details

Sequences are BLASTed against a global reference database, and the tool suggests locally occurring species which are most closely related (by taxonomy) to any of the best-matching BLAST hits (by bit score). Optionally, local sister taxonomic groups without reference sequences can be added to the local taxa suggestions by setting the include_missing argument to TRUE. If a local taxa list is not provided, then local taxa suggestions will be disabled, but all best-matching BLAST hits will still be returned. Alternatively, a reference database containing just the sequences of local species can be used, and local taxa suggestions can be disabled to return all best BLAST matches of local species. The reference database should be formatted with the format_reference_database function, and the local taxa lists can be prepared using the get_taxonomies.species_binomials and get_taxonomies.IUCN functions. Output field definitions are:

• Sequence_name: The query sequence name.

• Sequence: The query sequence.

• Best_match_references: Species binomials of all best-matching BLAST hits (by bit score) from the reference database.

• Best_match_E_value: The E-value associated with the best-matching BLAST hits.

• Best_match_bit_score: The bit score associated with the best-matching BLAST hits.

• Best_match_query_cover.mean: The mean query cover of all best-matching BLAST hits.

• Best_match_query_cover.SD: The standard deviation of query cover of all best-matching BLAST hits.

• Best_match_PID.mean: The mean percent identity of all best-matching BLAST hits.

• Best_match_PID.SD: The standard deviation of percent identity of all best-matching BLAST hits.

• Local_taxa (Field only present if a path to a local taxa list is provided): The finest taxonomic unit(s) which include both any species of the best-matching BLAST hits and any local species. If the species of any of the best-matching BLAST hits are local, then the finest taxonomic unit(s) are at the species level.

• Local_species (Field only present if a path to a local taxa list is provided): Species binomials of all local species which belong to the taxonomic unit(s) in the Local_taxa field.

• Local_taxa.include_missing (Field only present if both a path to a local taxa list is provided and the include_missing argument is set to TRUE): Local sister taxonomic groups without reference sequences are added to the local taxa suggestions from the Local_taxa field.

• Local_species.include_missing (Field only present if both a path to a local taxa list is provided and include_missing argument is set to TRUE): Species binomials of all local species which belong to the taxonomic unit(s) in the Local_taxa.include_missing field.

Value

No return value. Writes an output CSV file with fields defined in the details section.

References

A manuscript describing this taxonomic assignment method is in preparation.

See Also

format_reference_database for formatting reference databases.

get_taxonomies.species_binomials and get_taxonomies.IUCN for creating local taxa lists.

adjust_taxonomies for adjusting a taxonomy system.

Examples

# Get path to example query sequences FASTA file.
path_to_query_sequences<-system.file("extdata",
                                     "example_query_sequences.fasta",
                                     package="LocaTT",
                                     mustWork=TRUE)

# Get path to example BLAST reference database FASTA file.
path_to_BLAST_database<-system.file("extdata",
                                    "example_blast_database.fasta",
                                    package="LocaTT",
                                    mustWork=TRUE)

# Get path to example local taxa list CSV file.
path_to_local_taxa_list<-system.file("extdata",
                                     "example_local_taxa_list.csv",
                                     package="LocaTT",
                                     mustWork=TRUE)

# Create a temporary file path for the output CSV file.
path_to_output_file<-tempfile(fileext=".csv")

# Run the local taxa tool.
local_taxa_tool(path_to_query_sequences=path_to_query_sequences,
                path_to_BLAST_database=path_to_BLAST_database,
                path_to_output_file=path_to_output_file,
                path_to_local_taxa_list=path_to_local_taxa_list,
                include_missing=TRUE,
                full_names=TRUE,
                underscores=TRUE)

Description

Merges forward and reverse DNA sequence reads.

Usage

merge_pairs(forward_reads, reverse_reads, minimum_overlap = 10)

Arguments

Details

For each pair of forward and reverse DNA sequence reads, the reverse complement of the reverse read is internally derived using the reverse_complement function, and the read pair is merged into a single sequence if an overlap of at least the minimum length is found between the end of the forward read and the start of the reverse complement of the reverse read. If an overlap of the minimum length is not found, then an NA is returned for the merged read pair.

Value

A character vector of merged DNA sequence read pairs. NAs are returned for read pairs which could not be merged, which occurs when an overlap of at least the minimum length is not found between the end of the forward read and the start of the reverse complement of the reverse read.

See Also

truncate_and_merge_pairs for truncating and merging forward and reverse DNA sequence reads.

Examples

merge_pairs(forward_reads=c("CCTTACGAATCCTGT","TTCTCCACCCGCGGATA","CGCCCGGAGTCCCTGTAGTA"),
            reverse_reads=c("GACAAACAGGATTCG","CAATATCCGCGGGTG","TACTACAGGGACTCC"))

Description

Extract regression coefficients from multivariate logistic regression models.

Usage

mlcoef(fit, probs = c(0.025, 0.25, 0.5, 0.75, 0.975), dimnames)

Arguments

Details

Extracts regression coefficient estimates from a multivariate logistic regression model fit using the mlreg function. Summarizes estimates by the quantiles of their posterior distributions, and returns summaries in a 3-dimensional array. The dimensions of the 3D array represent the posterior quantiles, the predictor variables, and the response variables, respectively. Values at probs = 0.025 and 0.975 comprise 95% credible intervals. Values at probs = 0.25 and 0.75 comprise 50% credible intervals, and values at probs = 0.5 represent point estimates.

Value

Numeric 3-dimensional array of regression coefficient posterior quantiles.

See Also

mlreg for fitting multivariate logistic regression models.

mlcor for extracting residual correlations from multivariate logistic regression models.

mlformat for formatting output of multivariate logistic regression models.

Examples

# Define example data file path.
path<-system.file("extdata",
                  "example_mvlogistic_data.rds",
                  package="LocaTT",
                  mustWork=TRUE)

# Read in example regression data.
data<-readRDS(file=path)

# Extract regression coefficients.
out<-mlcoef(fit=data$fit)

Description

Extract residual correlations from multivariate logistic regression models.

Usage

mlcor(fit, probs = c(0.025, 0.25, 0.5, 0.75, 0.975), dimnames)

Arguments

Details

Extracts residual correlation estimates from a multivariate logistic regression model fit using the mlreg function. Summarizes estimates by the quantiles of their posterior distributions, and returns summaries in a 3-dimensional array. The dimensions of the 3D array represent the posterior quantiles (dimension 1) and the response variables (both dimensions 2 and 3). Values at probs = 0.025 and 0.975 comprise 95% credible intervals. Values at probs = 0.25 and 0.75 comprise 50% credible intervals, and values at probs = 0.5 represent point estimates.

Value

Numeric 3-dimensional array of residual correlation posterior quantiles.

See Also

mlreg for fitting multivariate logistic regression models.

mlcoef for extracting regression coefficients from multivariate logistic regression models.

mlformat for formatting output of multivariate logistic regression models.

Examples

# Define example data file path.
path<-system.file("extdata",
                  "example_mvlogistic_data.rds",
                  package="LocaTT",
                  mustWork=TRUE)

# Read in example regression data.
data<-readRDS(file=path)

# Extract residual correlations.
out<-mlcor(fit=data$fit)

Description

Format output of multivariate logistic regression models.

Usage

mlformat(fit, mode = "B", ci = 0.95, digits = 3, names.x, names.y)

Arguments

Details

Formats output of a multivariate logistic regression model fit using the mlreg function. When mode = "B" (the default), returns regression coefficient estimates. When mode = "R", returns residual correlation estimates (if mlreg is fit with multivariate = TRUE). Summarizes parameters by the quantiles of their posterior distributions, with a point estimate at the 50th percentile (i.e., the posterior median). Lower and upper limits are defined by the credible interval argument. At the default ci = 0.95, returns 95% credible intervals. When a credible interval does not overlap zero, the point estimate is appended with an asterisk.

Value

Numeric matrix of posterior summaries.

See Also

mlreg for fitting multivariate logistic regression models.

mlcoef for extracting regression coefficients from multivariate logistic regression models.

mlcor for extracting residual correlations from multivariate logistic regression models.

Examples

# Define example data file path.
path<-system.file("extdata",
                  "example_mvlogistic_data.rds",
                  package="LocaTT",
                  mustWork=TRUE)

# Read in example regression data.
data<-readRDS(file=path)

# Retrieve fitted regression model.
fit<-data$fit

# Retrieve predictor matrix.
X<-data$X

# Retrieve response matrix.
Y<-data$Y

# Extract regression coefficients.
B<-mlformat(fit=fit,mode="B",
            names.x=colnames(X),
            names.y=colnames(Y))

# Display regression coefficients.
print(B,quote=FALSE,right=TRUE)

# Extract residual correlations.
R<-mlformat(fit=fit,mode="R",
            names.y=colnames(Y))

# Display residual correlations.
print(R,quote=FALSE,right=TRUE)

Description

Generate predictions for multivariate logistic regression models.

Usage

mlpredict(X, fit, names)

Arguments

Details

Generates posterior predictions for multivariate logistic regression models fit with the mlreg function. Returns a list where each element contains a matrix of posterior predictions for the respective record of X. Field names for the element matrices can optionally be provided with the names argument.

Value

A list whose elements contain numeric matrices of posterior predictions. Within the list, one element is returned for each record of X. Element names are taken from the row names of X.

See Also

mlreg for fitting multivariate logistic regression models.

mlformat for formatting output of multivariate logistic regression models.

mlWAIC for computing widely applicable information criteria for multivariate logistic regression models.

Examples

# Define example data file path.
path<-system.file("extdata",
                  "example_mvlogistic_data.rds",
                  package="LocaTT",
                  mustWork=TRUE)

# Read in example regression data.
data<-readRDS(file=path)

# Predict with fitted multivariate logistic regression.
out<-mlpredict(X=data$X,fit=data$fit,names=colnames(data$Y))

Description

Fit a multivariate logistic regression model. Installation of the rstan package is required to use this function.

Usage

mlreg(
  Y,
  X,
  multivariate = TRUE,
  priors = c(B.mu = 0, B.sd = 1, lkj = 1),
  iter = 20000,
  thin = 20,
  control = list(adapt_delta = 0.99, max_treedepth = 20, stepsize = 0.01),
  ...
)

Arguments

Details

Fits a multivariate logistic regression model using the rstan interface to Stan (Carpenter et al. 2017). The multivariate logistic regression follows that of Ovaskainen et al. 2010, where the Bernoulli marginals are reparameterized as truncated continuous latent variables (Albert & Chib 1993). The latent variables z receive a positive constraint when y = 1 and a negative constraint when y = 0, where z is a linear combination of predictors with correlated standard logistic errors. Equivalently, the latent variables follow a multivariate logistic distribution with scale parameters fixed at one (O'Brien & Dunson 2004), constructed in Stan as a Gaussian copula with logistic marginals (Song 2000). A stanfit object of the fitted model is returned, which can be used with standard rstan functions to evaluate model convergence (e.g., posterior trace plots, R-hat convergence diagnostics, and effective sample sizes). By default, weakly informative priors are used on the regression coefficients (B) and residual correlation matrix (R).

Value

Returns a stanfit object of the fitted multivariate logistic regression model.

References

Albert JH, and Chib S. 1993. Bayesian analysis of binary and polychotomous response data. Journal of the American Statistical Association, 88(422): 669-679.

Carpenter B, Gelman A, Hoffman MD, Lee D, Goodrich B, Betancourt M, Brubaker M, Guo J, Li P, and Riddell A. 2017. Stan: A probabilistic programming language. Journal of Statistical Software, 76: 1-32. DOI: 10.18637/jss.v076.i01

O'Brien SM, and Dunson DB. 2004. Bayesian multivariate logistic regression. Biometrics, 60: 739-746. DOI: 10.1111/j.0006-341X.2004.00224.x

Ovaskainen O, Hottola J, and Siitonen J. 2010. Modeling species co-occurrence by multivariate logistic regression generates new hypotheses on fungal interactions. Ecology, 91(9): 2514-2521. DOI: 10.1890/10-0173.1

Song P. 2000. Multivariate dispersion models generated from Gaussian copula. Scandinavian Journal of Statistics, 27(2): 305-320. DOI: 10.1111/1467-9469.00191

See Also

mlformat for formatting output of multivariate logistic regression models.

mlpredict for generating predictions from multivariate logistic regression models.

mlWAIC for computing widely applicable information criteria for multivariate logistic regression models.

Examples

# Define example data file path.
path<-system.file("extdata",
                  "example_mvlogistic_data.rds",
                  package="LocaTT",
                  mustWork=TRUE)

# Read in example regression data.
data<-readRDS(file=path)

# Fit multivariate logistic regression.
out<-mlreg(Y=data$Y,X=data$X)

Description

Computes the widely applicable information criterion (WAIC) for multivariate logistic regression models. Serves as a wrapper for mlreg, mlpredict, djsdm, and waic for convenient WAIC calculations. Installation of the rstan package is required to use this function.

Usage

mlWAIC(
  Y,
  X,
  multivariate = TRUE,
  method = 2,
  priors = c(B.mu = 0, B.sd = 1, lkj = 1),
  iter = 20000,
  thin = 20,
  control = list(adapt_delta = 0.99, max_treedepth = 20, stepsize = 0.01),
  ...
)

Arguments

Details

For convenience, wraps the steps involved in WAIC calculations for Bayesian multivariate logistic regression models. Begins by fitting a Bayesian multivariate logistic regression model with the mlreg function, then generates resubstitution posterior predictions using the mlpredict function. The pointwise log-likelihood is calculated with the djsdm function given the response matrix and posterior predictions. WAIC is calculated from the pointwise log-likelihood using the waic function. Because djsdm does not consider residual correlations in density calculations, species interactions do not contribute to WAIC (i.e., response dimensions are independent).

Value

Returns numeric scalar of the widely applicable information criterion.

References

Carpenter B, Gelman A, Hoffman MD, Lee D, Goodrich B, Betancourt M, Brubaker M, Guo J, Li P, and Riddell A. 2017. Stan: A probabilistic programming language. Journal of Statistical Software, 76: 1-32. DOI: 10.18637/jss.v076.i01

Gelman A, Hwang J, and Vehtari A. 2014. Understanding predictive information criteria for Bayesian models. Statistics and Computing, 24(6): 997-1016. DOI: 10.1007/s11222-013-9416-2

O'Brien SM, and Dunson DB. 2004. Bayesian multivariate logistic regression. Biometrics, 60: 739-746. DOI: 10.1111/j.0006-341X.2004.00224.x

Ovaskainen O, Hottola J, and Siitonen J. 2010. Modeling species co-occurrence by multivariate logistic regression generates new hypotheses on fungal interactions. Ecology, 91(9): 2514-2521. DOI: 10.1890/10-0173.1

Watanabe S. 2010. Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning Research, 11(116): 3571-3594.

See Also

mlreg for fitting multivariate logistic regression models.

mlpredict for generating predictions from multivariate logistic regression models.

djsdm for probability mass function of a joint species distribution model.

waic for generic function to compute widely applicable information criterion.

Examples

# Define example data file path.
path<-system.file("extdata",
                  "example_mvlogistic_data.rds",
                  package="LocaTT",
                  mustWork=TRUE)

# Read in example regression data.
data<-readRDS(file=path)

# Compute WAIC for multivariate logistic regression.
out<-mlWAIC(Y=data$Y,X=data$X)

Description

Normalizes a vector or each record of a matrix into a simplex.

Usage

normalize(x)

Arguments

Details

Returns a vector or matrix whose elements (if vector) or records (if matrix) are computed as ⁠x/⁠sum(x). This normalizes a vector or matrix record into a set of proportions which sum to one (i.e., a simplex). If a matrix is provided for the x argument, then normalization is performed independently for each record.

Value

A numeric vector or matrix whose elements (if vector) or records (if matrix) sum to one.

See Also

softmax for the softmax function.

Examples

# Normalize vector.
normalize(x=c(3,1,5,7))

Description

Generates proportion plots for multiple groups.

Usage

proportion(
  x,
  r = 1,
  b = 0.025,
  v = 1000,
  w = 1,
  f = 0.5,
  c = "lightskyblue",
  s = FALSE,
  a = TRUE,
  m = 3,
  ...
)

Arguments

Details

Produces a pie-chart-like proportion plot with grouping structure. Each circle represents a group, and each sector represents a sample. When a = TRUE (the default), then the proportion of each sector filled with color represents the within-sample proportional abundance. Groups are sorted alphabetically (or inherit factor level ordering) and arranged from left to right and top to bottom. When s = FALSE (the default), then samples are sorted alphabetically and arranged in a clockwise orientation (from angle zero). When s = TRUE, then samples are sorted by decreasing proportional abundance. Samples are sorted independently for each group. This plot design is specialized for visualizing proportional abundance data.

Value

No return value.

References

A manuscript describing this plot design is in preparation.

See Also

singular.proportion for singular proportion plots.

detection for grouped detection plots.

Examples

set.seed(1234)
n.groups<-6
n.samples<-6
data<-list(g=rep(x=LETTERS[1:n.groups],each=n.samples),
           s=rep(x=letters[1:n.samples],times=n.groups),
           p=stats::rbeta(n=n.groups*n.samples,
                          shape1=1,shape2=1))
proportion(x=data)

Description

Reads FASTA files. Supports the reading of FASTA files with sequences wrapping multiple lines.

Usage

read.fasta(file)

Arguments

Value

A data frame with fields for sequence names and sequences.

See Also

write.fasta for writing FASTA files.

read.fastq for reading FASTQ files.

write.fastq for writing FASTQ files.

Examples

# Get path to example FASTA file.
path_to_fasta_file<-system.file("extdata",
                                "example_query_sequences.fasta",
                                package="LocaTT",
                                mustWork=TRUE)

# Read the example FASTA file.
read.fasta(file=path_to_fasta_file)

Description

Reads FASTQ files. Does not support the reading of FASTQ files with sequences or quality scores wrapping multiple lines.

Usage

read.fastq(file)

Arguments

Value

A data frame with fields for sequence names, sequences, comments, and quality scores.

See Also

write.fastq for writing FASTQ files.

read.fasta for reading FASTA files.

write.fasta for writing FASTA files.

Examples

# Get path to example FASTQ file.
path_to_fastq_file<-system.file("extdata",
                                "example_query_sequences.fastq",
                                package="LocaTT",
                                mustWork=TRUE)

# Read the example FASTQ file.
read.fastq(file=path_to_fastq_file)

Description

Gets the reverse complement of a DNA sequence. Ambiguous nucleotides are supported.

Usage

reverse_complement(sequence)

Arguments

Value

A string of the reverse complement of the DNA sequence.

Examples

reverse_complement(sequence="TTCTCCASCCGCGGATHTTG")

Description

Compute species richness from occupancy probabilities.

Usage

richness(psi)

Arguments

Details

Calculates species richness from occupancy probabilities. Given a vector of species occupancy probabilities, computes the expected number of species as the sum of the probabilities. If given a matrix of species occupancy probabilities (where each record represents a community), computes the expected number of species as the row sums.

Value

Numeric scalar or vector of species richness values.

See Also

diversity for computing Hill diversity from proportional abundances.

dissimilarity for computing Bray-Curtis dissimilarity from proportional abundances.

Examples

# Compute species richness.
richness(psi=c(0.506,0.825,0.135,0.683))

Description

Draws sector polygon.

Usage

sector(s, e, r, v = 1000, ...)

Arguments

Details

Draws a sector polygon given a start angle, end angle, and circle radius. The sector is drawn about the origin (i.e., x = 0, y = 0). Intended for use with template to generate detection and proportion plots.

Value

No return value.

See Also

circle for plotting circle polygons.

Examples

template(l=1)
sector(s=0,e=45,r=1)

Description

Generates a detection plot for a singular group.

Usage

singular.detection(
  x,
  r = 1,
  b = 0.025,
  v = 1000,
  w = 1,
  f = 0.5,
  c = "lightskyblue",
  t = "",
  ...
)

Arguments

Details

Produces a pie-chart-like detection plot without grouping structure. Each sector represents a sample, and each sub-sector represents a replicate. Filled replicates represent detections. Samples are sorted alphabetically and arranged in a clockwise orientation (from angle zero). This plot design is specialized for visualizing binary detection data.

Value

No return value.

References

A manuscript describing this plot design is in preparation.

See Also

detection for grouped detection plots.

proportion for grouped proportion plots.

Examples

set.seed(1234)
n.samples<-6
n.replicates<-3
data<-list(s=letters[1:n.samples],
           r=rep(x=n.replicates,times=n.samples),
           d=sample(x=0:n.replicates,size=n.samples,
                    replace=TRUE))
singular.detection(x=data)

Description

Generates a proportion plot for a singular group.

Usage

singular.proportion(
  x,
  r = 1,
  b = 0.025,
  v = 1000,
  w = 1,
  f = 0.5,
  c = "lightskyblue",
  t = "",
  s = FALSE,
  a = TRUE,
  ...
)

Arguments

Details

Produces a pie-chart-like proportion plot without grouping structure. Each sector represents a sample. When a = TRUE (the default), then the proportion of each sector filled with color represents the within-sample proportional abundance. When s = FALSE (the default), then samples are sorted alphabetically and arranged in a clockwise orientation (from angle zero). When s = TRUE, then samples are sorted by decreasing proportional abundance. This plot design is specialized for visualizing proportional abundance data.

Value

No return value.

References

A manuscript describing this plot design is in preparation.

See Also

proportion for grouped proportion plots.

singular.detection for singular detection plots.

Examples

set.seed(1234)
n.samples<-6
data<-list(s=letters[1:n.samples],
           p=stats::rbeta(n=n.samples,
                          shape1=1,shape2=1))
singular.proportion(x=data)

Description

Applies the softmax to a vector or each record of a matrix.

Usage

softmax(x)

Arguments

Details

Returns a vector or matrix whose elements (if vector) or records (if matrix) are computed as exp⁠(x)/⁠sum(exp⁠(x))⁠. The softmax converts a vector or matrix record into a set of proportions which sum to one. If a matrix is provided for the x argument, then the softmax is applied independently for each record.

Value

A numeric vector or matrix whose elements (if vector) or records (if matrix) sum to one.

See Also

normalize for vector or matrix normalization.

Examples

# Perform softmax on vector.
softmax(x=c(-0.25,0.75,1.5,0))

Description

Substitutes wildcard characters in a DNA sequence with their associated nucleotides surrounded by square brackets. The output is useful for matching in regular expressions.

Usage

substitute_wildcards(sequence)

Arguments

Value

A string of the DNA sequence in which wildcard characters are replaced with their associated nucleotides surrounded by square brackets.

Examples

substitute_wildcards(sequence="CAADATCCGCGGSTGGAGAA")

Description

For each base pair position, summarizes read length, Phred quality score, and the cumulative probability that all bases were called correctly.

Usage

summarize_quality_scores(
  forward_files,
  reverse_files,
  n.total = 10000,
  n.each = ceiling(n.total/length(forward_files)),
  seed = NULL,
  FUN = mean,
  ...
)

Arguments

Details

For each combination of base pair position and read direction, calculates summary statistics of read length, Phred quality score, and the cumulative probability that all bases were called correctly. The cumulative probability is calculated from the first base pair up to the current position. Quality scores are assumed to be encoded in Sanger format. Read pairs are selected by randomly sampling up to n.each read pairs from each pair of input FASTQ files. By default, n.each is derived from n.total, and n.total will be ignored if n.each is provided. By default, mean is used to compute the summary statistics, but the user may provide another summary function instead (e.g., median). Functions which return multiple summary statistics are also supported (e.g., summary and quantile). Arguments in ... are passed to the summary function.

Value

Returns a data frame containing summary statistics of read length and quality score at each base pair position. The returned data frame contains the following fields:

• Direction: The read direction (i.e., "Forward" or "Reverse").

• Position: The base pair position.

• Length: The summary statistic(s) of read lengths. If FUN returns multiple summary statistics, then a matrix of the summary statistics will be stored in this field, which can be accessed with $Length.

• Score: The summary statistic(s) of Phred quality scores. If FUN returns multiple summary statistics, then a matrix of the summary statistics will be stored in this field, which can be accessed with $Score.

• Probability: The summary statistic(s) of the cumulative probability that all bases were called correctly. If FUN returns multiple summary statistics, then a matrix of the summary statistics will be stored in this field, which can be accessed with $Probability.

See Also

decode_quality_scores for decoding quality scores.

Examples

# Get example forward FASTQ files.
forward_files<-system.file("extdata",
                           paste0("S0",1:3,"F.fastq"),
                           package="LocaTT",
                           mustWork=TRUE)

# Get example reverse FASTQ files.
reverse_files<-system.file("extdata",
                           paste0("S0",1:3,"R.fastq"),
                           package="LocaTT",
                           mustWork=TRUE)

# Summarize quality scores.
summarize_quality_scores(forward_files,reverse_files)

Description

Initiates a blank template plot.

Usage

template(l, b = 0.025)

Arguments

Details

Initiates a blank template plot with limits l and buffer b about the origin (i.e., x = 0, y = 0). l is used for axis limits in both the negative and positive directions. b extends the limits beyond l by a fixed proportion (i.e., l * (1 + b)). Intended for use with circle and sector.

Value

No return value.

See Also

circle for plotting circle polygons.

sector for plotting sector polygons.

Examples

template(l=1)
circle(r=1)

Description

Trims a target nucleotide sequence from the front or back of DNA sequences. Ambiguous nucleotides in the target nucleotide sequence are supported.

Usage

trim_sequences(
  sequences,
  target,
  anchor = "start",
  fixed = TRUE,
  required = TRUE,
  quality_scores
)

Arguments

Details

For each DNA sequence, the target nucleotide sequence is searched for at either the front or back of the DNA sequence, depending on the value of the anchor argument. If the target nucleotide sequence is found, then it is removed from the DNA sequence. If the required argument is set to TRUE, then DNA sequences in which the target nucleotide sequence was not found will be returned as NAs. If the required argument is set to FALSE, then untrimmed DNA sequences will be returned along with DNA sequences for which trimming was successful. Ambiguous nucleotides in the target nucleotide sequence are supported through the internal use of the substitute_wildcards function on the target nucleotide sequence, and a regular expression with a leading or ending anchor is used to search for the target nucleotide sequence in the DNA sequences. If the fixed argument is set to FALSE, then any number of characters are allowed between the start or end of the DNA sequences and the target nucleotide sequence. Trimming will fail for DNA sequences which contain ambiguous nucleotides (e.g., Ns) in their target nucleotide sequence region, resulting in NAs for those sequences if the required argument is set to TRUE.

Value

If quality scores are not provided, then a character vector of trimmed DNA sequences is returned. If quality scores are provided, then a list containing two elements is returned. The first element is a character vector of trimmed DNA sequences, and the second element is a character vector of quality scores which have been trimmed to their corresponding trimmed DNA sequences.

Examples

trim_sequences(sequences=c("ATATAGCGCG","TGCATATACG","ATCTATCACCGC"),
               target="ATMTA",
               anchor="start",
               fixed=TRUE,
               required=TRUE,
               quality_scores=c("989!.C;F@\"","A((#-#;,2F","HD8I/+67=1>?"))

Description

Removes DNA read pairs containing ambiguous nucleotides, truncates reads by length and quality score, and merges forward and reverse reads.

Usage

truncate_and_merge_pairs(
  forward_files,
  reverse_files,
  output_files,
  truncation_length = NA,
  threshold.quality_score = 3,
  threshold.probability = 0.5,
  minimum_overlap = 10,
  cores = 1,
  progress = FALSE
)

Arguments

Details

For each pair of input FASTQ files, removes DNA read pairs containing ambiguous nucleotides, truncates reads by length, quality score threshold, and probability threshold (in that order), and then merges forward and reverse reads. Merged reads are summarized by frequency of occurrence and written to a FASTA file. See contains_wildcards, truncate_sequences.length, truncate_sequences.quality_score, truncate_sequences.probability, and merge_pairs for methods. Quality scores are assumed to be encoded in Sanger format. Forward and reverse reads can be truncated by different thresholds (see truncation_length, threshold.quality_score, and threshold.probability arguments).

Multicore parallel processing is supported on Mac and Linux operating systems (not available on Windows). When cores > 1 (parallel processing enabled), warnings and errors are printed to the console in addition to being invisibly returned as a list (see the return value section), and errors produced while processing a pair of FASTQ files will not interrupt the processing of other FASTQ file pairs. When cores = 1, FASTQ file pairs are processed sequentially on a single core, and errors will prevent the processing of subsequent FASTQ file pairs (but warnings will not).

Value

If cores = 1, then no return value. Writes a FASTA file for each pair of input FASTQ files with DNA sequence counts stored in the header lines. If cores > 1, then also invisibly returns a list where each element contains warning or error messages associated with processing each pair of input FASTQ files. A NULL value in the returned list means that no warnings or errors were generated from processing the respective pair of FASTQ files.

References

A manuscript describing these methods is in preparation.

See Also

contains_wildcards for detecting ambiguous nucleotides in DNA sequences.

truncate_sequences.length for truncating DNA sequences to a specified length.

truncate_sequences.quality_score for truncating DNA sequences by Phred quality score.

truncate_sequences.probability for truncating DNA sequences by cumulative probability that all bases were called correctly.

merge_pairs for merging forward and reverse DNA sequence reads.

filter_sequences for filtering merged read pairs by PCR replicate.

Examples

# Get example forward FASTQ files.
forward_files<-system.file("extdata",
                           paste0("S0",1:3,"F.fastq"),
                           package="LocaTT",
                           mustWork=TRUE)

# Get example reverse FASTQ files.
reverse_files<-system.file("extdata",
                           paste0("S0",1:3,"R.fastq"),
                           package="LocaTT",
                           mustWork=TRUE)

# Create paths for temporary output files.
output_files<-tempfile(pattern=paste0("O",1:3),fileext=".fasta")

# Truncate and merge pairs.
truncate_and_merge_pairs(forward_files=forward_files,
                         reverse_files=reverse_files,
                         output_files=output_files)

Description

Truncates DNA sequences to a specified length.

Usage

truncate_sequences.length(sequences, length, quality_scores)

Arguments

Value

If quality scores are not provided, then a character vector of truncated DNA sequences is returned. If quality scores are provided, then a list containing two elements is returned. The first element is a character vector of truncated DNA sequences, and the second element is a character vector of quality scores which have been truncated to their corresponding truncated DNA sequences.

See Also

truncate_sequences.quality_score for truncating DNA sequences by Phred quality score.

truncate_sequences.probability for truncating DNA sequences by cumulative probability that all bases were called correctly.

truncate_and_merge_pairs for truncating and merging forward and reverse DNA sequence reads.

Examples

truncate_sequences.length(sequences=c("ATATAGCGCG","TGCCGATATA","ATCTATCACCGC"),
                          length=5,
                          quality_scores=c("989!.C;F@\"","A((#-#;,2F","HD8I/+67=1>?"))

Description

Calculates the cumulative probability that all bases were called correctly along each DNA sequence and truncates the DNA sequence immediately prior to the first occurrence of a probability being equal to or less than a specified value.

Usage

truncate_sequences.probability(sequences, quality_scores, threshold = 0.5)

Arguments

Value

A list containing two elements. The first element is a character vector of truncated DNA sequences, and the second element is a character vector of quality scores which have been truncated to their corresponding truncated DNA sequences.

See Also

truncate_sequences.length for truncating DNA sequences to a specified length.

truncate_sequences.quality_score for truncating DNA sequences by Phred quality score.

truncate_and_merge_pairs for truncating and merging forward and reverse DNA sequence reads.

Examples

truncate_sequences.probability(sequences=c("ATATAGCGCG","TGCCGATATA","ATCTATCACCGC"),
                               quality_scores=c("989!.C;F@\"","A((#-#;,2F","HD8I/+67=1>?"),
                               threshold=0.5)

Description

Truncates DNA sequences immediately prior to the first occurrence of a Phred quality score being equal to or less than a specified value.

Usage

truncate_sequences.quality_score(sequences, quality_scores, threshold = 3)

Arguments

Value

A list containing two elements. The first element is a character vector of truncated DNA sequences, and the second element is a character vector of quality scores which have been truncated to their corresponding truncated DNA sequences.

See Also

truncate_sequences.length for truncating DNA sequences to a specified length.

truncate_sequences.probability for truncating DNA sequences by cumulative probability that all bases were called correctly.

truncate_and_merge_pairs for truncating and merging forward and reverse DNA sequence reads.

Examples

truncate_sequences.quality_score(sequences=c("ATATAGCGCG","TGCCGATATA","ATCTATCACCGC"),
                                 quality_scores=c("989!.C;F@\"","A((#-#;,2F","HD8I/+67=1>?"),
                                 threshold=3)

Description

Generic function calculating widely applicable information criterion (WAIC) from the pointwise log-likelihood.

Usage

waic(loglik, method = 2)

Arguments

Details

Given the pointwise log-likelihood, calculates WAIC (Watanabe 2010) using the formulas described in Gelman et al. (2014). The expected log pointwise predictive density (elppd) is estimated as the log pointwise predictive density (lppd) adjusted by a bias correction (either pWAIC1 or pWAIC2). To reflect the deviance scale, WAIC is defined as the elppd times negative two. As recommended by Gelman et al. (2014), pWAIC2 is used as the default bias correction (method = 2). See Gelman et al. (2014) for details.

Value

Returns numeric scalar of the widely applicable information criterion.

References

Gelman A, Hwang J, and Vehtari A. 2014. Understanding predictive information criteria for Bayesian models. Statistics and Computing, 24(6): 997-1016. DOI: 10.1007/s11222-013-9416-2

Watanabe S. 2010. Asymptotic equivalence of Bayes cross validation and widely applicable information criterion in singular learning theory. Journal of Machine Learning Research, 11(116): 3571-3594.

See Also

dmWAIC for computing widely applicable information criteria for Dirichlet-multinomial regression models.

mlWAIC for computing widely applicable information criteria for multivariate logistic regression models.

Examples

# Define example data file path.
path<-system.file("extdata",
                  "example_regression_data.rds",
                  package="LocaTT",
                  mustWork=TRUE)

# Read in example regression data.
data<-readRDS(file=path)

# Compute WAIC from pointwise log-likelihood.
out<-waic(loglik=data$loglik)

Description

Writes FASTA files.

Usage

write.fasta(names, sequences, file)

Arguments

Value

No return value. Writes a FASTA file.

See Also

read.fasta for reading FASTA files.

write.fastq for writing FASTQ files.

read.fastq for reading FASTQ files.

Examples

# Get path to example sequences CSV file.
path_to_CSV_file<-system.file("extdata",
                              "example_query_sequences.csv",
                              package="LocaTT",
                              mustWork=TRUE)

# Read the example sequences CSV file.
df<-read.csv(file=path_to_CSV_file,stringsAsFactors=FALSE)

# Create a temporary file path for the FASTA file to write.
path_to_FASTA_file<-tempfile(fileext=".fasta")

# Write the example sequences as a FASTA file.
write.fasta(names=df$Name,
            sequences=df$Sequence,
            file=path_to_FASTA_file)

Description

Writes FASTQ files.

Usage

write.fastq(names, sequences, quality_scores, file, comments)

Arguments

Value

No return value. Writes a FASTQ file.

See Also

read.fastq for reading FASTQ files.

write.fasta for writing FASTA files.

read.fasta for reading FASTA files.

Examples

# Get path to example sequences CSV file.
path_to_CSV_file<-system.file("extdata",
                              "example_query_sequences.csv",
                              package="LocaTT",
                              mustWork=TRUE)

# Read the example sequences CSV file.
df<-read.csv(file=path_to_CSV_file,stringsAsFactors=FALSE)

# Create a temporary file path for the FASTQ file to write.
path_to_FASTQ_file<-tempfile(fileext=".fastq")

# Write the example sequences as a FASTQ file.
write.fastq(names=df$Name,
            sequences=df$Sequence,
            quality_scores=df$Quality_score,
            file=path_to_FASTQ_file,
            comments=df$Comment)