--- title: "Introduction to delimtools" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Introduction to delimtools} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) library(ape) library(delimtools) library(dplyr) library(here) ``` ## # Rationale There are at least three major steps when doing single-locus species delimitation analysis: * Preparing your data; * Running species delimitation analysis; * Visualizing results. The rationale of `delimtools` is to provide functions for all these steps. This guide introduces you to these functions and shows how to apply them to your datasets. In this tutorial, we will use the `geophagus` dataset provided in this package. This dataset contains `r length(geophagus)` cytochrome c oxidase sequences of the *Geophagus sensu stricto* species group downloaded from GenBank. Most of these sequences are from Ximenes et al. ([2021](https://doi.org/10.7717/peerj.12443)) and are documented in `?geophagus`. ## Preparing your data Usually, our main source of information for single-locus species delimitation analysis is a [FASTA](https://www.ncbi.nlm.nih.gov/genbank/fastaformat/) formated file containing nucleotide sequences. The most straightforward way to import a FASTA formatted file into R Environment is by using `ape::read.FASTA()`: ```{r eval=FALSE} path_to_file <- system.file("extdata/geophagus.fasta", package = "delimtools") geophagus <- ape::read.FASTA(file = path_to_file) ``` Since the `geophagus` dataset is provided along with `delimtools`, one can simply call it directly: ```{r} geophagus ``` ### Check identifiers across files Although the header of each FASTA file may contain both an ID and a Description, we recommend keeping the header simple and only using ID. The metadata (the Description) are easily stored in a tabular format file (e.g. `csv`, `tsv`, `xlsx`, etc.), and linked to the FASTA file via the ID. We provide such metadata in Darwin Core format. One can have a glimpse of our *Geophagus* sequence metadata by checking `geophagus_info`: ```{r} dplyr::glimpse(geophagus_info) ``` For downstream visualization, all FASTA IDs should have associated metadata; however, the metadata may have more entries than sequences in the FASTA file. The function `check_identifiers()` will check if all FASTA IDs have entries in the metadata, and if there are any duplicate IDs either in the FASTA file or in the metadata. Let's check if our FASTA file has associated metadata `check_identifiers()`: ```{r} check_identifiers(dna = geophagus, identifier = "gbAccession", data = geophagus_info) ``` Now, let's omit the first 10 sequences and check again: ```{r error= TRUE} check_identifiers(geophagus[11:354], "gbAccession", geophagus_info) ``` Now, let's omit the first 4 metadata entries and check again: ```{r error= TRUE} check_identifiers(geophagus, "gbAccession", geophagus_info[5:354,]) ``` Many other functions of `delimtools` rely on a tabular format file like `geophagus_info`, so please make sure to have one prior to any analysis. ### Cleaning DNA and Haplotype Collapsing Sometimes, our alignment file may contain non *ACTG* bases, like ambiguities (N, R, Y, etc), gaps (-), or missing data (?). Inclusion of these bases in delimitation analyses may have unintended consequences in delimitation analyses. We can use the function `clean_dna` to remove these ambiguities. `clean_dna` will remove both leading and trailing gaps from the alignment, therefore it will be necessary to realign the dataset after calling `clean_dna`. ```{r} geophagus_clean <- clean_dna(geophagus) # and display geophagus_clean ``` The presence of identical sequences (haplotypes) in an alignment will overestimate the number of lineages detected by tree-based species delimitation methods (e.g. bGMYC, GMYC, PTP, mPTP). Therefore it is best practice to collapse identical haplotypes into unique haplotypes. This is accomplished by calling the `hap_collapse` function. By default the function will also call `clean_dna`. ```{r} geophagus_collapsed <- hap_collapse(geophagus) # and display geophagus_collapsed ``` **Notice** that from the initial 354 DNA sequences, `hap_collapse` returned a total of 136 unique DNA sequences. This result may vary if setting `collapseSubstrings` to `FALSE`: ```{r} geophagus_collapsed <- hap_collapse(geophagus, collapseSubstrings = FALSE) # and display geophagus_collapsed ``` After setting `collapseSubstrings` to `FALSE`, we obtained 246 unique DNA sequences. These additional 108 sequences are in reality shorter but identical sequences which differ only by the presence of leading and trailing gaps. By setting both `clean_dna` and `collapseSubstrings` to `TRUE`, the default of `hap_collapse`, we ensure the removal of these leading and trailing gaps from sequences before collapsing these sequences into unique haplotypes. To visualize the relationships between these unique haplotypes and any duplicates, we may use the `haplotype_tbl` function: ```{r} haplotype_tbl(geophagus, verbose = FALSE) ``` We may also collapse other types of information by using `collapse_others`. For example we may take the collapsed haplotype table and onto it collapse our metadata to see, for example, which species and sampling localities are associated with each haplotype: ```{r} haps_df <- haplotype_tbl(geophagus, verbose = FALSE) collapse_others(data = geophagus_info, hap_tbl = haps_df, labels = "gbAccession", cols = c("scientificName", "country")) ``` ## Running species delimitation analysis After preparing your data (alignments and phylogenetic trees), the next step it to run species delimitation analyses and combine them into a single `tbl_df`. As of v0.3.0, `delimtools` includes executables of the ABGD, ASAP, GMYC, bGMYC, PTP and mPTP methods, therefore, it requires no outside dependencies to run these delimitation methods. Following are examples of species delimitation methods implemented in `delimtools`. The first three methods are distance-based methods and thus use the full aligned sequence matrix. Subsequent methods are tree-based methods; the phylogenetic trees need to be estimated externally, and should be based on unique haplotypes. For the convenience of this vignette, we provide the full data matrix (`geophagus`), a phylogram (`geophagus_raxml`) and a time-tree (`geophagus_beast`) as part of this package. ### ABGD The ABGD method is run by executing the `abgd` function and giving it an `ape::DNAbin` object which is a S3 R object containing an aligned FASTA file. It will also accept a path to an aligned FASTA. The result is then transformed into a table using the `abgd_tbl` function. ```{r} # geophagus - the `ape::DNAbin` object supplied with this package abgd_delim <- abgd(geophagus) abgd_df <- abgd_tbl(abgd_delim) # print the dataframe print(abgd_df) ``` ### ASAP The ASAP method is run by executing the `asap` function and giving it an `ape::DNAbin` object. It will also accept a path to an aligned FASTA. The result is then transformed into a table using the `asap_tbl` function. ```{r} # geophagus - the `ape::DNAbin` object supplied with this package asap_delim <- asap(geophagus) asap_df <- asap_tbl(asap_delim) # print the dataframe print(asap_df) ``` ### Local Minima We can run the Local Minima analysis by using the `spider` function `localMinima`. This function requires a genetic distance matrix as input for analysis. ```{r} dist_mat <- ape::dist.dna(geophagus, model = "raw", pairwise.deletion = TRUE) locmin <- spider::localMinima(dist_mat) ``` We can plot the `locmin` object to select the threshold. We will use the first value of the `locmin` object: ```{r fig.height= 3.5, fig.width= 6.9} plot(locmin) |> abline(v=locmin$localMinima[1],col="red") ``` Since `Local Minima` output is a numeric vector, we need to extract the results using `locmin_tbl`. ```{r} locmin_df <- locmin_tbl(dist_mat, threshold = locmin$localMinima[1]) print(locmin_df) ``` We can use patristic distances to generate species partitions as well. Using the `geophagus_beast` tree, we can apply a threshold to classify tips of the tree by using its branch lengths: ```{r fig.height= 3.5, fig.width= 6.9} dist_mat_tr <- ape::cophenetic.phylo(as.phylo(geophagus_beast)) |> as.dist() locmin_tr <- spider::localMinima(dist_mat_tr) plot(locmin_tr) |> abline(v=locmin_tr$localMinima[1],col="red") treedist_df <- locmin_tbl(dist_mat_tr, threshold = locmin_tr$localMinima[1], delimname = "treedist") print(treedist_df) ``` ## Fixed threshold Finally, we can also use `locmin_tbl` to create a species partition using a fixed threshold. For example, we can use a threshold of 2% to generate species partitions: ```{r} percent_df <- locmin_tbl(dist_mat, threshold = 0.02, delimname = "percent") print(percent_df) ``` ### GMYC The GMYC method is run by executing the `gmyc` function and giving it an `ape::phylo` object. The GMYC method requires an ultrametric tree file of class `phylo` as input for analysis. BEAST trees generated in Tree Annotator need to be read in using the `treeio:read.beast()` function, and then need to be converted into a `phylo` object by using `ape::as.phylo`. ```{r} # beast_tree - the Nexus formatted ultrametric phylogenetic tree object supplied with this package path_to_file <- system.file("extdata/geophagus_beast.nex", package = "delimtools") beast_tree <- treeio::read.beast(path_to_file) |> ape::as.phylo() ``` Run delimitation. ```{r} gmyc_delim <- gmyc(beast_tree) ``` Now, lets check the `summary`: ```{r} summary(gmyc_delim) ``` And let's transform the results into a table using the `gmyc_tbl` function. ```{r} gmyc_df <- gmyc_tbl(gmyc_delim) # print the dataframe print(gmyc_df) ``` ### bGMYC The bGMYC method is run by executing the `bgmyc` function and giving it an `ape::phylo` object. The bGMYC method requires an ultrametric tree file of class `phylo` as input for analysis. BEAST trees generated in Tree Annotator need to be read in using the `treeio:read.beast()` function, and then need to be converted into a `phylo` object by using `ape::as.phylo`. ```{r} # beast_tree already loaded and converted bgmyc_delim <- bgmyc(beast_tree) ``` The results then need to be transformed into a table using the `bgmyc_tbl` function. By default, `0.05` is used as default posterior probability threshold for clustering samples into species partitions. ```{r} bgmyc_df <- bgmyc_tbl(bgmyc_delim) print(bgmyc_df) ``` ### PTP and mPTP We can run both PTP and mPTP analysis by using the `mptp` function. Specifying the arguments ("multi" for mPTP; "single" for PTP) in the parameters determines which analysis is performed; the default is "multi" for mPTP. The `mptp` function accepts either phylo tree object or a Newick tree file. We can use either `ape::read.tree()` or `treeio::read.newick()` to read in our Newick file; it will be automatically stored as `phylo` object. ```{r} # raxml_tree - the Newick formatted phylogenetic tree object supplied with this package path_to_file <- system.file("extdata/geophagus_raxml.nwk", package = "delimtools") raxml_tree <- treeio::read.newick(path_to_file) ``` Run delimitation. ```{r} mptp_delim <- mptp(raxml_tree) mptp_df <- mptp_tbl(mptp_delim) print(mptp_df) ``` For convenience we have made a ptp() wrapper around mptp() that defaults to the "single" argument of the method parameter. Similarly, we also made a ptp_tbl() wrapper around mptp_tbl(). ```{r} # raxml_tree already loaded and converted ptp_delim <- ptp(raxml_tree) ptp_df <- ptp_tbl(ptp_delim) print(ptp_df) ``` To estimate a minimum branch length for your dataset, use `minbr_len`. ```{r} min_brlen(raxml_tree) ``` ### Morphology We can also turn species ranks or results from morphological analysis into a species partition by using `morph_tbl`. For the `geophagus` dataset, we will use the scientific names in `geophagus_info` to accomplish this. ```{r} morph_df <- morph_tbl(labels = geophagus_info$gbAccession, sppVector = geophagus_info$scientificName) |> dplyr::filter(labels %in% haps_df$labels) print(morph_df) ``` ## Joining Species Delimitation Outputs Since we have several species partitions, we need to merge them into a single data frame in order to visualize the results. However, since each method has its own rankings, we need to recode them prior to any further analysis. For this task, we can use `delim_join`. This function takes a `list` of species partitions to recode and then join. ```{r message= FALSE} all_delims <- list(abgd_df, asap_df, locmin_df, percent_df, treedist_df, bgmyc_df, gmyc_df, mptp_df, ptp_df, morph_df) all_delims_df <- delim_join(all_delims) ``` We can summarize the `all_delims_df` object by using the `report_delim` function. ```{r} report_delim(all_delims_df) ``` We can also check the match ratio statistics of Ahrens et al (2014) using `match_ratio`. ```{r} match_ratio(all_delims_df) |> dplyr::arrange(dplyr::desc(match_ratio)) |> print(n=Inf) ``` By checking the results, all comparisons including treedist have a low `match_ratio`. Since species partitions with very low `match_ratio` values can impact consensus calculations, we will remove treedist before visualizing our results. ```{r} all_delims_df <- dplyr::select(all_delims_df, -treedist) report_delim(all_delims_df) ``` ## Visualizing Results We can visualize the results by using `delim_autoplot`. This function will plot a phylogenetic tree alongside with the results contained in `all_delims_df`. We will use the `geophagus_beast` tree to visualize the results. ```{r warning=FALSE, message= FALSE, fig.height= 15, fig.width= 12} # geophagus_beast - the Nexus formatted ultrametric phylogenetic tree object used for visualization path_to_file <- system.file("extdata/geophagus_beast.nex", package = "delimtools") geophagus_tree <- treeio::read.beast(path_to_file) # customize tip labels of the tree tip_tab <- geophagus_info |> dplyr::filter(gbAccession %in% geophagus_beast@phylo$tip.label) |> dplyr::mutate(labs = glue::glue("{gbAccession} | {scientificName}")) |> dplyr::select(gbAccession, labs, scientificName) # create a customized color palette cols <- delim_brewer(delim = all_delims_df, package = "randomcoloR", seed = 42) # plot delim_autoplot(all_delims_df, geophagus_beast, consensus = TRUE, n_match = 5, tbl_labs = tip_tab, col_vec = cols, hexpand = 0.7, widths = c(0.5, 0.5)) ```