Package 'MetaScope'

Description

This object is a named vector of multiple options that can be chosen for functions that involve alignment with Rsubread, namely align_target() and filter_host(). Both functions take an object for the parameter settings, which are provided by align_details by default, or may be given by a user-created object containing the same information.

Usage

data(align_details)

Format

list

Details

The default options included in align_details are type = "dna", maxMismatches = 3, nsubreads = 10, phredOffset = 33, unique = FALSE, and nBestLocations = 16. Full descriptions of these parameters can be read by accessing ?Rsubread::align.

Examples

data("align_details")

Description

This is the main MetaScope target library mapping function, using Rsubread and multiple libraries. Aligns to each library separately, filters unmapped reads from each file, and then merges and sorts the .bam files from each library into one output file. If desired, output can be passed to 'filter_host()' to remove reads that also map to filter library genomes.

Usage

align_target(
  read1,
  read2 = NULL,
  lib_dir = NULL,
  libs,
  threads = 1,
  align_file = tools::file_path_sans_ext(read1),
  subread_options = align_details,
  quiet = TRUE
)

Arguments

Value

This function writes a merged and sorted .bam file after aligning to all reference libraries given, along with a summary report file, to the user's working directory. The function also outputs the new .bam filename.

Examples

#### Align example reads to an example reference library using Rsubread

## Create temporary directory
target_ref_temp <- tempfile()
dir.create(target_ref_temp)

tax <- "Ovine atadenovirus D"

## Create temporary taxonomizr accession
tmp_accession <- system.file("extdata", "example_accessions.sql", package = "MetaScope")

## Download genome
all_ref <- MetaScope::download_refseq(tax,
                                      reference = FALSE,
                                      representative = FALSE,
                                      compress = TRUE,
                                      out_dir = target_ref_temp,
                                      caching = TRUE,
                                      accession_path = tmp_accession)

## Create subread index
ind_out <- mk_subread_index(all_ref)

## Get path to example reads
readPath <- system.file("extdata", "reads.fastq",
                        package = "MetaScope")
## Copy the example reads to the temp directory
refPath <- file.path(target_ref_temp, "reads.fastq")
file.copy(from = readPath, to = refPath)

## Modify alignment parameters object
data("align_details")
align_details[["type"]] <- "rna"
align_details[["phredOffset"]] <- 50
# Just to get it to align - toy example!
align_details[["maxMismatches"]] <- 100

## Run alignment
target_map <- align_target(refPath,
                           libs = stringr::str_replace_all(tax, " ", "_"),
                           lib_dir = target_ref_temp,
                           subread_options = align_details)

## Remove temporary folder
unlink(target_ref_temp, recursive = TRUE)

Description

This is the main MetaScope target library mapping function, using Rbowtie2 and multiple libraries. Aligns to each library separately, filters unmapped reads from each file, and then merges and sorts the .bam files from each library into one output file. If desired, output can be passed to 'filter_host_bowtie()' to remove reads that also map to filter library genomes.

Usage

align_target_bowtie(
  read1,
  read2 = NULL,
  lib_dir,
  libs,
  align_dir,
  align_file,
  bowtie2_options = NULL,
  threads = 1,
  overwrite = FALSE,
  quiet = TRUE
)

Arguments

Details

The default parameters are the same that PathoScope 2.0 uses. "–very-sensitive-local -k 100 –score-min L,20,1.0"

If you experience any issues with reading the input files, make sure that the file(s) are not located in a read-only folder. This can be circumvented by copying files to a new location before running the function.

Value

Returns the path to where the output alignment file is stored.

Examples

#### Align example reads to an example reference library using Rbowtie2

## Create temporary directory to store file
target_ref_temp <- tempfile()
dir.create(target_ref_temp)

tmp_accession <- system.file("extdata", "example_accessions.sql", package = "MetaScope")

## Dowload reference genome
MetaScope::download_refseq("Morbillivirus hominis",
                           reference = FALSE,
                           representative = FALSE,
                           compress = TRUE,
                           out_dir = target_ref_temp,
                           caching = TRUE,
                           accession_path = tmp_accession
)

## Create temporary directory to store the indices
index_temp <- tempfile()
dir.create(index_temp)

## Create bowtie2 index
MetaScope::mk_bowtie_index(
  ref_dir = target_ref_temp,
  lib_dir = index_temp,
  lib_name = "target",
  overwrite = TRUE
)

## Create temporary directory for final file
output_temp <- tempfile()
dir.create(output_temp)

## Get path to example reads
readPath <- system.file("extdata", "virus_example.fastq",
                        package = "MetaScope")

## Align to target genomes
target_map <-
  MetaScope::align_target_bowtie(
    read1 = readPath,
    lib_dir = index_temp,
    libs = "target",
    align_dir = output_temp,
    align_file = "bowtie_target",
    overwrite = TRUE,
    bowtie2_options = "--very-sensitive-local"
  )

## Remove extra folders
unlink(target_ref_temp, recursive = TRUE)
unlink(index_temp, recursive = TRUE)
unlink(output_temp, recursive = TRUE)

Description

This function replaces the header from one .bam file with a header from a different .sam file. This function mimics the function of the 'reheader' function in samtools. It is not intended for use by users.

