metascope_id: Identify which genomes are represented in a processed sample

metascope_idR Documentation

Identify which genomes are represented in a processed sample

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,
  tmp_dir = dirname(input_file),
  convEM = 1/10000,
  maxitsEM = 25,
  update_bam = FALSE,
  num_species_plot = NULL,
  blast_fastas = FALSE,
  num_genomes = 100,
  num_reads = 50,
  quiet = TRUE
)

Arguments

input_file

The .bam or .csv.gz file of sample reads to be identified.

input_type

Extension of file input. Should be either "bam" or "csv.gz". Default is "csv.gz".

aligner

The aligner which was used to create the bam file. Default is "bowtie2" but can also be set to "subread" or "other".

db

Currently accepts one of c("ncbi", "silva", "other") Default is "ncbi", appropriate for samples aligned against indices compiled from NCBI whole genome databases. Alternatively, usage of an alternate database (like Greengenes2) should be specified with "other".

db_feature_table

If db = "other", a data.frame must be supplied with two columns, "Feature ID" matching the names of the alignment indices, and a second character column supplying the taxon identifying information.

accession_path

(character) Filepath to NCBI accessions SQL database. See taxonomzr::prepareDatabase().

tmp_dir

Path to a directory to which bam and updated bam files can be saved. Required.

convEM

The convergence parameter of the EM algorithm. Default set at 1/10000.

maxitsEM

The maximum number of EM iterations, regardless of whether the convEM is below the threshhold. Default set at 50. If set at 0, the algorithm skips the EM step and summarizes the .bam file 'as is'

update_bam

Whether to update BAM file with new read assignments. Default is FALSE. If TRUE, requires input_type = TRUE such that a BAM file is the input to the function.

num_species_plot

The number of genome coverage plots to be saved. Default is NULL, which saves coverage plots for the ten most highly abundant species.

blast_fastas

Logical, whether or not to output fasta files for MetaBlast. Default is FALSE.

num_genomes

Number of genomes to output fasta files for MetaBlast. Default is 100.

num_reads

Number of reads per genome per fasta file for MetaBlast. Default is 50.

quiet

Turns off most messages. Default is TRUE.

out_dir

The directory to which the .csv output file will be output. Defaults to dirname(input_file).

Value

This function returns a .csv file with annotated read counts to genomes with mapped reads. The function itself returns the output .csv file name. Depending on the parameters specified, can also output an updated BAM file, and fasta files for usage downstream with MetaBLAST.

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)

#### 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")

#### 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")

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


compbiomed/MetaScope documentation built on Jan. 4, 2025, 7:51 a.m.