## ----------------------------------------------------------------------------- library(single) pos_start <- 1 pos_end <- 100 refseq_fasta <- system.file("extdata", "ref_seq.fasta", package = "single") ref_seq <- Biostrings::subseq(Biostrings::readDNAStringSet(refseq_fasta), pos_start,pos_end) ## ----train model-------------------------------------------------------------- REF_READS <- system.file("extdata", "train_seqs_500.sorted.bam",package = "single") train <- single_train(bamfile=REF_READS, output="train", refseq_fasta=refseq_fasta, rates.matrix=mutation_rate, mean.n.mutations=5.4, pos_start=pos_start, pos_end=pos_end, verbose=FALSE, save_partial=FALSE, save_final= FALSE) print(head(train)) ## ----evaluate model----------------------------------------------------------- LIB_READS <- system.file("extdata","test_sequences.sorted.bam",package ="single") corrected_reads <- single_evaluate(bamfile = LIB_READS, single_fits = train, refseq_fasta=refseq_fasta, pos_start=pos_start,pos_end=pos_end, verbose=FALSE, gaps_weights = "minimum", save = FALSE, save_original_scores = FALSE) corrected_reads ## ----consensus---------------------------------------------------------------- BC_TABLE <- system.file("extdata", "Barcodes_table.txt",package = "single") consensus <- single_consensus_byBarcode(barcodes_table = BC_TABLE, sequences = corrected_reads, verbose = FALSE) consensus ## ----pileup------------------------------------------------------------------- counts_pnq <- pileup_by_QUAL(bam_file=REF_READS, pos_start=pos_start, pos_end=pos_end) head(counts_pnq) ## ----ppriorerrors------------------------------------------------------------- p_prior_errors <- p_prior_errors(counts_pnq=counts_pnq, save=FALSE) p_prior_errors ## ----ppriormutations---------------------------------------------------------- p_prior_mutations <- p_prior_mutations(rates.matrix = mutation_rate, mean.n.mut = 5,ref_seq = ref_seq,save = FALSE) head(p_prior_mutations) ## ----fits--------------------------------------------------------------------- fits <- fit_logregr(counts_pnq = counts_pnq,ref_seq=ref_seq, p_prior_errors = p_prior_errors, p_prior_mutations = p_prior_mutations, save=FALSE) head(fits) ## ----------------------------------------------------------------------------- evaluated_fits <- evaluate_fits(pos_range = c(1,5),q_range = c(0,10), data_fits = fits,ref_seq = ref_seq, save=FALSE,verbose = FALSE) head(evaluated_fits) ## ----sessionInfo, echo=FALSE-------------------------------------------------- sessionInfo()