### R code from vignette source 'ReQON.Rnw' ################################################### ### code chunk number 1: style ################################################### BiocStyle::latex() ################################################### ### code chunk number 2: ReQON.Rnw:46-54 ################################################### 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 3: ReQON.Rnw:67-88 ################################################### 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 4: ReQON.Rnw:109-117 ################################################### # 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 5: ReQON.Rnw:148-153 ################################################### # show the error counts by read position diagnostics$ReadPosErrors # plot ReadPosErrorPlot( diagnostics$ReadPosErrors, startpos = 0 ) ################################################### ### code chunk number 6: ReQON.Rnw:164-166 ################################################### # plot QualFreqPlot( diagnostics$QualFreqBefore, diagnostics$QualFreqAfter ) ################################################### ### code chunk number 7: ReQON.Rnw:188-200 ################################################### # 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