Usage

bam_reheader_R(
  head,
  old_bam,
  new_bam = paste(tools::file_path_sans_ext(old_bam), "h.bam", sep = "")
)

Arguments

Value

This function will return a new .bam file with a replaced header. The function also outputs the new .bam filename.

Description

Using the output from metascope_blast(), the blast_reassignment() function takes the results and alters the original metascope_id() output to reassign reads that were invalidated by the BLAST findings. Currently, the implementation of this function only reassigns reads to a taxon that was already found in the sample at a higher abundance.

Usage

blast_reassignment(
  metascope_blast_path,
  species_threshold,
  num_hits,
  blast_tmp_dir,
  out_dir,
  sample_name,
  reassign_validated = FALSE,
  reassign_to_validated = TRUE
)

Arguments

Value

Returns a data.frame with the reassigned taxa and read counts.

Description

This function calculates the best hit (genome with most blast read hits), uniqueness score (total number of genomes hit), species percentage hit (percentage of reads where MetaScope species also matched the blast hit species), genus percentage hit (percentage of reads where blast genus matched MetaScope aligned genus) and species contaminant score (percentage of reads that blasted to other species genomes) and genus contaminant score (percentage of reads that blasted to other genus genomes)

Usage

blast_result_metrics(blast_results_table_path, accession_path, db = NULL)

Arguments

Value

a vector with best_hit, uniqueness_score, species_percentage_hit genus_percentage_hit, species_contaminant_score, and genus_contaminant_score

Description

Reformat BLASTn results

Usage

blastn_results(
  results_table,
  bam_file,
  num_results = 10,
  num_reads_per_result = 100,
  hit_list = 10,
  num_threads = 1,
  db_path,
  out_path,
  db = NULL,
  sample_name = NULL,
  quiet = quiet,
  accession_path,
  fasta_dir = NULL,
  BPPARAM
)

Arguments

Value

Creates and exports num_results number of csv files with blast results from local blast

Description

blastn_single_result

Usage

blastn_single_result(
  results_table,
  bam_file,
  which_result,
  num_reads = 100,
  hit_list = 10,
  num_threads,
  db_path,
  quiet,
  accession_path,
  bam_seqs,
  out_path,
  sample_name,
  fasta_dir = NULL
)

Arguments

Value

Returns a dataframe of blast results for a metascope result

Description

This character string provides several Bowtie 2 options to provide an optimized alignment specifically optimized for 16S amplicon sequencing data. This object can be used with functions that use the Bowtie 2 aligner through the Rbowtie2 package, namely align_target_bowtie() and filter_host_bowtie. These settings can be substituted for default settings by passing to the bowtie2_options argument.

Usage

data(bt2_16S_params)

Format

list

Details

The default parameters listed in this object are "–local -R 2 -N 0 -L 25 -i S,1,0.75 -k 5 –score-min L,0,1.88"

Note that k is actually 10 and is doubled internally from 5. The score-min function was chosen such that the minimum alignment score allowed requires 98

Further delineation of Bowtie 2 parameters is provided in the Bowtie 2 manual.

Examples

data("bt2_16S_params")

Description

This character string provides several Bowtie 2 options to conduct an alignment useful for metagenomes, especially in the case where a genome may not be present in the reference database. This object can be used with functions that use the Bowtie 2 aligner through the Rbowtie2 package, namely align_target_bowtie() and filter_host_bowtie. These settings can be substituted for default settings by passing to the bowtie2_options argument.

Usage

data(bt2_missing_params)

Format

list

Details

The default parameters listed in this object are "–local -R 2 -N 0 -L 25 -i S,1,0.75 -k 5 –score-min L,0,1.4".

Further delineation of Bowtie 2 parameters is provided in the Bowtie 2 manual.

Examples

data("bt2_missing_params")

Description

This character string provides several Bowtie 2 options to provide a 95 alignment useful for metagenomes. This object can be used with functions that use the Bowtie 2 aligner through the Rbowtie2 package, namely align_target_bowtie() and filter_host_bowtie. These settings can be substituted for default settings by passing to the bowtie2_options argument.

Usage

data(bt2_regular_params)

Format

list

Details

The default parameters listed in this object are "–local -R 2 -N 0 -L 25 -i S,1,0.75 -k 5 –score-min L,0,1.7".

Further delineation of Bowtie 2 parameters is provided in the Bowtie 2 manual.

Examples

data("bt2_regular_params")

Description

This is an internal function that is not meant to be used outside of the package. It checks whether blastn exists on the system.

Usage

check_blastn_exists()

Details

Checks if blastn is installed

Value

Returns TRUE if blastn exists on the system, else FALSE.

Description

This is an internal function that is not meant to be used outside of the package. It checks whether samtools exists on the system.

Usage

check_samtools_exists()

Value

Returns TRUE if samtools exists on the system, else FALSE.

Description

This function generates a combined header from multiple .bam files from different reference libraries (e.g. a split bacterial library). It is not intended for use by users.

Usage

combined_header(bam_files, header_file = "header_tmp.sam")

Arguments

Value

This function will return a combined header from all the supplied .bam files.

Description

