To export or not export from phruta

Assembling a molecular dataset for particular target taxa with phruta can be performed almost entirely by saving objects into your workspace. This topic was covered in the introductory vignette to phruta (Using the phruta R package). However, in some situations, it might be desirable to save the outputs of different phruta functions to particular folders. This tutorial will cover that specific situation. Specifically, we will be reviewing how phruta can be used to export .fasta and .csv files that are generated in different steps of the pipeline.

From taxonomic names to sequence alignments, exporting data

In the introductory vignette to phruta (Using the phruta R package), we assembled a basic molecular dataset for inferring the phylogeny among three mammal genera. Let’s recreate the same tutorial, but in this case, exporting the intermediate files that are created after using several of the functions. Note that the structure of this tutorial closely follows that of “Using the phruta R package”. You should be able to follow this tutorial even without having reviewed the introductory vignette.

Please assume that we are interested in building a phylogenetic tree for the following three genera: Felis, Vulpes, and Phoca. All these three genera are classified within the Carnivora, a mammalian order. Both Felis and Vulpes are classified in different superfamilies within the Fissipedia. Finally, Phoca is part of another suborder, Pinnipedia. We’re going to root our tree with another mammal species, a Chinese Pangolin (Manis pentadactyla). Users can select additional target species and clades. However, for simplicity, we will run the analyses using three genera in the ingroup and a single outgroup species.

So far, we have decided the taxonomic make of our analyses in phruta. We will also need to determine the gene regions to be used in our analyses. Fortunately, mammals are extensively studied and a comprehensive list of potential gene regions to be analyzed is already available. For instance, we could use same gene regions sampled in Upham et al (2009). However, for this tutorial, we will simply try to find the gene regions are well sampled for the target taxa. I believe that figuring out the best sampled gene regions in genbank, instead of providing gene names, is potentially more valuable when working with poorly studied groups (e.g. invertebrates). Before we move on, please make sure that you you have set a working directory for this project. All the files will be saved to this directory. phruta will politely ask before writing files to your local directories.

Let’s start by loading phruta!

library(phruta)

Now, let’s look for the gene regions that are sampled for our target taxa. Again, this step is not always necessary. In some groups, gene-level sampling is very standard (e.g. COI, 12S). However, the structure of gene sampling sometimes becomes more blurry as you zoom out taxonomically. For instance, genes A and B can be extensively sampled in genus 1. However, genus 2 in the same family has mainly been studied using genes Y and Z. The idea here is that phruta will try to find those gene regions that are extensively sampled across species in the target taxa. We will use the gene.sampling.retrieve() function in phruta. The resulting data.frame, named gs.seqs in this example, will contain the list of full names for genes sampled in genbank for the target taxa.

gs.seqs <- gene.sampling.retrieve(organism = c("Felis", "Vulpes", "Phoca", "Manis_pentadactyla"), 
                                  speciesSampling = TRUE)

For the search terms used above, phruta was able to retrieve the names for 1594 gene regions. In the table below I summarize a few of those genes, with sampling frequency calculated at the level of species (see speciesSampling = TRUE argument above).

Thus, the gene.sampling.retrieve() function provides an estimate of the number of species in genbank that matches the taxonomic criteria and have sequences for a given gene region. Note that the estimates recovered by gene.sampling.retrieve() are only as good as the annotations that other researchers have provided for sequences deposited in genbank.

From here, we will generate a preliminary table summarizing accession numbers for the combination of taxa and gene regions that we’re interested in sampling. However, note that not all these accession numbers are expected to be in the final (curated) molecular dataset. For instance, several sequences might be dropped later after taxonomic information is curated. Now, we will assemble a species-level summary of accession numbers using the acc.table.retrieve() function. For simplicity, this tutorial will focus on sampling gene regions that are sampled in >30% of the species (targetGenes data.frame).

targetGenes <- gs.seqs[gs.seqs$PercentOfSampledSpecies > 30,]

acc.table <- acc.table.retrieve(
            clades  = c('Felis', 'Vulpes', 'Phoca'),
            species = 'Manis_pentadactyla' ,
            genes   = targetGenes$Gene,
            speciesLevel = TRUE
          )

The acc.table object is a data.frame that will be used below for downloading the relevant gene sequences. In this case, the dataset includes the following information:

Feel free to review this dataset, make changes, add new species, samples, etc. The integrity of this dataset is critical for the next steps so please take your time and review it carefully. For instance, let’s just make some minor changes to our dataset:

acc.table$Species <- sub("P.", "Phoca ", acc.table$Species, fixed = TRUE)
acc.table$Species <- sub("F.", "Felis ", acc.table$Species, fixed = TRUE)
acc.table$Species <- sub("V.", "Vulpes ", acc.table$Species, fixed = TRUE)
acc.table$Species <- sub("mitochondrial", "", acc.table$Species)
row.names(acc.table) <- NULL

Let’s check how the new table looks now…

