Nothing
R/local_taxa_tool.R
In LocaTT: Geographically-Conscious Taxonomic Assignment for Metabarcoding
Defines functions
local_taxa_tool
Documented in
local_taxa_tool
#' Perform Geographically-Conscious Taxonomic Assignment
#'
#' @description Performs taxonomic assignment of DNA metabarcoding sequences while considering geographic location.
#' @details Sequences are BLASTed against a global reference database, and the tool suggests locally occurring species which are most closely related (by taxonomy) to any of the best-matching BLAST hits (by bit score). Optionally, local sister taxonomic groups without reference sequences can be added to the local taxa suggestions by setting the `include_missing` argument to `TRUE`. If a local taxa list is not provided, then local taxa suggestions will be disabled, but all best-matching BLAST hits will still be returned. Alternatively, a reference database containing just the sequences of local species can be used, and local taxa suggestions can be disabled to return all best BLAST matches of local species. The reference database should be formatted with the [`format_reference_database`][format_reference_database()] function, and the local taxa lists can be prepared using the [`get_taxonomies.species_binomials`][get_taxonomies.species_binomials()] and [`get_taxonomies.IUCN`][get_taxonomies.IUCN()] functions. Output field definitions are:
#' * Sequence_name: The query sequence name.
#' * Sequence: The query sequence.
#' * Best_match_references: Species binomials of all best-matching BLAST hits (by bit score) from the reference database.
#' * Best_match_E_value: The E-value associated with the best-matching BLAST hits.
#' * Best_match_bit_score: The bit score associated with the best-matching BLAST hits.
#' * Best_match_query_cover.mean: The mean query cover of all best-matching BLAST hits.
#' * Best_match_query_cover.SD: The standard deviation of query cover of all best-matching BLAST hits.
#' * Best_match_PID.mean: The mean percent identity of all best-matching BLAST hits.
#' * Best_match_PID.SD: The standard deviation of percent identity of all best-matching BLAST hits.
#' * Local_taxa (Field only present if a path to a local taxa list is provided): The finest taxonomic unit(s) which include both any species of the best-matching BLAST hits and any local species. If the species of any of the best-matching BLAST hits are local, then the finest taxonomic unit(s) are at the species level.
#' * Local_species (Field only present if a path to a local taxa list is provided): Species binomials of all local species which belong to the taxonomic unit(s) in the Local_taxa field.
#' * Local_taxa.include_missing (Field only present if both a path to a local taxa list is provided and the `include_missing` argument is set to `TRUE`): Local sister taxonomic groups without reference sequences are added to the local taxa suggestions from the Local_taxa field.
#' * Local_species.include_missing (Field only present if both a path to a local taxa list is provided and `include_missing` argument is set to `TRUE`): Species binomials of all local species which belong to the taxonomic unit(s) in the Local_taxa.include_missing field.
#' @param path_to_query_sequences String specifying path to FASTA file containing sequences to classify. File path cannot contain spaces.
#' @param path_to_BLAST_database String specifying path to BLAST reference database in FASTA format. File path cannot contain spaces.
#' @param path_to_output_file String specifying path to output file of classified sequences in CSV format.
#' @param path_to_local_taxa_list String specifying path to list of local species in CSV format. The file should contain the following fields: 'Common_Name', 'Domain', 'Phylum', 'Class', 'Order', 'Family', 'Genus', 'Species'. There should be no 'NA's or blanks in the taxonomy fields. The species field should contain the binomial name without subspecies or other information below the species level. There should be no duplicate species (*i.e.*, multiple records with the same species binomial and taxonomy) in the local species list. If local taxa suggestions are not desired, set this variable to `NA` (the default).
#' @param include_missing Logical. If `TRUE`, then additional fields are included in the output CSV file in which local sister taxonomic groups without reference sequences are added to the local taxa suggestions. If `FALSE` (the default), then this is not performed.
#' @param blast_e_value Numeric. Maximum E-value of returned BLAST hits (lower E-values are associated with more 'significant' matches). The default is `1e-05`.
#' @param blast_max_target_seqs Numeric. Maximum number of BLAST target sequences returned per query sequence. Enough target sequences should be returned to ensure that all minimum E-value matches are returned for each query sequence. A warning will be produced if this value is not sufficient. The default is `2000`.
#' @param blast_task String specifying BLAST task specification. Use `'megablast'` (the default) to find very similar sequences (*e.g.*, intraspecies or closely related species). Use `'blastn-short'` for sequences shorter than 50 bases. See the blastn program help documentation for additional options and details.
#' @param full_names Logical. If `TRUE`, then full taxonomies are returned in the output CSV file. If `FALSE` (the default), then only the lowest taxonomic levels (*e.g.*, species binomials instead of the full species taxonomies) are returned in the output CSV file.
#' @param underscores Logical. If `TRUE`, then taxa names in the output CSV file use underscores instead of spaces. If `FALSE` (the default), then taxa names in the output CSV file use spaces.
#' @param separator String specifying the separator to use between taxa names in the output CSV file. The default is `', '`.
#' @param blastn_command String specifying path to the blastn program. The default (`'blastn'`) should work for standard BLAST installations. The user can provide a path to the blastn program for non-standard BLAST installations.
#' @param ... Accepts former argument names for backwards compatibility.
#' @returns No return value. Writes an output CSV file with fields defined in the details section.
#' @seealso
#' [`format_reference_database`][format_reference_database()] for formatting reference databases. \cr \cr
#' [`get_taxonomies.species_binomials`][get_taxonomies.species_binomials()] and [`get_taxonomies.IUCN`][get_taxonomies.IUCN()] for creating local taxa lists. \cr \cr
#' [`adjust_taxonomies`][adjust_taxonomies()] for adjusting a taxonomy system.
#' @references A manuscript describing this taxonomic assignment method is in preparation.
#' @examplesIf blast_command_found(blast_command="blastn")
#' # Get path to example query sequences FASTA file.
#' path_to_query_sequences<-system.file("extdata",
#' "example_query_sequences.fasta",
#' package="LocaTT",
#' mustWork=TRUE)
#'
#' # Get path to example BLAST reference database FASTA file.
#' path_to_BLAST_database<-system.file("extdata",
#' "example_blast_database.fasta",
#' package="LocaTT",
#' mustWork=TRUE)
#'
#' # Get path to example local taxa list CSV file.
#' path_to_local_taxa_list<-system.file("extdata",
#' "example_local_taxa_list.csv",
#' package="LocaTT",
#' mustWork=TRUE)
#'
#' # Create a temporary file path for the output CSV file.
#' path_to_output_file<-tempfile(fileext=".csv")
#'
#' # Run the local taxa tool.
#' local_taxa_tool(path_to_query_sequences=path_to_query_sequences,
#' path_to_BLAST_database=path_to_BLAST_database,
#' path_to_output_file=path_to_output_file,
#' path_to_local_taxa_list=path_to_local_taxa_list,
#' include_missing=TRUE,
#' full_names=TRUE,
#' underscores=TRUE)
#' @export
local_taxa_tool<-function(path_to_query_sequences,path_to_BLAST_database,path_to_output_file,path_to_local_taxa_list=NA,include_missing=FALSE,blast_e_value=1e-5,blast_max_target_seqs=2000,blast_task="megablast",full_names=FALSE,underscores=FALSE,separator=", ",blastn_command="blastn",...){
### Ensure backwards compatibility.
# Handle changes in argument names.
## Define data frame relating former to current argument names.
back.compat<-data.frame(arg.former=c("path_to_sequences_to_classify",
"path_to_list_of_local_taxa"),
arg.current=c("path_to_query_sequences",
"path_to_local_taxa_list"),
stringsAsFactors=FALSE)
## Collect extra arguments.
dots<-list(...)
## If there are extra arguments.
if(!missing(...)){
## Loop through each extra argument.
for(i in 1:length(dots)){
## Get the extra argument.
dot<-dots[i]
## If the extra argument is a former argument.
if(names(dot) %in% back.compat$arg.former){
## Get the argument translation record.
arg<-back.compat[back.compat$arg.former==names(dot),]
## Throw an error if the former argument name matches multiple current argument names.
if(nrow(arg) > 1) stop("Former argument name matches multiple current argument names.")
## If the current argument is missing.
if(do.call(missing,list(arg$arg.current))){
## Pass the former argument to the current argument.
assign(x=arg$arg.current,value=dot[[1]])
## Provide a warning that the former argument has been
## renamed to the current argument.
warning(paste0("Former argument '",names(dot),
"' has been renamed to current argument '",
arg$arg.current,"'."))
} else { ## If the current argument is not missing.
## If both old and current arguments are provided, throw an error.
stop(paste0("Multiple equivalent arguments found for current argument '",
arg$arg.current,"' (former argument '",arg$arg.former,"')."))
}
} else { # If the extra argument is not a former argument.
# Throw an error for the unused argument.
stop(paste0("unused argument (",names(dot)," = ",unname(unlist(dot)),")"))
}
}
}
### Check BLAST installation.
# blastn command.
## Throw an error if the blastn command is not a character string.
if(!is.character(blastn_command)) stop("blastn_command must be a character string.")
## Throw an error if the blastn command has multiple elements.
if(length(blastn_command) > 1) stop("blastn_command cannot have multiple elements.")
## Throw an error if the blastn command is NA.
if(is.na(blastn_command)) stop("blastn_command cannot be NA.")
## Throw an error if the blastn command cannot not be found.
if(!blast_command_found(blast_command=blastn_command)) stop("The blastn command could not be found. If using a non-standard installation of BLAST, set the path to the blastn command using the blastn_command argument.")
### Check arguments.
# Path to query sequences.
## Throw an error if the path to query sequences is not a character string.
if(!is.character(path_to_query_sequences)) stop("The path to query sequences must be a character string.")
## Throw an error if the path to query sequences has multiple elements.
if(length(path_to_query_sequences) > 1) stop("The path to query sequences cannot have multiple elements.")
## Throw an error if the path to query sequences is NA.
if(is.na(path_to_query_sequences)) stop("The path to query sequences cannot be NA.")
## Throw an error if the path to query sequences contains spaces.
if(grepl(pattern=" ",x=path_to_query_sequences)) stop("There cannot be spaces in the path to query sequences.")
# Path to BLAST database.
## Throw an error if the path to BLAST database is not a character string.
if(!is.character(path_to_BLAST_database)) stop("The path to BLAST database must be a character string.")
## Throw an error if the path to BLAST database has multiple elements.
if(length(path_to_BLAST_database) > 1) stop("The path to BLAST database cannot have multiple elements.")
## Throw an error if the path to BLAST database is NA.
if(is.na(path_to_BLAST_database)) stop("The path to BLAST database cannot be NA.")
## Throw an error if the path to BLAST database contains spaces.
if(grepl(pattern=" ",x=path_to_BLAST_database)) stop("There cannot be spaces in the path to BLAST database.")
# Path to output file.
## Throw an error if the path to output file is not a character string.
if(!is.character(path_to_output_file)) stop("The path to output file must be a character string.")
## Throw an error if the path to output file has multiple elements.
if(length(path_to_output_file) > 1) stop("The path to output file cannot have multiple elements.")
## Throw an error if the path to output file is NA.
if(is.na(path_to_output_file)) stop("The path to output file cannot be NA.")
# Path to local taxa list.
## Throw an error if the path to local taxa list is not a vector.
if(!is.vector(path_to_local_taxa_list)) stop("The path to local taxa list must be a vector.")
## Throw an error if the path to local taxa list has multiple elements.
if(length(path_to_local_taxa_list) > 1) stop("The path to local taxa list cannot have multiple elements.")
## Throw an error if the path to local taxa list is not a character string, and not NA.
if(!(is.na(path_to_local_taxa_list) | is.character(path_to_local_taxa_list))) stop("The path to local taxa list must be a character string or NA.")
# Include missing argument.
## Throw an error if the include missing argument is not logical.
if(!is.logical(include_missing)) stop("The include missing argument must be logical.")
## Throw an error if the include missing argument has multiple elements.
if(length(include_missing) > 1) stop("The include missing argument cannot have multiple elements.")
## Throw an error if the include missing argument is NA.
if(is.na(include_missing)) stop("The include missing argument cannot be NA.")
# BLAST E-value.
## Throw an error if the BLAST E-value argument is not numeric.
if(!is.numeric(blast_e_value)) stop("The BLAST E-value argument must be numeric.")
## Throw an error if the BLAST E-value argument has multiple elements.
if(length(blast_e_value) > 1) stop("The BLAST E-value argument cannot have multiple elements.")
## Throw an error if the BLAST E-value argument is NA.
if(is.na(blast_e_value)) stop("The BLAST E-value argument cannot be NA.")
# BLAST maximum target sequences.
## Throw an error if the BLAST maximum target sequences argument is not numeric.
if(!is.numeric(blast_max_target_seqs)) stop("The BLAST maximum target sequences argument must be numeric.")
## Throw an error if the BLAST maximum target sequences argument has multiple elements.
if(length(blast_max_target_seqs) > 1) stop("The BLAST maximum target sequences argument cannot have multiple elements.")
## Throw an error if the BLAST maximum target sequences argument is NA.
if(is.na(blast_max_target_seqs)) stop("The BLAST maximum target sequences argument cannot be NA.")
## Throw an error if the BLAST maximum target sequences argument is not an integer value.
if(round(x=blast_max_target_seqs,digits=0)!=blast_max_target_seqs) stop("The BLAST maximum target sequences argument must be an integer value.")
# BLAST task argument.
## Throw an error if the BLAST task argument is not a character string.
if(!is.character(blast_task)) stop("The BLAST task argument must be a character string.")
## Throw an error if the BLAST task argument has multiple elements.
if(length(blast_task) > 1) stop("The BLAST task argument cannot have multiple elements.")
## Throw an error if the BLAST task argument is NA.
if(is.na(blast_task)) stop("The BLAST task argument cannot be NA.")
# Full names argument.
## Throw an error if the full names argument is not logical.
if(!is.logical(full_names)) stop("The full names argument must be logical.")
## Throw an error if the full names argument has multiple elements.
if(length(full_names) > 1) stop("The full names argument cannot have multiple elements.")
## Throw an error if the full names argument is NA.
if(is.na(full_names)) stop("The full names argument cannot be NA.")
# Underscores argument.
## Throw an error if the underscores argument is not logical.
if(!is.logical(underscores)) stop("The underscores argument must be logical.")
## Throw an error if the underscores argument has multiple elements.
if(length(underscores) > 1) stop("The underscores argument cannot have multiple elements.")
## Throw an error if the underscores argument is NA.
if(is.na(underscores)) stop("The underscores argument cannot be NA.")
# Separator argument.
## Throw an error if the separator argument is not a character string.
if(!is.character(separator)) stop("The separator argument must be a character string.")
## Throw an error if the separator argument has multiple elements.
if(length(separator) > 1) stop("The separator argument cannot have multiple elements.")
## Throw an error if the separator argument is NA.
if(is.na(separator)) stop("The separator argument cannot be NA.")
# Throw a warning if include missing is enabled but a local taxa list is not provided.
if(include_missing & is.na(path_to_local_taxa_list)) warning("The include missing argument is ignored if a local taxa list is not provided.")
### Begin operations.
# Read in sequences to classify.
sequences_to_classify<-read.fasta(file=path_to_query_sequences)
colnames(sequences_to_classify)[1]<-"Sequence_name"
# Throw an error if there are NAs or blanks in the fasta file of query sequences.
if(any(is.na(sequences_to_classify) | sequences_to_classify=="")) stop("NAs or blanks are present in the sequence names or sequences of the fasta file containing the query sequences.")
# Throw an error if not all sequence names are unique in the fasta file of query sequences.
if(any(duplicated(sequences_to_classify$Sequence_name))) stop("Not all sequence names are unique in the fasta file containing the query sequences.")
# Perform the command line BLAST algorithm between the query and reference sequences.
blast_output<-system2(command=blastn_command,
args=c(paste0('-task ',blast_task),
paste0('-db ',path_to_BLAST_database),
paste0('-query ',path_to_query_sequences),
'-outfmt "6 qseqid sseqid evalue bitscore qcovs pident"',
paste0('-max_target_seqs ',blast_max_target_seqs),
paste0('-evalue ',blast_e_value)),
stdout=TRUE)
# Check that BLAST returned results for at least one sequence.
if(length(blast_output)==0) stop("BLAST did not return results for any sequences. Try increasing the user-defined minimum E-value parameter and run this function again.")
# Format the BLAST output.
## Get output as a data frame.
blast_output<-as.data.frame(do.call("rbind",strsplit(x=blast_output,split="\t")),
stringsAsFactors=FALSE)
## Provide field names.
colnames(blast_output)<-c("Sequence_name","Matches","E_value",
"Bit_score","Query_cover","PID")
## If spaces in taxa names are desired, replace underscores in matches with spaces.
if(!underscores){
blast_output$Matches<-gsub(pattern="_",replacement=" ",x=blast_output$Matches)
}
## Format the fields for E-value and after as numeric.
blast_output[,3:ncol(blast_output)]<-sapply(X=blast_output[,3:ncol(blast_output)],
FUN=as.numeric)
# Get best matches for each query sequence.
# We'll use the maximum bit score because the E-values
# can have underflow issues (i.e., for a given query sequence, there can
# be multiple hits with very low E-values, which underflow to zero, and
# these hits can have different bit scores). Using the maximum bit score
# is equivalent to using the minimum E-value.
## Get the maximum bit score match for each query sequence.
max_bit_score<-stats::aggregate(Bit_score~Sequence_name,data=blast_output,FUN=max)
## Rename field to maximum bit score.
colnames(max_bit_score)[2]<-"max_bit_score"
## Merge maximum bit score with the BLAST output.
blast_output<-merge(blast_output,max_bit_score)
## Keep just records which have the maximum bit score for each query sequence.
blast_output<-blast_output[blast_output$Bit_score==blast_output$max_bit_score,]
# Get number of best BLAST matches returned for each query sequence.
num_best_BLAST_matches_by_sequence<-table(blast_output$Sequence_name)
# Check whether each query sequence had the maximum number of returned best BLAST matches.
max_best_BLAST_matches_by_sequence<-num_best_BLAST_matches_by_sequence==blast_max_target_seqs
# If any of the unique sequences had no bit scores less than their maximum bit scores,
# then we may not have all of the best-matching taxa for these sequences!
if(any(max_best_BLAST_matches_by_sequence)){
# Get the names of unique sequences which had no bit scores less than their maximum
# bit scores.
sequences_at_max_target_seq_matches<-data.frame(
Sequence_name=names(which(max_best_BLAST_matches_by_sequence)),
stringsAsFactors=FALSE)
# Add the actual sequences to the data frame.
sequences_at_max_target_seq_matches<-merge(x=sequences_at_max_target_seq_matches,
sequences_to_classify,all.x=TRUE)
# Create a single string for each sequence.
sequences_at_max_target_seq_matches$Message<-paste0(
sequences_at_max_target_seq_matches$Sequence_name,
" (",sequences_at_max_target_seq_matches$Sequence,")")
# Issue a warning that blast_max_target_seqs is not high enough for these sequences.
warning(paste0("There are ",nrow(sequences_at_max_target_seq_matches),
" sequences for which a blast_max_target_seqs of ",blast_max_target_seqs,
" was insufficient. These sequences are the following: ",
paste(sequences_at_max_target_seq_matches$Message,collapse=", ")))
}
# Get average query cover and PID for the best-matching references for each sequence.
## Calculate means.
best_match_means<-stats::aggregate(.~Sequence_name,
data=blast_output[,c("Sequence_name","Query_cover","PID")],
FUN=mean)
## Add a mean label to the column names.
colnames(best_match_means)[2:3]<-paste0(colnames(best_match_means)[2:3],".mean")
# Get standard deviation of query cover and PID for the best-matching references for
# each sequence.
## Calculate standard deviation.
best_match_SDs<-stats::aggregate(.~Sequence_name,
data=blast_output[,c("Sequence_name","Query_cover","PID")],
FUN=stats::sd)
## Add a standard deviation label to the column names.
colnames(best_match_SDs)[2:3]<-paste0(colnames(best_match_SDs)[2:3],".SD")
# Merge the averages and standard deviations of query cover and PID for the
# best-matching references for each sequence.
best_match_means_and_SDs<-merge(best_match_means,best_match_SDs)
# Turn query to lower case.
colnames(best_match_means_and_SDs)<-gsub(pattern="^Query_",replacement="query_",
x=colnames(best_match_means_and_SDs))
# Append best match to the numeric fields.
colnames(best_match_means_and_SDs)[2:ncol(best_match_means_and_SDs)]<-paste0(
"Best_match_",colnames(best_match_means_and_SDs)[2:ncol(best_match_means_and_SDs)])
# Get just the sequence, matches, E-value, and bit-score fields in the BLAST output.
blast_output<-blast_output[,c("Sequence_name","Matches","E_value","Bit_score")]
# Define function for getting a given taxonomic level.
getTaxonomicLevel<-function(character_vector,desired_level_numeric){
# Split character vector names by semi-colon.
split_character_vector_names<-strsplit(x=character_vector,split=";")
# Define internal function for collapsing split names at the desired taxonomic level.
collapse_split_name<-function(split_name) paste(split_name[1:desired_level_numeric],
collapse=";")
# Apply the collapse names function across the list of split names.
collapsed_names<-sapply(X=split_character_vector_names,FUN=collapse_split_name)
# Return the collapsed names.
return(collapsed_names)
}
# If a local taxa list is provided.
if(!is.na(path_to_local_taxa_list)){
# Read in local taxa list.
local<-utils::read.csv(file=path_to_local_taxa_list,stringsAsFactors=FALSE)
# Check that the correct fields are present in the local taxa list.
if(!identical(colnames(local),c("Common_Name","Domain","Phylum","Class","Order",
"Family","Genus","Species"))){
stop('The field names in the local taxa list file should be: "Common_Name", "Domain", "Phylum", "Class", "Order", "Family", "Genus", "Species".')
}
# Check that there are no NAs in the taxonomies of the local taxa list.
if(any(is.na(local[,c("Domain","Phylum","Class","Order","Family","Genus","Species")]) | local[,c("Domain","Phylum","Class","Order","Family","Genus","Species")]=="")){
stop("There are NAs or blanks in the taxonomies of the local taxa list file. Please ensure that all taxonomy fields have entries for all records.")
}
# Check that there are no underscores in the taxonomy fields of the local taxa list.
if(any(t(apply(X=local[,c("Domain","Phylum","Class","Order","Family","Genus","Species")],
MARGIN=1,FUN=grepl,pattern="_")))){
stop("There cannot be underscores in the taxonomy fields of the local taxa list file.")
}
# If underscores are desired, replace spaces with underscores in the
# taxonomy fields of the local taxa list.
if(underscores){
local[,c("Domain","Phylum","Class","Order","Family","Genus","Species")]<-data.frame(
lapply(X=local[,c("Domain","Phylum","Class","Order",
"Family","Genus","Species")],
FUN=gsub,pattern=" ",replacement="_"))
}
# Collapse local taxa names by semi-colons.
local$Name<-apply(X=local[,c("Domain","Phylum","Class","Order",
"Family","Genus","Species")],
MARGIN=1,FUN=paste,collapse=";")
# Check that all species are unique.
if(any(duplicated(local$Name))) stop("There are duplicate species (i.e., there are multiple records with the same taxonomy) in the local taxa list file.")
# If include missing is enabled.
if(include_missing){
# Read in the reference database.
ref<-read.fasta(file=path_to_BLAST_database)
# If underscores are not desired.
if(!underscores){
# Replace underscres in reference names with spaces.
ref$Name<-gsub(pattern="_",replacement=" ",x=ref$Name)
}
}
# Get unique sequence names from BLAST output. Sequences not included
# in the BLAST output had no significant matches.
matched_sequences_to_classify<-unique(blast_output$Sequence_name)
# Get the number of unique sequences from BLAST output.
num_seq<-length(matched_sequences_to_classify)
# Create an empty data frame for storing sequence classifications.
sequence_classifications<-data.frame(NULL)
# Loop through each sequence to classify.
for(i in 1:num_seq){
# Get the sequence to classify.
sequence_to_classify<-matched_sequences_to_classify[i]
# Get the BLAST output for the sequence to classify.
sequence_to_classify_df<-blast_output[blast_output$Sequence_name==sequence_to_classify,]
# Get the maximum bit score and its associated minimum E-value.
# Using just the first values since all of them will always be the same.
## Minimum E-value.
min_E_value<-sequence_to_classify_df$E_value[1]
## Maximum bit-score.
bit_score<-sequence_to_classify_df$Bit_score[1]
# Get the unique taxa associated with the highest bit score.
best_matches<-sort(unique(sequence_to_classify_df$Matches))
# Get the species names for the best-matching reference sequences.
if(full_names){ # If full names are requested, get the full taxonomies.
best_match_references<-best_matches
} else { # If full names are not requested, get the species binomials.
best_match_references<-sapply(X=strsplit(x=best_matches,split=";"),FUN="[[",7)
}
# Set the while loop continuation variable.
inclusive_taxonomic_group_found<-"No"
# Set an initial current taxonomic level.
current_taxonomic_level<-7
# Simplify and check taxonomic names until matches are found between the taxonomies
# of the maximum bit score reference sequences and the sequence to classify.
while(inclusive_taxonomic_group_found=="No"){
# Get the current taxonomic level of the best-matches.
best_matches.current_level<-unique(getTaxonomicLevel(character_vector=best_matches,
desired_level_numeric=current_taxonomic_level))
# Get the current taxonomic level of the local species.
local.current_level<-getTaxonomicLevel(character_vector=local$Name,
desired_level_numeric=current_taxonomic_level)
# If any best matches share the current taxonomic level with local species.
if(any(best_matches.current_level %in% local.current_level)){
# Change the while loop continuation variable.
inclusive_taxonomic_group_found<-"Yes"
# Get the inclusive taxa.
inclusive_taxa<-sort(
best_matches.current_level[best_matches.current_level %in% local.current_level])
# Get potential local species.
potential_spp<-sort(local$Name[local.current_level %in% inclusive_taxa])
# If include missing is enabled.
if(include_missing){
# If the current taxonomic level is at phylum or below.
if(current_taxonomic_level > 1){
# Get the inclusive taxonomic level above the current one.
inclusive_taxa.upper_level<-unique(
getTaxonomicLevel(character_vector=inclusive_taxa,
desired_level_numeric=current_taxonomic_level-1))
# Get the taxonomic level above the current one for local species.
local.upper_level<-getTaxonomicLevel(character_vector=local$Name,
desired_level_numeric=current_taxonomic_level-1)
# Create a logical vector specifying whether the local taxonomic
# group is a sister group to the inclusive taxa.
sister_group<-(local.upper_level %in% inclusive_taxa.upper_level)
} else { # If the current taxonomic level is at domain.
# Set the sister group logical vector to TRUE, since all domains
# are sister groups.
sister_group<-TRUE
}
# If not all sister taxonomic groups are the inclusive taxa
# (i.e., if there are sister groups to the inclusive taxa).
if(!all(local.current_level[sister_group] %in% inclusive_taxa)){
# Get the current taxonomic level of the reference sequences.
ref.current_level<-getTaxonomicLevel(character_vector=ref$Name,
desired_level_numeric=current_taxonomic_level)
# Create a logical vector specifying whether the local taxonomic
# group is missing from the reference database.
missing_ref_seqs<-(!(local.current_level %in% ref.current_level))
# Create a logical vector specifying whether the local taxonomic
# group is both a sister group to the inclusive taxa and
# is missing from the reference database.
missing_sister_group<-sister_group & missing_ref_seqs
# If there are any local sister groups to the inclusive taxa
# which are missing from the reference database.
if(any(missing_sister_group)){
# Get the local taxonomic groups sister to the local taxa
# which are missing from the reference database.
local_absent.current_level<-unique(local.current_level[missing_sister_group])
# Get the local species belonging to the sister groups
# which are missing from the reference database.
local_absent.spp<-local$Name[missing_sister_group]
# Add the missing taxa to the inclusive taxa.
inclusive_taxa.include_missing<-sort(c(inclusive_taxa,
local_absent.current_level))
# Add the missing species to the potential species.
potential_spp.include_missing<-sort(c(potential_spp,local_absent.spp))
} else {
# If there are no local sister groups to the inclusive taxa
# which are missing from the reference database.
# Set the include missing taxa to the inclusive taxa.
inclusive_taxa.include_missing<-inclusive_taxa
# Set the include missing species to the potential species.
potential_spp.include_missing<-potential_spp
}
} else {
# If all sister taxonomic groups are the inclusive taxa
# (i.e., if there are no sister groups to the inclusive taxa),
# then there are no additional taxa to add.
# Set the include missing taxa to the inclusive taxa.
inclusive_taxa.include_missing<-inclusive_taxa
# Set the include missing species to the potential species.
potential_spp.include_missing<-potential_spp
}
}
# If full names are not requested, report simplified taxonomies.
if(!full_names){
# Potential local species.
potential_spp<-sapply(X=strsplit(x=potential_spp,split=";"),
FUN="[[",7)
# Inclusive taxa.
inclusive_taxa<-sapply(X=strsplit(x=inclusive_taxa,split=";"),
FUN="[[",current_taxonomic_level)
# If include missing is enabled.
if(include_missing){
# Potential local species (including missing).
potential_spp.include_missing<-sapply(X=strsplit(
x=potential_spp.include_missing,
split=";"),
FUN="[[",7)
# Inclusive taxa (including missing).
inclusive_taxa.include_missing<-sapply(X=strsplit(
x=inclusive_taxa.include_missing,
split=";"),
FUN="[[",current_taxonomic_level)
}
}
# Collect classification information in a data frame.
if(include_missing){
# If include missing is enabled.
sequence_information<-data.frame(Local_taxa=paste(inclusive_taxa,
collapse=separator),
Local_taxa.include_missing=paste(
inclusive_taxa.include_missing,
collapse=separator),
Local_species=paste(potential_spp,
collapse=separator),
Local_species.include_missing=paste(
potential_spp.include_missing,
collapse=separator),
Best_match_references=paste(
best_match_references,
collapse=separator),
Best_match_E_value=min_E_value,
Best_match_bit_score=bit_score,
Sequence_name=sequence_to_classify,
stringsAsFactors=FALSE)
} else {
# If include missing is not enabled.
sequence_information<-data.frame(Local_taxa=paste(inclusive_taxa,
collapse=separator),
Local_species=paste(potential_spp,
collapse=separator),
Best_match_references=paste(
best_match_references,collapse=separator),
Best_match_E_value=min_E_value,
Best_match_bit_score=bit_score,
Sequence_name=sequence_to_classify,
stringsAsFactors=FALSE)
}
} else {
# If no best matches share the current taxonomic level with local species.
# Move up a taxonomic level for the next iteration.
current_taxonomic_level<-current_taxonomic_level-1
# If the current taxonomic level counter reaches zero.
# Then the inclusive taxonomic group is root.
if(current_taxonomic_level==0){
# Change the while loop continuation variable.
inclusive_taxonomic_group_found<-"Yes"
# Collect classification information in a data frame.
if(include_missing){
# If include missing is enabled.
sequence_information<-data.frame(Local_taxa="root",
Local_taxa.include_missing="root",
Local_species=ifelse(full_names,
paste(sort(local$Name),collapse=separator),
paste(local$Species[order(local$Name)],
collapse=separator)),
Local_species.include_missing=ifelse(full_names,
paste(sort(local$Name),collapse=separator),
paste(local$Species[order(local$Name)],
collapse=separator)),
Best_match_references=paste(
best_match_references,collapse=separator),
Best_match_E_value=min_E_value,
Best_match_bit_score=bit_score,
Sequence_name=sequence_to_classify,
stringsAsFactors=FALSE)
} else {
# If include missing is not enabled.
sequence_information<-data.frame(Local_taxa="root",
Local_species=ifelse(full_names,
paste(sort(local$Name),collapse=separator),
paste(local$Species[order(local$Name)],
collapse=separator)),
Best_match_references=paste(
best_match_references,collapse=separator),
Best_match_E_value=min_E_value,
Best_match_bit_score=bit_score,
Sequence_name=sequence_to_classify,
stringsAsFactors=FALSE)
}
}
}
}
# Add sequence information to the storage data frame.
sequence_classifications<-rbind(sequence_classifications,sequence_information)
}
# Merge taxa suggestions with means and SDs of query cover and PID.
sequence_classifications<-merge(sequence_classifications,best_match_means_and_SDs)
# Add sequences to the classification data frame.
df<-merge(sequences_to_classify,sequence_classifications,all=TRUE)
# Keep certain fields in the classification data frame.
if(include_missing){
# If include missing is enabled.
df<-df[,c("Sequence_name","Sequence","Best_match_references","Best_match_E_value",
"Best_match_bit_score","Best_match_query_cover.mean",
"Best_match_query_cover.SD","Best_match_PID.mean","Best_match_PID.SD",
"Local_taxa","Local_species","Local_taxa.include_missing",
"Local_species.include_missing")]
} else {
# If include missing is not enabled.
df<-df[,c("Sequence_name","Sequence","Best_match_references","Best_match_E_value",
"Best_match_bit_score","Best_match_query_cover.mean",
"Best_match_query_cover.SD","Best_match_PID.mean","Best_match_PID.SD",
"Local_taxa","Local_species")]
}
} else { # If a local taxa list is not provided.
# Get unique best-matches for each sequence.
blast_output<-unique(blast_output)
# Order the records in the best-matches by sequence name and taxon name.
blast_output<-blast_output[order(blast_output$Sequence_name,blast_output$Matches),]
# Get the species names for the best-matching reference sequences.
if(full_names){ # If full names are requested, get the full taxonomies.
blast_output$Species_names<-blast_output$Matches
} else { # If full names are not requested, get the species binomials.
blast_output$Species_names<-sapply(X=strsplit(x=blast_output$Matches,split=";"),
FUN="[[",7)
}
# Get best-matching references.
best_matches<-stats::aggregate(Species_names~Sequence_name,data=blast_output,
FUN=paste,collapse=separator)
# Get the E-values of best-matches.
best_matches_E_value<-stats::aggregate(E_value~Sequence_name,data=blast_output,
FUN="[",1)
# Get the bit-scores of best-matches.
best_matches_Bit_score<-stats::aggregate(Bit_score~Sequence_name,data=blast_output,
FUN="[",1)
# Combine information on best-matching references.
best_matches<-merge(best_matches,best_matches_E_value)
best_matches<-merge(best_matches,best_matches_Bit_score)
best_matches<-merge(best_matches,best_match_means_and_SDs)
# Add sequences to the classification data frame.
df<-merge(sequences_to_classify,best_matches,all=TRUE)
# Rename certain fields.
colnames(df)[which(colnames(df)=="Species_names")]<-"Best_match_references"
colnames(df)[which(colnames(df)=="E_value")]<-"Best_match_E_value"
colnames(df)[which(colnames(df)=="Bit_score")]<-"Best_match_bit_score"
# Keep certain fields in the classification data frame.
df<-df[,c("Sequence_name","Sequence","Best_match_references","Best_match_E_value",
"Best_match_bit_score","Best_match_query_cover.mean",
"Best_match_query_cover.SD","Best_match_PID.mean","Best_match_PID.SD")]
}
# Set NA values in rows without an E-value to suggest that no significant similarity
# was found.
df[is.na(df$Best_match_E_value),3:ncol(df)]<-"No significant similarity found"
# Write out classified sequences.
utils::write.csv(x=df,file=path_to_output_file,row.names=FALSE)
}
Try the LocaTT package in your browser
Any scripts or data that you put into this service are public.
LocaTT documentation built on June 14, 2026, 1:06 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions local_taxa_tool
Documented in local_taxa_tool
#' Perform Geographically-Conscious Taxonomic Assignment
#'
#' @description Performs taxonomic assignment of DNA metabarcoding sequences while considering geographic location.
#' @details Sequences are BLASTed against a global reference database, and the tool suggests locally occurring species which are most closely related (by taxonomy) to any of the best-matching BLAST hits (by bit score). Optionally, local sister taxonomic groups without reference sequences can be added to the local taxa suggestions by setting the `include_missing` argument to `TRUE`. If a local taxa list is not provided, then local taxa suggestions will be disabled, but all best-matching BLAST hits will still be returned. Alternatively, a reference database containing just the sequences of local species can be used, and local taxa suggestions can be disabled to return all best BLAST matches of local species. The reference database should be formatted with the [`format_reference_database`][format_reference_database()] function, and the local taxa lists can be prepared using the [`get_taxonomies.species_binomials`][get_taxonomies.species_binomials()] and [`get_taxonomies.IUCN`][get_taxonomies.IUCN()] functions. Output field definitions are:
#' * Sequence_name: The query sequence name.
#' * Sequence: The query sequence.
#' * Best_match_references: Species binomials of all best-matching BLAST hits (by bit score) from the reference database.
#' * Best_match_E_value: The E-value associated with the best-matching BLAST hits.
#' * Best_match_bit_score: The bit score associated with the best-matching BLAST hits.
#' * Best_match_query_cover.mean: The mean query cover of all best-matching BLAST hits.
#' * Best_match_query_cover.SD: The standard deviation of query cover of all best-matching BLAST hits.
#' * Best_match_PID.mean: The mean percent identity of all best-matching BLAST hits.
#' * Best_match_PID.SD: The standard deviation of percent identity of all best-matching BLAST hits.
#' * Local_taxa (Field only present if a path to a local taxa list is provided): The finest taxonomic unit(s) which include both any species of the best-matching BLAST hits and any local species. If the species of any of the best-matching BLAST hits are local, then the finest taxonomic unit(s) are at the species level.
#' * Local_species (Field only present if a path to a local taxa list is provided): Species binomials of all local species which belong to the taxonomic unit(s) in the Local_taxa field.
#' * Local_taxa.include_missing (Field only present if both a path to a local taxa list is provided and the `include_missing` argument is set to `TRUE`): Local sister taxonomic groups without reference sequences are added to the local taxa suggestions from the Local_taxa field.
#' * Local_species.include_missing (Field only present if both a path to a local taxa list is provided and `include_missing` argument is set to `TRUE`): Species binomials of all local species which belong to the taxonomic unit(s) in the Local_taxa.include_missing field.
#' @param path_to_query_sequences String specifying path to FASTA file containing sequences to classify. File path cannot contain spaces.
#' @param path_to_BLAST_database String specifying path to BLAST reference database in FASTA format. File path cannot contain spaces.
#' @param path_to_output_file String specifying path to output file of classified sequences in CSV format.
#' @param path_to_local_taxa_list String specifying path to list of local species in CSV format. The file should contain the following fields: 'Common_Name', 'Domain', 'Phylum', 'Class', 'Order', 'Family', 'Genus', 'Species'. There should be no 'NA's or blanks in the taxonomy fields. The species field should contain the binomial name without subspecies or other information below the species level. There should be no duplicate species (*i.e.*, multiple records with the same species binomial and taxonomy) in the local species list. If local taxa suggestions are not desired, set this variable to `NA` (the default).
#' @param include_missing Logical. If `TRUE`, then additional fields are included in the output CSV file in which local sister taxonomic groups without reference sequences are added to the local taxa suggestions. If `FALSE` (the default), then this is not performed.
#' @param blast_e_value Numeric. Maximum E-value of returned BLAST hits (lower E-values are associated with more 'significant' matches). The default is `1e-05`.
#' @param blast_max_target_seqs Numeric. Maximum number of BLAST target sequences returned per query sequence. Enough target sequences should be returned to ensure that all minimum E-value matches are returned for each query sequence. A warning will be produced if this value is not sufficient. The default is `2000`.
#' @param blast_task String specifying BLAST task specification. Use `'megablast'` (the default) to find very similar sequences (*e.g.*, intraspecies or closely related species). Use `'blastn-short'` for sequences shorter than 50 bases. See the blastn program help documentation for additional options and details.
#' @param full_names Logical. If `TRUE`, then full taxonomies are returned in the output CSV file. If `FALSE` (the default), then only the lowest taxonomic levels (*e.g.*, species binomials instead of the full species taxonomies) are returned in the output CSV file.
#' @param underscores Logical. If `TRUE`, then taxa names in the output CSV file use underscores instead of spaces. If `FALSE` (the default), then taxa names in the output CSV file use spaces.
#' @param separator String specifying the separator to use between taxa names in the output CSV file. The default is `', '`.
#' @param blastn_command String specifying path to the blastn program. The default (`'blastn'`) should work for standard BLAST installations. The user can provide a path to the blastn program for non-standard BLAST installations.
#' @param ... Accepts former argument names for backwards compatibility.
#' @returns No return value. Writes an output CSV file with fields defined in the details section.
#' @seealso
#' [`format_reference_database`][format_reference_database()] for formatting reference databases. \cr \cr
#' [`get_taxonomies.species_binomials`][get_taxonomies.species_binomials()] and [`get_taxonomies.IUCN`][get_taxonomies.IUCN()] for creating local taxa lists. \cr \cr
#' [`adjust_taxonomies`][adjust_taxonomies()] for adjusting a taxonomy system.
#' @references A manuscript describing this taxonomic assignment method is in preparation.
#' @examplesIf blast_command_found(blast_command="blastn")
#' # Get path to example query sequences FASTA file.
#' path_to_query_sequences<-system.file("extdata",
#' "example_query_sequences.fasta",
#' package="LocaTT",
#' mustWork=TRUE)
#'
#' # Get path to example BLAST reference database FASTA file.
#' path_to_BLAST_database<-system.file("extdata",
#' "example_blast_database.fasta",
#' package="LocaTT",
#' mustWork=TRUE)
#'
#' # Get path to example local taxa list CSV file.
#' path_to_local_taxa_list<-system.file("extdata",
#' "example_local_taxa_list.csv",
#' package="LocaTT",
#' mustWork=TRUE)
#'
#' # Create a temporary file path for the output CSV file.
#' path_to_output_file<-tempfile(fileext=".csv")
#'
#' # Run the local taxa tool.
#' local_taxa_tool(path_to_query_sequences=path_to_query_sequences,
#' path_to_BLAST_database=path_to_BLAST_database,
#' path_to_output_file=path_to_output_file,
#' path_to_local_taxa_list=path_to_local_taxa_list,
#' include_missing=TRUE,
#' full_names=TRUE,
#' underscores=TRUE)
#' @export
local_taxa_tool<-function(path_to_query_sequences,path_to_BLAST_database,path_to_output_file,path_to_local_taxa_list=NA,include_missing=FALSE,blast_e_value=1e-5,blast_max_target_seqs=2000,blast_task="megablast",full_names=FALSE,underscores=FALSE,separator=", ",blastn_command="blastn",...){
### Ensure backwards compatibility.
# Handle changes in argument names.
## Define data frame relating former to current argument names.
back.compat<-data.frame(arg.former=c("path_to_sequences_to_classify",
"path_to_list_of_local_taxa"),
arg.current=c("path_to_query_sequences",
"path_to_local_taxa_list"),
stringsAsFactors=FALSE)
## Collect extra arguments.
dots<-list(...)
## If there are extra arguments.
if(!missing(...)){
## Loop through each extra argument.
for(i in 1:length(dots)){
## Get the extra argument.
dot<-dots[i]
## If the extra argument is a former argument.
if(names(dot) %in% back.compat$arg.former){
## Get the argument translation record.
arg<-back.compat[back.compat$arg.former==names(dot),]
## Throw an error if the former argument name matches multiple current argument names.
if(nrow(arg) > 1) stop("Former argument name matches multiple current argument names.")
## If the current argument is missing.
if(do.call(missing,list(arg$arg.current))){
## Pass the former argument to the current argument.
assign(x=arg$arg.current,value=dot[[1]])
## Provide a warning that the former argument has been
## renamed to the current argument.
warning(paste0("Former argument '",names(dot),
"' has been renamed to current argument '",
arg$arg.current,"'."))
} else { ## If the current argument is not missing.
## If both old and current arguments are provided, throw an error.
stop(paste0("Multiple equivalent arguments found for current argument '",
arg$arg.current,"' (former argument '",arg$arg.former,"')."))
}
} else { # If the extra argument is not a former argument.
# Throw an error for the unused argument.
stop(paste0("unused argument (",names(dot)," = ",unname(unlist(dot)),")"))
}
}
}
### Check BLAST installation.
# blastn command.
## Throw an error if the blastn command is not a character string.
if(!is.character(blastn_command)) stop("blastn_command must be a character string.")
## Throw an error if the blastn command has multiple elements.
if(length(blastn_command) > 1) stop("blastn_command cannot have multiple elements.")
## Throw an error if the blastn command is NA.
if(is.na(blastn_command)) stop("blastn_command cannot be NA.")
## Throw an error if the blastn command cannot not be found.
if(!blast_command_found(blast_command=blastn_command)) stop("The blastn command could not be found. If using a non-standard installation of BLAST, set the path to the blastn command using the blastn_command argument.")
### Check arguments.
# Path to query sequences.
## Throw an error if the path to query sequences is not a character string.
if(!is.character(path_to_query_sequences)) stop("The path to query sequences must be a character string.")
## Throw an error if the path to query sequences has multiple elements.
if(length(path_to_query_sequences) > 1) stop("The path to query sequences cannot have multiple elements.")
## Throw an error if the path to query sequences is NA.
if(is.na(path_to_query_sequences)) stop("The path to query sequences cannot be NA.")
## Throw an error if the path to query sequences contains spaces.
if(grepl(pattern=" ",x=path_to_query_sequences)) stop("There cannot be spaces in the path to query sequences.")
# Path to BLAST database.
## Throw an error if the path to BLAST database is not a character string.
if(!is.character(path_to_BLAST_database)) stop("The path to BLAST database must be a character string.")
## Throw an error if the path to BLAST database has multiple elements.
if(length(path_to_BLAST_database) > 1) stop("The path to BLAST database cannot have multiple elements.")
## Throw an error if the path to BLAST database is NA.
if(is.na(path_to_BLAST_database)) stop("The path to BLAST database cannot be NA.")
## Throw an error if the path to BLAST database contains spaces.
if(grepl(pattern=" ",x=path_to_BLAST_database)) stop("There cannot be spaces in the path to BLAST database.")
# Path to output file.
## Throw an error if the path to output file is not a character string.
if(!is.character(path_to_output_file)) stop("The path to output file must be a character string.")
## Throw an error if the path to output file has multiple elements.
if(length(path_to_output_file) > 1) stop("The path to output file cannot have multiple elements.")
## Throw an error if the path to output file is NA.
if(is.na(path_to_output_file)) stop("The path to output file cannot be NA.")
# Path to local taxa list.
## Throw an error if the path to local taxa list is not a vector.
if(!is.vector(path_to_local_taxa_list)) stop("The path to local taxa list must be a vector.")
## Throw an error if the path to local taxa list has multiple elements.
if(length(path_to_local_taxa_list) > 1) stop("The path to local taxa list cannot have multiple elements.")
## Throw an error if the path to local taxa list is not a character string, and not NA.
if(!(is.na(path_to_local_taxa_list) | is.character(path_to_local_taxa_list))) stop("The path to local taxa list must be a character string or NA.")
# Include missing argument.
## Throw an error if the include missing argument is not logical.
if(!is.logical(include_missing)) stop("The include missing argument must be logical.")
## Throw an error if the include missing argument has multiple elements.
if(length(include_missing) > 1) stop("The include missing argument cannot have multiple elements.")
## Throw an error if the include missing argument is NA.
if(is.na(include_missing)) stop("The include missing argument cannot be NA.")
# BLAST E-value.
## Throw an error if the BLAST E-value argument is not numeric.
if(!is.numeric(blast_e_value)) stop("The BLAST E-value argument must be numeric.")
## Throw an error if the BLAST E-value argument has multiple elements.
if(length(blast_e_value) > 1) stop("The BLAST E-value argument cannot have multiple elements.")
## Throw an error if the BLAST E-value argument is NA.
if(is.na(blast_e_value)) stop("The BLAST E-value argument cannot be NA.")
# BLAST maximum target sequences.
## Throw an error if the BLAST maximum target sequences argument is not numeric.
if(!is.numeric(blast_max_target_seqs)) stop("The BLAST maximum target sequences argument must be numeric.")
## Throw an error if the BLAST maximum target sequences argument has multiple elements.
if(length(blast_max_target_seqs) > 1) stop("The BLAST maximum target sequences argument cannot have multiple elements.")
## Throw an error if the BLAST maximum target sequences argument is NA.
if(is.na(blast_max_target_seqs)) stop("The BLAST maximum target sequences argument cannot be NA.")
## Throw an error if the BLAST maximum target sequences argument is not an integer value.
if(round(x=blast_max_target_seqs,digits=0)!=blast_max_target_seqs) stop("The BLAST maximum target sequences argument must be an integer value.")
# BLAST task argument.
## Throw an error if the BLAST task argument is not a character string.
if(!is.character(blast_task)) stop("The BLAST task argument must be a character string.")
## Throw an error if the BLAST task argument has multiple elements.
if(length(blast_task) > 1) stop("The BLAST task argument cannot have multiple elements.")
## Throw an error if the BLAST task argument is NA.
if(is.na(blast_task)) stop("The BLAST task argument cannot be NA.")
# Full names argument.
## Throw an error if the full names argument is not logical.
if(!is.logical(full_names)) stop("The full names argument must be logical.")
## Throw an error if the full names argument has multiple elements.
if(length(full_names) > 1) stop("The full names argument cannot have multiple elements.")
## Throw an error if the full names argument is NA.
if(is.na(full_names)) stop("The full names argument cannot be NA.")
# Underscores argument.
## Throw an error if the underscores argument is not logical.
if(!is.logical(underscores)) stop("The underscores argument must be logical.")
## Throw an error if the underscores argument has multiple elements.
if(length(underscores) > 1) stop("The underscores argument cannot have multiple elements.")
## Throw an error if the underscores argument is NA.
if(is.na(underscores)) stop("The underscores argument cannot be NA.")
# Separator argument.
## Throw an error if the separator argument is not a character string.
if(!is.character(separator)) stop("The separator argument must be a character string.")
## Throw an error if the separator argument has multiple elements.
if(length(separator) > 1) stop("The separator argument cannot have multiple elements.")
## Throw an error if the separator argument is NA.
if(is.na(separator)) stop("The separator argument cannot be NA.")
# Throw a warning if include missing is enabled but a local taxa list is not provided.
if(include_missing & is.na(path_to_local_taxa_list)) warning("The include missing argument is ignored if a local taxa list is not provided.")
### Begin operations.
# Read in sequences to classify.
sequences_to_classify<-read.fasta(file=path_to_query_sequences)
colnames(sequences_to_classify)[1]<-"Sequence_name"
# Throw an error if there are NAs or blanks in the fasta file of query sequences.
if(any(is.na(sequences_to_classify) | sequences_to_classify=="")) stop("NAs or blanks are present in the sequence names or sequences of the fasta file containing the query sequences.")
# Throw an error if not all sequence names are unique in the fasta file of query sequences.
if(any(duplicated(sequences_to_classify$Sequence_name))) stop("Not all sequence names are unique in the fasta file containing the query sequences.")
# Perform the command line BLAST algorithm between the query and reference sequences.
blast_output<-system2(command=blastn_command,
args=c(paste0('-task ',blast_task),
paste0('-db ',path_to_BLAST_database),
paste0('-query ',path_to_query_sequences),
'-outfmt "6 qseqid sseqid evalue bitscore qcovs pident"',
paste0('-max_target_seqs ',blast_max_target_seqs),
paste0('-evalue ',blast_e_value)),
stdout=TRUE)
# Check that BLAST returned results for at least one sequence.
if(length(blast_output)==0) stop("BLAST did not return results for any sequences. Try increasing the user-defined minimum E-value parameter and run this function again.")
# Format the BLAST output.
## Get output as a data frame.
blast_output<-as.data.frame(do.call("rbind",strsplit(x=blast_output,split="\t")),
stringsAsFactors=FALSE)
## Provide field names.
colnames(blast_output)<-c("Sequence_name","Matches","E_value",
"Bit_score","Query_cover","PID")
## If spaces in taxa names are desired, replace underscores in matches with spaces.
if(!underscores){
blast_output$Matches<-gsub(pattern="_",replacement=" ",x=blast_output$Matches)
}
## Format the fields for E-value and after as numeric.
blast_output[,3:ncol(blast_output)]<-sapply(X=blast_output[,3:ncol(blast_output)],
FUN=as.numeric)
# Get best matches for each query sequence.
# We'll use the maximum bit score because the E-values
# can have underflow issues (i.e., for a given query sequence, there can
# be multiple hits with very low E-values, which underflow to zero, and
# these hits can have different bit scores). Using the maximum bit score
# is equivalent to using the minimum E-value.
## Get the maximum bit score match for each query sequence.
max_bit_score<-stats::aggregate(Bit_score~Sequence_name,data=blast_output,FUN=max)
## Rename field to maximum bit score.
colnames(max_bit_score)[2]<-"max_bit_score"
## Merge maximum bit score with the BLAST output.
blast_output<-merge(blast_output,max_bit_score)
## Keep just records which have the maximum bit score for each query sequence.
blast_output<-blast_output[blast_output$Bit_score==blast_output$max_bit_score,]
# Get number of best BLAST matches returned for each query sequence.
num_best_BLAST_matches_by_sequence<-table(blast_output$Sequence_name)
# Check whether each query sequence had the maximum number of returned best BLAST matches.
max_best_BLAST_matches_by_sequence<-num_best_BLAST_matches_by_sequence==blast_max_target_seqs
# If any of the unique sequences had no bit scores less than their maximum bit scores,
# then we may not have all of the best-matching taxa for these sequences!
if(any(max_best_BLAST_matches_by_sequence)){
# Get the names of unique sequences which had no bit scores less than their maximum
# bit scores.
sequences_at_max_target_seq_matches<-data.frame(
Sequence_name=names(which(max_best_BLAST_matches_by_sequence)),
stringsAsFactors=FALSE)
# Add the actual sequences to the data frame.
sequences_at_max_target_seq_matches<-merge(x=sequences_at_max_target_seq_matches,
sequences_to_classify,all.x=TRUE)
# Create a single string for each sequence.
sequences_at_max_target_seq_matches$Message<-paste0(
sequences_at_max_target_seq_matches$Sequence_name,
" (",sequences_at_max_target_seq_matches$Sequence,")")
# Issue a warning that blast_max_target_seqs is not high enough for these sequences.
warning(paste0("There are ",nrow(sequences_at_max_target_seq_matches),
" sequences for which a blast_max_target_seqs of ",blast_max_target_seqs,
" was insufficient. These sequences are the following: ",
paste(sequences_at_max_target_seq_matches$Message,collapse=", ")))
}
# Get average query cover and PID for the best-matching references for each sequence.
## Calculate means.
best_match_means<-stats::aggregate(.~Sequence_name,
data=blast_output[,c("Sequence_name","Query_cover","PID")],
FUN=mean)
## Add a mean label to the column names.
colnames(best_match_means)[2:3]<-paste0(colnames(best_match_means)[2:3],".mean")
# Get standard deviation of query cover and PID for the best-matching references for
# each sequence.
## Calculate standard deviation.
best_match_SDs<-stats::aggregate(.~Sequence_name,
data=blast_output[,c("Sequence_name","Query_cover","PID")],
FUN=stats::sd)
## Add a standard deviation label to the column names.
colnames(best_match_SDs)[2:3]<-paste0(colnames(best_match_SDs)[2:3],".SD")
# Merge the averages and standard deviations of query cover and PID for the
# best-matching references for each sequence.
best_match_means_and_SDs<-merge(best_match_means,best_match_SDs)
# Turn query to lower case.
colnames(best_match_means_and_SDs)<-gsub(pattern="^Query_",replacement="query_",
x=colnames(best_match_means_and_SDs))
# Append best match to the numeric fields.
colnames(best_match_means_and_SDs)[2:ncol(best_match_means_and_SDs)]<-paste0(
"Best_match_",colnames(best_match_means_and_SDs)[2:ncol(best_match_means_and_SDs)])
# Get just the sequence, matches, E-value, and bit-score fields in the BLAST output.
blast_output<-blast_output[,c("Sequence_name","Matches","E_value","Bit_score")]
# Define function for getting a given taxonomic level.
getTaxonomicLevel<-function(character_vector,desired_level_numeric){
# Split character vector names by semi-colon.
split_character_vector_names<-strsplit(x=character_vector,split=";")
# Define internal function for collapsing split names at the desired taxonomic level.
collapse_split_name<-function(split_name) paste(split_name[1:desired_level_numeric],
collapse=";")
# Apply the collapse names function across the list of split names.
collapsed_names<-sapply(X=split_character_vector_names,FUN=collapse_split_name)
# Return the collapsed names.
return(collapsed_names)
}
# If a local taxa list is provided.
if(!is.na(path_to_local_taxa_list)){
# Read in local taxa list.
local<-utils::read.csv(file=path_to_local_taxa_list,stringsAsFactors=FALSE)
# Check that the correct fields are present in the local taxa list.
if(!identical(colnames(local),c("Common_Name","Domain","Phylum","Class","Order",
"Family","Genus","Species"))){
stop('The field names in the local taxa list file should be: "Common_Name", "Domain", "Phylum", "Class", "Order", "Family", "Genus", "Species".')
}
# Check that there are no NAs in the taxonomies of the local taxa list.
if(any(is.na(local[,c("Domain","Phylum","Class","Order","Family","Genus","Species")]) | local[,c("Domain","Phylum","Class","Order","Family","Genus","Species")]=="")){
stop("There are NAs or blanks in the taxonomies of the local taxa list file. Please ensure that all taxonomy fields have entries for all records.")
}
# Check that there are no underscores in the taxonomy fields of the local taxa list.
if(any(t(apply(X=local[,c("Domain","Phylum","Class","Order","Family","Genus","Species")],
MARGIN=1,FUN=grepl,pattern="_")))){
stop("There cannot be underscores in the taxonomy fields of the local taxa list file.")
}
# If underscores are desired, replace spaces with underscores in the
# taxonomy fields of the local taxa list.
if(underscores){
local[,c("Domain","Phylum","Class","Order","Family","Genus","Species")]<-data.frame(
lapply(X=local[,c("Domain","Phylum","Class","Order",
"Family","Genus","Species")],
FUN=gsub,pattern=" ",replacement="_"))
}
# Collapse local taxa names by semi-colons.
local$Name<-apply(X=local[,c("Domain","Phylum","Class","Order",
"Family","Genus","Species")],
MARGIN=1,FUN=paste,collapse=";")
# Check that all species are unique.
if(any(duplicated(local$Name))) stop("There are duplicate species (i.e., there are multiple records with the same taxonomy) in the local taxa list file.")
# If include missing is enabled.
if(include_missing){
# Read in the reference database.
ref<-read.fasta(file=path_to_BLAST_database)
# If underscores are not desired.
if(!underscores){
# Replace underscres in reference names with spaces.
ref$Name<-gsub(pattern="_",replacement=" ",x=ref$Name)
}
}
# Get unique sequence names from BLAST output. Sequences not included
# in the BLAST output had no significant matches.
matched_sequences_to_classify<-unique(blast_output$Sequence_name)
# Get the number of unique sequences from BLAST output.
num_seq<-length(matched_sequences_to_classify)
# Create an empty data frame for storing sequence classifications.
sequence_classifications<-data.frame(NULL)
# Loop through each sequence to classify.
for(i in 1:num_seq){
# Get the sequence to classify.
sequence_to_classify<-matched_sequences_to_classify[i]
# Get the BLAST output for the sequence to classify.
sequence_to_classify_df<-blast_output[blast_output$Sequence_name==sequence_to_classify,]
# Get the maximum bit score and its associated minimum E-value.
# Using just the first values since all of them will always be the same.
## Minimum E-value.
min_E_value<-sequence_to_classify_df$E_value[1]
## Maximum bit-score.
bit_score<-sequence_to_classify_df$Bit_score[1]
# Get the unique taxa associated with the highest bit score.
best_matches<-sort(unique(sequence_to_classify_df$Matches))
# Get the species names for the best-matching reference sequences.
if(full_names){ # If full names are requested, get the full taxonomies.
best_match_references<-best_matches
} else { # If full names are not requested, get the species binomials.
best_match_references<-sapply(X=strsplit(x=best_matches,split=";"),FUN="[[",7)
}
# Set the while loop continuation variable.
inclusive_taxonomic_group_found<-"No"
# Set an initial current taxonomic level.
current_taxonomic_level<-7
# Simplify and check taxonomic names until matches are found between the taxonomies
# of the maximum bit score reference sequences and the sequence to classify.
while(inclusive_taxonomic_group_found=="No"){
# Get the current taxonomic level of the best-matches.
best_matches.current_level<-unique(getTaxonomicLevel(character_vector=best_matches,
desired_level_numeric=current_taxonomic_level))
# Get the current taxonomic level of the local species.
local.current_level<-getTaxonomicLevel(character_vector=local$Name,
desired_level_numeric=current_taxonomic_level)
# If any best matches share the current taxonomic level with local species.
if(any(best_matches.current_level %in% local.current_level)){
# Change the while loop continuation variable.
inclusive_taxonomic_group_found<-"Yes"
# Get the inclusive taxa.
inclusive_taxa<-sort(
best_matches.current_level[best_matches.current_level %in% local.current_level])
# Get potential local species.
potential_spp<-sort(local$Name[local.current_level %in% inclusive_taxa])
# If include missing is enabled.
if(include_missing){
# If the current taxonomic level is at phylum or below.
if(current_taxonomic_level > 1){
# Get the inclusive taxonomic level above the current one.
inclusive_taxa.upper_level<-unique(
getTaxonomicLevel(character_vector=inclusive_taxa,
desired_level_numeric=current_taxonomic_level-1))
# Get the taxonomic level above the current one for local species.
local.upper_level<-getTaxonomicLevel(character_vector=local$Name,
desired_level_numeric=current_taxonomic_level-1)
# Create a logical vector specifying whether the local taxonomic
# group is a sister group to the inclusive taxa.
sister_group<-(local.upper_level %in% inclusive_taxa.upper_level)
} else { # If the current taxonomic level is at domain.
# Set the sister group logical vector to TRUE, since all domains
# are sister groups.
sister_group<-TRUE
}
# If not all sister taxonomic groups are the inclusive taxa
# (i.e., if there are sister groups to the inclusive taxa).
if(!all(local.current_level[sister_group] %in% inclusive_taxa)){
# Get the current taxonomic level of the reference sequences.
ref.current_level<-getTaxonomicLevel(character_vector=ref$Name,
desired_level_numeric=current_taxonomic_level)
# Create a logical vector specifying whether the local taxonomic
# group is missing from the reference database.
missing_ref_seqs<-(!(local.current_level %in% ref.current_level))
# Create a logical vector specifying whether the local taxonomic
# group is both a sister group to the inclusive taxa and
# is missing from the reference database.
missing_sister_group<-sister_group & missing_ref_seqs
# If there are any local sister groups to the inclusive taxa
# which are missing from the reference database.
if(any(missing_sister_group)){
# Get the local taxonomic groups sister to the local taxa
# which are missing from the reference database.
local_absent.current_level<-unique(local.current_level[missing_sister_group])
# Get the local species belonging to the sister groups
# which are missing from the reference database.
local_absent.spp<-local$Name[missing_sister_group]
# Add the missing taxa to the inclusive taxa.
inclusive_taxa.include_missing<-sort(c(inclusive_taxa,
local_absent.current_level))
# Add the missing species to the potential species.
potential_spp.include_missing<-sort(c(potential_spp,local_absent.spp))
} else {
# If there are no local sister groups to the inclusive taxa
# which are missing from the reference database.
# Set the include missing taxa to the inclusive taxa.
inclusive_taxa.include_missing<-inclusive_taxa
# Set the include missing species to the potential species.
potential_spp.include_missing<-potential_spp
}
} else {
# If all sister taxonomic groups are the inclusive taxa
# (i.e., if there are no sister groups to the inclusive taxa),
# then there are no additional taxa to add.
# Set the include missing taxa to the inclusive taxa.
inclusive_taxa.include_missing<-inclusive_taxa
# Set the include missing species to the potential species.
potential_spp.include_missing<-potential_spp
}
}
# If full names are not requested, report simplified taxonomies.
if(!full_names){
# Potential local species.
potential_spp<-sapply(X=strsplit(x=potential_spp,split=";"),
FUN="[[",7)
# Inclusive taxa.
inclusive_taxa<-sapply(X=strsplit(x=inclusive_taxa,split=";"),
FUN="[[",current_taxonomic_level)
# If include missing is enabled.
if(include_missing){
# Potential local species (including missing).
potential_spp.include_missing<-sapply(X=strsplit(
x=potential_spp.include_missing,
split=";"),
FUN="[[",7)
# Inclusive taxa (including missing).
inclusive_taxa.include_missing<-sapply(X=strsplit(
x=inclusive_taxa.include_missing,
split=";"),
FUN="[[",current_taxonomic_level)
}
}
# Collect classification information in a data frame.
if(include_missing){
# If include missing is enabled.
sequence_information<-data.frame(Local_taxa=paste(inclusive_taxa,
collapse=separator),
Local_taxa.include_missing=paste(
inclusive_taxa.include_missing,
collapse=separator),
Local_species=paste(potential_spp,
collapse=separator),
Local_species.include_missing=paste(
potential_spp.include_missing,
collapse=separator),
Best_match_references=paste(
best_match_references,
collapse=separator),
Best_match_E_value=min_E_value,
Best_match_bit_score=bit_score,
Sequence_name=sequence_to_classify,
stringsAsFactors=FALSE)
} else {
# If include missing is not enabled.
sequence_information<-data.frame(Local_taxa=paste(inclusive_taxa,
collapse=separator),
Local_species=paste(potential_spp,
collapse=separator),
Best_match_references=paste(
best_match_references,collapse=separator),
Best_match_E_value=min_E_value,
Best_match_bit_score=bit_score,
Sequence_name=sequence_to_classify,
stringsAsFactors=FALSE)
}
} else {
# If no best matches share the current taxonomic level with local species.
# Move up a taxonomic level for the next iteration.
current_taxonomic_level<-current_taxonomic_level-1
# If the current taxonomic level counter reaches zero.
# Then the inclusive taxonomic group is root.
if(current_taxonomic_level==0){
# Change the while loop continuation variable.
inclusive_taxonomic_group_found<-"Yes"
# Collect classification information in a data frame.
if(include_missing){
# If include missing is enabled.
sequence_information<-data.frame(Local_taxa="root",
Local_taxa.include_missing="root",
Local_species=ifelse(full_names,
paste(sort(local$Name),collapse=separator),
paste(local$Species[order(local$Name)],
collapse=separator)),
Local_species.include_missing=ifelse(full_names,
paste(sort(local$Name),collapse=separator),
paste(local$Species[order(local$Name)],
collapse=separator)),
Best_match_references=paste(
best_match_references,collapse=separator),
Best_match_E_value=min_E_value,
Best_match_bit_score=bit_score,
Sequence_name=sequence_to_classify,
stringsAsFactors=FALSE)
} else {
# If include missing is not enabled.
sequence_information<-data.frame(Local_taxa="root",
Local_species=ifelse(full_names,
paste(sort(local$Name),collapse=separator),
paste(local$Species[order(local$Name)],
collapse=separator)),
Best_match_references=paste(
best_match_references,collapse=separator),
Best_match_E_value=min_E_value,
Best_match_bit_score=bit_score,
Sequence_name=sequence_to_classify,
stringsAsFactors=FALSE)
}
}
}
}
# Add sequence information to the storage data frame.
sequence_classifications<-rbind(sequence_classifications,sequence_information)
}
# Merge taxa suggestions with means and SDs of query cover and PID.
sequence_classifications<-merge(sequence_classifications,best_match_means_and_SDs)
# Add sequences to the classification data frame.
df<-merge(sequences_to_classify,sequence_classifications,all=TRUE)
# Keep certain fields in the classification data frame.
if(include_missing){
# If include missing is enabled.
df<-df[,c("Sequence_name","Sequence","Best_match_references","Best_match_E_value",
"Best_match_bit_score","Best_match_query_cover.mean",
"Best_match_query_cover.SD","Best_match_PID.mean","Best_match_PID.SD",
"Local_taxa","Local_species","Local_taxa.include_missing",
"Local_species.include_missing")]
} else {
# If include missing is not enabled.
df<-df[,c("Sequence_name","Sequence","Best_match_references","Best_match_E_value",
"Best_match_bit_score","Best_match_query_cover.mean",
"Best_match_query_cover.SD","Best_match_PID.mean","Best_match_PID.SD",
"Local_taxa","Local_species")]
}
} else { # If a local taxa list is not provided.
# Get unique best-matches for each sequence.
blast_output<-unique(blast_output)
# Order the records in the best-matches by sequence name and taxon name.
blast_output<-blast_output[order(blast_output$Sequence_name,blast_output$Matches),]
# Get the species names for the best-matching reference sequences.
if(full_names){ # If full names are requested, get the full taxonomies.
blast_output$Species_names<-blast_output$Matches
} else { # If full names are not requested, get the species binomials.
blast_output$Species_names<-sapply(X=strsplit(x=blast_output$Matches,split=";"),
FUN="[[",7)
}
# Get best-matching references.
best_matches<-stats::aggregate(Species_names~Sequence_name,data=blast_output,
FUN=paste,collapse=separator)
# Get the E-values of best-matches.
best_matches_E_value<-stats::aggregate(E_value~Sequence_name,data=blast_output,
FUN="[",1)
# Get the bit-scores of best-matches.
best_matches_Bit_score<-stats::aggregate(Bit_score~Sequence_name,data=blast_output,
FUN="[",1)
# Combine information on best-matching references.
best_matches<-merge(best_matches,best_matches_E_value)
best_matches<-merge(best_matches,best_matches_Bit_score)
best_matches<-merge(best_matches,best_match_means_and_SDs)
# Add sequences to the classification data frame.
df<-merge(sequences_to_classify,best_matches,all=TRUE)
# Rename certain fields.
colnames(df)[which(colnames(df)=="Species_names")]<-"Best_match_references"
colnames(df)[which(colnames(df)=="E_value")]<-"Best_match_E_value"
colnames(df)[which(colnames(df)=="Bit_score")]<-"Best_match_bit_score"
# Keep certain fields in the classification data frame.
df<-df[,c("Sequence_name","Sequence","Best_match_references","Best_match_E_value",
"Best_match_bit_score","Best_match_query_cover.mean",
"Best_match_query_cover.SD","Best_match_PID.mean","Best_match_PID.SD")]
}
# Set NA values in rows without an E-value to suggest that no significant similarity
# was found.
df[is.na(df$Best_match_E_value),3:ncol(df)]<-"No significant similarity found"
# Write out classified sequences.
utils::write.csv(x=df,file=path_to_output_file,row.names=FALSE)
}
Try the LocaTT package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.