Upon completion of the MetaScope pipeline, users can analyze and visualize abundances in their samples using the animalcules package. This function allows interoperability of metascope_id output with both animalcules and QIIME. After running this function, the user should save the returned MAE to an RDS file using a function like saveRDS to upload the output into the animalcules package.

Usage

convert_animalcules(
  meta_counts,
  annot_path,
  which_annot_col,
  end_string = ".metascope_id.csv",
  qiime_biom_out = FALSE,
  path_to_write = ".",
  accession_path = NULL
)

Arguments

Value

Returns a MultiAssay Experiment file of combined sample counts data and/or biom file and mapping file for analysis with QIIME. The MultiAssay Experiment will have a counts assay ("MGX").

Examples

tempfolder <- tempfile()
dir.create(tempfolder)

# Create three different samples
samp_names <- c("X123", "X456", "X789")
all_files <- file.path(tempfolder,
                       paste0(samp_names, ".csv"))

create_IDcsv <- function (out_file) {
  final_taxids <- c("273036", "418127", "11234")
  final_genomes <- c(
    "Staphylococcus aureus RF122, complete sequence",
    "Staphylococcus aureus subsp. aureus Mu3, complete sequence",
    "Measles virus, complete genome")
  best_hit <- sample(seq(100, 1050), 3)
  proportion <- best_hit/sum(best_hit) |> round(2)
  EMreads <- best_hit + round(runif(3), 1)
  EMprop <- proportion + 0.003
  dplyr::tibble(TaxonomyID = final_taxids,
                Genome = final_genomes,
                read_count = best_hit, Proportion = proportion,
                EMreads = EMreads, EMProportion = EMprop) |>
    dplyr::arrange(dplyr::desc(.data$read_count)) |>
    utils::write.csv(file = out_file, row.names = FALSE)
  message("Done!")
  return(out_file)
}
out_files <- vapply(all_files, create_IDcsv, FUN.VALUE = character(1))

# Create annotation data for samples
annot_dat <- file.path(tempfolder, "annot.csv")
dplyr::tibble(Sample = samp_names, RSV = c("pos", "neg", "pos"),
              month = c("March", "July", "Aug"),
              yrsold = c(0.5, 0.6, 0.2)) |>
  utils::write.csv(file = annot_dat,
                   row.names = FALSE)

# Create temporary taxonomizr accession
tmp_accession <- system.file("extdata", "example_accessions.sql", package = "MetaScope")

# Convert samples to MAE
outMAE <- convert_animalcules(meta_counts = out_files,
                              annot_path = annot_dat,
                              which_annot_col = "Sample",
                              end_string = ".metascope_id.csv",
                              qiime_biom_out = FALSE,
                              accession_path = tmp_accession)

unlink(tempfolder, recursive = TRUE)

Description

This function serves as a legacy integration method for usage with PathoScope 2.0 outputs. Upon completion of the PathoScope 2.0 pipeline, users can analyze and visualize abundances in their samples using the animalcules package. After running this function, the user should save the returned MAE to an RDS file using a function like saveRDS to upload the output into the animalcules package.

Usage

convert_animalcules_patho(
  patho_counts,
  annot_path,
  which_annot_col,
  end_string = "-sam-report.tsv"
)

Arguments

Value

Returns a MultiAssay Experiment file of combined sample counts data. The MultiAssay Experiment will have a counts assay ("MGX").

Description

Upon completion of the MetaScope pipeline, users can analyze and visualize abundances in their samples using the animalcules package. This function allows interoperability of metascope_id output with both animalcules and QIIME. After running this function, the user should save the returned MAE to an RDS file using a function like saveRDS to upload the output into the animalcules package. NOTE: This function is for outputs that were generated with the SILVA 13_8 database.

Usage

convert_animalcules_silva(
  meta_counts,
  annot_path,
  which_annot_col,
  end_string = ".metascope_id.csv",
  qiime_biom_out = FALSE,
  path_to_write = ".",
  caching = TRUE
)

Arguments

Value

Returns a MultiAssay Experiment file of combined sample counts data and/or saved biom file and mapping file for analysis with QIIME. The MultiAssayExperiment will have a counts assay ("MGX").

Examples

tempfolder <- tempfile()
dir.create(tempfolder)

# Create three different samples
samp_names <- c("X123", "X456", "X789")
all_files <- file.path(tempfolder,
                       paste0(samp_names, ".csv"))

create_IDcsv <- function (out_file) {
  final_taxids <- c("AY846380.1.2583", "AY909584.1.2313", "HG531388.1.1375")
  final_genomes <- rep("Genome name", 3)
  best_hit <- sample(seq(100, 1050), 3)
  proportion <- best_hit/sum(best_hit) |> round(2)
  EMreads <- best_hit + round(runif(3), 1)
  EMprop <- proportion + 0.003
  dplyr::tibble("TaxonomyID" = final_taxids,
                "Genome" = final_genomes,
                "read_count" = best_hit, "Proportion" = proportion,
                "EMreads" = EMreads, "EMProportion" = EMprop) |>
    dplyr::arrange(dplyr::desc(.data$read_count)) |>
    utils::write.csv(file = out_file, row.names = FALSE)
  message("Done!")
  return(out_file)
}
out_files <- vapply(all_files, create_IDcsv, FUN.VALUE = character(1))