Now, since we’re going to retrieve sequences from genbank using an existing preliminary accession numbers table, we will use the sq.retrieve.indirect() function in phruta. I’m going to spend some time in here to explain the differences between the two versions of sq.retrieve.* in phruta. The one that we’re using in this tutorial, sq.retrieve.indirect(), retrieves sequences “indirectly” because it follows the initial step of generating a table summarizing accession numbers (see the acc.table.retrieve() function above). I present the information in this vignette using sq.retrieve.indirect() instead of sq.retrieve.direct() because the first function is way more flexible and allows for correcting issues prior to download any sequence. For instance, you can add new sequences, species, populations to the resulting data.frame from acc.table.retrieve(). Additionally, you could even manually assemble your own dataset of accession numbers to be retrieved using sq.retrieve.indirect(). Instead, sq.retrieve.direct() does its best to directly (i.e. without potential input from the user) retrieve sequences for a target set of taxa and set of gene regions. In short, you should be able to catch errors using sq.retrieve.indirect() but mistakes will be harder to spot and fix if you’re using sq.retrieve.direct(). Note that the functionality of sq.retrieve.direct() is outlined in the “Using phruta with defined target genes” vignette.

We still need to retrieve all the sequences from the accessions table that was generated avobe using acc.table. The sq.retrieve.indirect() function will write all the resulting fasta files into a newly created folder 0.Sequences located in our working directory (please check the download.sqs = TRUE argument).

sq.retrieve.indirect(acc.table, download.sqs = TRUE)

Next, we’re going to make sure that we include only sequences that are reliable and from species that we are actually interested in analyzing. For this, we will be using the sq.curate() function. We need to provide a list of taxonomic names to filter out incorrect sequences (filterTaxonomicCriteria argument). For simplicity, our criteria can be the genera that we’re interested in analyzing. Note that the outgroup’s name should also be included in the list. If the taxonomic information for a sequence retrieved from genbank does not match with any of these strings, this species will be dropped. You will have to specify whether sampling is for animals or plants (kingdom argument). Finally, you might have already noticed that the same gene regions can have different names. For instance, sometimes our searches retrieve both “cytochrome oxidase subunit 1” and “cytochrome c oxidase subunit I” as widely sampled genes for the target species. In that case, we can combine the sequences in these two files into a single file name COI. To merge gene files, you will have to provide a named list to the mergeGeneFiles argument of the sq.curate function. This named list (tb.merged below) will have a length that corresponds to the number of final files that should be constructed.

tb.merged <- list('COI' = c("cytochrome oxidase subunit 1", "cytochrome c oxidase subunit I"))
sq.curate(filterTaxonomicCriteria = 'Felis|Vulpes|Phoca|Manis',
          mergeGeneFiles = tb.merged,
          kingdom = 'animals', 
          folder = '0.Sequences',
          removeOutliers = FALSE)

Running the line of code above will create the a folder 1.CuratedSequences containing (1) the curated sequences with original names, (2) the curated sequences with species-level names (renamed_* prefix), (3) a table of accession numbers (0.AccessionTable.csv), and (4) a summary of the taxonomic information for all the species sampled in the files (1.Taxonomy.csv). We’ll use the renamed_* and 1.Taxonomy.csv files in the next steps. Let’s take a look at the sampling per gene region in the 0.AccessionTable.csv table.

We’ll now align the sequences that we just curated. For this, we just use sq.aln() with default parameters. We need to indicate that we’re interested in aligning only the "renamed" fasta files in our 1.CuratedSequences folder.

sq.aln(folder = '1.CuratedSequences', FilePatterns = "renamed")

The resulting multiple sequence alignments will be saved to the 2.Alignments folder. In that new folder, we will have two types of files: (1) raw alignments (same file names as in 1.CuratedSequences) and (2) alignments with ambiguous sites removed (Masked_* prefix). Masked alignments are only created if the mask argument in sq.aln is set to TRUE. In that case, one additional .csv file is created for each of the alignments (0.Masked.Information_*). Each of these datasets list the number of sites in the masked alignment that (1) are not gaps (NonGaps column), (2) if the sequence was removed due to the elevated number of gaps (removedPerGaps; controlled using the threshold argument in sq.aln), or (3) if it was removed directly in the masking step (removedMasking). Note that, for some gene regions, making can fail. In that case, only the original alignment file is saved to the 2.Alignments folder.

Note that we could use these resulting alignments directly to infer our phylogenies. We cover these steps within phruta in another vignette: “Phylogenetics with the phruta R package”. For now, let’s wrap up and plot one of our (cool) alignments. Let’s first check the raw alignments!

A figure showing raw alignments

Now, the masked alignments…!!

A figure showing curated alignments

And we’re done for now!! Thanks for following this tutorial…:)

In total, this vignette took 11 minutes to render in my local machine. You can now try to run phruta using your favorite groups organisms! Don’t forget to check the other tutorials and get in touch if you find any issues…Buena suerte!