Using phruta with a defined set of target genes

In this vignette, we will be retrieving sequences from genbank using phruta. We will need a list of taxa and the names of particular gene regions that need to be sampled. If you have no idea on what genes you should sample - phruta’ve got your back. Just check the “Introduction to the `phruta` R package” and “To export or not export `phruta` outputs” vignettes.

Now, given that we already have some target genes in mind, we can avoid using gene.sampling.retrieve(). We won’t generate an accession table using acc.table.retrieve() at the beginning, mostly because the sampling that we’re using in phruta this time. Therefore, sq.retrieve.direct() is less flexible than sq.retrieve.indirect(). As noted before, sq.retrieve.direct() does its best to directly (i.e. without input from the user) retrieve sequences for a target set of taxa and set of gene regions. You are more likely to catch errors using sq.retrieve.indirect(). However, mistakes will be harder to spot and fix when using sq.retrieve.direct().

Let’s first start by loading phruta!

library(phruta)

Let’s download gene sequences for the taxa in Felis, Vulpes, and Phoca. Tjhis species-level sampling corresponds with the one in other vignettes. We will be sampling only two gene regions using sq.retrieve.direct(). Note that we’re going to be using the list of species for the target genera from the gbif taxonomic backbone (hence db = "gbif"). There are more databases available for this (please check taxize::downstream()).

sq.retrieve.direct(
  clades = c("Felis", "Vulpes", "Phoca"),
  species = "Manis_pentadactyla",
  genes = c("ADORA3", "CYTB"),
  db = "gbif"
)

The code above will create a folder 0.Sequences in your working directory. Note that sq.retrieve.direct() does not return objects to the environment. The rest of this tutorial follow the same structure as the “To export or not export `phruta` outputs” vignette.

For instance, we can use the sq.curate() function to remove highly divergent sequences. This function will also remove observations from species that are not within our target taxonomic groups.

sq.curate(filterTaxonomicCriteria = 'Felis|Vulpes|Phoca|Manis',
                         kingdom = 'animals', 
                         folder = "0.Sequences",
                         removeOutliers = FALSE,
                         minSeqs = 2)

Let’s review the current sampling…

Now, we’ll align the sequences that we just curated. For this, we use sq.aln() with under default parameters.

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

The resulting multiple sequence alignments will be saved to sqs.aln() object, a list. For each of the gene regions, we will have access to the original alignment (Aln.Original), the masked one (Aln.Masked), and information on the masking process. Let’s first check the raw alignments…

Now, the masked alignments…

A figure showing curated alignments

And we’re done for now!! You can compare the sampling between the main three tutorials using sq.retrieve.direct() and sq.retrieve.indirect(). In total, this vignette took 4 minutes to render in my local machine.