# Create annotation data for samples
annot_dat <- file.path(tempfolder, "annot.csv")
dplyr::tibble(Sample = samp_names, RSV = c("pos", "neg", "pos"),
              month = c("March", "July", "Aug"),
              yrsold = c(0.5, 0.6, 0.2)) |>
  utils::write.csv(file = annot_dat,
                   row.names = FALSE)

# Convert samples to MAE
outMAE <- convert_animalcules_silva(meta_counts = out_files,
                                    annot_path = annot_dat,
                                    which_annot_col = "Sample",
                                    end_string = ".metascope_id.csv",
                                    qiime_biom_out = FALSE,
                                    caching = TRUE)

unlink(tempfolder, recursive = TRUE)

Description

The 'CIGAR' (Compact Idiosyncratic Gapped Alignment Report) string is how the SAM/BAM format represents spliced alignments. This function will accept a CIGAR string for a single read and a single character indicating the operation to be parsed in the string. An operation is a type of column that appears in the alignment, e.g. a match or gap. The integer following the operator specifies a number of consecutive operations. The count_matches() function will identify all occurrences of the operator in the string input, add them, and return an integer number representing the total number of operations for the read that was summarized by the input CIGAR string.

Usage

count_matches(x, char = "M")

Arguments

Details

This function is best used on a vector of CIGAR strings using an apply function (see examples).

Value

an integer number representing the total number of alignment operations for the read that was summarized by the input CIGAR string.

Examples

# A single cigar string: 3M + 3M + 5M
cigar1 <- "3M1I3M1D5M"
count_matches(cigar1, char = "M")

# Parse with operator "P": 2P
cigar2 <- "4M1I2P9M"
count_matches(cigar2, char = "P")

# Apply to multiple strings: 1I + 1I + 5I
cigar3 <- c("3M1I3M1D5M", "4M1I1P9M", "76M13M5I")
vapply(cigar3, count_matches, char = "I",
       FUN.VALUE = numeric(1))

Description

This is a necessary step for all samples utilizing NCBI and SILVA databases in the MetaScope pipeline. As specified by the user, download_accessions will automatically download the NCBI accessions database, the SILVA taxonomy database, and or the NCBI Blast 16S database and prepare consolidated databases for downstream use with the MetaID and MetaBLAST modules. This package relies on the taxonomizr package.

Usage

download_accessions(
  ind_dir,
  tmp_dir = file.path(ind_dir, "tmp"),
  remove_tmp_dir = TRUE,
  NCBI_accessions_database = TRUE,
  NCBI_accessions_name = "accessionTaxa",
  silva_taxonomy_database = TRUE,
  silva_taxonomy_name = "all_silva_headers",
  blast_16S_database = TRUE,
  blast_16S_name = "16S_ribosomal_RNA"
)

Arguments

Value

Exports database(s) with names and to location specified by the user.

Examples

## Not run: 
  download_accessions(
    ind_dir = "C:/Users/JohnSmith/Research",
    tmp_dir = file.path(ind_dir, "tmp"),
    remove_tmp_dir = TRUE,
    NCBI_accessions_database = TRUE,
    NCBI_accessions_name = "accessionTaxa.sql",
    silva_taxonomy_database = TRUE,
    silva_taxonomy_name = "all_silva_headers.rds")

## End(Not run)

Description

This function will automatically download RefSeq genome libraries in a fasta format from the specified taxon. The function will first download the summary report at: ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/**kingdom**/assembly_summary.txt, and then use this file to download the genome(s) and combine them in a single compressed or uncompressed .fasta file.

Usage

download_refseq(
  taxon,
  reference = TRUE,
  representative = FALSE,
  compress = TRUE,
  patho_out = FALSE,
  out_dir = NULL,
  caching = FALSE,
  quiet = TRUE,
  accession_path = NULL
)

Arguments

Details

When selecting the taxon to be downloaded, if you receive an error saying Your input is not a valid taxon, please take a look at the taxonomy_table object, which can be accessed with the command MetaScope:::taxonomy_table). Only taxa with exact spelling as they appear at any level of the table will be acknowledged.

Value

