### R code from vignette source 'vignettes/ReQON/inst/doc/ReQON.Rnw' ################################################### ### code chunk number 1: ReQON.Rnw:36-44 ################################################### library( ReQON ) library( seqbias ) library( Rsamtools ) reads_fn <- system.file( "extra/example.bam", package = "seqbias" ) # sort BAM sorted <- sortBam( reads_fn, tempfile() ) # index BAM indexBam( sorted ) ################################################### ### code chunk number 2: ReQON.Rnw:57-78 ################################################### ref_fn <- system.file( "extra/example.fa", package = "seqbias" ) ref_f <- FaFile( ref_fn ) open.FaFile( ref_f ) seqs <- scanFa( ref_f ) len <- length( seqs[[1]] ) ref <- matrix( nrow = len, ncol = 3 ) # chromosome/sequence name in column 1 ref[,1] <- rep( "seq1", len ) # sequence position in column 2 ref[,2] <- c( 1:len ) # reference base in column 3 str <- toString( subseq( seqs[[1]], 1, len ) ) s <- strsplit( str, NULL ) ref[,3] <- s[[1]] # look at reference file ref[1:5,] # save reference file write.table( ref, file = "ref_seq.txt", sep = "\t", quote = FALSE, row.names = FALSE, col.names = FALSE ) ################################################### ### code chunk number 3: ReQON.Rnw:99-107 ################################################### # set training region reg <- paste( "seq1:1-", len, sep = "" ) diagnostics <- ReQON( sorted, "Recalibrated_example.bam", region = reg, RefSeq = "ref_seq.txt", nerr = 20, nraf = 0.25, plotname = "Recalibrated_example_plots.pdf" ) # delete reference sequence file unlink( "ref_seq.txt" ) ################################################### ### code chunk number 4: ReQON.Rnw:138-143 ################################################### # show the error counts by read position diagnostics$ReadPosErrors # plot ReadPosErrorPlot( diagnostics$ReadPosErrors, startpos = 0 ) ################################################### ### code chunk number 5: ReQON.Rnw:154-156 ################################################### # plot QualFreqPlot( diagnostics$QualFreqBefore, diagnostics$QualFreqAfter ) ################################################### ### code chunk number 6: ReQON.Rnw:178-190 ################################################### # report FWSE from ReQON output diagnostics$FWSE # plot before recalibration and report FWSE f1 <- FWSEplot(diagnostics$ErrRatesBefore, diagnostics$QualFreqBefore, col = "blue", main_title = "Reported vs. Empirical Quality: Before") f1 # plot after recalibration and report FWSE f2 <- FWSEplot(diagnostics$ErrRatesAfter, diagnostics$QualFreqAfter, col = "red", main_title = "Reported vs. Empirical Quality: After") f2