indexing dartR objects correctly...

Description

indexing dartR objects correctly...

Usage

## S4 method for signature 'dartR,ANY,ANY,ANY'
x[i, j, ..., pop = NULL, treatOther = TRUE, quiet = TRUE, drop = FALSE]

Arguments

adjust cbind for dartR

Description

cbind is a bit lazy and does not take care for the metadata (so data in the other slot is lost). You can get most of the loci metadata back using gl.compliance.check.

Usage

## S3 method for class 'dartR'
cbind(...)

Arguments

Value

A genlight object

Examples

t1 <- platypus.gl
class(t1) <- "dartR"
t2 <- cbind(t1[,1:10],t1[,11:20])

Estimates expected Heterozygosity

Description

Estimates expected Heterozygosity

Usage

gl.He(gl)

Arguments

Value

A simple vector whit Ho for each loci

Author(s)

Bernd Gruber (Post to https://groups.google.com/d/forum/dartr)

Estimates observed Heterozygosity

Description

Estimates observed Heterozygosity

Usage

gl.Ho(gl)

Arguments

Value

A simple vector whit Ho for each loci

Author(s)

Bernd Gruber (bugs? Post to https://groups.google.com/d/forum/dartr)

Adds metadata into a genlight object

Description

This function adds the metadata information to the slot ind.metrics and populates population and coordinates information slots if the they are found in the metadata.

Usage

gl.add.indmetrics(x, ind.metafile, verbose = NULL)

Arguments

Details

The ind.metadata file needs to have very specific headings. First a column with a heading named 'id'. Here the ids must match the ids in the genlight object, e.g. indNames(your_genlight). The following column headings are optional:

• 'pop' - specifies the population membership of each individual.

• 'lat' - latitude coordinates (in decimal degrees WGS1984 format).

• 'lon' - longitude coordinates (in decimal degrees WGS1984 format).

Additional columns with individual metadata can be imported (e.g. age, sex, etc).

Value

A genlight object with metadata information for each individual.

Author(s)

Custodian: Luis Mijangos – Post to https://groups.google.com/d/forum/dartr

Examples

dartfile <- system.file('extdata','testset_SNPs_2Row.csv', package='dartR.data')
metadata <- system.file('extdata','testset_metadata.csv', package='dartR.data')
gl <- gl.read.dart(dartfile, probar=TRUE)
gl <- gl.add.indmetrics(gl, ind.metafile = metadata)

Calculates allele frequency of the first and second allele for each locus A very simple function to report allele frequencies

Description

Calculates allele frequency of the first and second allele for each locus A very simple function to report allele frequencies

Usage

gl.alf(x)

Arguments

Value

A simple data.frame with ref (reference allele), alt (alternate allele).

Author(s)

Bernd Gruber (bugs? Post to https://groups.google.com/d/forum/dartr)

See Also

Other utilities: utils.check.datatype(), utils.collapse.matrix(), utils.dart2genlight(), utils.dist.binary(), utils.flag.start(), utils.hamming(), utils.het.pop(), utils.impute, utils.is.fixed(), utils.jackknife(), utils.n.var.invariant(), utils.plot.save(), utils.read.fasta(), utils.read.ped(), utils.recalc.avgpic(), utils.recalc.callrate(), utils.recalc.freqhets(), utils.recalc.freqhomref(), utils.recalc.freqhomsnp(), utils.recalc.maf(), utils.reset.flags(), utils.transpose(), utils.vcfr2genlight.polyploid()

Examples

#for the first 10 loci only
#Deprecated:
gl.alf(possums.gl[,1:10])
barplot(t(as.matrix(gl.alf(possums.gl[,1:10]))))
#Current:
gl.allele.freq(possums.gl[,1:10],simple=TRUE)
barplot(t(as.matrix(gl.allele.freq(possums.gl[,1:10],simple=TRUE))))

Generates percentage allele frequencies by locus and population

Description

This is a support script, to take SNP data or SilicoDArT presence/absence data grouped into populations in a genlight object {adegenet} and generate a table of allele frequencies for each population and locus

Usage

gl.allele.freq(x, percent = FALSE, by = "pop", simple = FALSE, verbose = NULL)

Arguments

Value

A matrix with allele (SNP data) or presence/absence frequencies (Tag P/A data) broken down by population and locus

Author(s)

Custodian: Arthur Georges (Post to https://groups.google.com/d/forum/dartr)

See Also

Other unmatched report: gl.report.basics(), gl.report.diversity(), gl.report.excess.het(), gl.report.heterozygosity(), gl.report.polyploid_heterozygosity()

Examples

gl.allele.freq(testset.gl,percent=FALSE,by='pop')
gl.allele.freq(testset.gl,percent=FALSE,by="loc")
gl.allele.freq(testset.gl,percent=FALSE,by="popxloc")
gl.allele.freq(testset.gl,simple=TRUE)

Performs AMOVA using genlight data

Description

This script performs an AMOVA based on the genetic distance matrix from stamppNeisD() [package StAMPP] using the amova() function from the package PEGAS for exploring within and between population variation. For detailed information use their help pages: ?pegas::amova, ?StAMPP::stamppAmova. Be aware due to a conflict of the amova functions from various packages I had to 'hack' StAMPP::stamppAmova to avoid a namespace conflict.

Usage

gl.amova(x, distance = NULL, permutations = 100, verbose = NULL)

Arguments

Value

An object of class 'amova' which is a list with a table of sums of square deviations (SSD), mean square deviations (MSD), and the number of degrees of freedom, and a vector of variance components.

Author(s)

Bernd Gruber (bugs? Post to https://groups.google.com/d/forum/dartr)

Examples

#permutations should be higher, here set to 1 because of speed
out <- gl.amova(bandicoot.gl, permutations=1)

Checks the current global verbosity

Description

The verbosity can be set in one of two ways – (a) explicitly by the user by passing a value using the parameter verbose in a function, or (b) by setting the verbosity globally as part of the r environment (gl.set.verbosity).

Usage

gl.check.verbosity(x = NULL)

Arguments

Value

The verbosity, in variable verbose

Author(s)

Bernd Gruber (Post to https://groups.google.com/d/forum/dartr)

See Also

Other environment: gl.check.wd(), gl.print.history(), gl.set.wd(), theme_dartR()

Examples

gl.check.verbosity()

Checks the global working directory

Description

The working directory can be set in one of two ways – (a) explicitly by the user by passing a value using the parameter plot.dir in a function, or (b) by setting the working directory globally as part of the r environment (gl.setwd). The default is in acccordance to CRAN set to tempdir().

Usage

gl.check.wd(wd = NULL, verbose = NULL)

Arguments

Value

the working directory

Author(s)

Custodian: Bernd Gruber (Post to https://groups.google.com/d/forum/dartr)

See Also

Other environment: gl.check.verbosity(), gl.print.history(), gl.set.wd(), theme_dartR()

Examples

gl.check.wd()

Returns a list of colors for use in plots

Description

Creates a vector of colors in hex (e.g. "#3B9AB2" "#78B7C5") based on user selected category (parameter type).

• "2" [two colors]

• "2c" [two colors contrast]

• "3" [three colors]

• 4" [four colors]

• "pal" [need to be specify the palette type and the number of colors]

A palette of colors can be specified via "div" [divergent], "dis" [discrete], "con" [convergent], "vir" [viridis]. Be aware a palette needs the number of colors specified as well. It returns a function and therefore the number of colors needs to be a part of the function call. Check the examples to see how this works.

Usage

gl.colors(type = 2, verbose = NULL)

Arguments

Value

returns colors as a vector

Author(s)

Custodian: Bernd Gruber – Post to https://groups.google.com/d/forum/dartr

See Also

Other graphics: gl.map.interactive(), gl.plot.heatmap(), gl.report.ld.map(), gl.select.colors(), gl.select.shapes(), gl.smearplot(), gl.tree.nj()

Examples

gl.colors(2)
gl.colors("2")
gl.colors("2c")
#five discrete colors
gl.colors(type="dis")(5)
#seven divergent colors
gl.colors("div")(7)

Checks a genlight object to see if it complies with dartR expectations and amends it to comply if necessary @family environment

Description

This function will check to see that the genlight object conforms to expectation in regard to dartR requirements (see details), and if it does not, will rectify it.

Usage

gl.compliance.check(x, verbose = NULL)

Arguments

Details

A genlight object used by dartR has a number of requirements that allow functions within the package to operate correctly. The genlight object comprises:

1. The SNP genotypes or Tag Presence/Absence data (SilicoDArT);

2. An associated dataframe (gl@other$loc.metrics) containing the locus metrics (e.g. Call Rate, Repeatability, etc);

3. An associated dataframe (gl@other$ind.metrics) containing the individual/sample metrics (e.g. sex, latitude (=lat), longitude(=lon), etc);

4. A specimen identity field (indNames(gl)) with the unique labels applied to each individual/sample;

5. A population assignment (popNames) for each individual/specimen;

6. Flags that indicate whether or not calculable locus metrics have been updated.

Value

A genlight object that conforms to the expectations of dartR

Author(s)

Custodian: Luis Mijangos – Post to https://groups.google.com/d/forum/dartr

Examples

x <- gl.compliance.check(testset.gl)
x <- gl.compliance.check(testset.gs)

Defines a new population in a genlight object for specified individuals

Description

The script reassigns existing individuals to a new population and removes their existing population assignment. The script returns a genlight object with the new population assignment.

Usage

gl.define.pop(x, ind.list, new, verbose = NULL)

Arguments

Value

A genlight object with the redefined population structure.

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

Other data manipulation: gl.drop.ind(), gl.drop.loc(), gl.drop.pop(), gl.edit.recode.pop(), gl.impute(), gl.join(), gl.keep.ind(), gl.keep.loc(), gl.keep.pop(), gl.make.recode.ind(), gl.merge.pop(), gl.reassign.pop(), gl.recode.ind(), gl.recode.pop(), gl.rename.pop(), gl.sample(), gl.sim.genotypes(), gl.sort(), gl.subsample.ind(), gl.subsample.loc()

Examples

popNames(testset.gl)
gl <- gl.define.pop(testset.gl, ind.list=c('AA019073','AA004859'), 
new='newguys')
popNames(gl)
indNames(gl)[pop(gl)=='newguys']

Provides descriptive stats and plots to diagnose potential problems with Hardy-Weinberg proportions @family matched report

Description

Different causes may be responsible for lack of Hardy-Weinberg proportions. This function helps diagnose potential problems.

Usage

gl.diagnostics.hwe(
  x,
  alpha_val = 0.05,
  bins = 20,
  stdErr = TRUE,
  colors.hist = gl.colors(2),
  colors.barplot = gl.colors("2c"),
  plot.theme = theme_dartR(),
  n.cores = "auto",
  plot.file = NULL,
  plot.dir = NULL,
  verbose = NULL
)

Arguments

Details

This function initially runs gl.report.hwe and reports the ternary plots. The remaining outputs follow the recommendations from Waples (2015) paper and De Meeûs 2018. These include:

1. A histogram with the distribution of p-values of the HWE tests. The distribution should be roughly uniform across equal-sized bins.

2. A bar plot with observed and expected (null expectation) number of significant HWE tests for the same locus in multiple populations (that is, the x-axis shows whether a locus results significant in 1, 2, ..., n populations. The y axis is the count of these occurrences. The zero value on x-axis shows the number of non-significant tests). If HWE tests are significant by chance alone, observed and expected number of HWE tests should have roughly a similar distribution.

3. A scatter plot with a linear regression between Fst and Fis, averaged across subpopulations. De Meeûs 2018 suggests that in the case of Null alleles, a strong positive relationship is expected (together with the Fis standard error much larger than the Fst standard error, see below). Note, this is not the scatter plot that Waples 2015 presents in his paper. In the lower right corner of the plot, the Pearson correlation coefficient is reported.

4. The Fis and Fst (averaged over loci and subpopulations) standard errors are also printed on screen and reported in the returned list (if stdErr=TRUE). These are computed with the Jackknife method over loci (See De Meeûs 2007 for details on how this is computed) and it may take some time for these computations to complete. De Meeûs 2018 suggests that under a global significant heterozygosity deficit: - if the correlation between Fis and Fst is strongly positive, and StdErrFis >> StdErrFst, Null alleles are likely to be the cause. - if the correlation between Fis and Fst is ~0 or mildly positive, and StdErrFis > StdErrFst, Wahlund may be the cause. - if the correlation between Fis and Fst is ~0, and StdErrFis ~ StdErrFst, selfing or sib mating could to be the cause. It is important to realise that these statistics only suggest a pattern (pointers). Their absence is not conclusive evidence of the absence of the problem, as their presence does not confirm the cause of the problem.

5. A table where the number of observed and expected significant HWE tests are reported by each population, indicating whether these are due to heterozygosity excess or deficiency. These can be used to have a clue of potential problems (e.g. deficiency might be due to a Wahlund effect, presence of null alleles or non-random sampling; excess might be due to sex linkage or different selection between sexes, demographic changes or small Ne. See Table 1 in Wapples 2015). The last two columns of the table generated by this function report chisquare values and their associated p-values. Chisquare is computed following Fisher's procedure for a global test (Fisher 1970). This basically tests whether there is at least one test that is truly significant in the series of tests conducted (De Meeûs et al 2009).

Value

A list with the table with the summary of the HWE tests and (if stdErr=TRUE) a named vector with the StdErrFis and StdErrFst.

Author(s)

Custodian: Carlo Pacioni – Post to https://groups.google.com/d/forum/dartr

References

• de Meeûs, T., McCoy, K.D., Prugnolle, F., Chevillon, C., Durand, P., Hurtrez-Boussès, S., Renaud, F., 2007. Population genetics and molecular epidemiology or how to “débusquer la bête”. Infection, Genetics and Evolution 7, 308-332.

• De Meeûs, T., Guégan, J.-F., Teriokhin, A.T., 2009. MultiTest V.1.2, a program to binomially combine independent tests and performance comparison with other related methods on proportional data. BMC Bioinformatics 10, 443-443.

• De Meeûs, T., 2018. Revisiting FIS, FST, Wahlund Effects, and Null Alleles. Journal of Heredity 109, 446-456.

• Fisher, R., 1970. Statistical methods for research workers Edinburgh: Oliver and Boyd.

• Waples, R. S. (2015). Testing for Hardy–Weinberg proportions: have we lost the plot?. Journal of heredity, 106(1), 1-19.

See Also

gl.report.hwe

Examples

require("dartR.data")
res <- gl.diagnostics.hwe(x = gl.filter.allna(platypus.gl[,1:50]), 
stdErr=FALSE, n.cores=1)

Calculates a distance matrix for individuals defined in a genlight object

Description

Calculates various distances between individuals based on allele frequencies or presence-absence data

Usage

gl.dist.ind(
  x,
  method = NULL,
  scale = FALSE,
  swap = FALSE,
  type = "dist",
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = NULL,
  plot.file = NULL,
  plot.dir = NULL,
  verbose = NULL
)

Arguments

Details

The distance measure for SNP genotypes can be one of:

• Euclidean Distance [method = "Euclidean"]

• Scaled Euclidean Distance [method='Euclidean", scale=TRUE]

• Simple Mismatch Distance [method="Simple"]

• Absolute Mismatch Distance [method="Absolute"]

• Czekanowski (Manhattan) Distance [method="Manhattan"]

The distance measure for Sequence Tag Presence/Absence data (binary) can be one of:

• Euclidean Distance [method = "Euclidean"]

• Scaled Euclidean Distance [method='Euclidean", scale=TRUE]

• Simple Matching Distance [method="Simple"]

• Jaccard Distance [method="Jaccard"]

• Bray-Curtis Distance [method="Bray-Curtis"]

Refer to the documentation of functions in https://doi.org/10.1101/2023.03.22.533737 for algorithms and definitions.

Value

An object of class 'matrix' or dist' giving distances between individuals

Author(s)

Author(s): Custodian: Arthur Georges – Post to #' https://groups.google.com/d/forum/dartr

See Also

Other distance: gl.dist.pop(), gl.fdsim(), utils.dist.ind.snp()

Examples

D <- gl.dist.ind(testset.gl[1:20,], method='manhattan')
D <- gl.dist.ind(testset.gs[1:20,], method='Jaccard',swap=TRUE)

D <- gl.dist.ind(testset.gl[1:20,], method='euclidean',scale=TRUE)

Generates a distance matrix from a SNP genlight object taking into account a substitution model

Description

Generates a distance matrix for individuals or populations in a genlight object using one of a selection of substitution models.

Usage

gl.dist.phylo(
  xx,
  subst.model = "F81",
  min.tag.len = NULL,
  pairwise.missing = TRUE,
  by.pop = TRUE,
  verbose = NULL
)

Arguments

Details

The script takes a genlight object as input, creates a set of sequences from the trimmed sequence tags for each individual, calculates distances between the individuals and then optionally averages those distances between the populations defined in the genlight object (typically OTUs).

min.tag.length : Sequence tags can vary considerably in length, which results in large numbers of Ns in alignments. This can have an impact of distance measures depending on how missing values are managed. To minimize this effect, you might elect to filter on tag length using this parameter.

subst.model : Use this parameter to specify the substitution model, selecting from the list used by package ape.

• raw: This is simply the proportion or the number of sites that differ between each pair of sequences. This may be useful to draw “saturation plots”. The options variance and gamma have no effect, but pairwise.deletion can.

• TS, TV: These are the numbers of transitions and transversions, respectively.

• JC69: This model was developed by Jukes and Cantor (1969). It assumes that all substitutions (i.e. a change of a base by another one) have the same probability. This probability is the same for all sites along the DNA sequence. This last assumption can be relaxed by assuming that the substition rate varies among site following a gamma distribution which parameter must be given by the user. By default, no gamma correction is applied. Another assumption is that the base frequencies are balanced and thus equal to 0.25.

• K80: The distance derived by Kimura (1980), sometimes referred to as “Kimura's 2-parameters distance”, has the same underlying assumptions than the Jukes–Cantor distance except that two kinds of substitutions are considered: transitions (A <-> G, C <-> T), and transversions (A <-> C, A <-> T, C <-> G, G <-> T). They are assumed to have different probabilities. A transition is the substitution of a purine (C, T) by another one, or the substitution of a pyrimidine (A, G) by another one. A transversion is the substitution of a purine by a pyrimidine, or vice-versa. Both transition and transversion rates are the same for all sites along the DNA sequence. Jin and Nei (1990) modified the Kimura model to allow for variation among sites following a gamma distribution. Like for the Jukes–Cantor model, the gamma parameter must be given by the user. By default, no gamma correction is applied.

• F81: Felsenstein (1981) generalized the Jukes–Cantor model by relaxing the assumption of equal base frequencies. The formulae used in this function were taken from McGuire et al. (1999).

• K81: Kimura (1981) generalized his model (Kimura 1980) by assuming different rates for two kinds of transversions: A <-> C and G <-> T on one side, and A <-> T and C <-> G on the other. This is what Kimura called his “three substitution types model” (3ST), and is sometimes referred to as “Kimura's 3-parameters distance”.

• F84: This model generalizes K80 by relaxing the assumption of equal base frequencies. It was first introduced by Felsenstein in 1984 in Phylip, and is fully described by Felsenstein and Churchill (1996). The formulae used in this function were taken from McGuire et al. (1999).

• BH87: Barry and Hartigan (1987) developed a distance based on the observed proportions of changes among the four bases. This distance is not symmetric.

• T92: Tamura (1992) generalized the Kimura model by relaxing the assumption of equal base frequencies. This is done by taking into account the bias in G+C content in the sequences. The substitution rates are assumed to be the same for all sites along the DNA sequence.

• TN93: Tamura and Nei (1993) developed a model which assumes distinct rates for both kinds of transition (A <-> G versus C <-> T), and transversions. The base frequencies are not assumed to be equal and are estimated from the data. A gamma correction of the inter-site variation in substitution rates is possible.

• GG95: Galtier and Gouy (1995) introduced a model where the G+C content may change through time. Different rates are assumed for transitons and transversions.

• logdet: The Log-Det distance, developed by Lockhart et al. (1994), is related to BH87. However, this distance is symmetric. Formulae from Gu and Li (1996) are used. dist.logdet in phangorn uses a different implementation that gives substantially different distances for low-diverging sequences.

• paralin: Lake (1994) developed the paralinear distance which can be viewed as another variant of the Barry–Hartigan distance.

• pairwise.missing : If TRUE, then missing values in the sequence (NNNs) will be accommodated in the calculations pair of taxa at a time; otherwise, the deletion of data at positions in the sequence will be global (deleted if any missing data at the position in any individual).

Value

The distance matrix as an object of class dist.

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

Other phylogeny: gl.tree.fitch()

Examples

## Not run: 
tmp <- gl.filter.monomorphs(testset.gl)
gl.dist.phylo(xx=tmp,subst.model="F80")

## End(Not run)

Calculates a distance matrix for populations with SNP or Silicodart genotypes in a genlight object

Description

This script calculates various distances between populations based on allele frequencies (SNP genotypes) or frequency of presences in PA (SilicoDArT) data

Usage

gl.dist.pop(
  x,
  as.pop = NULL,
  method = "euclidean",
  scale = FALSE,
  type = "dist",
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = NULL,
  plot.file = NULL,
  plot.dir = NULL,
  verbose = NULL
)

Arguments

Details

For SNP data, the distance measure can be one of 'euclidean', 'fixed-diff', 'reynolds', 'nei' and 'chord'. For SilicoDArT data, the distance measure can be one of 'Refer to the documentation of functions in https://doi.org/10.1101/2023.03.22.533737 for algorithms and definitions.

Value

An object of class 'dist' giving distances between populations

Author(s)

author(s): Arthur Georges. Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

Other distance: gl.dist.ind(), gl.fdsim(), utils.dist.ind.snp()

Examples

 
# SNP genotypes
D <- gl.dist.pop(possums.gl, method='euclidean')
D <- gl.dist.pop(possums.gl, method='euclidean',scale=TRUE)
D <- gl.dist.pop(possums.gl, method='nei')
D <- gl.dist.pop(possums.gl, method='reynolds')
D <- gl.dist.pop(possums.gl, method='chord')
D <- gl.dist.pop(possums.gl, method='fixed-diff')
#Presence-Absence data [only 10 individuals due to speed]
D <- gl.dist.pop(testset.gs[1:10,], method='euclidean')

Removes specified individuals from a dartR genlight object

Description

This function deletes individuals and their associated metadata. Monomorphic loci and loci that are scored all NA are optionally deleted (mono.rm=TRUE). The script also optionally recalculates locus metatdata statistics to accommodate the deletion of individuals from the dataset (recalc=TRUE).

The script returns a dartR genlight object with the retained individuals and the recalculated locus metadata. The script works with both genlight objects containing SNP genotypes and Tag P/A data (SilicoDArT).

Usage

gl.drop.ind(x, ind.list, recalc = FALSE, mono.rm = FALSE, verbose = NULL)

Arguments

Value

A reduced dartR genlight object

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

gl.keep.ind to keep rather than drop specified individuals

Other data manipulation: gl.define.pop(), gl.drop.loc(), gl.drop.pop(), gl.edit.recode.pop(), gl.impute(), gl.join(), gl.keep.ind(), gl.keep.loc(), gl.keep.pop(), gl.make.recode.ind(), gl.merge.pop(), gl.reassign.pop(), gl.recode.ind(), gl.recode.pop(), gl.rename.pop(), gl.sample(), gl.sim.genotypes(), gl.sort(), gl.subsample.ind(), gl.subsample.loc()

Examples

 # SNP data
   gl2 <- gl.drop.ind(testset.gl,
   ind.list=c('AA019073','AA004859'))
 # Tag P/A data
   gs2 <- gl.drop.ind(testset.gs,
   ind.list=c('AA020656','AA19077','AA004859'))
   gs2 <- gl.drop.ind(testset.gs, ind.list=c('AA020656'
   ,'AA19077','AA004859'),mono.rm=TRUE, recalc=TRUE)
   

Removes specified loci from a dartR genlight object

Description

This function deletes individuals and their associated metadata. The script returns a dartR genlight object with the retained loci. The script works with both genlight objects containing SNP genotypes and Tag P/A data (SilicoDArT).

Usage

gl.drop.loc(x, loc.list = NULL, first = NULL, last = NULL, verbose = NULL)

Arguments

Value

A reduced dartR genlight object

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

gl.keep.loc to keep rather than drop specified loci

Other data manipulation: gl.define.pop(), gl.drop.ind(), gl.drop.pop(), gl.edit.recode.pop(), gl.impute(), gl.join(), gl.keep.ind(), gl.keep.loc(), gl.keep.pop(), gl.make.recode.ind(), gl.merge.pop(), gl.reassign.pop(), gl.recode.ind(), gl.recode.pop(), gl.rename.pop(), gl.sample(), gl.sim.genotypes(), gl.sort(), gl.subsample.ind(), gl.subsample.loc()

Examples

# SNP data
  gl2 <- gl.drop.loc(testset.gl, loc.list=c('100051468|42-A/T', '100049816-51-A/G'),verbose=3)
# Tag P/A data
  gs2 <- gl.drop.loc(testset.gs, loc.list=c('20134188','19249144'),verbose=3)
  

Removes specified populations from a dartR genlight object

Description

Individuals are assigned to populations based on associated specimen metadata stored in the dartR genlight object. This function deletes all individuals in the nominated populations (pop.list). Monomorphic loci and loci that are scored all NA are optionally deleted (mono.rm=TRUE). The script also optionally recalculates locus metatdata statistics to accommodate the deletion of individuals from the dataset (recalc=TRUE). The script returns a dartR genlight object with the retained populations and the recalculated locus metadata. The script works with both genlight objects containing SNP genotypes and Tag P/A data (SilicoDArT).

Usage

gl.drop.pop(
  x,
  pop.list,
  as.pop = NULL,
  recalc = FALSE,
  mono.rm = FALSE,
  verbose = NULL
)

Arguments

Value

A reduced dartR genlight object

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

gl.keep.pop to keep rather than drop specified populations

Other data manipulation: gl.define.pop(), gl.drop.ind(), gl.drop.loc(), gl.edit.recode.pop(), gl.impute(), gl.join(), gl.keep.ind(), gl.keep.loc(), gl.keep.pop(), gl.make.recode.ind(), gl.merge.pop(), gl.reassign.pop(), gl.recode.ind(), gl.recode.pop(), gl.rename.pop(), gl.sample(), gl.sim.genotypes(), gl.sort(), gl.subsample.ind(), gl.subsample.loc()

Examples

 # SNP data
   gl2 <- gl.drop.pop(testset.gl,
   pop.list=c('EmsubRopeMata','EmvicVictJasp'),verbose=3)
   gl2 <- gl.drop.pop(testset.gl, pop.list=c('EmsubRopeMata','EmvicVictJasp'),
   mono.rm=TRUE,recalc=TRUE)
   gl2 <- gl.drop.pop(testset.gl,as.pop='sex',pop.list=c('Male','Unknown'),verbose=3)
 # Tag P/A data
   gs2 <- gl.drop.pop(testset.gs, pop.list=c('EmsubRopeMata','EmvicVictJasp'))

Creates or edits individual (=specimen) names, creates a recode_ind file and applies the changes to a genlight object @family data manipulation

Description

A function to edit names of individual in a dartR genlight object, or to create a reassignment table taking the individual labels from a genlight object, or to edit existing individual labels in an existing recode_ind file. The amended recode table is then applied to the genlight object.

Usage

gl.edit.recode.ind(
  x,
  out.recode.file = NULL,
  outpath = NULL,
  recalc = FALSE,
  mono.rm = FALSE,
  verbose = NULL
)

Arguments

Details

Renaming individuals may be required when there have been errors in labeling arising in the passage of samples to sequencing. There may be occasions where renaming individuals is required for preparation of figures. This function will input an existing recode table for editing and optionally save it as a new table, or if the name of an input table is not supplied, will generate a table using the individual labels in the parent genlight object. When caution needs to be exercised because of the potential for breaking the 'chain of evidence' associated with the samples, recoding individuals using a recode table (csv) can provide a durable record of the changes. For SNP genotype data, the function, having deleted individuals, optionally identifies resultant monomorphic loci or loci with all values missing and deletes them. The script also optionally recalculates the locus metadata as appropriate. The optional deletion of monomorphic loci and the optional recalculation of locus statistics is not available for Tag P/A data (SilicoDArT). Use outpath=getwd() when calling this function to direct output files to your working directory. The function returns a dartR genlight object with the new population assignments and the recalculated locus metadata.

Value

An object of class ('genlight') with the revised individual labels.

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

gl.recode.ind, gl.drop.ind, gl.keep.ind

Examples

#this is an interactive example
if(interactive()){
gl <- gl.edit.recode.ind(testset.gl)
gl <- gl.edit.recode.ind(testset.gl, out.recode.file='ind.recode.table.csv')
}

Creates or edits and applies a population re-assignment table

Description

A function to edit population assignments in a dartR genlight object, or to create a reassignment table taking the population assignments from a genlight object, or to edit existing population assignments in a pop.recode.table. The amended recode table is then applied to the genlight object.

Usage

gl.edit.recode.pop(
  x,
  pop.recode = NULL,
  out.recode.file = NULL,
  outpath = NULL,
  recalc = FALSE,
  mono.rm = FALSE,
  verbose = NULL
)

Arguments

Value

A genlight object with the revised population assignments

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

gl.recode.pop, gl.drop.pop, gl.keep.pop, gl.merge.pop, gl.reassign.pop

Other data manipulation: gl.define.pop(), gl.drop.ind(), gl.drop.loc(), gl.drop.pop(), gl.impute(), gl.join(), gl.keep.ind(), gl.keep.loc(), gl.keep.pop(), gl.make.recode.ind(), gl.merge.pop(), gl.reassign.pop(), gl.recode.ind(), gl.recode.pop(), gl.rename.pop(), gl.sample(), gl.sim.genotypes(), gl.sort(), gl.subsample.ind(), gl.subsample.loc()

Examples

#this is an interactive example
if(interactive()){
gl <- gl.edit.recode.pop(testset.gl)
gs <- gl.edit.recode.pop(testset.gs)
}

# See also -------------------

Estimates the rate of false positives in a fixed difference analysis

Description

This function takes two populations and generates allele frequency profiles for them. It then samples an allele frequency for each, at random, and estimates a sampling distribution for those two allele frequencies. Drawing two samples from those sampling distributions, it calculates whether or not they represent a fixed difference. This is applied to all loci, and the number of fixed differences so generated are counted, as an expectation. The script distinguished between true fixed differences (with a tolerance of delta), and false positives. The simulation is repeated a given number of times (default=1000) to provide an expectation of the number of false positives, given the observed allele frequency profiles and the sample sizes. The probability of the observed count of fixed differences is greater than the expected number of false positives is calculated.

Usage

gl.fdsim(
  x,
  poppair,
  obs = NULL,
  sympatric = FALSE,
  reps = 1000,
  delta = 0.02,
  verbose = NULL
)

Arguments

Value

A list containing the following square matrices [[1]] observed fixed differences; [[2]] mean expected number of false positives for each comparison; [[3]] standard deviation of the no. of false positives for each comparison; [[4]] probability the observed fixed differences arose by chance for each comparison.

Author(s)

Custodian: Arthur Georges (Post to https://groups.google.com/d/forum/dartr)

See Also

Other distance: gl.dist.ind(), gl.dist.pop(), utils.dist.ind.snp()

Examples

fd <- gl.fdsim(testset.gl[,1:100],poppair=c('EmsubRopeMata','EmmacBurnBara'),
sympatric=TRUE,verbose=3)

Filters loci that are all NA across individuals and/or populations with all NA across loci

Description

This script deletes loci or individuals with all calls missing (NA), from a genlight object.

A DArT dataset will not have loci for which the calls are scored all as missing (NA) for a particular individual, but such loci can arise rarely when populations or individuals are deleted. Similarly, a DArT dataset will not have individuals for which the calls are scored all as missing (NA) across all loci, but such individuals may sneak in to the dataset when loci are deleted. Retaining individual or loci with all NAs can cause issues for several functions.

Also, on occasions an analysis will require that there are some loci scored in each population. Setting by.pop=TRUE will result in removal of loci when they are all missing in any one population.

Note that loci that are missing for all individuals in a population are not imputed with method 'frequency' or 'HW'. Consider using the function gl.filter.allna with by.pop=TRUE.

Usage

gl.filter.allna(x, by.pop = FALSE, recalc = FALSE, verbose = NULL)

Arguments

Value

A genlight object having removed individuals that are scored NA across all loci, or loci that are scored NA across all individuals.

Author(s)

Author(s): Arthur Georges. Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

Other filter functions: gl.filter.hwe(), gl.report.allna()

Examples

# SNP data
  result <- gl.filter.allna(testset.gl, verbose=3)
# Tag P/A data
  result <- gl.filter.allna(testset.gs, verbose=3)

Filters loci or specimens in a genlight {adegenet} object based on call rate

Description

SNP datasets generated by DArT have missing values primarily arising from failure to call a SNP because of a mutation at one or both of the restriction enzyme recognition sites. The script gl.filter.callrate() will filter out the loci with call rates below a specified threshold. Tag Presence/Absence datasets (SilicoDArT) have missing values where it is not possible to determine reliably if there the sequence tag can be called at a particular locus.

Usage

gl.filter.callrate(
  x,
  method = "loc",
  threshold = 0.95,
  mono.rm = FALSE,
  recalc = FALSE,
  recursive = FALSE,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = NULL,
  plot.file = NULL,
  plot.dir = NULL,
  bins = 25,
  verbose = NULL
)

Arguments

Details

Because this filter operates on call rate, this function recalculates Call Rate, if necessary, before filtering. If individuals are removed using method='ind', then the call rate stored in the genlight object is, optionally, recalculated after filtering. Note that when filtering individuals on call rate, the initial call rate is calculated and compared against the threshold. After filtering, if mono.rm=TRUE, the removal of monomorphic loci will alter the call rates. Some individuals with a call rate initially greater than the nominated threshold, and so retained, may come to have a call rate lower than the threshold. If this is a problem, repeated iterations of this function will resolve the issue. This is done by setting mono.rm=TRUE and recursive=TRUE, or it can be done manually. Callrate is summarized by locus or by individual to allow sensible decisions on thresholds for filtering taking into consideration consequential loss of data. The summary is in the form of a tabulation and plots. Plot themes can be obtained from

• https://ggplot2.tidyverse.org/reference/ggtheme.html and

• https://yutannihilation.github.io/allYourFigureAreBelongToUs/ggthemes/

Resultant ggplot(s) and the tabulation(s) are saved to the session's temporary directory.

Value

The reduced genlight or genind object, plus a summary

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

gl.report.callrate

Other matched filter: gl.filter.hamming(), gl.filter.ld(), gl.filter.locmetric(), gl.filter.maf(), gl.filter.monomorphs(), gl.filter.overshoot(), gl.filter.pa(), gl.filter.secondaries()

Examples

 
# SNP data
  result <- gl.filter.callrate(testset.gl[1:10], method='loc', threshold=0.8,
   verbose=3)
  result <- gl.filter.callrate(testset.gl[1:10], method='ind', threshold=0.8,
   verbose=3)
  result <- gl.filter.callrate(testset.gl[1:10], method='pop', threshold=0.8,
   verbose=3)
# Tag P/A data
  result <- gl.filter.callrate(testset.gs[1:10], method='loc', 
  threshold=0.95, verbose=3)
  result <- gl.filter.callrate(testset.gs[1:10], method='ind', 
  threshold=0.8, verbose=3)
  result <- gl.filter.callrate(testset.gs[1:10], method='pop', 
  threshold=0.8, verbose=3)
  
  res <- gl.filter.callrate(platypus.gl)

Filters excessively-heterozygous loci from a genlight object

Description

Calculates excess of heterozygosity in a genlight object and remove those loci

Usage

gl.filter.excess.het(
  x,
  Yates = FALSE,
  mono.rm = FALSE,
  recalc = FALSE,
  verbose = NULL
)

Arguments

Value

Returns unaltered genlight object

Author(s)

Author(s): Jesús Castrejón-Figueroa, Diana A Robledo-Ruiz (Custodian: Ching Ching Lau) – Post to https://groups.google.com/d/forum/dartr

References

• https://github.com/drobledoruiz/conservation_genomics/tree/main/filter.excess.het

• Robledo‐Ruiz, D. A., Austin, L., Amos, J. N., Castrejón‐Figueroa, J., Harley, D. K., Magrath, M. J., Sunnucks, P. & Pavlova, A. (2023). Easy‐to‐use R functions to separate reduced‐representation genomic datasets into sex‐linked and autosomal loci, and conduct sex assignment. Molecular Ecology Resources.

See Also

gl.filter.callrate

Other matched report: gl.report.allna(), gl.report.callrate(), gl.report.hamming(), gl.report.locmetric(), gl.report.maf(), gl.report.overshoot(), gl.report.pa(), gl.report.rdepth(), gl.report.reproducibility(), gl.report.secondaries(), gl.report.taglength()

Examples

filtered.gl <- gl.filter.excess.het(x = LBP, Yates = TRUE)
# Use below function to output information of the loci with Yates's continuity correction specified 
filtered.table <- gl.report.excess.het(x = LBP, Yates = TRUE)

Filters loci based on factor loadings for a PCA or PCoA

Description

Extracts the factor loadings from a glPCA object (generated by gl.pcoa) and filters loci based on a user specified threshold for the ABSOLUTE value of the factor loadings.

Usage

gl.filter.factorloadings(
  x,
  pca,
  axis = 1,
  threshold,
  retain = FALSE,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = NULL,
  plot.file = NULL,
  plot.dir = NULL,
  bins = 25,
  verbose = NULL,
  ...
)

Arguments

Details

The function extracts the factor loadings for a given axis from a PCA object generated by gl.pcoa and then filters loci on the basis of a user specified threshold. The threshold value is decided using gl.report.factorloadings. The function can be used to filter out loci that load high with a particular axis or alternatively if retain=TRUE, to retain loci that load high on a specified axis.

Note that this function also removes monomorphic loci because PCA is performed only on polymorphic loci.

A color vector can be obtained with gl.select.colors() and then passed to the function with the plot.colors parameter.

Themes can be obtained from in

• https://ggplot2.tidyverse.org/reference/ggtheme.html and

• https://yutannihilation.github.io/allYourFigureAreBelongToUs/ggthemes/

If a plot.file is given, the ggplot arising from this function is saved as an "RDS" binary file using saveRDS(); can be reloaded with readRDS(). A file name must be specified for the plot to be saved. If a plot directory (plot.dir) is specified, the ggplot binary is saved to that directory; otherwise to the tempdir().

Value

The unchanged genlight object

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

Examples

pca <- gl.pcoa(testset.gl)
gl.report.factorloadings(pca = pca)
gl2 <- gl.filter.factorloadings(pca=pca,x=testset.gl,threshold=0.2)

Filters loci based on pairwise Hamming distance between sequence tags

Description

Hamming distance is calculated as the number of base differences between two sequences which can be expressed as a count or a proportion. Typically, it is calculated between two sequences of equal length. In the context of DArT trimmed sequences, which differ in length but which are anchored to the left by the restriction enzyme recognition sequence, it is sensible to compare the two trimmed sequences starting from immediately after the common recognition sequence and terminating at the last base of the shorter sequence.

Usage

gl.filter.hamming(
  x,
  threshold = 0.2,
  rs = 5,
  tag.length = 69,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = NULL,
  plot.file = NULL,
  plot.dir = NULL,
  pb = FALSE,
  verbose = NULL
)

Arguments

Details

Hamming distance can be computed by exploiting the fact that the dot product of two binary vectors x and (1-y) counts the corresponding elements that are different between x and y. This approach can also be used for vectors that contain more than two possible values at each position (e.g. A, C, T or G). If a pair of DNA sequences are of differing length, the longer is truncated. The algorithm is that of Johann de Jong https://johanndejong.wordpress.com/2015/10/02/faster-hamming-distance-in-r-2/ as implemented in utils.hamming. Only one of two loci are retained if their Hamming distance is less that a specified percentage. 5 base differences out of 100 bases is a 20

Value

A genlight object filtered on Hamming distance.

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

Other matched filter: gl.filter.callrate(), gl.filter.ld(), gl.filter.locmetric(), gl.filter.maf(), gl.filter.monomorphs(), gl.filter.overshoot(), gl.filter.pa(), gl.filter.secondaries()

Examples

# SNP data
test <- gl.subsample.loc(platypus.gl,n=50)
result <- gl.filter.hamming(test, threshold=0.6, verbose=3)

Filters individuals with average heterozygosity greater than a specified upper threshold or less than a specified lower threshold @family matched filter

Description

Calculates the observed heterozygosity for each individual in a genlight object and filters individuals based on specified threshold values. Use gl.report.heterozygosity to determine the appropriate thresholds.

Usage

gl.filter.heterozygosity(x, t.upper = 0.7, t.lower = 0, verbose = NULL)

Arguments

Value

The filtered genlight object.

Author(s)

Custodian: Luis Mijangos – Post to https://groups.google.com/d/forum/dartr

Examples

 result <- gl.filter.heterozygosity(testset.gl,t.upper=0.06,verbose=3)
 tmp <- gl.report.heterozygosity(result,method='ind')
 

Filters loci that show significant departure from Hardy-Weinberg Equilibrium @family matched filter

Description

This function filters out loci showing significant departure from H-W proportions based on observed frequencies of reference homozygotes, heterozygotes and alternate homozygotes. Loci are filtered out if they show HWE departure either in any one population (n.pop.threshold =1) or in at least X number of populations (n.pop.threshold > 1).

Usage

gl.filter.hwe(
  x,
  subset = "each",
  n.pop.threshold = 1,
  test.type = "Exact",
  mult.comp.adj = FALSE,
  mult.comp.adj.method = "BY",
  alpha = 0.05,
  pvalue.type = "midp",
  cc.val = 0.5,
  n.min = 5,
  verbose = NULL
)

Arguments

Details

Several factors can cause deviations from Hardy-Weinberg equilibrium including: mutation, finite population size, selection, population structure, age structure, assortative mating, sex linkage, nonrandom sampling and genotyping errors. Refer to Waples (2015). Note that tests for departure from H-W equilibrium are only valid if there is no population substructure (assuming random mating) and have sufficient power only when there is sufficient sample size (n individuals > 15). Populations can be defined in three ways:

• Merging all populations in the dataset using subset = 'all'.

• Within each population separately using: subset = 'each'.

• Within selected populations using for example: subset = c('pop1','pop2').

Two different statistical methods to test for deviations from Hardy Weinberg proportions:

• The classical chi-square test (test.type='ChiSquare') based on the function HWChisq of the R package HardyWeinberg. By default a continuity correction is applied (cc.val=0.5). The continuity correction can be turned off (by specifying cc.val=0), for example when extreme allele frequencies occur continuity correction can lead to excessive Type I error rates.

• The exact test (test.type='Exact') based on the exact calculations contained in the function HWExactStats of the R package HardyWeinberg as described by Wigginton et al. (2005). The exact test is recommended in most cases. Three different methods to estimate p-values (pvalue.type) in the Exact test can be used:

  • 'dost' p-value is computed as twice the tail area of a one-sided test.

  • 'selome' p-value is computed as the sum of the probabilities of all samples less or equally likely as the current sample.

  • 'midp', p-value is computed as half the probability of the current sample + the probabilities of all samples that are more extreme.

  The standard exact p-value is overly conservative, in particular for small minor allele frequencies. The mid p-value ameliorates this problem by bringing the rejection rate closer to the nominal level, at the price of occasionally exceeding the nominal level (Graffelman & Moreno, 2013).

Correction for multiple tests can be applied using the following methods based on the function p.adjust:

• 'holm' is also known as the sequential Bonferroni technique (Rice, 1989). This method has a greater statistical power than the standard Bonferroni test, however this method becomes very stringent when many tests are performed and many real deviations from the null hypothesis can go undetected (Waples, 2015).

• 'hochberg' based on Hochberg, 1988.

• 'hommel' based on Hommel, 1988. This method is more powerful than Hochberg's, but the difference is usually small.

• 'bonferroni' in which p-values are multiplied by the number of tests. This method is very stringent and therefore has reduced power to detect multiple departures from the null hypothesis.

• 'BH' based on Benjamini & Hochberg, 1995.

• 'BY' based on Benjamini & Yekutieli, 2001.

The first four methods are designed to give strong control of the family-wise error rate. The last two methods control the false discovery rate (FDR), the expected proportion of false discoveries among the rejected hypotheses. The false discovery rate is a less stringent condition than the family-wise error rate, so these methods are more powerful than the others, especially when number of tests is large. The number of tests on which the adjustment for multiple comparisons is the number of populations times the number of loci. From v2.1 gl.filter.hwe takes the argument n.pop.threshold. if n.pop.threshold > 1 loci will be removed only if they are concurrently significant (after adjustment if applied) out of hwe in >= n.pop.threshold > 1.

Value

A genlight object with the loci departing significantly from H-W proportions removed.

Author(s)

Custodian: Luis Mijangos – Post to https://groups.google.com/d/forum/dartr

References

• Benjamini, Y., and Yekutieli, D. (2001). The control of the false discovery rate in multiple testing under dependency. Annals of Statistics, 29, 1165–1188.

• Graffelman, J. (2015). Exploring Diallelic Genetic Markers: The Hardy Weinberg Package. Journal of Statistical Software 64:1-23.

• Graffelman, J. & Morales-Camarena, J. (2008). Graphical tests for Hardy-Weinberg equilibrium based on the ternary plot. Human Heredity 65:77-84.

• Graffelman, J., & Moreno, V. (2013). The mid p-value in exact tests for Hardy-Weinberg equilibrium. Statistical applications in genetics and molecular biology, 12(4), 433-448.

• Hochberg, Y. (1988). A sharper Bonferroni procedure for multiple tests of significance. Biometrika, 75, 800–803.

• Hommel, G. (1988). A stagewise rejective multiple test procedure based on a modified Bonferroni test. Biometrika, 75, 383–386.

• Rice, W. R. (1989). Analyzing tables of statistical tests. Evolution, 43(1), 223-225.

• Waples, R. S. (2015). Testing for Hardy–Weinberg proportions: have we lost the plot?. Journal of heredity, 106(1), 1-19.

• Wigginton, J.E., Cutler, D.J., & Abecasis, G.R. (2005). A Note on Exact Tests of Hardy-Weinberg Equilibrium. American Journal of Human Genetics 76:887-893.

See Also

gl.report.hwe

Other filter functions: gl.filter.allna(), gl.report.allna()

Examples

result <- gl.filter.hwe(x = bandicoot.gl)

Filters loci based on linkage disequilibrium (LD)

Description

This function uses the statistic set in the parameter stat.keep from function gl.report.ld.map to choose the SNP to keep when two SNPs are in LD. When a SNP is selected to be filtered out in each pairwise comparison, the function stores its name in a list. In subsequent pairwise comparisons, if the SNP is already in the list, the other SNP will be kept.

Usage

gl.filter.ld(
  x,
  ld.report,
  threshold = 0.2,
  pop.limit = ceiling(nPop(x)/2),
  verbose = NULL
)

Arguments

Value

The reduced genlight object.

Author(s)

Custodian: Luis Mijangos – Post to https://groups.google.com/d/forum/dartr

See Also

gl.report.ld.map

Other matched filter: gl.filter.callrate(), gl.filter.hamming(), gl.filter.locmetric(), gl.filter.maf(), gl.filter.monomorphs(), gl.filter.overshoot(), gl.filter.pa(), gl.filter.secondaries()

Filters loci on the basis of numeric information stored in other$loc.metrics in a genlight {adegenet} object

Description

This script uses any field with numeric values stored in $other$loc.metrics to filter loci. The loci to keep can be within the upper and lower thresholds ('within') or outside of the upper and lower thresholds ('outside').

Usage

gl.filter.locmetric(x, metric, upper, lower, keep = "within", verbose = NULL)

Arguments

Details

The fields that are included in dartR, and a short description, are found below. Optionally, the user can also set his/her own filter by adding a vector into $other$loc.metrics as shown in the example.

1. SnpPosition - position (zero is position 1) in the sequence tag of the defined SNP variant base.

2. CallRate - proportion of samples for which the genotype call is non-missing (that is, not '-' ).

3. OneRatioRef - proportion of samples for which the genotype score is 0.

4. OneRatioSnp - proportion of samples for which the genotype score is 2.

5. FreqHomRef - proportion of samples homozygous for the Reference allele.

6. FreqHomSnp - proportion of samples homozygous for the Alternate (SNP) allele.

7. FreqHets - proportion of samples which score as heterozygous, that is, scored as 1.

8. PICRef - polymorphism information content (PIC) for the Reference allele.

9. PICSnp - polymorphism information content (PIC) for the SNP.

10. AvgPIC - average of the polymorphism information content (PIC) of the Reference and SNP alleles.

11. AvgCountRef - sum of the tag read counts for all samples, divided by the number of samples with non-zero tag read counts, for the Reference allele row.

12. AvgCountSnp - sum of the tag read counts for all samples, divided by the number of samples with non-zero tag read counts, for the Alternate (SNP) allele row.

13. RepAvg - proportion of technical replicate assay pairs for which the marker score is consistent.

Value

The reduced genlight dataset.

Author(s)

Luis Mijangos – Post to https://groups.google.com/d/forum/dartr

See Also

Other matched filter: gl.filter.callrate(), gl.filter.hamming(), gl.filter.ld(), gl.filter.maf(), gl.filter.monomorphs(), gl.filter.overshoot(), gl.filter.pa(), gl.filter.secondaries()

Examples

# adding dummy data
test <- testset.gl
test$other$loc.metrics$test <- 1:nLoc(test)
result <- gl.filter.locmetric(x=test, metric= 'test', upper=255,
lower=200, keep= 'within', verbose=3)

Filters loci on the basis of minor allele frequency (MAF) or minor allele count (MAC)

Description

This script calculates the minor allele frequency for each locus and updates the locus metadata for FreqHomRef, FreqHomSnp, FreqHets and MAF (if it exists). It then uses the updated metadata for MAF to filter loci.

Usage

gl.filter.maf(
  x,
  threshold = 0.01,
  by.pop = FALSE,
  pop.limit = ceiling(nPop(x)/2),
  ind.limit = 10,
  recalc = FALSE,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = NULL,
  plot.file = NULL,
  plot.dir = NULL,
  bins = 25,
  verbose = NULL
)

Arguments

Details

Careful consideration needs to be given to the settings to be used for this function. When the filter is applied globally (i.e. by.pop=FALSE) but the data include multiple population, there is the risk to remove markers because the allele frequencies is low (at global level) but the allele frequencies for the same markers may be high within some of the populations (especially if the per-population sample size is small). Similarly, not always it is a sensible choice to run this function using by.pop=TRUE because allele that are rare in a population may be very common in other, but the (possible) allele frequencies will depend on the sample size within each population. Where the purpose of filtering for MAF is to remove possible spurious alleles (i.e. sequencing errors), it is perhaps better to filter based on the number of times an allele is observed (MAC, Minimum Allele Count), under the assumption that if an allele is observed > MAC, it is fairly rare to be an error.

From v2.1 The threshold can take values > 1. In this case, these are interpreted as a threshold for MAC.

Value

The reduced genlight dataset

Author(s)

Custodian: Luis Mijangos – Post to https://groups.google.com/d/forum/dartr

See Also

Other matched filter: gl.filter.callrate(), gl.filter.hamming(), gl.filter.ld(), gl.filter.locmetric(), gl.filter.monomorphs(), gl.filter.overshoot(), gl.filter.pa(), gl.filter.secondaries()

Examples

result <- gl.filter.maf(platypus.gl, threshold = 0.05, verbose = 3)
#result <- gl.filter.maf(platypus.gl, by.pop = TRUE, threshold = 0.05, verbose = 3)

Filters monomorphic loci, including those with all NAs

Description

This script deletes monomorphic loci from a genlight {adegenet} object A DArT dataset will not have monomorphic loci, but they can arise, along with loci that are scored all NA, when populations or individuals are deleted. Retaining monomorphic loci unnecessarily increases the size of the dataset and will affect some calculations. Note that for SNP data, NAs likely represent null alleles; in tag presence/absence data, NAs represent missing values (presence/absence could not be reliably scored)

Usage

gl.filter.monomorphs(x, verbose = NULL)

Arguments

Value

A genlight object with monomorphic (and all NA) loci removed.

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

Other matched filter: gl.filter.callrate(), gl.filter.hamming(), gl.filter.ld(), gl.filter.locmetric(), gl.filter.maf(), gl.filter.overshoot(), gl.filter.pa(), gl.filter.secondaries()

Examples

# SNP data
  result <- gl.filter.monomorphs(testset.gl, verbose=3)
# Tag P/A data
  result <- gl.filter.monomorphs(testset.gs, verbose=3)

Filters loci for which the SNP has been trimmed from the sequence tag along with the adaptor

Description

This function checks the position of the SNP within the trimmed sequence tag and identifies those for which the SNP position is outside the trimmed sequence tag. This can happen, rarely, when the sequence containing the SNP resembles the adaptor. The SNP genotype can still be used in most analyses, but functions like gl2fasta() will present challenges if the SNP has been trimmed from the sequence tag. Not fatal, but should apply this filter before gl.filter.secondaries, for obvious reasons.

Usage

gl.filter.overshoot(x, verbose = NULL)

Arguments

Value

A new genlight object with the recalcitrant loci deleted

Author(s)

Author(s): Arthur Georges; Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

Other matched filter: gl.filter.callrate(), gl.filter.hamming(), gl.filter.ld(), gl.filter.locmetric(), gl.filter.maf(), gl.filter.monomorphs(), gl.filter.pa(), gl.filter.secondaries()

Examples

result <- gl.filter.overshoot(testset.gl, verbose=3)

Filters loci that contain private (and fixed alleles) between two populations

Description

This script is meant to be used prior to gl.nhybrids to maximise the information content of the SNPs used to identify hybrids (currently newhybrids does allow only 200 SNPs). The idea is to use first all loci that have fixed alleles between the potential source populations and then 'fill up' to 200 loci using loci that have private alleles between those. The functions filters for those loci (if invers is set to TRUE, the opposite is returned (all loci that are not fixed and have no private alleles - not sure why yet, but maybe useful.)

Usage

gl.filter.pa(x, pop1, pop2, invers = FALSE, verbose = NULL)

Arguments

Value

The reduced genlight dataset, containing now only fixed and private alleles.

Author(s)

Authors: Bernd Gruber & Ella Kelly (University of Melbourne); Custodian: Luis Mijangos – Post to https://groups.google.com/d/forum/dartr

See Also

Other matched filter: gl.filter.callrate(), gl.filter.hamming(), gl.filter.ld(), gl.filter.locmetric(), gl.filter.maf(), gl.filter.monomorphs(), gl.filter.overshoot(), gl.filter.secondaries()

Examples

result <- gl.filter.pa(testset.gl, pop1=pop(testset.gl)[1], 
pop2=pop(testset.gl)[2],verbose=3)

Filters loci based on counts of sequence tags scored at a locus (read depth) @family matched filter

Description

SNP datasets generated by DArT report AvgCountRef and AvgCountSnp as counts of sequence tags for the reference and alternate alleles respectively. These can be used to back calculate Read Depth. Fragment presence/absence datasets as provided by DArT (SilicoDArT) provide Average Read Depth and Standard Deviation of Read Depth as standard columns in their report. Filtering on Read Depth using the companion script gl.filter.rdepth can be on the basis of loci with exceptionally low counts, or loci with exceptionally high counts.

Usage

gl.filter.rdepth(
  x,
  lower = 5,
  upper = 1000,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = NULL,
  plot.file = NULL,
  plot.dir = NULL,
  verbose = NULL
)

Arguments

Details

For examples of themes, see:

• https://ggplot2.tidyverse.org/reference/ggtheme.html and

• https://yutannihilation.github.io/allYourFigureAreBelongToUs/ggthemes/

Value

Returns a genlight object retaining loci with a Read Depth in the range specified by the lower and upper threshold.

Author(s)

Custodian: Arthur Georges (Post to https://groups.google.com/d/forum/dartr)

See Also

gl.filter.rdepth

Examples

 
# SNP data
  gl.report.rdepth(testset.gl)
  result <- gl.filter.rdepth(testset.gl, lower=8, upper=50, verbose=3)
# Tag P/A data
  result <- gl.filter.rdepth(testset.gs, lower=8, upper=50, verbose=3)
  
  res <- gl.filter.rdepth(platypus.gl)

Filters loci in a genlight {adegenet} object based on average repeatability of alleles at a locus @family matched filter

Description

SNP datasets generated by DArT have an index, RepAvg, generated by reproducing the data independently for 30 of alleles that give a repeatable result, averaged over both alleles for each locus. SilicoDArT datasets generated by DArT have a similar index, Reproducibility. For these fragment presence/absence data, repeatability is the percentage of scores that are repeated in the technical replicate dataset.

Usage

gl.filter.reproducibility(
  x,
  threshold = 0.99,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = NULL,
  plot.file = NULL,
  plot.dir = NULL,
  verbose = NULL
)

Arguments

Value

Returns a genlight object retaining loci with repeatability (Repavg or Reproducibility) greater than the specified threshold.

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

gl.report.reproducibility

Examples

 
# SNP data
  gl.report.reproducibility(testset.gl)
  result <- gl.filter.reproducibility(testset.gl, threshold=0.99, verbose=3)
# Tag P/A data
  gl.report.reproducibility(testset.gs)
  result <- gl.filter.reproducibility(testset.gs, threshold=0.99)
  
  res <- gl.filter.reproducibility(testset.gl)
  

Filters loci that represent secondary SNPs in a genlight object

Description

SNP datasets generated by DArT include fragments with more than one SNP and record them separately with the same CloneID (=AlleleID). These multiple SNP loci within a fragment (secondaries) are likely to be linked, and so you may wish to remove secondaries. This script filters out all but the first sequence tag with the same CloneID after ordering the genlight object on based on repeatability, avgPIC in that order (method='best') or at random (method='random'). The filter has not been implemented for tag presence/absence data.

Usage

gl.filter.secondaries(x, method = "random", verbose = NULL)

Arguments

Value

The genlight object, with the secondary SNP loci removed.

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

Other matched filter: gl.filter.callrate(), gl.filter.hamming(), gl.filter.ld(), gl.filter.locmetric(), gl.filter.maf(), gl.filter.monomorphs(), gl.filter.overshoot(), gl.filter.pa()

Examples

gl.report.secondaries(testset.gl)
result <- gl.filter.secondaries(testset.gl)

Filters loci in a genlight {adegenet} object based on sequence tag length @family matched filter

Description

SNP datasets generated by DArT typically have sequence tag lengths ranging from 20 to 69 base pairs.

Usage

gl.filter.taglength(x, lower = 20, upper = 69, verbose = NULL)

Arguments

Value

Returns a genlight object retaining loci with a sequence tag length in the range specified by the lower and upper threshold.

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

Examples

 
# SNP data
  gl.report.taglength(testset.gl)
  result <- gl.filter.taglength(testset.gl,lower=60)
  gl.report.taglength(result)
# Tag P/A data
  gl.report.taglength(testset.gs)
  result <- gl.filter.taglength(testset.gs,lower=60)
  gl.report.taglength(result)
  
  

Generates a matrix of fixed differences and associated statistics for populations taken pairwise

Description

This script takes SNP data or sequence tag P/A data grouped into populations in a genlight object (DArTSeq) and generates a matrix of fixed differences between populations taken pairwise

Usage

gl.fixed.diff(
  x,
  tloc = 0,
  test = FALSE,
  delta = 0.02,
  alpha = 0.05,
  reps = 1000,
  mono.rm = TRUE,
  pb = FALSE,
  verbose = NULL
)

Arguments

Details

A fixed difference at a locus occurs when two populations share no alleles or where all members of one population has a sequence tag scored, and all members of the other population has the sequence tag absent. The challenge with this approach is that when sample sizes are finite, fixed differences will occur through sampling error, compounded when many loci are examined. Simulations suggest that sample sizes of n1=5 and n2=5 are adequate to reduce the probability of [experiment-wide] type 1 error to negligible levels [ploidy=2]. A warning is issued if comparison between two populations involves sample sizes less than 5, taking into account allele drop-out. Optionally, if test=TRUE, the script will test the fixed differences between final OTUs for statistical significance, using simulation, and then further amalgamate populations that for which there are no significant fixed differences at a specified level of significance (alpha). To avoid conflation of true fixed differences with false positives in the simulations, it is necessary to decide a threshold value (delta) for extreme true allele frequencies that will be considered fixed for practical purposes. That is, fixed differences in the sample set will be considered to be positives (not false positives) if they arise from true allele frequencies of less than 1-delta in one or both populations. The parameter delta is typically set to be small (e.g. delta = 0.02). NOTE: The above test will only be calculated if tloc=0, that is, for analyses of absolute fixed differences. The test applies in comparisons of allopatric populations only. For sympatric populations, use gl.pval.sympatry(). An absolute fixed difference is as defined above. However, one might wish to score fixed differences at some lower level of allele frequency difference, say where percent allele frequencies are 95,5 and 5,95 rather than 100:0 and 0:100. This adjustment can be done with the tloc parameter. For example, tloc=0.05 means that SNP allele frequencies of 95,5 and 5,95 percent will be regarded as fixed when comparing two populations at a locus.

Value

A list of Class 'fd' containing the gl object and square matrices, as follows:

1. $gl – the output genlight object;

2. $fd – raw fixed differences;

3. $pcfd – percent fixed differences;

4. $nobs – mean no. of individuals used in each comparison;

5. $nloc – total number of loci used in each comparison;

6. $expfpos – if test=TRUE, the expected count of false positives for each comparison [by simulation];

7. $sdfpos – if test=TRUE, the standard deviation of the count of false positives for each comparison [by simulation];

8. $pval – if test=TRUE, the significance of the count of fixed differences [by simulation])

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

utils.is.fixed

Examples

fd <- gl.fixed.diff(testset.gl, tloc=0, verbose=3 )
fd <- gl.fixed.diff(testset.gl, tloc=0, test=TRUE, delta=0.02, reps=100, verbose=3 )

Calculates a pairwise Fst values for populations in a genlight object This script calculates pairwise Fst values based on the implementation in the StAMPP package (?stamppFst). It allows to run bootstrap to estimate probability of Fst values to be different from zero. For detailed information please check the help pages (?stamppFst).

Description

Calculates a pairwise Fst values for populations in a genlight object This script calculates pairwise Fst values based on the implementation in the StAMPP package (?stamppFst). It allows to run bootstrap to estimate probability of Fst values to be different from zero. For detailed information please check the help pages (?stamppFst).

Usage

gl.fst.pop(x, nboots = 1, percent = 95, nclusters = 1, verbose = NULL)

Arguments

Value

A matrix of distances between populations (class dist), if nboots =1, otherwise a list with Fsts (in a matrix), Pvalues (a matrix of pvalues), Bootstraps results (data frame of all runs). Hint: Use as.matrix(as.dist(fsts)) if you want to have a squared matrix with symmetric entries returned, instead of a dist object.

Author(s)

Bernd Gruber (bugs? Post to https://groups.google.com/d/forum/dartr)

Examples

test <- gl.filter.callrate(platypus.gl,threshold = 1)
test <- gl.filter.monomorphs(test)
out <- gl.fst.pop(test, nboots=1)

Performs Hardy-Weinberg tests over loci and populations

Description

Hardy-Weinberg tests are performed for each loci in each of the populations as defined by the pop slot in a genlight object.

Usage

gl.hwe.pop(
  x,
  alpha_val = 0.05,
  plot.out = TRUE,
  plot_theme = theme_dartR(),
  plot_colors = c("gray90", "deeppink"),
  HWformat = FALSE,
  verbose = NULL
)

Arguments

Details

This function employs the HardyWeinberg package, which needs to be installed. The function that is used is HWExactStats, but there are several other great functions implemented in the package regarding HWE. Therefore, this function can return the data in the format expected by the HWE package expects, via HWformat=TRUE and then use this to run other functions of the package. This functions performs a HWE test for every population (rows) and loci (columns) and returns a true false matrix. True is reported if the p-value of an HWE-test for a particular loci and population was below the specified threshold (alpha_val, default=0.05). The thinking behind this approach is that loci that are not in HWE in several populations have most likely to be treated (e.g. filtered if loci under selection are of interest). If plot=TRUE a barplot on the loci and the sum of deviation over all population is returned. Loci that deviate in the majority of populations can be identified via colSums on the resulting matrix. Plot themes can be obtained from

• https://ggplot2.tidyverse.org/reference/ggtheme.html and

• https://yutannihilation.github.io/allYourFigureAreBelongToUs/ggthemes/

Resultant ggplots and the tabulation are saved to the session's temporary directory.

Value

The function returns a list with up to three components:

• 'HWE' is the matrix over loci and populations

• 'plot' is a plot (ggplot) which shows the significant results for population and loci (can be amended further using ggplot syntax)

• 'HWEformat=TRUE' the 'HWformat' entails SNP data for each population in 'HardyWeinberg'-format to be used with other functions of the package (e.g HWPerm or HWExactPrevious).

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

Examples

out <- gl.hwe.pop(bandicoot.gl[,1:33], alpha_val=0.05, plot.out=TRUE, HWformat=FALSE)

Imputes missing data

Description

This function imputes genotypes on a population-by-population basis, where populations can be considered panmictic, or imputes the state for presence-absence data.

Usage

gl.impute(
  x,
  method = "neighbour",
  fill.residual = TRUE,
  parallel = FALSE,
  verbose = NULL
)

Arguments

Details

We recommend that imputation be performed on sampling locations, before any aggregation. Imputation is achieved by replacing missing values using either of four methods:

• If "frequency", genotypes scored as missing at a locus in an individual are imputed using the average allele frequencies at that locus in the population from which the individual was drawn.

• If "HW", genotypes scored as missing at a locus in an individual are imputed by sampling at random assuming Hardy-Weinberg equilibrium. Applies only to genotype data.

• If "neighbour", substitute the missing values for the focal individual with the values taken from the nearest neighbour. Repeat with next nearest and so on until all missing values are replaced.

• if "random", missing data are substituted by random values (0, 1 or 2).

The nearest neighbour is the one at the smallest Euclidean distancefrom the focal individual The advantage of this approach is that it works regardless of how many individuals are in the population to which the focal individual belongs, and the displacement of the individual is haphazard as opposed to: (a) Drawing the individual toward the population centroid (HW and Frequency). (b) Drawing the individual toward the global centroid (glPCA). Note that loci that are missing for all individuals in a population are not imputed with method 'frequency' or 'HW' and can give unpredictable results for particular individuals using 'neighbour'. Consider using the function gl.filter.allna with by.pop=TRUE to remove them first.

Value

A genlight object with the missing data imputed.

Author(s)

Custodian: Luis Mijangos (Post to https://groups.google.com/d/forum/dartr)

See Also

Other data manipulation: gl.define.pop(), gl.drop.ind(), gl.drop.loc(), gl.drop.pop(), gl.edit.recode.pop(), gl.join(), gl.keep.ind(), gl.keep.loc(), gl.keep.pop(), gl.make.recode.ind(), gl.merge.pop(), gl.reassign.pop(), gl.recode.ind(), gl.recode.pop(), gl.rename.pop(), gl.sample(), gl.sim.genotypes(), gl.sort(), gl.subsample.ind(), gl.subsample.loc()

Examples

 
require("dartR.data")
# SNP genotype data
gl <- gl.filter.callrate(platypus.gl,threshold=0.95)
gl <- gl.filter.allna(gl)
gl <- gl.impute(gl,method="neighbour")
# Sequence Tag presence-absence data
gs <- gl.filter.callrate(testset.gs,threshold=0.95)
gl <- gl.filter.allna(gl)
gs <- gl.impute(gs, method="neighbour")

gs <- gl.impute(platypus.gl,method ="random")

Combines two dartR genlight objects

Description

This function combines two genlight objects and their associated metadata. The history associated with the two genlight objects is cleared from the new genlight object. The individuals/samples must be the same in each genlight object. The function is typically used to combine datasets from the same service where the files have been split because of size limitations. The data is read in from multiple csv files, then the resultant genlight objects are combined. This function works with both SNP and Tag P/A data.

Usage

gl.join(x1, x2, method = "sidebyside", verbose = NULL)

Arguments

Details

This script joins two genlight objects together along with the associated metadata. if method='sidebyside' (the default), the individuals in the two genlight objects must be the same and in the same order. The loci are combined.

If method='end2end', the loci in the two genlight objects must be the same and in the same order. The data for the two sets of individuals are combined. Note that if two individuals have the same names, they will be made unique.#'

Value

A new genlight object

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

Other data manipulation: gl.define.pop(), gl.drop.ind(), gl.drop.loc(), gl.drop.pop(), gl.edit.recode.pop(), gl.impute(), gl.keep.ind(), gl.keep.loc(), gl.keep.pop(), gl.make.recode.ind(), gl.merge.pop(), gl.reassign.pop(), gl.recode.ind(), gl.recode.pop(), gl.rename.pop(), gl.sample(), gl.sim.genotypes(), gl.sort(), gl.subsample.ind(), gl.subsample.loc()

Examples

x1 <- testset.gl[,1:100]
x1@other$loc.metrics <-  testset.gl@other$loc.metrics[1:100,]
nLoc(x1)
x2 <- testset.gl[,101:150]
x2@other$loc.metrics <-  testset.gl@other$loc.metrics[101:150,]
nLoc(x2)
gl <- gl.join(x1, x2, verbose=2)
nLoc(gl)

x1 <- testset.gl[,1:100]
x1@other$loc.metrics <-  testset.gl@other$loc.metrics[1:100,]
nLoc(x1)
x2 <- testset.gl[,1:100]
x2@other$loc.metrics <-  testset.gl@other$loc.metrics[1:100,]
nLoc(x2)
gl <- gl.join(x1, x2, method="end2end", verbose=2)
nInd(gl)

Removes all but the specified individuals from a dartR genlight object

Description

This script deletes all individuals apart from those listed (ind.list). Monomorphic loci and loci that are scored all NA are optionally deleted (mono.rm=TRUE). The script also optionally recalculates locus metatdata statistics to accommodate the deletion of individuals from the dataset (recalc=TRUE). The script returns a dartR genlight object with the retained individuals and the recalculated locus metadata. The script works with both genlight objects containing SNP genotypes and Tag P/A data (SilicoDArT).

Usage

gl.keep.ind(x, ind.list, recalc = FALSE, mono.rm = FALSE, verbose = NULL)

Arguments

Value

A reduced dartR genlight object

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

gl.drop.pop to drop rather than keep specified populations

Other data manipulation: gl.define.pop(), gl.drop.ind(), gl.drop.loc(), gl.drop.pop(), gl.edit.recode.pop(), gl.impute(), gl.join(), gl.keep.loc(), gl.keep.pop(), gl.make.recode.ind(), gl.merge.pop(), gl.reassign.pop(), gl.recode.ind(), gl.recode.pop(), gl.rename.pop(), gl.sample(), gl.sim.genotypes(), gl.sort(), gl.subsample.ind(), gl.subsample.loc()

Examples

  # SNP data
    gl2 <- gl.keep.ind(testset.gl, ind.list=c('AA019073','AA004859'))
  # Tag P/A data
   gs2 <- gl.keep.ind(testset.gs, ind.list=c('AA020656','AA19077','AA004859'))
   

Removes all but the specified loci from a genlight object

Description

This function deletes loci that are not specified to keep, and their associated metadata. The script returns a dartR genlight object with the retained loci. The script works with both genlight objects containing SNP genotypes and Tag P/A data (SilicoDArT).

Usage

gl.keep.loc(x, loc.list = NULL, first = NULL, last = NULL, verbose = NULL)

Arguments

Value

A genlight object with the reduced data

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

gl.drop.loc to drop rather than keep specified loci

Other data manipulation: gl.define.pop(), gl.drop.ind(), gl.drop.loc(), gl.drop.pop(), gl.edit.recode.pop(), gl.impute(), gl.join(), gl.keep.ind(), gl.keep.pop(), gl.make.recode.ind(), gl.merge.pop(), gl.reassign.pop(), gl.recode.ind(), gl.recode.pop(), gl.rename.pop(), gl.sample(), gl.sim.genotypes(), gl.sort(), gl.subsample.ind(), gl.subsample.loc()

Examples

# SNP data
  gl2 <- gl.keep.loc(testset.gl, loc.list=c('100051468|42-A/T', '100049816-51-A/G'))
# Tag P/A data
  gs2 <- gl.keep.loc(testset.gs, loc.list=c('20134188','19249144'))
  

Removes all but the specified populations from a dartR genlight object

Description

Individuals are assigned to populations based on associated specimen metadata stored in the dartR genlight object. This script deletes all individuals apart from those in listed populations (pop.list). Monomorphic loci and loci that are scored all NA are optionally deleted (mono.rm=TRUE). The script also optionally recalculates locus metatdata statistics to accommodate the deletion of individuals from the dataset (recalc=TRUE). The script returns a dartR genlight object with the retained populations and the recalculated locus metadata. The script works with both genlight objects containing SNP genotypes and Tag P/A data (SilicoDArT).

Usage

gl.keep.pop(
  x,
  pop.list,
  as.pop = NULL,
  recalc = FALSE,
  mono.rm = FALSE,
  verbose = NULL
)

Arguments

Value

A reduced dartR genlight object

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

gl.drop.pop to drop rather than keep specified populations

Other data manipulation: gl.define.pop(), gl.drop.ind(), gl.drop.loc(), gl.drop.pop(), gl.edit.recode.pop(), gl.impute(), gl.join(), gl.keep.ind(), gl.keep.loc(), gl.make.recode.ind(), gl.merge.pop(), gl.reassign.pop(), gl.recode.ind(), gl.recode.pop(), gl.rename.pop(), gl.sample(), gl.sim.genotypes(), gl.sort(), gl.subsample.ind(), gl.subsample.loc()

Examples

 # SNP data
   gl2 <- gl.keep.pop(testset.gl, pop.list=c('EmsubRopeMata', 'EmvicVictJasp'))
   gl2 <- gl.keep.pop(testset.gl, pop.list=c('EmsubRopeMata', 'EmvicVictJasp'),
   mono.rm=TRUE,recalc=TRUE)
   gl2 <- gl.keep.pop(testset.gl, pop.list=c('Female'),as.pop='sex')
 # Tag P/A data
   gs2 <- gl.keep.pop(testset.gs, pop.list=c('EmsubRopeMata','EmvicVictJasp'))
   

Loads an object from compressed binary format produced by gl.save()

Description

This is a wrapper for readRDS() The function loads the object from the current workspace, checks if it is a dartR genlight object, converts it if it is not, and returns the gl object. A compliance check can be requested.

Usage

gl.load(file, compliance = FALSE, verbose = NULL)

Arguments

Value

The loaded object

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

gl.save

Other io: gl.read.csv(), gl.read.dart(), gl.read.fasta(), gl.read.silicodart(), gl.save(), gl.write.csv(), utils.read.dart()

Assigns individuals to populations with an associated probability

Description

Uses Mahalanobis Distances between individuals and group centroids to calculate probability of group membership

Usage

gl.mahal.assign(
  x,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = NULL,
  plot.file = NULL,
  plot.dir = NULL,
  verbose = NULL
)

Arguments

Details

The function generates 200 simulated individuals for each group (=population in the genlight object) drawing from the observed allele frequencies for each group. The group centroids and covariance matricies are calculated for these simulated groups. The covariance matrix is inverted using package MASS::ginv to overcome the singularities that would otherwise arise with typical SNP data. Mahanobilis Distances are calculated using stats::mahanalobis for each individual in the dataset and associated Chi Square probabilities of group membership are calculated for each individual in the original genlight object. The resultant table can be used for decisions on group membership. A special group (=population in the genlight object) called 'unknowns' can be used to specifically identify individuals with unknown group membership.

A color vector can be obtained with gl.select.colors() and then passed to the function with the plot.colors parameter. Themes can be obtained from in

• https://ggplot2.tidyverse.org/reference/ggtheme.html and

• https://yutannihilation.github.io/allYourFigureAreBelongToUs/ggthemes/

If a plot.file is given, the ggplot arising from this function is saved as an "RDS" binary file using saveRDS(); can be reloaded with readRDS(). A file name must be specified for the plot to be saved. If a plot directory (plot.dir) is specified, the ggplot binary is saved to that directory; otherwise to the tempdir().

Value

The unchanged genlight object

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

Other matched reports: gl.report.bases(), gl.report.factorloadings(), gl.report.fstat(), gl.report.monomorphs()

Creates a proforma recode_ind file for reassigning individual (=specimen) names

Description

Renaming individuals may be required when there have been errors in labeling arising in the process from sample to sequencing files. There may be occasions where renaming individuals is required for preparation of figures.

Usage

gl.make.recode.ind(
  x,
  out.recode.file = "default_recode_ind.csv",
  outpath = NULL,
  verbose = NULL
)

Arguments

Details

This function facilitates the construction of a recode table by producing a proforma file with current individual (=specimen) names in two identical columns. Edit the second column to reassign individual names. Use keyword 'Delete' to delete an individual. When caution needs to be exercised because of the potential for breaking the 'chain of evidence' associated with the samples, recoding individuals using a recode table (csv) can provide a clear record of the changes. Use outpath=getwd() or when calling this function to direct output files to your working directory. The function works with both genlight objects containing SNP genotypes and Tag P/A data (SilicoDArT). Apply the recoding using gl.recode.ind().

Value

A vector containing the new individual names.

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

Other data manipulation: gl.define.pop(), gl.drop.ind(), gl.drop.loc(), gl.drop.pop(), gl.edit.recode.pop(), gl.impute(), gl.join(), gl.keep.ind(), gl.keep.loc(), gl.keep.pop(), gl.merge.pop(), gl.reassign.pop(), gl.recode.ind(), gl.recode.pop(), gl.rename.pop(), gl.sample(), gl.sim.genotypes(), gl.sort(), gl.subsample.ind(), gl.subsample.loc()

Examples

result <- gl.make.recode.ind(testset.gl, out.recode.file ='Emmac_recode_ind.csv',outpath=tempdir())

Creates a proforma recode_pop_table file for reassigning population names @family data manipulation

Description

Renaming populations may be required when there have been errors in assignment arising in the process from sample to sequence files or when one wishes to amalgamate populations, or delete populations. Recoding populations can also be done with a recode table (csv).

Usage

gl.make.recode.pop(
  x,
  out.recode.file = "recode_pop_table.csv",
  outpath = NULL,
  verbose = NULL
)

Arguments

Details

This function facilitates the construction of a recode table by producing a proforma file with current population names in two identical columns. Edit the second column to reassign populations. Use keyword 'Delete' to delete a population. When caution needs to be exercised because of the potential for breaking the 'chain of evidence' associated with the samples, recoding individuals using a recode table (csv) can provide a clear record of the changes. Use outpath=getwd() or when calling this function to direct output files to your working directory. The function works with both genlight objects containing SNP genotypes and Tag P/A data (SilicoDArT). Apply the recoding using gl.recode.pop().

Value

A vector containing the new population names.

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

Examples

result <- gl.make.recode.pop(testset.gl,out.recode.file='test.csv',outpath=tempdir(),verbose=2)

Creates an interactive map (based on latlon) from a genlight object

Description

Creates an interactive map (based on latlon) from a genlight object

Usage

gl.map.interactive(
  x,
  matrix = NULL,
  standard = TRUE,
  symmetric = TRUE,
  pop.labels = TRUE,
  pop.labels.cex = 12,
  ind.circles = TRUE,
  ind.circle.cols = rainbow,
  ind.circle.cex = 10,
  ind.circle.transparency = 0.8,
  palette.links = NULL,
  legend.title = NULL,
  provider = "Esri.NatGeoWorldMap",
  raster.image = NULL,
  raster.opacity = 0.5,
  raster.colors = (scales::viridis_pal(option = "D"))(255),
  verbose = NULL
)

Arguments

Details

A wrapper around the leaflet package. For possible background maps check as specified via the provider: http://leaflet-extras.github.io/leaflet-providers/preview/index.html

The palette.links argument can be any of the following: A character vector of RGB or named colors. Examples: palette(), c("#000000", "#0000FF", "#FFFFFF"), topo.colors(10)

The name of an RColorBrewer palette, e.g. "BuPu" or "Greens".

The full name of a viridis palette: "viridis", "magma", "inferno", or "plasma".

A function that receives a single value between 0 and 1 and returns a color. Examples: colorRamp(c("#000000", "#FFFFFF"), interpolate = "spline").

Value

plots a map

Author(s)

Bernd Gruber – Post to https://groups.google.com/d/forum/dartr

See Also

Other graphics: gl.colors(), gl.plot.heatmap(), gl.report.ld.map(), gl.select.colors(), gl.select.shapes(), gl.smearplot(), gl.tree.nj()

Examples

require("dartR.data")
gl.map.interactive(bandicoot.gl)
cols <- c("red","blue","yellow")
gl.map.interactive(platypus.gl, ind.circle.cols=cols, ind.circle.cex=10, 
ind.circle.transparency=0.5)

Merges two or more populations in a dartR genlight object into one population

Description

Individuals are assigned to populations based on the specimen metadata data file (csv) used with gl.read.dart(). This function assigns individuals from two nominated populations into a new single population. It can also be used to rename populations. The function works with both SNP and Tag P/A (silicoDArT) data. The function returns a genlight object with the new population assignments.

Usage

gl.merge.pop(x, old = NULL, new = NULL, verbose = NULL)

Arguments

Value

A genlight object with the new population assignments.

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

Other data manipulation: gl.define.pop(), gl.drop.ind(), gl.drop.loc(), gl.drop.pop(), gl.edit.recode.pop(), gl.impute(), gl.join(), gl.keep.ind(), gl.keep.loc(), gl.keep.pop(), gl.make.recode.ind(), gl.reassign.pop(), gl.recode.ind(), gl.recode.pop(), gl.rename.pop(), gl.sample(), gl.sim.genotypes(), gl.sort(), gl.subsample.ind(), gl.subsample.loc()

Examples

   gl <- gl.merge.pop(testset.gl, old=c('EmsubRopeMata','EmvicVictJasp'), new='Outgroup')
   

Ordination applied to genotypes in a genlight object (PCA), in an fd object, or to a distance matrix (PCoA)

Description

This function takes the genotypes for individuals and undertakes a Pearson Principal Component analysis (PCA) on SNP or Sequence tag P/A (SilicoDArT) data; it undertakes a Gower Principal Coordinate analysis (PCoA) if supplied with a distance matrix.

Usage

gl.pcoa(
  x,
  nfactors = 5,
  pc.select = "broken-stick",
  correction = NULL,
  mono.rm = TRUE,
  parallel = FALSE,
  n.cores = 1,
  plot.out = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = gl.colors(2),
  plot.file = NULL,
  plot.dir = NULL,
  verbose = NULL
)

Arguments

Details

The function is essentially a wrapper for glPca {adegenet} or pcoa {ape} with default settings apart from those specified as parameters in this function. Sources of stress in the visual representation While, technically, any distance matrix can be represented in an ordinated space, the representation will not typically be exact.There are three major sources of stress in a reduced-representation of distances or dissimilarities among entities using PCA or PCoA. By far the greatest source comes from the decision to select only the top two or three axes from the ordinated set of axes derived from the PCA or PCoA. The representation of the entities such a heavily reduced space will not faithfully represent the distances in the input distance matrix simply because of the loss of information in deeper informative dimensions. For this reason, it is not sensible to be too precious about managing the other two sources of stress in the visual representation. The measure of distance between entities in a PCA is the Pearson Correlation Coefficient, essentially a standardized Euclidean distance. This is both a metric distance and a Euclidean distance and so the distances in the final ordination are faithful to those between entities in the dataset. Note that missing values are imputed in PCA, and that this can be a source of disparity between the distances between the entities in the dataset and the distances in the ordinated space.

In PCoA, the second source of stress is the choice of distance measure or dissimilarity measure. While pretty much any distance or dissimilarity matrix can be represented in an ordinated space, the distances between entities can be faithfully represented in that space (that is, without stress) only if the distances are metric. Furthermore, for distances between entities to be faithfully represented in a rigid Cartesian space, the distance measure needs to be Euclidean. If this is not the case, the distances between the entities in the ordinated visualized space will not #' exactly represent the distances in the input matrix (stress will be non-zero). This source of stress will be evident as negative eigenvalues in the deeper dimensions.

A third source of stress arises from having a sparse dataset, one with missing values. This affects both PCA and PCoA. If the original data matrix is not fully populated, that is, if there are missing values, then even a Euclidean distance matrix will not necessarily be 'positive definite'. It follows that some of the eigenvalues may be negative, even though the distance metric is Euclidean. This issue is exacerbated when the number of loci greatly exceeds the number of individuals, as is typically the case when working with SNP data. The impact of missing values can be minimized by stringently filtering on Call Rate, albeit with loss of data. An alternative is given in a paper 'Honey, I shrunk the sample covariance matrix' and more recently by Ledoit and Wolf (2018), but their approach has not been implemented here.

Options for imputing missing values while minimizing distortion are provided in the function gl.impute(). The good news is that, unless the sum of the negative eigenvalues, arising from a non-Euclidean distance measure or from missing values, approaches those of the final PCA or PCoA axes to be displayed, the distortion is probably of no practical consequence and certainly not comparable to the stress arising from selecting only two or three final dimensions out of several informative dimensions for the visual representation. Function's output Two diagnostic plots are produced. The first is a Scree Plot, showing the percentage variation explained by each of the PCA or PCoA axes, for those axes are considered informative. The scree plot informs a decision on the number of dimensions to be retained in the visual summaries. Various approaches are available for identifying which axes are informative (in terms of containing biologically significant variation) and which are noise axes. The simplest method is to consider only those axes that explain more variance than the original variables on average as being informative (pc.select="Kaiser-Guttman"). A second method (the default) is the broken-stick method (pc.select="broken-stick"). A third method is the Tracy-Widom statistical approach (pc.select="Tracy-Widom").

Once you have the informative axes identified, a judgement call is made as to how many dimensions to retain and present as results. This requires a decision on how much information on structure in the data is to be discarded. Retaining at least those axes that explain 10 The second graph is for diagnostic purposes only. It shows the distribution of eigenvalues for the remaining uninformative (noise) axes, including those with negative eigenvalues. If a plot.file is given, the ggplot arising from this function is saved as an binary file using gl.save(); can be reloaded with gl.load(). A file name must be specified for the plot to be saved. If a plot directory (plot.dir) is specified, the ggplot binary is saved to that directory; otherwise to the tempdir(). Action is recommended (verbose >= 2) if the negative eigenvalues are dominant, their sum approaching in magnitude the eigenvalues for axes selected for the final visual solution. Output is a glPca object conforming to adegenet::glPca but with only the following retained.

• $call - The call that generated the PCA/PCoA

• $eig - Eigenvalues – All eigenvalues (positive, null, negative).

• $scores - Scores (coefficients) for each individual

• $loadings - Loadings of each SNP for each principal component

Examples of other themes that can be used can be consulted in

• https://ggplot2.tidyverse.org/reference/ggtheme.html and

• https://yutannihilation.github.io/allYourFigureAreBelongToUs/ggthemes/

PCA was developed by Pearson (1901) and Hotelling (1933), whilst the best modern reference is Jolliffe (2002). PCoA was developed by Gower (1966) while the best modern reference is Legendre & Legendre (1998).

Value

An object of class pcoa containing the eigenvalues and factor scores

Author(s)

Author(s): Arthur Georges and Jesus Castrejon. Custodian: Arthur Georges (Post to https://groups.google.com/d/forum/dartr)

References

• Cailliez, F. (1983) The analytical solution of the additive constant problem. Psychometrika, 48, 305-308.

• Gower, J. C. (1966) Some distance properties of latent root and vector methods used in multivariate analysis. Biometrika, 53, 325-338.

• Hotelling, H., 1933. Analysis of a complex of statistical variables into Principal Components. Journal of Educational Psychology 24:417-441, 498-520.

• Jolliffe, I. (2002) Principal Component Analysis. 2nd Edition, Springer, New York.

• Ledoit, O. and Wolf, M. (2018). Analytical nonlinear shrinkage of large-dimensional covariance matrices. University of Zurich, Department of Economics, Working Paper No. 264, Revised version. Available at SSRN: https://ssrn.com/abstract=3047302 or http://dx.doi.org/10.2139/ssrn.3047302

• Legendre, P. and Legendre, L. (1998). Numerical Ecology, Volume 24, 2nd Edition. Elsevier Science, NY.

• Lingoes, J. C. (1971) Some boundary conditions for a monotone analysis of symmetric matrices. Psychometrika, 36, 195-203.

• Pearson, K. (1901). On lines and planes of closest fit to systems of points in space. Philosophical Magazine. Series 6, vol. 2, no. 11, pp. 559-572.

See Also

gl.pcoa.plot

Examples

# PCA (using SNP genlight object)
gl <- possums.gl
pca <- gl.pcoa(possums.gl[1:50,],verbose=2)
gl.pcoa.plot(pca,gl)

gs <- testset.gs
levels(pop(gs))<-c(rep('Coast',5),rep('Cooper',3),rep('Coast',5),
rep('MDB',8),rep('Coast',6),'Em.subglobosa','Em.victoriae')
# PCA (using SilicoDArT genlight object)
pca <- gl.pcoa(gs)
gl.pcoa.plot(pca,gs)
# Using a distance matrix
D <- gl.dist.ind(testset.gs, method='jaccard')
pcoa <- gl.pcoa(D,correction="cailliez")
gl.pcoa.plot(pcoa,gs)

Bivariate or trivariate plot of the results of an ordination generated using gl.pcoa()

Description

This script takes output from the ordination generated by gl.pcoa() and plots the individuals classified by population.

Usage

gl.pcoa.plot(
  glPca,
  x,
  scale = FALSE,
  ellipse = FALSE,
  plevel = 0.95,
  pop.labels = "pop",
  interactive = FALSE,
  as.pop = NULL,
  hadjust = 1.5,
  vadjust = 1,
  xaxis = 1,
  yaxis = 2,
  zaxis = NULL,
  pt.size = 2,
  pt.colors = NULL,
  pt.shapes = NULL,
  label.size = 1,
  axis.label.size = 1.5,
  plot.file = NULL,
  plot.dir = NULL,
  verbose = NULL
)

Arguments

Details

The factor scores are taken from the output of gl.pcoa() and the population assignments are taken from from the original data file. In the bivariate plots, the specimens are shown optionally with adjacent labels and enclosing ellipses. Population labels on the plot are shuffled so as not to overlap (using package {directlabels}). This can be a bit clunky, as the labels may be some distance from the points to which they refer, but it provides the opportunity for moving labels around using graphics software (e.g. Adobe Illustrator). 3D plotting is activated by specifying a zaxis. Any pair or trio of axes can be specified from the ordination, provided they are within the range of the nfactors value provided to gl.pcoa(). In the 2D plots, axes can be scaled to represent the proportion of variation explained. In any case, the proportion of variation explained by each axis is provided in the axis label. Colors and shapes of the points can be altered by passing a vector of shapes and/or a vector of colors. These vectors can be created with gl.select.shapes() and gl.select.colors() and passed to this script using the pt.shapes and pt.colors parameters. Points displayed in the ordination can be identified if the option interactive=TRUE is chosen, in which case the resultant plot is ggplotly() friendly. Identification of points is by moving the mouse over them. Refer to the plotly package for further information. The interactive option is automatically enabled for 3D plotting.

If a plot.file is given, the ggplot arising from this function is saved as an "RDS" binary file using saveRDS(); can be reloaded with readRDS(). A file name must be specified for the plot to be saved. If a plot directory (plot.dir) is specified, the ggplot binary is saved to that directory; otherwise to the tempdir().

Value

returns no value (i.e. NULL)

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

gl.pcoa

Examples

 test <- gl.pcoa(platypus.gl)
 gl.pcoa.plot(glPca = test, x = platypus.gl)

# SET UP DATASET
gl <- testset.gl
levels(pop(gl))<-c(rep('Coast',5),rep('Cooper',3),rep('Coast',5),
rep('MDB',8),rep('Coast',7),'Em.subglobosa','Em.victoriae')
# RUN PCA
pca<-gl.pcoa(gl,nfactors=5)
# VARIOUS EXAMPLES
gl.pcoa.plot(pca, gl, ellipse=TRUE, plevel=0.95, pop.labels='pop', 
axis.label.size=1, hadjust=1.5,vadjust=1)
gl.pcoa.plot(pca, gl, ellipse=TRUE, plevel=0.99, pop.labels='legend', 
axis.label.size=1)
gl.pcoa.plot(pca, gl, ellipse=TRUE, plevel=0.99, pop.labels='legend', 
axis.label.size=1.5,scale=TRUE)
gl.pcoa.plot(pca, gl, ellipse=TRUE, axis.label.size=1.2, xaxis=1, yaxis=3, 
scale=TRUE)
gl.pcoa.plot(pca, gl, pop.labels='none',scale=TRUE)
gl.pcoa.plot(pca, gl, axis.label.size=1.2, interactive=TRUE)
gl.pcoa.plot(pca, gl, ellipse=TRUE, plevel=0.99, xaxis=1, yaxis=2, zaxis=3)
# COLOR AND SHAPE ADJUSTMENTS
shp <- gl.select.shapes(select=c(16,17,17,0,2))
col <- gl.select.colors(library='brewer',palette='Spectral',ncolors=11,
select=c(1,9,3,11,11))
gl.pcoa.plot(pca, gl, ellipse=TRUE, plevel=0.95, pop.labels='pop', 
pt.colors=col, pt.shapes=shp, axis.label.size=1, hadjust=1.5,vadjust=1)
gl.pcoa.plot(pca, gl, ellipse=TRUE, plevel=0.99, pop.labels='legend',
 pt.colors=col, pt.shapes=shp, axis.label.size=1)
# DISTANCE MATRIX
 D <- gl.dist.ind(gl)
 pco <- gl.pcoa(D)
 gl.pcoa.plot(pco,gl,ellipse=TRUE)
 

Represents a distance matrix as a heatmap

Description

The script plots a heat map to represent the distances in the distance or dissimilarity matrix. This function is a wrapper for heatmap.2 (package gplots).

Usage

gl.plot.heatmap(D, palette.divergent = gl.colors("div"), verbose = NULL, ...)

Arguments

Value

returns no value (i.e. NULL)

Author(s)

Custodian: Luis Mijangos – Post to https://groups.google.com/d/forum/dartr)

See Also

Other graphics: gl.colors(), gl.map.interactive(), gl.report.ld.map(), gl.select.colors(), gl.select.shapes(), gl.smearplot(), gl.tree.nj()

Examples

   gl <- testset.gl[1:10,]
   D <- dist(as.matrix(gl),upper=TRUE,diag=TRUE)
   gl.plot.heatmap(D)
   D2 <- gl.dist.pop(possums.gl)
   gl.plot.heatmap(D2)
   D3 <- gl.fixed.diff(testset.gl)
   gl.plot.heatmap(D3)
   
   if ((requireNamespace("gplots", quietly = TRUE))) {
   D2 <- gl.dist.pop(possums.gl)
   gl.plot.heatmap(D2)
   }
   

Prints history of a genlight object

Description

Prints history of a genlight object

Usage

gl.print.history(x = NULL, history = NULL)

Arguments

Value

Prints a table with all history records. Currently the style cannot be changed.

Author(s)

Bernd Gruber (bugs? Post to https://groups.google.com/d/forum/dartr)

See Also

Other environment: gl.check.verbosity(), gl.check.wd(), gl.set.wd(), theme_dartR()

Examples

dartfile <- system.file('extdata','testset_SNPs_2Row.csv', package='dartR.data')
metadata <- system.file('extdata','testset_metadata.csv', package='dartR.data')
gl <- gl.read.dart(dartfile, ind.metafile = metadata, probar=FALSE)
gl2 <- gl.filter.callrate(gl, method='loc', threshold=0.9)
gl3 <- gl.filter.callrate(gl2, method='ind', threshold=0.95)
gl.print.history(gl3)

Calculates a similarity (distance) matrix for individuals on the proportion of shared alleles @family distance

Description

This script calculates an individual based distance matrix. It uses an C++ implementation, so package Rcpp needs to be installed and it is therefore really fast (once it has compiled the function after the first run).

Usage

gl.propShared(x)

Arguments

Value

A similarity matrix

Author(s)

Bernd Gruber (Post to https://groups.google.com/d/forum/dartr)

Examples

#takes some time at the first run of the function...

res <- gl.propShared(bandicoot.gl)
res[1:5,1:7] #show only a small part of the matrix

Randomly changes the allocation of 0's and 2's in a genlight object

Description

This function samples randomly half of the SNPs and re-codes, in the sampled SNP's, 0's by 2's.

This function samples randomly half of the SNPs and re-codes, in the sampled SNP's, 0's by 2's.

Usage

gl.randomize.snps(
  x,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = NULL,
  plot.file = NULL,
  plot.dir = NULL,
  verbose = NULL
)

gl.randomize.snps(
  x,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = NULL,
  plot.file = NULL,
  plot.dir = NULL,
  verbose = NULL
)

Arguments

Details

DArT calls the most common allele as the reference allele. In a genlight object, homozygous for the reference allele are coded with a '0' and homozygous for the alternative allele are coded with a '2'. This causes some distortions in visuals from time to time. If plot.display = TRUE, two smear plots (pre-randomisation and post-randomisation) are presented using a random subset of individuals (10) and loci (100) to provide an overview of the changes. Resultant ggplots are saved to the session's temporary directory.

DArT calls the most common allele as the reference allele. In a genlight object, homozygous for the reference allele are coded with a '0' and homozygous for the alternative allele are coded with a '2'. This causes some distortions in visuals from time to time. If plot.display = TRUE, two smear plots (pre-randomisation and post-randomisation) are presented using a random subset of individuals (10) and loci (100) to provide an overview of the changes. Resultant ggplots are saved to the session's temporary directory.

Value

Returns a genlight object with half of the loci re-coded.

Returns a genlight object with half of the loci re-coded.

Author(s)

Custodian: Luis Mijangos – Post to https://groups.google.com/d/forum/dartr

Examples

require("dartR.data")
res <- gl.randomize.snps(platypus.gl[1:5,1:5],verbose = 5)
gl <- gl.filter.monomorphs(testset.gl)
res <- gl.randomize.snps(gl,verbose = 5)

Reads PLINK data file into a genlight object

Description

This function imports PLINK data into a genlight object and append available metadata.

Usage

gl.read.PLINK(
  filename,
  ind.metafile = NULL,
  loc.metafile = NULL,
  plink.cmd = "plink",
  plink.path = "path",
  plink.flags = NULL,
  verbose = NULL
)

Arguments

Details

This function handles .ped or .bed file (with the associate files - e.g. .fam, .bim). However, if a .ped file is provided, PLINK needs to be installed and it is used to convert the .ped into a .bed, which is then converted into a genlight.

Additional metadata can be included passing .csv files. These will be appended to the existing metadata present in the PLINK files.

The locus metadata needs to be in a csv file with headings, with a mandatory column headed AlleleID corresponding exactly to the locus identity labels provided with the SNP data

Value

A genlight object with the SNP data and associated metadata included.

Author(s)

Custodian: Carlo Pacioni – Post to https://groups.google.com/d/forum/dartr

Reads SNP data from a csv file into a genlight object

Description

This script takes SNP genotypes from a csv file, combines them with individual and locus metrics and creates a genlight object. The SNP data need to be in one of two forms. SNPs can be coded 0 for homozygous reference, 2 for homozygous alternate, 1 for heterozygous, and NA for missing values; or the SNP data can be coded A/A, A/C, C/T, G/A etc, and -/- for missing data. In this format, the reference allele is the most frequent allele, as used by DArT. Other formats will throw an error. The SNP data need to be individuals as rows, labeled, and loci as columns, also labeled. If the orientation is individuals as columns and loci by rows, then set transpose=TRUE. The individual metrics need to be in a csv file, with headings, with a mandatory id column corresponding exactly to the individual identity labels provided with the SNP data and in the same order. The locus metadata needs to be in a csv file with headings, with a mandatory column headed AlleleID corresponding exactly to the locus identity labels provided with the SNP data and in the same order. Note that the locus metadata will be complemented by calculable statistics corresponding to those that would be provided by Diversity Arrays Technology (e.g. CallRate).

Usage

gl.read.csv(
  filename,
  transpose = FALSE,
  ind.metafile = NULL,
  loc.metafile = NULL,
  verbose = NULL
)

Arguments

Value

A genlight object with the SNP data and associated metadata included.

Author(s)

Custodian: Luis Mijangos – Post to https://groups.google.com/d/forum/dartr

See Also

Other io: gl.load(), gl.read.dart(), gl.read.fasta(), gl.read.silicodart(), gl.save(), gl.write.csv(), utils.read.dart()

Imports DArT data into dartR and converts it into a dartR genlight object

Description

This function is a wrapper function that allows you to convert your DArT file into a genlight object of class dartR.

Usage

gl.read.dart(
  filename,
  ind.metafile = NULL,
  recalc = TRUE,
  mono.rm = FALSE,
  nas = "-",
  topskip = NULL,
  lastmetric = NULL,
  covfilename = NULL,
  service.row = 1,
  plate.row = 3,
  probar = FALSE,
  verbose = NULL
)

Arguments

Details

The function will determine automatically if the data are in Diversity Arrays one-row csv format or two-row csv format. The first row of data is determined from the number of rows with an * in the first column. This can be alternatively specified with the topskip parameter. The DArT service code is added to the ind.metrics of the genlight object. The row containing the service code for each individual can be specified with the service.row parameter. #'The DArT plate well is added to the ind.metrics of the genlight object. The row containing the plate well for each individual can be specified with the plate.row parameter. If individuals have been deleted from the input file manually, then the locus metrics supplied by DArT will no longer be correct and some loci may be monomorphic. To accommodate this, set mono.rm and recalc to TRUE.

Value

A dartR genlight object that contains individual and locus metrics [if data were provided] and locus metrics [from a DArT report].

Author(s)

Custodian: Bernd Gruber (Post to https://groups.google.com/d/forum/dartr)

See Also

utils.read.dart

Other io: gl.load(), gl.read.csv(), gl.read.fasta(), gl.read.silicodart(), gl.save(), gl.write.csv(), utils.read.dart()

Examples

dartfile <- system.file('extdata','testset_SNPs_2Row.csv', package='dartR.data')
metadata <- system.file('extdata','testset_metadata.csv', package='dartR.data')
gl <- gl.read.dart(dartfile, ind.metafile = metadata, probar=TRUE)

Reads FASTA files and converts them to genlight object

Description

The following IUPAC Ambiguity Codes are taken as heterozygotes:

• M is heterozygote for AC and CA

• R is heterozygote for AG and GA

• W is heterozygote for AT and TA

• S is heterozygote for CG and GC

• Y is heterozygote for CT and TC

• K is heterozygote for GT and TG

The following IUPAC Ambiguity Codes are taken as missing data:

• V

• H

• D

• B

• N

The function can deal with missing data in individuals, e.g. when FASTA files have different number of individuals due to missing data. The allele with the highest frequency is taken as the reference allele. SNPs with more than two alleles are skipped.

Usage

gl.read.fasta(fasta.files, parallel = FALSE, n.cores = NULL, verbose = NULL)

Arguments

Details

Ambiguity characters are often used to code heterozygotes. However, using heterozygotes as ambiguity characters may bias many estimates. See more information in the link below: https://evodify.com/heterozygotes-ambiguity-characters/

Value

A genlight object.

Author(s)

Custodian: Luis Mijangos – Post to https://groups.google.com/d/forum/dartr

See Also

Other io: gl.load(), gl.read.csv(), gl.read.dart(), gl.read.silicodart(), gl.save(), gl.write.csv(), utils.read.dart()

Imports presence/absence data from SilicoDArT to genlight {agegenet} format (ploidy=1)

Description

DaRT provide the data as a matrix of entities (individual animals) across the top and attributes (P/A of sequenced fragment) down the side in a format that is unique to DArT. This program reads the data in to adegenet format for consistency with other programming activity. The script may require modification as DArT modify their data formats from time to time.

Usage

gl.read.silicodart(
  filename,
  ind.metafile = NULL,
  nas = "-",
  topskip = NULL,
  lastmetric = "Reproducibility",
  probar = TRUE,
  verbose = NULL
)

Arguments

Details

gl.read.silicodart() opens the data file (csv comma delimited) and skips the first n=topskip lines. The script assumes that the next line contains the entity labels (specimen ids) followed immediately by the SNP data for the first locus. It reads the presence/absence data into a matrix of 1s and 0s, and inputs the locus metadata and specimen metadata. The locus metadata comprises a series of columns of values for each locus including the essential columns of CloneID and the desirable variables Reproducibility and PIC. Refer to documentation provide by DArT for an explanation of these columns. The specimen metadata provides the opportunity to reassign specimens to populations, and to add other data relevant to the specimen. The key variables are id (specimen identity which must be the same and in the same order as the SilicoDArT file, each unique), pop (population assignment), lat (latitude, optional) and lon (longitude, optional). id, pop, lat, lon are the column headers in the csv file. Other optional columns can be added. The data matrix, locus names (forced to be unique), locus metadata, specimen names, specimen metadata are combined into a genind object. Refer to the documentation for {adegenet} for further details.

Value

An object of class genlight with ploidy set to 1, containing the presence/absence data, and locus and individual metadata.

Author(s)

Custodian: Bernd Gruber – Post to https://groups.google.com/d/forum/dartr

See Also

gl.read.dart

Other io: gl.load(), gl.read.csv(), gl.read.dart(), gl.read.fasta(), gl.save(), gl.write.csv(), utils.read.dart()

Examples

silicodartfile <- system.file('extdata','testset_SilicoDArT.csv', package='dartR.data')
metadata <- system.file('extdata',ind.metafile ='testset_metadata_silicodart.csv',
package='dartR.data')
testset.gs <- gl.read.silicodart(filename = silicodartfile, ind.metafile = metadata)

Converts a vcf file into a genlight object

Description

This function needs package vcfR, please install it.

Usage

gl.read.vcf(vcffile, ind.metafile = NULL, mode = "genotype", verbose = NULL)

Arguments

Details

The ind.metadata file needs to have very specific headings. First a heading called id. Here the ids have to match the ids in the dartR object. The following column headings are optional. pop: specifies the population membership of each individual. lat and lon specify spatial coordinates (in decimal degrees WGS1984 format). Additional columns with individual metadata can be imported (e.g. age, gender). Note also that this function checks to see if there are input of mode, missing input of mode will issue the user with an error. "Dosage" mode of this function assign ploidy levels as maximum copy number of alternate alleles. Please carefully check the data if "dosage" mode is used.

Value

A genlight object.

Author(s)

Bernd Gruber, Ching Ching Lau (Post to https://groups.google.com/d/forum/dartr)

Examples

## Not run: 
# read in vcf and convert to format as DArT data
obj <- gl.read.vcf(system.file('extdata/test.vcf', package='dartR'), 
                   ind.metafile = "metafile.csv")
# read in vcf and convert to format as dosage
obj <- gl.read.vcf(system.file('extdata/test.vcf', package='dartR'), 
                   ind.metafile = "metafile.csv", mode="dosage")

## End(Not run)

Assigns an individual metric as pop in a genlight {adegenet} object

Description

Individuals are assigned to populations based on the individual/sample/specimen metrics file (csv) used with gl.read.dart(). One might want to define the population structure in accordance with another classification, such as using an individual metric (e.g. sex, male or female). This script discards the current population assignments and replaces them with new population assignments defined by a specified individual metric. The function returns a genlight object with the new population assignments. Note that the original population assignments are lost.

Usage

gl.reassign.pop(x, as.pop, verbose = NULL)

Arguments

Value

A genlight object with the reassigned populations.

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

Other data manipulation: gl.define.pop(), gl.drop.ind(), gl.drop.loc(), gl.drop.pop(), gl.edit.recode.pop(), gl.impute(), gl.join(), gl.keep.ind(), gl.keep.loc(), gl.keep.pop(), gl.make.recode.ind(), gl.merge.pop(), gl.recode.ind(), gl.recode.pop(), gl.rename.pop(), gl.sample(), gl.sim.genotypes(), gl.sort(), gl.subsample.ind(), gl.subsample.loc()

Examples

# SNP data
   popNames(testset.gl)
   gl <- gl.reassign.pop(testset.gl, as.pop='sex',verbose=3)
   popNames(gl)
# Tag P/A data
   popNames(testset.gs)
   gs <- gl.reassign.pop(testset.gs, as.pop='sex',verbose=3)
   popNames(gs)
   

Recalculates locus metrics when individuals or populations are deleted from a genlight {adegenet} object @family environment

Description

When individuals,or populations, are deleted from a genlight object, the locus metrics no longer apply. For example, the Call Rate may be different considering the subset of individuals, compared with the full set. This script recalculates those affected locus metrics, namely, avgPIC, CallRate, freqHets, freqHomRef, freqHomSnp, OneRatioRef, OneRatioSnp, PICRef and PICSnp. Metrics that remain unaltered are RepAvg and TrimmedSeq as they are unaffected by the removal of individuals. The script optionally removes resultant monomorphic loci or loci with all values missing and deletes them (using gl.filter.monomorphs.r). The script returns a genlight object with the recalculated locus metadata.

Usage

gl.recalc.metrics(x, mono.rm = FALSE, verbose = NULL)

Arguments

Value

A genlight object with the recalculated locus metadata.

Author(s)

Custodian: Luis Mijangos (Post to https://groups.google.com/d/forum/dartr)

See Also

gl.filter.monomorphs

Examples

  gl <- gl.recalc.metrics(testset.gl, verbose=2)
  

Recodes individual (=specimen = sample) labels in a genlight object

Description

This function recodes individual labels and/or deletes individuals from a DaRT genlight SNP file based on a lookup table provided as a csv file.

Usage

gl.recode.ind(x, ind.recode, recalc = FALSE, mono.rm = FALSE, verbose = NULL)

Arguments

Details

Renaming individuals may be required when there have been errors in labeling arising in the process from sample to sequence files. There may be occasions where renaming individuals is required for preparation of figures. When caution needs to be exercised because of the potential for breaking the 'chain of evidence' associated with the samples, recoding individuals using a recode table (csv) can provide a durable record of the changes. The function works with genlight objects containing SNP genotypes and Tag P/A data (SilicoDArT). For SNP genotype data, the function, having deleted individuals, optionally identifies resultant monomorphic loci or loci with all values missing and deletes them. The script also optionally recalculates the locus metadata as appropriate. The optional deletion of monomorphic loci and the optional recalculation of locus statistics is not available for Tag P/A data (SilicoDArT). The script returns a dartR genlight object with the new individual names and the recalculated locus metadata.

Value

A genlight or genind object with the recoded and reduced data.

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

gl.filter.monomorphs for filtering monomorphs, gl.recalc.metrics for recalculating locus metrics, gl.recode.pop for recoding populations

Other data manipulation: gl.define.pop(), gl.drop.ind(), gl.drop.loc(), gl.drop.pop(), gl.edit.recode.pop(), gl.impute(), gl.join(), gl.keep.ind(), gl.keep.loc(), gl.keep.pop(), gl.make.recode.ind(), gl.merge.pop(), gl.reassign.pop(), gl.recode.pop(), gl.rename.pop(), gl.sample(), gl.sim.genotypes(), gl.sort(), gl.subsample.ind(), gl.subsample.loc()

Examples

  file <- system.file('extdata','testset_ind_recode.csv', package='dartR.data')
  gl <- gl.recode.ind(testset.gl, ind.recode=file, verbose=3)
  

Recodes population assignments in a genlight object

Description

This function recodes population assignments and/or deletes populations from a DaRT genlight object based on information provided in a csv population recode file.

Usage

gl.recode.pop(x, pop.recode, recalc = FALSE, mono.rm = FALSE, verbose = NULL)

Arguments

Details

Individuals are assigned to populations based on the specimen metadata data file (csv) used with gl.read.dart(). Recoding can be used to amalgamate populations or to selectively delete or retain populations. When caution needs to be exercised because of the potential for breaking the 'chain of evidence' associated with the samples, recoding individuals using a recode table (csv) can provide a durable record of the changes. The population recode file contains a list of populations taken from the genlight object as the first column of the csv file, and the new population assignments are located in the second column of the csv file. The keyword 'Delete' used as a new population assignment will result in the associated specimen being dropped from the dataset. The function works with genlight objects containing SNP genotypes and Tag P/A data (SilicoDArT). For SNP genotype data, the function, having deleted populations, optionally identifies resultant monomorphic loci or loci with all values missing and deletes them. The script also optionally recalculates the locus metadata as appropriate. The optional deletion of monomorphic loci and the optional recalculation of locus statistics is not available for Tag P/A data (SilicoDArT).

Value

A genlight object with the recoded and reduced data.

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

gl.filter.monomorphs

gl.recode.pop

Other data manipulation: gl.define.pop(), gl.drop.ind(), gl.drop.loc(), gl.drop.pop(), gl.edit.recode.pop(), gl.impute(), gl.join(), gl.keep.ind(), gl.keep.loc(), gl.keep.pop(), gl.make.recode.ind(), gl.merge.pop(), gl.reassign.pop(), gl.recode.ind(), gl.rename.pop(), gl.sample(), gl.sim.genotypes(), gl.sort(), gl.subsample.ind(), gl.subsample.loc()

Examples

  mfile <- system.file('extdata', 'testset_pop_recode.csv', package='dartR.data')
  nPop(testset.gl)
  gl <- gl.recode.pop(testset.gl, pop.recode=mfile, verbose=3)
 

Renames a population in a genlight object

Description

Individuals are assigned to populations based on the specimen metadata data file (csv) used with gl.read.dart(). This script renames a nominated population. The script returns a genlight object with the new population name.

Usage

gl.rename.pop(x, old = NULL, new = NULL, verbose = NULL)

Arguments

Value

A genlight object with the new population name.

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

Other data manipulation: gl.define.pop(), gl.drop.ind(), gl.drop.loc(), gl.drop.pop(), gl.edit.recode.pop(), gl.impute(), gl.join(), gl.keep.ind(), gl.keep.loc(), gl.keep.pop(), gl.make.recode.ind(), gl.merge.pop(), gl.reassign.pop(), gl.recode.ind(), gl.recode.pop(), gl.sample(), gl.sim.genotypes(), gl.sort(), gl.subsample.ind(), gl.subsample.loc()

Examples

   gl <- gl.rename.pop(testset.gl, old='EmsubRopeMata', new='Outgroup')
   

Report allelic richness per population from a genlight object

Description

This function needs package adegenet, please install it.

Usage

gl.report.allelerich(
  x,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.dir = NULL,
  plot.file = NULL,
  error.bar = "SD",
  verbose = 2
)

Arguments

Details

details

• raw allele count

• allelic richness

Value

A dataframe containing richness per site, richness per population, raw reference allele count, raw alternate allele count.

Author(s)

Ching Ching Lau (Post to https://groups.google.com/d/forum/dartr)

References

• El Mousadik, A., & Petit, R. J. (1996). High level of genetic differentiation for allelic richness among populations of the argan tree [Argania spinosa (L.) Skeels] endemic to Morocco. Theoretical and applied genetics, 92, 832-839.

Examples

 gl.report.allelerich(possums.gl)

Reports loci that are all NA across individuals and/or populations with all NA across loci

Description

This script reports loci or individuals with all calls missing (NA), from a genlight object.

A DArT dataset will not have loci for which the calls are scored all as missing (NA) for a particular individual, but such loci can arise rarely when populations or individuals are deleted. Similarly, a DArT dataset will not have individuals for which the calls are scored all as missing (NA) across all loci, but such individuals may sneak in to the dataset when loci are deleted. Retaining individual or loci with all NAs can cause issues for several functions.

Also, on occasions an analysis will require that there are some loci scored in each population. Setting by.pop=TRUE will result in removal of loci when they are all missing in any one population.

Note that loci that are missing for all individuals in a population are not imputed with method 'frequency' or 'HW'. Consider using the function gl.filter.allna with by.pop=TRUE.

Usage

gl.report.allna(x, by.pop = FALSE, verbose = NULL)

Arguments

Value

gl.report.allna

Author(s)

Author(s): Arthur Georges. Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

Other matched report: gl.filter.excess.het(), gl.report.callrate(), gl.report.hamming(), gl.report.locmetric(), gl.report.maf(), gl.report.overshoot(), gl.report.pa(), gl.report.rdepth(), gl.report.reproducibility(), gl.report.secondaries(), gl.report.taglength()

Other filter functions: gl.filter.allna(), gl.filter.hwe()

Examples

# SNP data
  result <- gl.report.allna(testset.gl, verbose=3)
# Tag P/A data
  result <- gl.report.allna(testset.gs, verbose=3)

Reports summary of base pair frequencies

Description

Calculates the frequencies of the four DNA nucleotide bases: adenine (A), cytosine (C), 'guanine (G) and thymine (T), and the frequency of transitions (Ts) and transversions (Tv) in a DArT genlight object.

Usage

gl.report.bases(
  x,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = NULL,
  plot.file = NULL,
  plot.dir = NULL,
  verbose = NULL,
  ...
)

Arguments

Details

The function checks first if trimmed sequences are included in the locus metadata (@other$loc.metrics$TrimmedSequence), and if so, tallies up the numbers of A, T, G and C bases. Only the reference state at the SNP locus is counted. Counts of transitions (Ts) and transversions (Tv) assume that there is no directionality, that is C->T is the same as T->C, because the reference state is arbitrary. For presence/absence data (SilicoDArT), it is not possible to count transversions or transitions or transversions/transitions ratio because the SNP data are not available, only a single sequence tag per locus. Only base frequencies are provided. A color vector can be obtained with gl.select.colors() and then passed to the function with the plot.colors parameter. If a plot.file is given, the ggplot arising from this function is saved as an "RDS" binary file using saveRDS(); can be reloaded with readRDS(). A file name must be specified for the plot to be saved. If a plot directory (plot.dir) is specified, the ggplot binary is saved to that directory; otherwise to the tempdir(). To avoid issues from inadvertent use of this function in an assignment statement, the function returns the genlight object unaltered.

Value

The unchanged genlight object

Author(s)

Author(s); Arthur Georges. Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

Other matched reports: gl.mahal.assign(), gl.report.factorloadings(), gl.report.fstat(), gl.report.monomorphs()

Examples

# SNP data
  out <- gl.report.bases(testset.gl)
  col <- gl.select.colors(select=c(6,1),palette=rainbow, verbose=0)
  out <- gl.report.bases(testset.gl,plot.colors=col)
# Tag P/A data
  out <- gl.report.bases(testset.gs)

Basic statistics for a genlight object

Description

Calculates basic statistics for a genlight object.

Usage

gl.report.basics(x, verbose = NULL)

Arguments

Author(s)

Custodian: Arthur Georges (Post to https://groups.google.com/d/forum/dartr)

See Also

Other unmatched report: gl.allele.freq(), gl.report.diversity(), gl.report.excess.het(), gl.report.heterozygosity(), gl.report.polyploid_heterozygosity()

Examples

gl.report.basics(platypus.gl)

Reports summary of Call Rate for loci or individuals

Description

SNP datasets generated by DArT have missing values primarily arising from failure to call a SNP because of a mutation at one or both of the restriction enzyme recognition sites. P/A datasets (SilicoDArT) have missing values because it was not possible to call whether a sequence tag was amplified or not.

Usage

gl.report.callrate(
  x,
  method = "loc",
  ind.to.list = 20,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = NULL,
  plot.dir = NULL,
  plot.file = NULL,
  bins = 50,
  verbose = NULL,
  ...
)

Arguments

Details

This function expects a genlight object, containing either SNP data or SilicoDArT (=presence/absence data).

Callrate is summarized by locus or by individual to allow sensible decisions on thresholds for filtering taking into consideration consequential loss of data. The summary is in the form of a tabulation and plots.

The table of quantiles is useful for deciding a threshold for subsequent filtering as it provides an indication of the percentages of loci that will be retained and lost.

In the case of method='ind', a list of individuals to be deleted is provided. To manage the screen output, this list is limited to ind.to.list individuals (or nInd(x)) whichever is the smaller.

To avoid issues from inadvertent use of this function in an assignment statement, the function returns the genlight object unaltered.

A color vector can be obtained with gl.select.colors() and then passed to the function with the plot.colors parameter.

If a plot.file is given, the ggplot arising from this function is saved as an "RDS" binary file using saveRDS(); can be reloaded with readRDS(). A file name must be specified for the plot to be saved.

If a plot directory (plot.dir) is specified, the ggplot binary is saved to that directory; otherwise to the tempdir().

Value

Returns unaltered genlight object

Author(s)

Author(s): Arthur Georges. Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

gl.filter.callrate

Other matched report: gl.filter.excess.het(), gl.report.allna(), gl.report.hamming(), gl.report.locmetric(), gl.report.maf(), gl.report.overshoot(), gl.report.pa(), gl.report.rdepth(), gl.report.reproducibility(), gl.report.secondaries(), gl.report.taglength()

Examples

 
# SNP data
  test.gl <- testset.gl[1:20,]
  gl.report.callrate(test.gl)
  gl.report.callrate(test.gl,method='ind')
  gl.report.callrate(test.gl,method='ind',plot.file="test")
  gl.report.callrate(test.gl,method='loc',by.pop=TRUE)
  gl.report.callrate(test.gl,method='loc',by.pop=TRUE,plot.file="test")
# Tag P/A data
  test.gs <- testset.gs[1:20,]
  gl.report.callrate(test.gs)
  gl.report.callrate(test.gs,method='ind')
  
  test.gl <- testset.gl[1:20,]
  gl.report.callrate(test.gl)
  

Calculates diversity indexes for SNPs

Description

This script takes a genlight object and calculates alpha and beta diversity for q = 0:2. Formulas are taken from Sherwin et al. 2017. The paper describes nicely the relationship between the different q levels and how they relate to population genetic processes such as dispersal and selection. The citation below also includes a link to a 3-minute video that explains, q, D and H.

Usage

gl.report.diversity(
  x,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors.pop = gl.colors("dis"),
  plot.dir = NULL,
  plot.file = NULL,
  table = "DH",
  verbose = NULL
)

Arguments

Details

For all indexes, the entropies (H) and corresponding effective numbers, i.e. Hill numbers (D), which reflect the number of needed entities to get the observed values, are calculated. In a nutshell, the alpha indexes between the different q-values should be similar if there is no deviation from expected allele frequencies and occurrences (e.g. all loci in HWE & equilibrium). If there is a deviation of an index, this links to a process causing it, such as dispersal, selection or strong drift. For a detailed explanation of all the indexes, we recommend resorting to the literature provided below. Error bars are +/- 1 standard deviation.

Function's output

If the function's parameter "table" = "DH" (the default value) is used, the output of the function is 20 tables.

The first two show the number of loci used. The name of each of the rest of the tables starts with three terms separated by underscores.

The first term refers to the q value (0 to 2). The q values identify different ways of summarising diversity (H): q=0 is simply the number of alleles per locus, with no information about their relative proportions; q=2 is the expected heterozygosity, ie the chance of drawing two different alleles at random from the population; q=1 is the Shannon measure of ‘surprise, relating to how likely it is that the next allele drawn will be one that has not been seen before (Sherwin et al 2017, 2021, and associated video).

The second term refers to whether it is the diversity measure (H) or its transformation to Hill numbers (D) The D value tells you how many equally-frequent alleles there would need to be to give the corresponding H-value (in the actual population) The D-values are all in units of numbers of alleles, so they can be plotted against the q-value to get a rich representation of the diversity (Box 1, Fig II in Sherwin et al 2017, 2021, and associated video).

The third term refers to whether the diversity is calculated within populations (alpha) or between populations (beta).

In the case of alpha diversity tables, standard deviations have their own table, which finishes with a fourth term: "sd".

In the case of beta diversity tables, standard deviations are in the upper triangle of the matrix and diversity values are in the lower triangle of the matrix.

Plotting

Plot colours can be set with gl.select.colors().

If plot.file is specified, plots are saved to the directory specified by the user, or the global default working directory set by gl.set.wd() or to the tempdir().

Examples of other themes that can be used can be consulted in

• https://ggplot2.tidyverse.org/reference/ggtheme.html and

• https://yutannihilation.github.io/allYourFigureAreBelongToUs/ggthemes/

Value

A list of entropy indexes for each level of q and equivalent numbers for alpha and beta diversity.

Author(s)

Bernd Gruber (Post to https://groups.google.com/d/forum/dartr), Contributors: William B. Sherwin, Alexander Sentinella

References

Sherwin, W.B., Chao, A., Johst, L., Smouse, P.E. (2017, 2021). Information Theory Broadens the Spectrum of Molecular Ecology and Evolution. TREE 32(12) 948-963. doi:10.1016/j.tree.2017.09.12 AND TREE 36:955-6 doi.org/10.1016/j.tree.2021.07.005 AND 3-Minute video: ars.els-cdn.com/content/image/1-s2.0-S0169534717302550-mmc2.mp4

See Also

Other unmatched report: gl.allele.freq(), gl.report.basics(), gl.report.excess.het(), gl.report.heterozygosity(), gl.report.polyploid_heterozygosity()

Examples

div <- gl.report.diversity(bandicoot.gl, table=FALSE)
div$zero_H_alpha
div$two_H_beta
names(div)

Report loci with excess of heterozygosity

Description

Calculates excess of heterozygosity in a genlight object

Usage

gl.report.excess.het(
  x,
  Yates = FALSE,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = NULL,
  plot.file = NULL,
  plot.dir = NULL,
  verbose = NULL
)

Arguments

Details

Loci with observed heterozygosity larger than 0.5 and expected heterozygosity would be indicated as excess (p-value <= 0.05). You can remove the loci with excess of heterozygosity from genlight object using gl.filter.excess.het

Function's output

If a plot.file is given, the ggplot arising from this function is saved as an "RDS" binary file using saveRDS(); can be reloaded with readRDS(). A file name must be specified for the plot to be saved.

If a plot directory (plot.dir) is specified, the ggplot binary is saved to that directory; otherwise to the tempdir().

Examples of other themes that can be used can be consulted in:

• https://ggplot2.tidyverse.org/reference/ggtheme.html and

• https://yutannihilation.github.io/allYourFigureAreBelongToUs/ggthemes/

A color vector can be obtained with gl.select.colors() and then passed to the function with the plot.colors parameter.

Value

1. Table with information of excessively-heterozygous loci
2. Two plots of heterozygosity of the loci before and after filtering (i.e. without excessively heterozygous loci).
3. A vector with the names of loci to be remove by gl.filter.excess.het

Author(s)

Jesús Castrejón-Figueroa, Diana A Robledo-Ruiz (Custodian: Ching Ching Lau) – Post to https://groups.google.com/d/forum/dartr

References

• https://github.com/drobledoruiz/conservation_genomics/tree/main/filter.excess.het

• Robledo‐Ruiz, D. A., Austin, L., Amos, J. N., Castrejón‐Figueroa, J., Harley, D. K., Magrath, M. J., Sunnucks, P. & Pavlova, A. (2023). Easy‐to‐use R functions to separate reduced‐representation genomic datasets into sex‐linked and autosomal loci, and conduct sex assignment. Molecular Ecology Resources.

See Also

gl.filter.excess.het

Other unmatched report: gl.allele.freq(), gl.report.basics(), gl.report.diversity(), gl.report.heterozygosity(), gl.report.polyploid_heterozygosity()

Examples

filtered.table <- gl.report.excess.het(x = LBP, Yates = TRUE)

Reports factor loadings for a PCA or PCoA

Description

Extracts the factor loadings from a glPCA object (generated by gl.pcoa) and plots their distribuion.

Usage

gl.report.factorloadings(
  pca,
  axis = 1,
  n.display = 15,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = NULL,
  plot.file = NULL,
  plot.dir = NULL,
  bins = 25,
  verbose = NULL,
  ...
)

Arguments

Details

The function extracts the factor loadings for a given axis from a PCA object generated by gl.pcoa and plots their magnitudes. Useful for identifying loci that load high for a given axis.

A color vector can be obtained with gl.select.colors() and then passed to the function with the plot.colors parameter.

Themes can be obtained from in

• https://ggplot2.tidyverse.org/reference/ggtheme.html and

• https://yutannihilation.github.io/allYourFigureAreBelongToUs/ggthemes/

If a plot.file is given, the ggplot arising from this function is saved as an "RDS" binary file using saveRDS(); can be reloaded with readRDS(). A file name must be specified for the plot to be saved. If a plot directory (plot.dir) is specified, the ggplot binary is saved to that directory; otherwise to the tempdir().

Value

The unchanged genlight object

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

Other matched reports: gl.mahal.assign(), gl.report.bases(), gl.report.fstat(), gl.report.monomorphs()

Examples

pca <- gl.pcoa(testset.gl)
gl.report.factorloadings(pca = pca)

Reports various statistics of genetic differentiation between populations with confidence intervals

Description

This function calculates four genetic differentiation between populations statistics (see the "Details" section for further information).

• Fst - Measure of the degree of genetic differentiation of subpopulations (Nei, 1987).

• Fstp - Unbiased (i.e. corrected for sampling error, see explanation below) Fst (Nei, 1987).

• Dest - Jost’s D (Jost, 2008).

• Gst_H - Gst standardized by the maximum level that it can obtain for the observed amount of genetic variation (Hedrick 2005).

Sampling errors arise because allele frequencies in our samples differ from those in the subpopulations from which they were taken (Holsinger, 2012).

Confidence Intervals are obtained using bootstrapping.

Usage

gl.report.fstat(
  x,
  nboots = 0,
  conf = 0.95,
  CI.type = "bca",
  ncpus = 1,
  plot.stat = "Fstp",
  plot.display = TRUE,
  palette.divergent = gl.colors("div"),
  font.size = 0.5,
  plot.dir = NULL,
  plot.file = NULL,
  verbose = NULL,
  ...
)

Arguments

Details

Even though Fst and its relatives can predict evolutionary processes (Holsinger & Weir, 2009), they are not true measures of genetic differentiation in the sense that they are dependent on the diversity within populations (Meirmans & Hedrick, 2011), the number of populations analysed (Alcala & Rosenberg, 2017) and are not monotonic (Sherwin et al., 2017). Recent approaches have been developed to accommodate these mathematical restrictions (G'ST; "Gst_H"; Hedrick, 2005, and Jost's D; "Dest"; Jost, 2008). More recently, novel approaches based on information theory (Mutual Information; Sherwin et al., 2017) and allele frequencies (Allele Frequency Difference; Berner, 2019) have distinct properties that make them valuable resources to interpret genetic differentiation between populations.

Note that each measure of genetic differentiation has advantages and drawbacks, and the decision of using a particular measure is usually based on the research question.

Statistics calculated

The equations used to calculate the statistics are shown below.

• Ho - Unbiased estimate of observed heterozygosity across subpopulations (Nei, 1987, pp. 164, eq. 7.38) is calculated as:

  

  where Pkii represents the proportion of homozygote ii for allele i in individual k and s represents the number of subpopulations.

• Hs - Unbiased estimate of the expected heterozygosity under Hardy-Weinberg equilibrium across subpopulations (Nei, 1987, pp. 164, eq. 7.39) is calculated as:

  

  where ñ is the harmonic mean of nk (the number of individuals in each subpopulation), pki is the proportion (sometimes misleadingly called frequency) of allele i in subpopulation k.

• Ht - Heterozygosity for the total population (Nei, 1987, pp. 164, eq. 7.40) is calculated as:

  

• Dst - The average allele frequency differentiation between populations (Nei, 1987, pp. 163) is calculated as:

  

• Htp - Unbiased estimate of Heterozygosity for the total population (Nei, 1987, pp. 165) is calculated as:

  

• Dstp - Unbiased estimate of the average allele frequency differentiation between populations (Nei, 1987, pp. 165) is calculated as:

  

• Fst - Measure of the extent of genetic differentiation of subpopulations (Nei, 1987, pp. 162, eq. 7.34) is calculated as:

  

• Fstp - Unbiased measure of the extent of genetic differentiation of subpopulations (Nei, 1987, pp. 163, eq. 7.36) is calculated as:

  

• Dest - Jost’s D (Jost, 2008, eq. 12) is calculated as:

  

• Gst-max - The maximum level that Gst can obtain for the observed amount of genetic variation (Hedrick 2005, eq. 4a) is calculated as:

  

• Gst-H - Gst standardized by the maximum level that it can obtain for the observed amount of genetic variation (Hedrick 2005, eq. 4b) is calculated as:

  

Confidence Intervals

The uncertainty of a parameter, in this case the mean of the statistic, can be summarised by a confidence interval (CI) which includes the true parameter value with a specified probability (i.e. confidence level; the parameter "conf" in this function).

In this function, CI are obtained using Bootstrap which is an inference method that samples with replacement the data (i.e. loci) and calculates the statistics every time.

This function uses the function boot (package boot) to perform the bootstrap replicates and the function boot.ci (package boot) to perform the calculations for the CI.

Four different types of nonparametric CI can be calculated (parameter "CI.type" in this function):

• First order normal approximation interval ("norm").

• Basic bootstrap interval ("basic").

• Bootstrap percentile interval ("perc").

• Adjusted bootstrap percentile interval ("bca").

The studentized bootstrap interval ("stud") was not included in the CI types because it is computationally intensive, it may produce estimates outside the range of plausible values and it has been found to be erratic in practice, see for example the "Studentized (t) Intervals" section in:

https://www.r-bloggers.com/2019/09/understanding-bootstrap-confidence-interval-output-from-the-r-boot-package/

Nice tutorials about the different types of CI can be found in:

https://www.datacamp.com/tutorial/bootstrap-r

and

https://www.r-bloggers.com/2019/09/understanding-bootstrap-confidence-interval-output-from-the-r-boot-package/

Efron and Tibshirani (1993, p. 162) and Davison and Hinkley (1997, p. 194) suggest that the number of bootstrap replicates should be between 1000 and 2000.

It is important to note that unreliable confidence intervals will be obtained if too few number of bootstrap replicates are used. Therefore, the function boot.ci will throw warnings and errors if bootstrap replicates are too few. Consider increasing then number of bootstrap replicates to at least 200.

The "bca" interval is often cited as the best for theoretical reasons, however it may produce unstable results if the bootstrap distribution is skewed or has extreme values. For example, you might get the warning "extreme order statistics used as endpoints" or the error "estimated adjustment 'a' is NA". In this case, you may want to use more bootstrap replicates or a different method or check your data for outliers.

The error "estimated adjustment 'w' is infinite" means that the estimated adjustment ‘w’ for the "bca" interval is infinite, which can happen when the empirical influence values are zero or very close to zero. This can be caused by various reasons, such as:

The number of bootstrap replicates is too small, the statistic of interest is constant or nearly constant across the bootstrap samples, the data contains outliers or extreme values.

You can try some possible solutions, such as:

Increasing the number of bootstrap replicates, using a different type of bootstrap confidence interval or removing or transforming the outliers or extreme values.

Plotting

The plot can be customised by including any parameter(s) from the function heatmap.2 (package gplots).

For the color palette you could try for example:

> library(viridis)

> res <- gl.report.fstat(platypus.gl, palette.divergent = viridis)

If a plot.file is given, the plot arising from this function is saved as an "RDS" binary file using the function saveRDS (package base); can be reloaded with function readRDS (package base). A file name must be specified for the plot to be saved.

If a plot directory (plot.dir) is specified, the gplot binary is saved to that directory; otherwise to the tempdir().

Your plot might not shown in full because your 'Plots' pane is too small (in RStudio). Increase the size of the 'Plots' pane before running the function. Alternatively, use the parameter 'plot.file' to save the plot to a file.

Parallelisation

If the parameter ncpus > 1, parallelisation is enabled. In Windows, parallel computing employs a "socket" approach that starts new copies of R on each core. POSIX systems, on the other hand (Mac, Linux, Unix, and BSD), utilise a "forking" approach that replicates the whole current version of R and transfers it to a new core.

Opening and terminating R sessions in each core involves a significant amount of processing time, therefore parallelisation in Windows machines is only quicker than not using parallelisation when nboots > 1000-2000.

Value

Two lists, the first list contains matrices with genetic statistics taken pairwise by population, the second list contains tables with the genetic statistics for each pair of populations. If nboots > 0, tables with the four statistics calculated with Low Confidence Intervals (LCI) and High Confidence Intervals (HCI).

Author(s)

Custodian: Luis Mijangos – Post to https://groups.google.com/d/forum/dartr

References

• Alcala, N., & Rosenberg, N. A. (2017). Mathematical constraints on FST: Biallelic markers in arbitrarily many populations. Genetics (206), 1581-1600.

• Berner, D. (2019). Allele frequency difference AFD–an intuitive alternative to FST for quantifying genetic population differentiation. Genes, 10(4), 308.

• Davison AC, Hinkley DV (1997). Bootstrap Methods and their Application. Cambridge University Press: Cambridge.

• Efron, B. (1979). Bootstrap methods: Another look at the jackknife. Annals of Statistics 7, 1–26.

• Efron B, Tibshirani RJ (1993). An Introduction to the Bootstrap. Chapman and Hall: London.

• Hedrick, P. W. (2005). A standardized genetic differentiation measure. Evolution, 59(8), 1633-1638.

• Holsinger, K. E. (2012). Lecture notes in population genetics.

• Holsinger, K. E., & Weir, B. S. (2009). Genetics in geographically structured populations: defining, estimating and interpreting FST. Nature Reviews Genetics, 10(9), 639- 650.

• Jost, L. (2008). GST and its relatives do not measure differentiation. Molecular Ecology, 17(18), 4015-4026.

• Meirmans, P. G., & Hedrick, P. W. (2011). Assessing population structure: FST and related measures. Molecular Ecology Resources, 11(1), 5-18.

• Nei, M. (1987). Molecular evolutionary genetics: Columbia University Press.

• Sherwin, W. B., Chao, A., Jost, L., & Smouse, P. E. (2017). Information theory broadens the spectrum of molecular ecology and evolution. Trends in Ecology & Evolution, 32(12), 948-963.

See Also

Other matched reports: gl.mahal.assign(), gl.report.bases(), gl.report.factorloadings(), gl.report.monomorphs()

Examples

res <- gl.report.fstat(platypus.gl)

Calculates the pairwise Hamming distance between DArT trimmed DNA sequences

Description

Hamming distance is calculated as the number of base differences between two sequences which can be expressed as a count or a proportion. Typically, it is calculated between two sequences of equal length. In the context of DArT trimmed sequences, which differ in length but which are anchored to the left by the restriction enzyme recognition sequence, it is sensible to compare the two trimmed sequences starting from immediately after the common recognition sequence and terminating at the last base of the shorter sequence.

Usage

gl.report.hamming(
  x,
  rs = 5,
  threshold = 3,
  tag.length = 69,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = NULL,
  plot.dir = NULL,
  plot.file = NULL,
  probar = FALSE,
  verbose = NULL
)

Arguments

Details

The function gl.filter.hamming will filter out one of two loci if their Hamming distance is less than a specified percentage Hamming distance can be computed by exploiting the fact that the dot product of two binary vectors x and (1-y) counts the corresponding elements that are different between x and y. This approach can also be used for vectors that contain more than two possible values at each position (e.g. A, C, T or G). If a pair of DNA sequences are of differing length, the longer is truncated. The algorithm is that of Johann de Jong https://johanndejong.wordpress.com/2015/10/02/faster-hamming-distance-in-r-2/ as implemented in utils.hamming If plot.file is specified, plots are saved to the directory specified by the user, or the global default working directory set by gl.set.wd() or to the tempdir(). Examples of other themes that can be used can be consulted in

• https://ggplot2.tidyverse.org/reference/ggtheme.html and

• https://yutannihilation.github.io/allYourFigureAreBelongToUs/ggthemes/

Value

Returns unaltered genlight object

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

gl.filter.hamming

Other matched report: gl.filter.excess.het(), gl.report.allna(), gl.report.callrate(), gl.report.locmetric(), gl.report.maf(), gl.report.overshoot(), gl.report.pa(), gl.report.rdepth(), gl.report.reproducibility(), gl.report.secondaries(), gl.report.taglength()

Examples

 
gl.report.hamming(testset.gl[,1:100])
gl.report.hamming(testset.gs[,1:100])

#' # SNP data
test <- platypus.gl
test <- gl.subsample.loc(platypus.gl,n=50)
result <- gl.report.hamming(test, verbose=3)
result <- gl.report.hamming(test, plot.file="ttest", verbose=3)

Reports observed, expected and unbiased heterozygosities and FIS (inbreeding coefficient) by population or by individual from SNP data

Description

Calculates the observed, expected and unbiased expected (i.e. corrected for sample size) heterozygosities and FIS (inbreeding coefficient) for each population or the observed heterozygosity for each individual in a genlight object.

Usage

gl.report.heterozygosity(
  x,
  method = "pop",
  n.invariant = 0,
  subsample.pop = FALSE,
  n.limit = 10,
  nboots = 0,
  boot.method = "ind",
  conf = 0.95,
  CI.type = "bca",
  ncpus = 1,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors.pop = gl.colors("dis"),
  plot.colors.ind = gl.colors(2),
  plot.file = NULL,
  plot.dir = NULL,
  error.bar = "SD",
  verbose = NULL
)

Arguments

Details

Observed heterozygosity for a population takes the proportion of heterozygous loci for each individual and averages it over all individuals in that population. The calculations take into account missing values.

Expected heterozygosity for a population takes the expected proportion of heterozygotes, that is, expected under Hardy-Weinberg equilibrium, for each locus, then averages this across the loci for an average estimate for the population.

The unbiased expected heterozygosity is calculated using the correction for sample size following equation 2 from Nei 1978.

Accuracy of all heterozygosity estimates is affected by small sample sizes, and so is their comparison between populations or repeated analysis. Expected heterozygosities are less affected because their calculations are based on allele frequencies while observed heterozygosities are strongly susceptible to sampling effects when the sample size is small.

Observed heterozygosity for individuals is calculated as the proportion of loci that are heterozygous for that individual.

Finally, the loci that are invariant across all individuals in the dataset (that is, across populations), is typically unknown. This can render estimates of heterozygosity analysis specific, and so it is not valid to compare such estimates across species or even across different analyses (see Schimdt et al 2021). This is a similar problem faced by microsatellites. If you have an estimate of the number of invariant sequence tags (loci) in your data, such as provided by gl.report.secondaries, you can specify it with the n.invariant parameter to standardize your estimates of heterozygosity. This is called autosomal heterozygosities by Schimddt et al (2021).

NOTE: It is important to realise that estimation of adjusted (autosomal) heterozygosity requires that secondaries not to be removed.

Heterozygosities and FIS (inbreeding coefficient) are calculated by locus within each population using the following equations, and then averaged across all loci:

• Observed heterozygosity (Ho) = number of heterozygotes / n_Ind, where n_Ind is the number of individuals without missing data for that locus.

• Observed heterozygosity adjusted (Ho.adj) <- Ho * n_Loc / (n_Loc + n.invariant), where n_Loc is the number of loci that do not have all missing data and n.invariant is an estimate of the number of invariant loci to adjust heterozygosity.

• Expected heterozygosity (He) = 1 - (p^2 + q^2), where p is the frequency of the reference allele and q is the frequency of the alternative allele.

• Expected heterozygosity adjusted (He.adj) = He * n_Loc / (n_Loc + n.invariant)

• Unbiased expected heterozygosity (uHe) = He * (2 * n_Ind / (2 * n_Ind - 1))

• Inbreeding coefficient (FIS) = 1 - Ho / uHe

Function's output

Output for method='pop' is an ordered barchart of observed heterozygosity, unbiased expected heterozygosity and FIS (Inbreeding coefficient) across populations together with a table of mean observed and expected heterozygosities and FIS by population and their respective standard deviations (SD). In the output, it is also reported by population: the number of loci used to estimate heterozygosity (n.Loc), the number of polymorphic loci (polyLoc), the number of monomorphic loci (monoLoc) and loci with all missing data (all_NALoc).

Output for method='ind' is a histogram and a boxplot of heterozygosity across individuals.

If a plot.file is given, the ggplot arising from this function is saved as an "RDS" binary file using saveRDS(); can be reloaded with readRDS(). A file name must be specified for the plot to be saved.

If a plot directory (plot.dir) is specified, the ggplot binary is saved to that directory; otherwise to the tempdir().

Examples of other themes that can be used can be consulted in:

• https://ggplot2.tidyverse.org/reference/ggtheme.html and

• https://yutannihilation.github.io/allYourFigureAreBelongToUs/ggthemes/

Subsampling populations

To test the effect of five population sample sizes (n = 10, 5, 4, 3, 2) on observed heterozygosity estimates, the function subsamples individuals, without replacement. The subsampling is repeated 10 times for each sample size n. This approach is not an implementation of Schmidt et al (2021). Please refer to this paper for additional complexities in estimating heterozygosity using SNP data.

Error bars

The best method for presenting or assessing genetic statistics depends on the type of data you have and the specific questions you're trying to answer. Here's a brief overview of when you might use each method:

1. Confidence Intervals ("CI"):

- Usage: Often used to convey the precision of an estimate.

- Advantage: Confidence intervals give a range in which the true parameter (like a population mean) is likely to fall, given the data and a specified probability (like 95%).

- In Context: For genetic statistics, if you're estimating a parameter, a 95% CI gives you a range in which you're 95% confident the true parameter lies.

2. Standard Deviation ("SD"):

- Usage: Describes the amount of variation from the average in a set of data.

- Advantage: Allows for an understanding of the spread of individual data points around the mean.

- In Context: If you're looking at the distribution of a quantitative trait (like height) in a population with a particular genotype, the SD can describe how much individual heights vary around the average height.

3. Standard Error ("SE"):

- Usage: Describes the precision of the sample mean as an estimate of the population mean.

- Advantage: Smaller than the SD in large samples; it takes into account both the SD and the sample size.

- In Context: If you want to know how accurately your sample mean represents the population mean, you'd look at the SE.

Recommendation:

- If you're trying to convey the precision of an estimate, confidence intervals are very useful.

- For understanding variability within a sample, standard deviation is key.

- To see how well a sample mean might estimate a population mean, consider the standard error.

In practice, geneticists often use a combination of these methods to analyze and present their data, depending on their research questions and the nature of the data.

Confidence Intervals

The uncertainty of a parameter, in this case the mean of the statistic, can be summarised by a confidence interval (CI) which includes the true parameter value with a specified probability (i.e. confidence level; the parameter "conf" in this function).

In this function, CI are obtained using Bootstrap which is an inference method that samples with replacement the data (i.e. loci) and calculates the statistics every time.

This function uses the function boot (package boot) to perform the bootstrap replicates and the function boot.ci (package boot) to perform the calculations for the CI.

Four different types of nonparametric CI can be calculated (parameter "CI.type" in this function):

• First order normal approximation interval ("norm").

• Basic bootstrap interval ("basic").

• Bootstrap percentile interval ("perc").

• Adjusted bootstrap percentile interval ("bca").

The studentized bootstrap interval ("stud") was not included in the CI types because it is computationally intensive, it may produce estimates outside the range of plausible values and it has been found to be erratic in practice, see for example the "Studentized (t) Intervals" section in:

https://www.r-bloggers.com/2019/09/understanding-bootstrap-confidence-interval-output-from-the-r-boot-package/

Nice tutorials about the different types of CI can be found in:

https://www.datacamp.com/tutorial/bootstrap-r

and

https://www.r-bloggers.com/2019/09/understanding-bootstrap-confidence-interval-output-from-the-r-boot-package/

Efron and Tibshirani (1993, p. 162) and Davison and Hinkley (1997, p. 194) suggest that the number of bootstrap replicates should be between 1000 and 2000.

It is important to note that unreliable confidence intervals will be obtained if too few number of bootstrap replicates are used. Therefore, the function boot.ci will throw warnings and errors if bootstrap replicates are too few. Consider increasing the number of bootstrap replicates to at least 200.

The "bca" interval is often cited as the best for theoretical reasons, however it may produce unstable results if the bootstrap distribution is skewed or has extreme values. For example, you might get the warning "extreme order statistics used as endpoints" or the error "estimated adjustment 'a' is NA". In this case, you may want to use more bootstrap replicates or a different method or check your data for outliers.

The error "estimated adjustment 'w' is infinite" means that the estimated adjustment ‘w’ for the "bca" interval is infinite, which can happen when the empirical influence values are zero or very close to zero. This can be caused by various reasons, such as:

The number of bootstrap replicates is too small, the statistic of interest is constant or nearly constant across the bootstrap samples, the data contains outliers or extreme values.

You can try some possible solutions, such as:

Increasing the number of bootstrap replicates, using a different type of bootstrap confidence interval or removing or transforming the outliers or extreme values.

Parallelisation

If the parameter ncpus > 1, parallelisation is enabled. In Windows, parallel computing employs a "socket" approach that starts new copies of R on each core. POSIX systems, on the other hand (Mac, Linux, Unix, and BSD), utilise a "forking" approach that replicates the whole current version of R and transfers it to a new core.

Opening and terminating R sessions in each core involves a significant amount of processing time, therefore parallelisation in Windows machines is only quicker than not using parallelisation when nboots > 1000-2000.

Value

A dataframe containing population labels, heterozygosities, FIS, their standard deviations and sample sizes.

Author(s)

Custodian: Luis Mijangos (Post to https://groups.google.com/d/forum/dartr)

References

• Nei, M. (1978). Estimation of average heterozygosity and genetic distance from a small number of individuals. Genetics, 89(3), 583-590.

• Schmidt, T. L., Jasper, M. E., Weeks, A. R., & Hoffmann, A. A. (2021). Unbiased population heterozygosity estimates from genome‐wide sequence data. Methods in Ecology and Evolution, 12(10), 1888-1898.

See Also

gl.filter.heterozygosity

Other unmatched report: gl.allele.freq(), gl.report.basics(), gl.report.diversity(), gl.report.excess.het(), gl.report.polyploid_heterozygosity()

Examples

 
require("dartR.data")
df <- gl.report.heterozygosity(platypus.gl)
df <- gl.report.heterozygosity(platypus.gl,method='ind')
n.inv <- gl.report.secondaries(platypus.gl)
gl.report.heterozygosity(platypus.gl, n.invariant = n.inv[7, 2])
gl.report.heterozygosity(platypus.gl, subsample.pop = TRUE)

df <- gl.report.heterozygosity(platypus.gl)

Reports departure from Hardy-Weinberg proportions

Description

Calculates the probabilities of agreement with H-W proportions based on observed frequencies of reference homozygotes, heterozygotes and alternate homozygotes.

Usage

gl.report.hwe(
  x,
  subset = "each",
  method_sig = "Exact",
  multi_comp = FALSE,
  multi_comp_method = "BY",
  alpha_val = 0.05,
  pvalue_type = "midp",
  cc_val = 0.5,
  sig_only = TRUE,
  min_sample_size = 5,
  plot.out = TRUE,
  plot_colors = gl.colors("2c"),
  max_plots = 4,
  save2tmp = FALSE,
  verbose = NULL
)

Arguments

Details

There are several factors that can cause deviations from Hardy-Weinberg proportions including: mutation, finite population size, selection, population structure, age structure, assortative mating, sex linkage, nonrandom sampling and genotyping errors. Therefore, testing for Hardy-Weinberg proportions should be a process that involves a careful evaluation of the results, a good place to start is Waples (2015). Note that tests for H-W proportions are only valid if there is no population substructure (assuming random mating) and have sufficient power only when there is sufficient sample size (n individuals > 15). Populations can be defined in three ways:

• Merging all populations in the dataset using subset = 'all'.

• Within each population separately using: subset = 'each'.

• Within selected populations using for example: subset = c('pop1','pop2').

Two different statistical methods to test for deviations from Hardy Weinberg proportions:

• The classical chi-square test (method_sig='ChiSquare') based on the function HWChisq of the R package HardyWeinberg. By default a continuity correction is applied (cc_val=0.5). The continuity correction can be turned off (by specifying cc_val=0), for example in cases of extreme allele frequencies in which the continuity correction can lead to excessive type 1 error rates.

• The exact test (method_sig='Exact') based on the exact calculations contained in the function HWExactStats of the R package HardyWeinberg, and described in Wigginton et al. (2005). The exact test is recommended in most cases (Wigginton et al., 2005). Three different methods to estimate p-values (pvalue_type) in the Exact test can be used:

  • 'dost' p-value is computed as twice the tail area of a one-sided test.

  • 'selome' p-value is computed as the sum of the probabilities of all samples less or equally likely as the current sample.

  • 'midp', p-value is computed as half the probability of the current sample + the probabilities of all samples that are more extreme.

  The standard exact p-value is overly conservative, in particular for small minor allele frequencies. The mid p-value ameliorates this problem by bringing the rejection rate closer to the nominal level, at the price of occasionally exceeding the nominal level (Graffelman & Moreno, 2013).

Correction for multiple tests can be applied using the following methods based on the function p.adjust:

• 'holm' is also known as the sequential Bonferroni technique (Rice, 1989). This method has a greater statistical power than the standard Bonferroni test, however this method becomes very stringent when many tests are performed and many real deviations from the null hypothesis can go undetected (Waples, 2015).

• 'hochberg' based on Hochberg, 1988.

• 'hommel' based on Hommel, 1988. This method is more powerful than Hochberg's, but the difference is usually small.

• 'bonferroni' in which p-values are multiplied by the number of tests. This method is very stringent and therefore has reduced power to detect multiple departures from the null hypothesis.

• 'BH' based on Benjamini & Hochberg, 1995.

• 'BY' based on Benjamini & Yekutieli, 2001.

The first four methods are designed to give strong control of the family-wise error rate. The last two methods control the false discovery rate (FDR), the expected proportion of false discoveries among the rejected hypotheses. The false discovery rate is a less stringent condition than the family-wise error rate, so these methods are more powerful than the others, especially when number of tests is large. The number of tests on which the adjustment for multiple comparisons is the number of populations times the number of loci. Ternary plots 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. See Graffelman & Morales-Camarena (2008) for further details. Ternary plots are based on the function HWTernaryPlot from 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.

Value

A dataframe containing loci, counts of reference SNP homozygotes, heterozygotes and alternate SNP homozygotes; probability of departure from H-W proportions, per locus significance with and without correction for multiple comparisons and the number of population where the same locus is significantly out of HWE.

Author(s)

Custodian: Luis Mijangos – Post to https://groups.google.com/d/forum/dartr

References

• Benjamini, Y., and Yekutieli, D. (2001). The control of the false discovery rate in multiple testing under dependency. Annals of Statistics, 29, 1165–1188.

• Graffelman, J. (2015). Exploring Diallelic Genetic Markers: The Hardy Weinberg Package. Journal of Statistical Software 64:1-23.

• Graffelman, J. & Morales-Camarena, J. (2008). Graphical tests for Hardy-Weinberg equilibrium based on the ternary plot. Human Heredity 65:77-84.

• Graffelman, J., & Moreno, V. (2013). The mid p-value in exact tests for Hardy-Weinberg equilibrium. Statistical applications in genetics and molecular biology, 12(4), 433-448.

• Hochberg, Y. (1988). A sharper Bonferroni procedure for multiple tests of significance. Biometrika, 75, 800–803.

• Hommel, G. (1988). A stagewise rejective multiple test procedure based on a modified Bonferroni test. Biometrika, 75, 383–386.

• Rice, W. R. (1989). Analyzing tables of statistical tests. Evolution, 43(1), 223-225.

• Waples, R. S. (2015). Testing for Hardy–Weinberg proportions: have we lost the plot?. Journal of heredity, 106(1), 1-19.

• Wigginton, J.E., Cutler, D.J., & Abecasis, G.R. (2005). A Note on Exact Tests of Hardy-Weinberg Equilibrium. American Journal of Human Genetics 76:887-893.

See Also

gl.filter.hwe

Calculates pairwise population based Linkage Disequilibrium across all loci using the specified number of cores @family matched report

Description

This function is implemented in a parallel fashion to speed up the process. There is also the ability to restart the function if crashed by specifying the chunk file names or restarting the function exactly in the same way as in the first run. This is implemented because sometimes, due to connectivity loss between cores, the function may crash half way. Before running the function, it is advisable to use the function gl.filter.allna to remove loci with all missing data.

Usage

gl.report.ld(
  x,
  name = NULL,
  save = TRUE,
  outpath = tempdir(),
  nchunks = 2,
  ncores = 1,
  chunkname = NULL,
  probar = FALSE,
  verbose = NULL
)

Arguments

Value

Returns calculation of pairwise LD across all loci between subpopulations. This functions uses if specified many cores on your computer to speed up. And if save is used can restart (if save=TRUE is used) with the same command starting where it crashed. The final output is a data frame that holds all statistics of pairwise LD between loci. (See ?LD in package genetics for details).

Author(s)

Bernd Gruber (Post to https://groups.google.com/d/forum/dartr)

Calculates pairwise linkage disequilibrium by population

Description

This function calculates pairwise linkage disequilibrium (LD) by population using the function ld from package snpStats. If SNPs are not mapped to a reference genome, the parameter ld.max.pairwise should be set as NULL (the default). In this case, the function will assign the same chromosome ("1") to all the SNPs in the dataset and assign a sequence from 1 to n loci as the position of each SNP. The function will then calculate LD for all possible SNP pair combinations. If SNPs are mapped to a reference genome, the parameter ld.max.pairwise should be filled out (i.e. not NULL). In this case, the information for SNP's position should be stored in the genlight accessor "@position" and the SNP's chromosome name in the accessor "@chromosome" (see examples). The function will then calculate LD within each chromosome and for all possible SNP pair combinations within a distance of ld.max.pairwise.

Usage

gl.report.ld.map(
  x,
  ld.max.pairwise = 1e+06,
  maf = 0.05,
  ld.stat = "R.squared",
  ind.limit = 10,
  stat.keep = "AvgPIC",
  ld.threshold.pops = 0.2,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.file = NULL,
  plot.dir = NULL,
  histogram.colors = NULL,
  boxplot.colors = NULL,
  bins = 50,
  verbose = NULL
)

Arguments

Details

This function reports LD between SNP pairs by population. The function gl.filter.ld filters out the SNPs in LD using as input the results of gl.report.ld.map. The actual number of SNPs to be filtered out depends on the parameters set in the function gl.filter.ld. Boxplots of LD by population and a histogram showing LD frequency are presented.

Value

A dataframe with information for each SNP pair in LD.

Author(s)

Custodian: Luis Mijangos – Post to https://groups.google.com/d/forum/dartr

See Also

gl.filter.ld

Other graphics: gl.colors(), gl.map.interactive(), gl.plot.heatmap(), gl.select.colors(), gl.select.shapes(), gl.smearplot(), gl.tree.nj()

Examples

require("dartR.data")
x <- platypus.gl
x <- gl.filter.callrate(x,threshold = 1)
x <- gl.filter.monomorphs(x)
x$position <- x$other$loc.metrics$ChromPos_Platypus_Chrom_NCBIv1
x$chromosome <- as.factor(x$other$loc.metrics$Chrom_Platypus_Chrom_NCBIv1)
ld_res <- gl.report.ld.map(x,ld.max.pairwise = 10000000)

Reports summary of the slot $other$loc.metrics

Description

This function reports summary statistics (mean, minimum, average, quantiles), histograms and boxplots for any loc.metric with numeric values (stored in $other$loc.metrics) to assist the decision of choosing thresholds for the filter function gl.filter.locmetric.

Usage

gl.report.locmetric(
  x,
  metric,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = NULL,
  plot.dir = NULL,
  plot.file = NULL,
  verbose = NULL
)

Arguments

Details

The function gl.filter.locmetric will filter out the loci with a locmetric value below a specified threshold. The fields that are included in dartR, and a short description, are found below. Optionally, the user can also set his/her own field by adding a vector into $other$loc.metrics as shown in the example. You can check the names of all available loc.metrics via: names(gl$other$loc.metrics).

• SnpPosition - position (zero is position 1) in the sequence tag of the defined SNP variant base.

• CallRate - proportion of samples for which the genotype call is non-missing (that is, not '-' ).

• OneRatioRef - proportion of samples for which the genotype score is 0.

• OneRatioSnp - proportion of samples for which the genotype score is 2.

• FreqHomRef - proportion of samples homozygous for the Reference allele.

• FreqHomSnp - proportion of samples homozygous for the Alternate (SNP) allele.

• FreqHets - proportion of samples which score as heterozygous, that is, scored as 1.

• PICRef - polymorphism information content (PIC) for the Reference allele.

• PICSnp - polymorphism information content (PIC) for the SNP.

• AvgPIC - average of the polymorphism information content (PIC) of the reference and SNP alleles.

• AvgCountRef - sum of the tag read counts for all samples, divided by the number of samples with non-zero tag read counts, for the Reference allele row.

• AvgCountSnp - sum of the tag read counts for all samples, divided by the number of samples with non-zero tag read counts, for the Alternate (SNP) allele row.

• RepAvg - proportion of technical replicate assay pairs for which the marker score is consistent.

• rdepth - read depth.

Function's output The minimum, maximum, mean and a tabulation of quantiles of the locmetric values against thresholds rate are provided. Output also includes a boxplot and a histogram. Quantiles are partitions of a finite set of values into q subsets of (nearly) equal sizes. In this function q = 20. Quantiles are useful measures because they are less susceptible to long-tailed distributions and outliers. Plot colours can be set with gl.select.colors(). If plot.file is specified, plots are saved to the directory specified by the user, or the global default working directory set by gl.set.wd() or to the tempdir(). Examples of other themes that can be used can be consulted in:

• https://ggplot2.tidyverse.org/reference/ggtheme.html and

• https://yutannihilation.github.io/allYourFigureAreBelongToUs/ggthemes/

Value

An unaltered genlight object.

Author(s)

Luis Mijangos (Post to https://groups.google.com/d/forum/dartr)

See Also

gl.filter.locmetric

Other matched report: gl.filter.excess.het(), gl.report.allna(), gl.report.callrate(), gl.report.hamming(), gl.report.maf(), gl.report.overshoot(), gl.report.pa(), gl.report.rdepth(), gl.report.reproducibility(), gl.report.secondaries(), gl.report.taglength()

Examples

# SNP data
out <- gl.report.locmetric(testset.gl,metric='SnpPosition')
# Tag P/A data
out <- gl.report.locmetric(testset.gs,metric='AvgReadDepth')

Reports minor allele frequency (MAF) for each locus in a SNP dataset

Description

This script provides summary histograms of MAF for each population and an overall histogram to assist the decision of choosing thresholds for the filter function gl.filter.maf

Usage

gl.report.maf(
  x,
  as.pop = NULL,
  maf.limit = 0.5,
  ind.limit = 5,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = NULL,
  plot.dir = NULL,
  plot.file = NULL,
  bins = 25,
  verbose = NULL
)

Arguments

Details

The function gl.filter.maf will filter out the loci with MAF below a specified threshold.

Function's output

The minimum, maximum, mean and a tabulation of MAF quantiles against thresholds rate are provided. Output also includes a boxplot and a histogram.

This function reports the MAF for each of several quantiles. Quantiles are partitions of a finite set of values into q subsets of (nearly) equal sizes. In this function q = 20. Quantiles are useful measures because they are less susceptible to long-tailed distributions and outliers.

Plot colours can be set with gl.select.colors().

If plot.file is specified, plots are saved to the directory specified by the user, or the global default working directory set by gl.set.wd() or to the tempdir().

Examples of other themes that can be used can be consulted in

• https://ggplot2.tidyverse.org/reference/ggtheme.html and

• https://yutannihilation.github.io/allYourFigureAreBelongToUs/ggthemes/

Value

An unaltered genlight object

Author(s)

Custodian: Arthur Georges (Post to https://groups.google.com/d/forum/dartr)

See Also

gl.filter.maf

Other matched report: gl.filter.excess.het(), gl.report.allna(), gl.report.callrate(), gl.report.hamming(), gl.report.locmetric(), gl.report.overshoot(), gl.report.pa(), gl.report.rdepth(), gl.report.reproducibility(), gl.report.secondaries(), gl.report.taglength()

Examples

gl <- gl.filter.allna(platypus.gl)
gl.report.maf(gl)

Reports monomorphic loci

Description

This script reports the number of monomorphic loci and those with all NAs in a genlight {adegenet} object

Usage

gl.report.monomorphs(x, verbose = NULL)

Arguments

Details

A DArT dataset will not have monomorphic loci, but they can arise, along with loci that are scored all NA, when populations or individuals are deleted. Retaining monomorphic loci unnecessarily increases the size of the dataset and will affect some calculations. Note that for SNP data, NAs likely represent null alleles; in tag presence/absence data, NAs represent missing values (presence/absence could not be reliably scored)

Value

An unaltered genlight object

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

gl.filter.monomorphs

Other matched reports: gl.mahal.assign(), gl.report.bases(), gl.report.factorloadings(), gl.report.fstat()

Examples

# SNP data
  gl.report.monomorphs(testset.gl)
# SilicoDArT data
  gl.report.monomorphs(testset.gs)

Reports loci for which the SNP has been trimmed from the sequence tag along with the adaptor

Description

This function checks the position of the SNP within the trimmed sequence tag and identifies those for which the SNP position is outside the trimmed sequence tag. This can happen, rarely, when the sequence containing the SNP resembles the adaptor.

Usage

gl.report.overshoot(x, verbose = NULL)

Arguments

Details

The SNP genotype can still be used in most analyses, but functions like gl2fasta() will present challenges if the SNP has been trimmed from the sequence tag. Resultant ggplot(s) and the tabulation(s) are saved to the session's temporary directory.

Value

An unaltered genlight object

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

gl.filter.overshoot

Other matched report: gl.filter.excess.het(), gl.report.allna(), gl.report.callrate(), gl.report.hamming(), gl.report.locmetric(), gl.report.maf(), gl.report.pa(), gl.report.rdepth(), gl.report.reproducibility(), gl.report.secondaries(), gl.report.taglength()

Examples

gl.report.overshoot(testset.gl)

Reports private alleles (and fixed alleles) per pair of populations

Description

This function reports private alleles in one population compared with a second population, for all populations taken pairwise. It also reports a count of fixed allelic differences and the mean absolute allele frequency differences (AFD) between pairs of populations.

Usage

gl.report.pa(
  x,
  x2 = NULL,
  method = "pairwise",
  loc.names = FALSE,
  test.asym = FALSE,
  test.asym.boot = 100,
  plot.display = FALSE,
  matrix.pa = FALSE,
  plot.font = 14,
  map.interactive = FALSE,
  provider = "Esri.NatGeoWorldMap",
  palette.discrete = NULL,
  plot.file = NULL,
  plot.dir = NULL,
  verbose = NULL
)

Arguments

Details

Note that the number of paired alleles between two populations is not a symmetric dissimilarity measure.

If no x2 is provided, the function uses the pop(gl) hierarchy to determine pairs of populations, otherwise it runs a single comparison between x and x2.

Hint: in case you want to run comparisons between individuals (assuming individual names are unique), you can simply redefine your population names with your individual names, as below:

pop(gl) <- indNames(gl)

Definition of fixed and private alleles

The table below shows the possible cases of allele frequencies between two populations (0 = homozygote for Allele 1, x = both Alleles are present, 1 = homozygote for Allele 2).

• p: cases where there is a private allele in pop1 compared to pop2 (but not viceversa)

• f: cases where there is a fixed allele in pop1 (and pop2, as those cases are symmetric)

The absolute allele frequency difference (AFD) in this function is a simple differentiation metric displaying intuitive properties which provides a valuable alternative to FST. For details about its properties and how it is calculated see Berner (2019).

The Bootstrap test for significant differences of private alleles uses a bootstrap simulation by shuffling individuals between a pair of populations and drawing with replacement. For each bootstrap the ratio of private alleles is compared to the actual ratio and recorded how often it is larger than the simulated one. If number of individuals are different between population bootstrap is done using the smaller number of samples in both populations.

The function also reports an estimation of the lower bound of the number of undetected private alleles using the Good-Turing frequency formula, originally developed for cryptography, which estimates in an ecological context the true frequencies of rare species in a single assemblage based on an incomplete sample of individuals. The approach is described in Chao et al. (2017). For this function, the equation 2c is used. This estimate is reported in the output table as Chao1 and Chao2.

In this function a Sankey Diagram is used to visualize patterns of private alleles between populations. This diagram allows to display flows (private alleles) between nodes (populations). Their links are represented with arcs that have a width proportional to the importance of the flow (number of private alleles).

if save2temp=TRUE, resultant plot(s) and the tabulation(s) are saved to the session's temporary directory.

Value

A data.frame. Each row shows, for each pair of populations the number of individuals in each population, the number of loci with fixed differences (same for both populations) in pop1 (compared to pop2) and viceversa. Same for private alleles and finally the absolute mean allele frequency difference between loci (AFD). If loc.names = TRUE, loci names with private alleles and fixed differences are reported in a list in addition to the dataframe.

Author(s)

Custodian: Bernd Gruber – Post to https://groups.google.com/d/forum/dartr

References

• Berner, D. (2019). Allele frequency difference AFD – an intuitive alternative to FST for quantifying genetic population differentiation. Genes, 10(4), 308.

• Chao, Anne, et al. "Deciphering the enigma of undetected species, phylogenetic, and functional diversity based on Good-Turing theory." Ecology 98.11 (2017): 2914-2929.

See Also

Other matched report: gl.filter.excess.het(), gl.report.allna(), gl.report.callrate(), gl.report.hamming(), gl.report.locmetric(), gl.report.maf(), gl.report.overshoot(), gl.report.rdepth(), gl.report.reproducibility(), gl.report.secondaries(), gl.report.taglength()

Other report functions: gl.report.replicates()

Examples

out <- gl.report.pa(platypus.gl)

out <- gl.report.pa(platypus.gl)

Reports observed, expected and unbiased heterozygosities and FIS (inbreeding coefficient) by population or by individual from SNP data

Description

Calculates the observed, expected and unbiased expected (i.e. corrected for sample size) heterozygosities and FIS (inbreeding coefficient) for each population or the observed heterozygosity for each individual in a genlight object.

Usage

gl.report.polyploid_heterozygosity(
  x,
  method = "pop",
  n.invariant = 0,
  subsample.pop = FALSE,
  n.limit = 10,
  nboots = 0,
  conf = 0.95,
  CI.type = "bca",
  ncpus = 1,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors.pop = gl.colors("dis"),
  plot.colors.ind = gl.colors(2),
  plot.file = NULL,
  plot.dir = NULL,
  error.bar = "SD",
  verbose = NULL
)

Arguments

Details

Observed heterozygosity for a population takes the proportion of heterozygous loci for each individual and averages it over all individuals in that population. The calculations take into account missing values.

Expected heterozygosity for a population takes the expected proportion of heterozygotes, that is, expected under Hardy-Weinberg equilibrium, for each locus, then averages this across the loci for an average estimate for the population.

The unbiased expected heterozygosity is calculated using the correction for sample size following equation 2 from Nei 1978.

Accuracy of all heterozygosity estimates is affected by small sample sizes, and so is their comparison between populations or repeated analysis. Expected heterozygosities are less affected because their calculations are based on allele frequencies while observed heterozygosities are strongly susceptible to sampling effects when the sample size is small.

Observed heterozygosity for individuals is calculated as the proportion of heterozygous gametes that could be produced by that individual.

Finally, the loci that are invariant across all individuals in the dataset (that is, across populations), is typically unknown. This can render estimates of heterozygosity analysis specific, and so it is not valid to compare such estimates across species or even across different analyses (see Schimdt et al 2021). This is a similar problem faced by microsatellites. If you have an estimate of the number of invariant sequence tags (loci) in your data, such as provided by gl.report.secondaries, you can specify it with the n.invariant parameter to standardize your estimates of heterozygosity. This is called autosomal heterozygosities by Schimddt et al (2021).

NOTE: It is important to realise that estimation of adjusted (autosomal) heterozygosity requires that secondaries not to be removed.

Heterozygosities and FIS (inbreeding coefficient) are calculated by locus within each population using the following equations, and then averaged across all loci:

• Observed heterozygosity (Ho) = number of heterozygous gametes / all combinations of gametes, where n_Ind is the number of individuals without missing data for that locus.

• Observed heterozygosity adjusted (Ho.adj) <- Ho * n_Loc / (n_Loc + n.invariant), where n_Loc is the number of loci that do not have all missing data and n.invariant is an estimate of the number of invariant loci to adjust heterozygosity.

• Expected heterozygosity (He) = 1 - (p^2 + q^2), where p is the frequency of the reference allele and q is the frequency of the alternative allele.

• Expected heterozygosity adjusted (He.adj) = He * n_Loc / (n_Loc + n.invariant)

• Unbiased expected heterozygosity (uHe) = He * (2 * n_Ind / (2 * n_Ind - 1))

• Inbreeding coefficient (FIS) = 1 - Ho / uHe

Function's output

Output for method='pop' is an ordered barchart of observed heterozygosity, unbiased expected heterozygosity and FIS (Inbreeding coefficient) across populations together with a table of mean observed and expected heterozygosities and FIS by population and their respective standard deviations (SD). In the output, it is also reported by population: the number of loci used to estimate heterozygosity (n.Loc), the number of polymorphic loci (polyLoc), the number of monomorphic loci (monoLoc) and loci with all missing data (all_NALoc).

Output for method='ind' is a histogram and a boxplot of heterozygosity across individuals.

If a plot.file is given, the ggplot arising from this function is saved as an "RDS" binary file using saveRDS(); can be reloaded with readRDS(). A file name must be specified for the plot to be saved.

If a plot directory (plot.dir) is specified, the ggplot binary is saved to that directory; otherwise to the tempdir().

Examples of other themes that can be used can be consulted in:

• https://ggplot2.tidyverse.org/reference/ggtheme.html and

• https://yutannihilation.github.io/allYourFigureAreBelongToUs/ggthemes/

Subsampling populations

To test the effect of five population sample sizes (n = 10, 5, 4, 3, 2) on observed heterozygosity estimates, the function subsamples individuals, without replacement. The subsampling is repeated 10 times for each sample size n. This approach is an implementation of Schmidt et al (2021).

Error bars

The best method for presenting or assessing genetic statistics depends on the type of data you have and the specific questions you're trying to answer. Here's a brief overview of when you might use each method:

1. Confidence Intervals ("CI"):

- Usage: Often used to convey the precision of an estimate.

- Advantage: Confidence intervals give a range in which the true parameter (like a population mean) is likely to fall, given the data and a specified probability (like 95%).

- In Context: For genetic statistics, if you're estimating a parameter, a 95% CI gives you a range in which you're 95% confident the true parameter lies.

2. Standard Deviation ("SD"):

- Usage: Describes the amount of variation from the average in a set of data.

- Advantage: Allows for an understanding of the spread of individual data points around the mean.

- In Context: If you're looking at the distribution of a quantitative trait (like height) in a population with a particular genotype, the SD can describe how much individual heights vary around the average height.

3. Standard Error ("SE"):

- Usage: Describes the precision of the sample mean as an estimate of the population mean.

- Advantage: Smaller than the SD in large samples; it takes into account both the SD and the sample size.

- In Context: If you want to know how accurately your sample mean represents the population mean, you'd look at the SE.

Recommendation:

- If you're trying to convey the precision of an estimate, confidence intervals are very useful.

- For understanding variability within a sample, standard deviation is key.

- To see how well a sample mean might estimate a population mean, consider the standard error.

analyze and present their data, depending on their research questions and the nature of the data.

Confidence Intervals

The uncertainty of a parameter, in this case the mean of the statistic, can be summarised by a confidence interval (CI) which includes the true parameter value with a specified probability (i.e. confidence level; the parameter "conf" in this function).

In this function, CI are obtained using Bootstrap which is an inference method that samples with replacement the data (i.e. loci) and calculates the statistics every time.

This function uses the function boot (package boot) to perform the bootstrap replicates and the function boot.ci (package boot) to perform the calculations for the CI.

Four different types of nonparametric CI can be calculated (parameter "CI.type" in this function):

• First order normal approximation interval ("norm").

• Basic bootstrap interval ("basic").

• Bootstrap percentile interval ("perc").

• Adjusted bootstrap percentile interval ("bca").

The studentized bootstrap interval ("stud") was not included in the CI types because it is computationally intensive, it may produce estimates outside the range of plausible values and it has been found to be erratic in practice, see for example the "Studentized (t) Intervals" section in:

https://www.r-bloggers.com/2019/09/understanding-bootstrap-confidence-interval-output-from-the-r-boot-package/

Nice tutorials about the different types of CI can be found in:

https://www.datacamp.com/tutorial/bootstrap-r

and

https://www.r-bloggers.com/2019/09/understanding-bootstrap-confidence-interval-output-from-the-r-boot-package/

Efron and Tibshirani (1993, p. 162) and Davison and Hinkley (1997, p. 194) suggest that the number of bootstrap replicates should be between 1000 and 2000.

It is important to note that unreliable confidence intervals will be obtained if too few number of bootstrap replicates are used. Therefore, the function boot.ci will throw warnings and errors if bootstrap replicates are too few. Consider increasing the number of bootstrap replicates to at least 200.

The "bca" interval is often cited as the best for theoretical reasons, however it may produce unstable results if the bootstrap distribution is skewed or has extreme values. For example, you might get the warning "extreme order statistics used as endpoints" or the error "estimated adjustment 'a' is NA". In this case, you may want to use more bootstrap replicates or a different method or check your data for outliers.

The error "estimated adjustment 'w' is infinite" means that the estimated adjustment ‘w’ for the "bca" interval is infinite, which can happen when the empirical influence values are zero or very close to zero. This can be caused by various reasons, such as:

The number of bootstrap replicates is too small, the statistic of interest is constant or nearly constant across the bootstrap samples, the data contains outliers or extreme values.

You can try some possible solutions, such as:

Increasing the number of bootstrap replicates, using a different type of bootstrap confidence interval or removing or transforming the outliers or extreme values.

Parallelisation

If the parameter ncpus > 1, parallelisation is enabled. In Windows, parallel computing employs a "socket" approach that starts new copies of R on each core. POSIX systems, on the other hand (Mac, Linux, Unix, and BSD), utilise a "forking" approach that replicates the whole current version of R and transfers it to a new core.

Opening and terminating R sessions in each core involves a significant amount of processing time, therefore parallelisation in Windows machines is only quicker than not using parallelisation when nboots > 1000-2000.

Value

A dataframe containing population labels, heterozygosities, FIS, their standard deviations and sample sizes.

Author(s)

Custodian: Ching Ching Lau (Post to https://groups.google.com/d/forum/dartr)

References

• Moody, M. E., Mueller, L. D., & Soltis, D. E. (1993). Genetic variation and random drift in autotetraploid populations. Genetics, 134(2), 649-657.

See Also

gl.filter.heterozygosity

Other unmatched report: gl.allele.freq(), gl.report.basics(), gl.report.diversity(), gl.report.excess.het(), gl.report.heterozygosity()

Reports summary of Read Depth for each locus

Description

SNP datasets generated by DArT report AvgCountRef and AvgCountSnp as counts of sequence tags for the reference and alternate alleles respectively. These can be used to back calculate Read Depth. Fragment presence/absence datasets as provided by DArT (SilicoDArT) provide Average Read Depth and Standard Deviation of Read Depth as standard columns in their report. This function reports the read depth by locus for each of several quantiles.

Usage

gl.report.rdepth(
  x,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = NULL,
  plot.dir = NULL,
  plot.file = NULL,
  verbose = NULL
)

Arguments

Details

The function displays a table of minimum, maximum, mean and quantiles for read depth against possible thresholds that might subsequently be specified in gl.filter.rdepth. If plot.display=TRUE, display also includes a boxplot and a histogram to guide in the selection of a threshold for filtering on read depth. Plot colours can be set with gl.select.colors(). If plot.file is specified, plots are saved to the directory specified by the user, or the global default working directory set by gl.set.wd() or to the tempdir(). For examples of themes, see

• https://ggplot2.tidyverse.org/reference/ggtheme.html and

• https://yutannihilation.github.io/allYourFigureAreBelongToUs/ggthemes/

Value

An unaltered genlight object

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

gl.filter.rdepth

Other matched report: gl.filter.excess.het(), gl.report.allna(), gl.report.callrate(), gl.report.hamming(), gl.report.locmetric(), gl.report.maf(), gl.report.overshoot(), gl.report.pa(), gl.report.reproducibility(), gl.report.secondaries(), gl.report.taglength()

Examples

# SNP data
df <- gl.report.rdepth(testset.gl)
df <- gl.report.rdepth(testset.gs)

Identify replicated individuals

Description

Identify replicated individuals

Usage

gl.report.replicates(
  x,
  loc_threshold = 100,
  perc_geno = 0.95,
  plot.out = TRUE,
  plot_theme = theme_dartR(),
  plot_colors = c("#2171B5", "#6BAED6"),
  bins = 100,
  verbose = NULL
)

Arguments

Details

This function uses an C++ implementation, so package Rcpp needs to be installed and it is therefore fast (once it has compiled the function after the first run).

Ideally, in a large dataset with related and unrelated individuals and several replicated individuals, such as in a capture/mark/recapture study, the first histogram should have four "peaks". The first peak should represent unrelated individuals, the second peak should correspond to second-degree relationships (such as cousins), the third peak should represent first-degree relationships (like parent/offspring and full siblings), and the fourth peak should represent replicated individuals.

In order to ensure that replicated individuals are properly identified, it's important to have a clear separation between the third and fourth peaks in the second histogram. This means that there should be bins with zero counts between these two peaks.

Value

A list with three elements:

• table.rep: A dataframe with pairwise results of percentage of same genotypes between two individuals, the number of loci used in the comparison and the missing data for each individual.

• ind.list.drop: A vector of replicated individuals to be dropped. Replicated individual with the least missing data is reported.

• ind.list.rep: A list of of each individual that has replicates in the dataset, the name of the replicates and the percentage of the same genotype.

Author(s)

Custodian: Luis Mijangos – Post to https://groups.google.com/d/forum/dartr

See Also

Other report functions: gl.report.pa()

Examples

res_rep <- gl.report.replicates(platypus.gl, loc_threshold = 500, 
perc_geno = 0.85)

Reports summary of RepAvg (repeatability averaged over both alleles for each locus) or reproducibility (repeatability of the scores for fragment presence/absence)

Description

SNP datasets generated by DArT have an index, RepAvg, generated by reproducing the data independently for 30 of alleles that give a repeatable result, averaged over both alleles for each locus. In the case of fragment presence/absence data (SilicoDArT), repeatability is the percentage of scores that are repeated in the technical replicate dataset.

Usage

gl.report.reproducibility(
  x,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = NULL,
  plot.dir = NULL,
  plot.file = NULL,
  verbose = NULL
)

Arguments

Details

The function displays a table of minimum, maximum, mean and quantiles for repeatbility against possible thresholds that might subsequently be specified in gl.filter.reproducibility. If plot.display=TRUE, display also includes a boxplot and a histogram to guide in the selection of a threshold for filtering on repeatability. If plot.file is specified, plots are saved to the directory specified by the user, or the global default working directory set by gl.set.wd() or to the tempdir(). For examples of themes, see:

• https://ggplot2.tidyverse.org/reference/ggtheme.html and

• https://yutannihilation.github.io/allYourFigureAreBelongToUs/ggthemes/

Value

An unaltered genlight object

Author(s)

Custodian: Arthur Georges – Post to https://groups.google.com/d/forum/dartr

See Also

gl.filter.reproducibility

Other matched report: gl.filter.excess.het(), gl.report.allna(), gl.report.callrate(), gl.report.hamming(), gl.report.locmetric(), gl.report.maf(), gl.report.overshoot(), gl.report.pa(), gl.report.rdepth(), gl.report.secondaries(), gl.report.taglength()

Examples

 
# SNP data
  out <- gl.report.reproducibility(testset.gl)
  
# Tag P/A data
  out <- gl.report.reproducibility(testset.gs)

Reports loci containing secondary SNPs in sequence tags and calculates number of invariant sites

Description

SNP datasets generated by DArT include fragments with more than one SNP (that is, with secondaries). They are recorded separately with the same CloneID (=AlleleID). These multiple SNP loci within a fragment are likely to be linked, and so you may wish to remove secondaries. This function reports statistics associated with secondaries, and the consequences of filtering them out, and provides three plots. The first is a boxplot, the second is a barplot of the frequency of secondaries per sequence tag, and the third is the Poisson expectation for those frequencies including an estimate of the zero class (no. of sequence tags with no SNP scored).

Usage

gl.report.secondaries(
  x,
  nsim = 1000,
  taglength = 69,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.colors = NULL,
  plot.dir = NULL,
  plot.file = NULL,
  verbose = NULL
)

Arguments

Details

The function gl.filter.secondaries will filter out the loci with secondaries retaining only one sequence tag. Heterozygosity as estimated by the function gl.report.heterozygosity is in a sense relative, because it is calculated against a background of only those loci that are polymorphic somewhere in the dataset. To allow intercompatibility across studies and species, any measure of heterozygosity needs to accommodate loci that are invariant (autosomal heterozygosity. See Schmidt et al 2021). However, the number of invariant loci are unknown given the SNPs are detected as single point mutational variants and invariant sequences are discarded, and because of the particular additional filtering pre-analysis. Modelling the counts of SNPs per sequence tag as a Poisson distribution in this script allows estimate of the zero class, that is, the number of invariant loci. This is reported, and the veracity of the estimate can be assessed by the correspondence of the observed frequencies against those under Poisson expectation in the associated graphs. The number of invariant loci can then be optionally provided to the function gl.report.heterozygosity via the parameter n.invariants. In case the calculations for the Poisson expectation of the number of invariant sequence tags fail to converge, try to rerun the analysis with a larger nsim values. This function now also calculates the number of invariant sites (i.e. nucleotides) of the sequence tags (if TrimmedSequence is present in x$other$loc.metrics) or estimate these by assuming that the average length of the sequence tags is 69 nucleotides. Based on the Poisson expectation of the number of invariant sequence tags, it also estimates the number of invariant sites for these to eventually provide an estimate of the total number of invariant sites. Note, previous version of dartR would only return an estimate of the number of invariant sequence tags (not sites). If plot.file is specified, plots are saved to the directory specified by the user, or the global default working directory set by gl.set.wd() or to the tempdir(). Examples of other themes that can be used can be consulted in:

• https://ggplot2.tidyverse.org/reference/ggtheme.html and

• https://yutannihilation.github.io/allYourFigureAreBelongToUs/ggthemes/

• n.total.tags Number of sequence tags in total

• n.SNPs.secondaries Number of secondary SNP loci that would be removed on filtering

• n.invariant.tags Estimated number of invariant sequence tags

• n.tags.secondaries Number of sequence tags with secondaries

• n.inv.gen Number of invariant sites in sequenced tags

• mean.len.tag Mean length of sequence tags

• n.invariant Total Number of invariant sites (including invariant sequence tags)

• k Lambda: mean of the Poisson distribution of number of SNPs in the sequence tags

Value

A data.frame with the list of parameter values

Author(s)

Custodian: Arthur Georges (Post to https://groups.google.com/d/forum/dartr)

References

Schmidt, T.L., Jasper, M.-E., Weeks, A.R., Hoffmann, A.A., 2021. Unbiased population heterozygosity estimates from genome-wide sequence data. Methods in Ecology and Evolution n/a.

See Also

gl.filter.secondaries,gl.report.heterozygosity, utils.n.var.invariant

Other matched report: gl.filter.excess.het(), gl.report.allna(), gl.report.callrate(), gl.report.hamming(), gl.report.locmetric(), gl.report.maf(), gl.report.overshoot(), gl.report.pa(), gl.report.rdepth(), gl.report.reproducibility(), gl.report.taglength()

Examples

require("dartR.data")
test <- gl.filter.callrate(platypus.gl,threshold = 1)
n.inv <- gl.report.secondaries(test)
gl.report.heterozygosity(test, n.invariant = n.inv[7, 2])

Report SNP diversity from a genlight object, with reference to Ma, Z., Li, L., & Zhang, Y. P. (2020). Defining individual-level genetic diversity and similarity profiles. Scientific reports, 10(1), 5805.

Description

This function needs package adegenet, please install it.

Usage

gl.report.shannon(
  x,
  plot.display = TRUE,
  plot.theme = theme_dartR(),
  plot.dir = NULL,
  plot.file = NULL,
  level = "alpha",
  order = 5,
  verbose = 2
)

Arguments