Returns a .fasta or .fasta.gz file of the desired RefSeq genomes. This file is named after the kingdom selected and saved to the current directory (e.g. 'bacteria.fasta.gz'). This function also has the option to return a .fasta file formatted for PathoScope as well (e.g. bacteria.pathoscope.fasta.gz') if path_out = TRUE.

Examples

#### Download RefSeq genomes

## Create temporary taxonomizr accession
tmp_accession <- system.file("extdata", "example_accessions.sql", package = "MetaScope")

## Download all RefSeq reference Bovismacovirus genus genomes
download_refseq('Bovismacovirus', reference = FALSE, representative = FALSE,
                out_dir = NULL, compress = TRUE, patho_out = FALSE,
                caching = TRUE, accession_path = tmp_accession)

Description

This function will automatically download the 16S rRNA RefSeq libraries from https://ftp.ncbi.nlm.nih.gov/refseq/TargetedLoci/ and combine them into a single .fna file

Usage

download_refseq_16S(out_dir, combined_name = "refseq_16S.fna")

Arguments

Value

Returns a character string of the path to the combined 16S rRNA .fna file

Note

This function requires the suggested packages RCurl and R.utils. If they are not installed, you will need to install them manually.

Examples

## Not run: 
#### Download 16S rRNA Genomes

download_refseq_16S(out_dir = "out_dir",
  combined_name = "refseq_16S.fna")
  
## End(Not run)

Description

Helper function for demultiplexing sequencing reads, designed in a way to allow for parallelization across barcodes (parallel extraction of reads by barcode). This function takes a specific barcode (numeric index) from lists of sample names/barcodes, a Biostrings::DNAStringSet of barcodes by sequence header, and a Biostrings::QualityScaledXStringSet of reads corresponding to the barcodes. Based on the barcode index given, it extracts all reads for the indexed barcode and writes all the reads from that barcode to a separate .fastq file.

Usage

extract_reads(
  barcodeIndex,
  barcodes,
  sampleNames,
  index,
  reads,
  location = "./demultiplex_fastq",
  rcBarcodes = TRUE,
  hDist = 0,
  quiet = TRUE
)

Arguments

Value

Writes a single .fastq file that contains all reads whose index matches the barcode specified. This file will be written to the location directory, and will be named based on the specified sampleName and barcode, e.g. './demultiplex_fastq/SampleName1_GGAATTATCGGT.fastq.gz'

Examples

## Create temporary directory
ref_temp <- tempfile()
dir.create(ref_temp)

## Load example barcode, index, and read data into R session
barcodePath <- system.file("extdata", "barcodes.txt", package = "MetaScope")
bcFile <- read.table(barcodePath, sep = "\t", header = TRUE)

indexPath <- system.file("extdata", "virus_example_index.fastq",
package = "MetaScope")
inds <- Biostrings::readDNAStringSet(indexPath, format = "fastq")

readPath <- system.file("extdata", "virus_example.fastq",
                        package = "MetaScope")
reads <- Biostrings::readQualityScaledDNAStringSet(readPath)

## Extract reads from the first barcode
results <- extract_reads(1, bcFile[, 2], bcFile[, 1], inds, reads,
                        rcBarcodes = FALSE, location = ref_temp)

## Extract reads from multiple barcodes
more_results <- lapply(1:6, extract_reads, bcFile[, 2], bcFile[, 1], inds,
                       reads, rcBarcodes = FALSE, location = ref_temp)

## Remove temporary directory
unlink(ref_temp, recursive = TRUE)

Description

After aligning your sample to a target library with align_target(), use filter_host() to remove unwelcome host contamination using filter reference libraries. This function takes as input the name of the .bam file produced via align_target(), and produces a sorted .bam file with any reads that match the filter libraries removed. This resulting .bam file may be used upstream for further analysis. This function uses Rsubread. For the Rbowtie2 equivalent of this function, see filter_host_bowtie.

Usage

filter_host(
  reads_bam,
  lib_dir = NULL,
  libs,
  make_bam = FALSE,
  output = paste(tools::file_path_sans_ext(reads_bam), "filtered", sep = "."),
  subread_options = align_details,
  YS = 1e+05,
  threads = 1,
  quiet = TRUE
)

Arguments

Details

A compressed .csv can be created to produce a smaller output file that is created more efficiently and is still compatible with metascope_id().

Value

The name of a filtered, sorted .bam file written to the user's current working directory. Or, if make_bam = FALSE, a .csv.gz file containing a data frame of only requisite information to run metascope_id().

Examples

#### Filter reads from bam file that align to any of the filter libraries

## Assuming a bam file has been created previously with align_target()

## Create temporary directory
filter_ref_temp <- tempfile()
dir.create(filter_ref_temp)

## Create temporary taxonomizr accession
tmp_accession <- system.file("extdata", "example_accessions.sql", package = "MetaScope")


## Download filter genome
all_species <- c("Staphylococcus aureus subsp. aureus str. Newman")
all_ref <- vapply(all_species, MetaScope::download_refseq,
                  reference = FALSE, representative = FALSE, compress = TRUE,
                  out_dir = filter_ref_temp, caching = FALSE,
                  accession_path = tmp_accession,
                  FUN.VALUE = character(1))
ind_out <- vapply(all_ref, mk_subread_index, FUN.VALUE = character(1))

## Get path to example reads
readPath <- system.file("extdata", "subread_target.bam",
                        package = "MetaScope")

## Copy the example reads to the temp directory
refPath <- file.path(filter_ref_temp, "subread_target.bam")
file.copy(from = readPath, to = refPath)
utils::data("align_details")
align_details[["type"]] <- "rna"
align_details[["phredOffset"]] <- 10
# Just to get it to align - toy example!
align_details[["maxMismatches"]] <- 10

## Align and filter reads
filtered_map <- filter_host(
  refPath, lib_dir = filter_ref_temp,
  libs = stringr::str_replace_all(all_species, " ", "_"),
  threads = 1, subread_options = align_details)

## Remove temporary directory
unlink(filter_ref_temp, recursive = TRUE)

Description

After a sample is aligned to a target library with align_target_bowtie(), we may use filter_host_bowtie() to remove unwelcome host contamination using filter reference libraries. This function takes as input the name of the .bam file produced via align_target_bowtie(), and produces a sorted .bam or .csv.gz file with any reads that match the filter libraries removed. This resulting .bam file may be used downstream for further analysis. This function uses Rbowtie2 For the Rsubread equivalent of this function, see filter_host.

Usage

filter_host_bowtie(
  reads_bam,
  lib_dir,
  libs,
  make_bam = FALSE,
  output = paste(tools::file_path_sans_ext(reads_bam), "filtered", sep = "."),
  bowtie2_options = NULL,
  YS = 1e+05,
  threads = 1,
  overwrite = FALSE,
  quiet = TRUE
)

Arguments

Details

A compressed .csv can be created to produce a smaller output file that is created more efficiently and is still compatible with metascope_id().

The default parameters are the same that PathoScope 2.0 uses. "–very-sensitive-local -k 100 –score-min L,20,1.0"

Value

The name of a filtered, sorted .bam file written to the user's current working directory. Or, if make_bam = FALSE, a .csv.gz file containing a data frame of only requisite information to run metascope_id().

Examples

#### Filter reads from bam file that align to any of the filter libraries

## Assuming a bam file has already been created with align_target_bowtie()
# Create temporary filter library
filter_ref_temp <- tempfile()
dir.create(filter_ref_temp)

## Create temporary taxonomizr accession
tmp_accession <- system.file("extdata", "example_accessions.sql", package = "MetaScope")

## Download reference genome
MetaScope::download_refseq("Orthoebolavirus zairense",
                           reference = FALSE,
                           representative = FALSE,
                           compress = TRUE,
                           out_dir = filter_ref_temp,
                           caching = TRUE,
                           accession_path = tmp_accession)

## Create temp directory to store the indices
index_temp <- tempfile()
dir.create(index_temp)

## Create filter index
MetaScope::mk_bowtie_index(
  ref_dir = filter_ref_temp,
  lib_dir = index_temp,
  lib_name = "filter",
  overwrite = TRUE
)

## Create temporary folder to hold final output file
output_temp <- tempfile()
dir.create(output_temp)

## Get path to example bam
bamPath <- system.file("extdata", "bowtie_target.bam",
                       package = "MetaScope")
target_copied <- file.path(output_temp, "bowtie_target.bam")
file.copy(bamPath, target_copied)

## Align and filter reads
filter_out <-
  filter_host_bowtie(
    reads_bam = target_copied,
    lib_dir = index_temp,
    libs = "filter",
    threads = 1
  )

## Remove temporary directories
unlink(filter_ref_temp, recursive = TRUE)
unlink(index_temp, recursive = TRUE)
unlink(output_temp, recursive = TRUE)

Description

This function will remove all unmapped reads or lines in a .bam file (warning: overwrites the original file!). This function is needed because combining multiple .bam files from different microbial libraries may lead to some reads that mapped to one library and have unmapped entries from another library. This will remove any unmapped entries and leave all reference mapped lines in the .bam file.

Usage

filter_unmapped_reads(bamfile)

Arguments

Details

It is not intended for direct use.

Value

This function will overwrite the existing .bam file with a new .bam file in the same location that has only mapped lines. The function itself returns the output .bam file name.

Description

This function will utilize an organism classification table to obtain all children species and/or strains with available NCBI reference sequences given a parent taxon and its rank.

Usage

get_children(input_taxon, input_rank, tax_dat = NULL)

Arguments

Value

Returns a vector of all the child species and/or strains of the input taxon.

Examples

## Get all child species and strains in bacteria superkingdom
get_children('Bacteria','superkingdom')

## Get all child species and strains in fungi kingdom
get_children('Fungi', 'kingdom')

## Get all child species in primate order
get_children('Primates', 'order')

Description

Returns fasta sequences from a bam file with a given taxonomy ID

Usage

get_seqs(id, bam_file, n = 10, bam_seqs)

Arguments

Value

Biostrings format sequences

Description

Used to create plots of genome coverage for any number of accession numbers

Usage

locations(
  which_taxid,
  which_genome,
  accessions,
  taxids,
  reads,
  out_base,
  out_dir
)

Arguments

Value

A plot of the read coverage for a given genome

Description

This function merges .bam files. It first used the combined_header function to generate a combined header for all the files, reheaders the files, and then merges and sorts the .bam files. It is similar to the 'samtools merge' function, but it allows the .bam files to have different headers. It is not intended for direct use.

Usage

merge_bam_files(
  bam_files,
  destination,
  head_file = paste(destination, "_header.sam", sep = ""),
  quiet = TRUE
)

Arguments

Value

This function merges .bam files and combines them into a single file. The function also outputs the new .bam filename.

Description

Function for demultiplexing sequencing reads arranged in a common format provided by sequencers (such as Illumina) generally for 16S data. This function takes a matrix of sample names/barcodes, a .fastq file of barcodes by sequence header, and a .fastq file of reads corresponding to the barcodes. Based on the barcodes given, the function extracts all reads for the indexed barcode and writes all the reads from that barcode to separate .fastq files.

Usage

meta_demultiplex(
  barcodeFile,
  indexFile,
  readFile,
  rcBarcodes = TRUE,
  location = NULL,
  threads = 1,
  hammingDist = 0,
  quiet = TRUE
)

Arguments

Value

Returns multiple .fastq files that contain all reads whose index matches the barcodes given. These files will be written to the location directory, and will be named based on the given sampleNames and barcodes, e.g. './demultiplex_fastq/SampleName1_GGAATTATCGGT.fastq.gz'

Examples

## Get barcode, index, and read data locations
barcodePath <- system.file("extdata", "barcodes.txt", package = "MetaScope")
indexPath <- system.file("extdata", "virus_example_index.fastq",
                         package = "MetaScope")
readPath <- system.file("extdata", "virus_example.fastq",
                         package = "MetaScope")

## Demultiplex
demult <- meta_demultiplex(barcodePath, indexPath, readPath, rcBarcodes = FALSE,
                      hammingDist = 2)
demult

Description

This function allows the user to check a subset of identified reads against NCBI BLAST and the nucleotide database to confirm or contradict results provided by MetaScope. It aligns the top 'metascope_id()' results to NCBI BLAST database. It REQUIRES that command-line BLAST and a separate nucleotide database have already been installed on the host machine. It returns a csv file updated with BLAST result metrics.

Usage

metascope_blast(
  metascope_id_path,
  bam_file_path = list.files(tmp_dir, ".updated.bam$", full.names = TRUE)[1],
  tmp_dir,
  out_dir,
  sample_name,
  fasta_dir = NULL,
  num_results = 10,
  num_reads = 100,
  hit_list = 10,
  num_threads = 1,
  db_path,
  quiet = FALSE,
  db = NULL,
  accession_path = NULL
)

Arguments

Details

This function assumes that you used the NCBI nucleotide database to process samples, with a download date of 2021 or later. This is necessary for compatibility with the bam file headers.

This is highly computationally intensive and should be ran with multiple cores, submitted as a multi-threaded computing job if possible.

Note, if best_hit_strain is FALSE, then no strain was observed more often among the BLAST results.

Value

This function writes an updated csv file with metrics.

Examples

## Not run: 
### Create temporary directory
file_temp <- tempfile()
dir.create(file_temp)

bamPath <- system.file("extdata", "bowtie_target.bam",
                       package = "MetaScope")
file.copy(bamPath, file_temp)

metascope_id(file.path(file_temp, "bowtie_target.bam"), aligner = "bowtie2",
             input_type = "bam", out_dir = file_temp, num_species_plot = 0,
             update_bam = TRUE)

## Run metascope blast
### Get export name and metascope id results
out_base <- bamPath |> basename() |> tools::file_path_sans_ext() |>
  tools::file_path_sans_ext()
metascope_id_path <- file.path(file_temp, paste0(out_base, ".metascope_id.csv"))

# NOTE: change db_path to the location where your BLAST database is stored!
db <- "/restricted/projectnb/pathoscope/data/blastdb/nt/nt"

tmp_accession <- system.file("extdata", "example_accessions.sql", package = "MetaScope"

metascope_blast(metascope_id_path,
                bam_file_path = file.path(file_temp, "bowtie_target.bam"),
                tmp_dir = file_temp,
                out_dir = file_temp,
                sample_name = out_base,
                db_path = db,
                num_results = 10,
                num_reads = 5,
                hit_list = 10,
                num_threads = 3,
                db = "ncbi",
                quiet = FALSE,
                fasta_dir = NULL,
                accession_path = tmp_accession)

## Remove temporary directory
unlink(file_temp, recursive = TRUE)

## End(Not run)

Description

This function will read in a .bam or .csv.gz file, annotate the taxonomy and genome names, reduce the mapping ambiguity using a mixture model, and output a .csv file with the results. Currently, it assumes that the genome library/.bam files use NCBI accession names for reference names (rnames in .bam file).

Usage

metascope_id(
  input_file,
  input_type = "csv.gz",
  aligner = "bowtie2",
  db = "ncbi",
  db_feature_table = NULL,
  accession_path = NULL,
  priors_df = NULL,
  tmp_dir = dirname(input_file),
  out_dir = dirname(input_file),
  convEM = 1/10000,
  maxitsEM = 25,
  update_bam = FALSE,
  num_species_plot = NULL,
  group_by_taxa = "species",
  quiet = TRUE
)

Arguments

Value

This function exports a .csv file with annotated read counts to genomes with mapped reads to the location returned by the function. Depending on the parameters specified, can also output an updated BAM file, and fasta files for additional analysis downstream.

Examples

#### Align reads to reference library and then apply metascope_id()
## Assuming filtered bam files already exist

## Create temporary directory
file_temp <- tempfile()
dir.create(file_temp)

## Get temporary accessions database
tmp_accession <- system.file("extdata", "example_accessions.sql", package = "MetaScope")

#### Subread aligned bam file

## Create object with path to filtered subread csv.gz file
filt_file <- "subread_target.filtered.csv.gz"
bamPath <- system.file("extdata", filt_file, package = "MetaScope")
file.copy(bamPath, file_temp)

## Run metascope id with the aligner option set to subread
metascope_id(input_file = file.path(file_temp, filt_file),
             aligner = "subread", num_species_plot = 0,
             input_type = "csv.gz", accession_path = tmp_accession)

#### Bowtie 2 aligned .csv.gz file

## Create object with path to filtered bowtie2 bam file
bowtie_file <- "bowtie_target.filtered.csv.gz"
bamPath <- system.file("extdata", bowtie_file, package = "MetaScope")
file.copy(bamPath, file_temp)

## Run metascope id with the aligner option set to bowtie2
metascope_id(file.path(file_temp, bowtie_file), aligner = "bowtie2",
             num_species_plot = 0, input_type = "csv.gz",
             accession_path = tmp_accession)

## Remove temporary directory
unlink(file_temp, recursive = TRUE)

Description

This function can be use to call the bowtie2-inspect wrapper which wraps the bowtie2-inspect-s and bowtie2-inspect-l binaries.

Usage

metascope_inspect(bt2_base, ...)

Arguments

Details

All additional arguments in ... are intepreted as parameters to be passed on to bowtie2-inspect wrapper. All of them Character or Numeric scalar. You can put all additional arguments in one Character (e.g. "–across 60 –names") with white space separation, or put them in different Character (e.g. "–across","60","–names"). See the output of metascope_inspect_usage() for details about available parameters.

References

Langmead B, Wilks C, Antonescu V, Charles R. Scaling read aligners to hundreds of threads on general-purpose processors. Bioinformatics. 2018 Jul 18. doi: 10.1093/bioinformatics/bty648.

Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nature Methods. 2012 Mar 4;9(4):357-9. doi: 10.1038/nmeth.1923.

Examples

## Not run: 
metascope_inspect(bt2_base)

## End(Not run)

Description

Calling metascope_inspect_usage() prints the available arguments that can be passed to the ... argument of the metascope_inspect() function of the package. Note that some arguments are invalid if they are already handled as explicit function arguments.

Usage

metascope_inspect_usage()

Value

Information about available arguments that can be passed to metascope_inspect()

References

Langmead B, Wilks C, Antonescu V, Charles R. Scaling read aligners to hundreds of threads on general-purpose processors. Bioinformatics. 2018 Jul 18. doi: 10.1093/bioinformatics/bty648.

Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nature Methods. 2012 Mar 4;9(4):357-9. doi: 10.1038/nmeth.1923.

Examples

metascope_inspect_usage()

Description

This function is a wrapper for the Rbowtie2::bowtie2_build function. It will create either small (.bt2) or large Bowtie2 indexes (.bt2l) depending on the combined size of the reference fasta files.

Usage

mk_bowtie_index(
  ref_dir,
  lib_dir,
  lib_name,
  bowtie2_build_options,
  threads = 1,
  overwrite = FALSE
)

Arguments

Value

Creates the Bowtie2 indexes of the supplied reference .fasta files. Returns the path to the directory containing these files.

Examples

#### Create a bowtie index from the example reference library

## Create a temporary directory to store the reference library
ref_temp <- tempfile()
dir.create(ref_temp)

tmp_accession <- system.file("extdata", "example_accessions.sql", package = "MetaScope")

## Download reference genome
download_refseq('Bovismacovirus', reference = FALSE, representative = FALSE,
                out_dir = ref_temp, compress = TRUE, patho_out = FALSE,
                caching = TRUE, accession_path = tmp_accession)

## Create the reference library index files in the current directory
mk_bowtie_index(ref_dir = ref_temp, lib_dir = ref_temp,
                lib_name = "target", threads = 1, overwrite = FALSE)

## Remove temporary directory
unlink(ref_temp, recursive = TRUE)

Description

This function is a wrapper for the Rsubread::buildindex function. It will generate one or more Subread indexes from a .fasta file. If the library is too large (default >4GB) it will automatically be split into multiple indexes, with _1, _2, etc at the end of the ref_lib basename.

Usage

mk_subread_index(ref_lib, split = 4, mem = 8000, quiet = TRUE)

Arguments

Value

Creates one or more Subread indexes for the supplied reference .fasta file. If multiple indexes are created, the libraries will be named the ref_lib basename + "_1", "_2", etc. The function returns the names of the folders holding these files.

Examples

#### Create a subread index from the example reference library

## Create a temporary directory to store the reference library
ref_temp <- tempfile()
dir.create(ref_temp)
tmp_accession <- system.file("extdata", "example_accessions.sql", package = "MetaScope")

## Download reference genome
out_fasta <- download_refseq('Orthoebolavirus zairense', reference = FALSE,
                             representative = FALSE, out_dir = ref_temp,
                             compress = TRUE, patho_out = FALSE,
                             caching = TRUE, accession_path = tmp_accession)

## Make subread index of reference library
mk_subread_index(out_fasta)
unlink(ref_temp)

Description

Using the filter_host() function, we align our sequencing sample to all filter libraries of interest. The remove_matches() function allows for removal of any target reads that are also aligned to filter libraries.

Usage

remove_matches(
  reads_bam,
  read_names,
  output,
  YS,
  threads,
  aligner,
  make_bam,
  quiet
)

Arguments

Details

This function is not intended for direct use.

Value

Depending on input make_bam, either the name of a filtered, sorted .bam file written to the user's current working directory, or an RDS file containing a data frame of only requisite information to run metascope_id().

Description

Converts NCBI taxonomy ID to scientific name

Usage

taxid_to_name(taxids, accession_path)

Arguments

Value

Returns a dataframe of blast results for a metascope result