single_evaluate: Evaluate SINGLE model

In rocioespci/single: Accurate consensus sequence from nanopore reads of a gene library

Evaluate SINGLE model

Description

Main function to evaluate a gene library using a SINGLE model.

Usage

single_evaluate(
  bamfile,
  single_fits,
  refseq_fasta,
  pos_start = NULL,
  pos_end = NULL,
  gaps_weights,
  save = FALSE,
  output_file,
  verbose = FALSE,
  save_original_scores = FALSE
)

Arguments

bamfile	File containing the counts per position returned by samtools mpileup
single_fits	Results of the SINGLE model as returned by single_train(). It can be either the output data.frame or the saved file.
refseq_fasta	Fasta file containing reference sequence
pos_start	Numeric. Position to start analyzing, counting starts from 1 and it refers to reference used for minimap2 alignment.
pos_end	Numeric. Position to stop analyzing, counting starts from 1 and it refers to reference used for minimap2 alignment.
gaps_weights	One of "minimum","none","mean". How to assign qscores to deletions.
save	Logical. Should data be saved in a output_file?
output_file	File name for output, if save=TRUE.
verbose	Logical
save_original_scores	Logical. Should original Qscores be saved? If TRUE, they are stored in file whose name finishes in _original.fastq

Details

Before running single_evaluate_function you have to align your INPUT data to a REFERENCE using minimap2 and count the nucleotides per position using samtools using these lines:

minimap2 -ax map-ont --sam-hit-only REFERENCE.fasta INPUT.fastq >ALIGNMENT.sam

samtools view -S -b ALIGNMENT.sam > ALIGNMENT.bam

samtools sort ALIGNMENT.bam -o ALIGNMENT.sorted.bam

samtools mpileup -Q 0 ALIGNMENT.sorted.bam > COUNTS.txt

Value

Creates file output_prefix_corrected.txt with the Qscores re-scaled by SINGLE. Columns are SeqID position nucleotide isWT original_quality p_SINGLe

Examples

refseq_fasta = system.file("extdata", "ref_seq_10bases.fasta", package = "single")
train_file <- system.file("extdata", "train_example.txt", package = "single")
train <- read.table(train_file, header=TRUE)
test_reads_example <- system.file("extdata", "test_sequences.sorted.bam",
   package = "single")
corrected_reads <- single_evaluate(bamfile = test_reads_example,
                 single_fits = train,refseq_fasta = refseq_fasta,
                 pos_start=1,pos_end=10,gaps_weights = "minimum")
corrected_reads

rocioespci/single documentation built on April 18, 2023, 8:48 p.m.