1 Introduction

The Rbowtie package provides an R wrapper around the popular bowtie (Langmead et al. 2009) short read aligner and around SpliceMap (Au et al. 2010) a de novo splice junction discovery and alignment tool, which makes use of the bowtie software package.

The package is used by the QuasR (Gaidatzis et al. 2015) bioconductor package to _qu_antify and _a_nnotate _s_hort _r_eads. We recommend to use the QuasR package instead of using Rbowtie directly. The QuasR package provides a simpler interface than Rbowtie and covers the whole analysis workflow of typical ultra-high throughput sequencing experiments, starting from the raw sequence reads, over pre-processing and alignment, up to quantification.

2 Preliminaries

2.1 Citing Rbowtie

If you use Rbowtie (Hahne, Lerch, and Stadler 2012) in your work, you can cite it as follows:

citation("Rbowtie")
## 
## The Rbowtie package contains code from two separate software
## projects. If using bowtie only, it can be cited as Langmead et al.
## (2009). If using also SpliceMap, it can be cited in addition Au et
## al. (2010). The Rbowtie package can be cited using the third
## reference below:
## 
##   Langmead B, Trapnell C, Pop M, Salzberg SL. Ultrafast and
##   memory-efficient alignment of short DNA sequences to the human
##   genome. Genome Biology 10(3):R25 (2009).
## 
##   Au KF, Jiang H, Lin L, Xing Y, Wong WH. Detection of splice
##   junctions from paired-end RNA-seq data by SpliceMap. Nucleic Acids
##   Research, 38(14):4570-8 (2010).
## 
##   Hahne F, Lerch A, Stadler MB. Rbowtie: An R wrapper for bowtie and
##   SpliceMap short read aligners. (unpublished)
## 
## This free open-source software implements academic research by the
## authors and co-workers. If you use it, please support the project by
## citing the appropriate journal articles.
## To see these entries in BibTeX format, use 'print(<citation>,
## bibtex=TRUE)', 'toBibtex(.)', or set
## 'options(citation.bibtex.max=999)'.

2.2 Installation

Rbowtie is a package for the R computing environment and it is assumed that you have already installed R. See the R project at (http://www.r-project.org). To install the latest version of Rbowtie, you will need to be using the latest version of R. Rbowtie is part of the Bioconductor project at (http://www.bioconductor.org). To get Rbowtie together with its dependencies you can use

if (!require("BiocManager"))
    install.packages("BiocManager")
BiocManager::install("Rbowtie")

2.3 Loading of Rbowtie

In order to run the code examples in this vignette, the Rbowtie library need to be loaded.

library(Rbowtie)

2.4 How to get help

Most questions about Rbowtie will hopefully be answered by the documentation or references. If you’ve run into a question which isn’t addressed by the documentation, or you’ve found a conflict between the documentation and software itself, then there is an active support community which can offer help.

The authors of the package (maintainer: Michael Stadler [email protected]) always appreciate receiving reports of bugs in the package functions or in the documentation. The same goes for well-considered suggestions for improvements.

Any other questions or problems concerning Rbowtie should be posted to the Bioconductor support site (https://support.bioconductor.org). Users posting to the support site for the first time should read the helpful posting guide at (https://support.bioconductor.org/info/faq/). Note that each function in Rbowtie has it’s own help page, e.g. help("bowtie"). Posting etiquette requires that you read the relevant help page carefully before posting a problem to the site.

3 Example usage for individual Rbowtie functions

Please refer to the Rbowtie reference manual or the function documentation (e.g. using ?bowtie) for a complete description of Rbowtie functions. The descriptions provided below are meant to give and overview over all functions and summarize the purpose of each one.

3.1 Build the reference index with bowtie_build

To be able to align short reads to a genome, an index has to be build first using the function bowtie_build. Information about arguments can be found with the help of the bowtie_build_usage function or in the manual page ?bowtie_build.

bowtie_build_usage()
##  [1] "Usage: bowtie2-build-s [options]* <reference_in> <ebwt_outfile_base>"            
##  [2] "    reference_in            comma-separated list of files with ref sequences"    
##  [3] "    ebwt_outfile_base       write Ebwt data to files with this dir/basename"     
##  [4] "Options:"                                                                        
##  [5] "    -f                      reference files are Fasta (default)"                 
##  [6] "    -c                      reference sequences given on cmd line (as <seq_in>)" 
##  [7] "    -C/--color              build a colorspace index"                            
##  [8] "    -a/--noauto             disable automatic -p/--bmax/--dcv memory-fitting"    
##  [9] "    -p/--packed             use packed strings internally; slower, uses less mem"
## [10] "    --bmax <int>            max bucket sz for blockwise suffix-array builder"    
## [11] "    --bmaxdivn <int>        max bucket sz as divisor of ref len (default: 4)"    
## [12] "    --dcv <int>             diff-cover period for blockwise (default: 1024)"     
## [13] "    --nodc                  disable diff-cover (algorithm becomes quadratic)"    
## [14] "    -r/--noref              don't build .3/.4.ebwt (packed reference) portion"   
## [15] "    -3/--justref            just build .3/.4.ebwt (packed reference) portion"    
## [16] "    -o/--offrate <int>      SA is sampled every 2^offRate BWT chars (default: 5)"
## [17] "    -t/--ftabchars <int>    # of chars consumed in initial lookup (default: 10)" 
## [18] "    --threads <int>         # of threads"                                        
## [19] "    --ntoa                  convert Ns in reference to As"                       
## [20] "    --seed <int>            seed for random number generator"                    
## [21] "    -q/--quiet              verbose output (for debugging)"                      
## [22] "    -h/--help               print detailed description of tool and its options"  
## [23] "    --usage                 print this usage message"                            
## [24] "    --version               print version information and quit"

refFiles below is a vector with filenames of the reference sequence in FASTA format, and indexDir specifies an output directory for the index files that will be generated when calling bowtie_build:

refFiles <- dir(system.file(package="Rbowtie", "samples", "refs"), full=TRUE)
indexDir <- file.path(tempdir(), "refsIndex")

tmp <- bowtie_build(references=refFiles, outdir=indexDir, prefix="index", force=TRUE)
head(tmp)
## [1] "Settings:"                                                 
## [2] "  Output files: \"/tmp/Rtmp31IA1d/refsIndex/index.*.ebwt\""
## [3] "  Line rate: 6 (line is 64 bytes)"                         
## [4] "  Lines per side: 1 (side is 64 bytes)"                    
## [5] "  Offset rate: 5 (one in 32)"                              
## [6] "  FTable chars: 10"

3.2 Create alignment with bowtie

Information about the arguments supported by the bowtie function can be obtained with the help of the bowtie_usage function or in the manual page ?bowtie.

bowtie_usage()
##  [1] "Usage: "                                                                                        
##  [2] "bowtie-align-s [options]* <ebwt> {-1 <m1> -2 <m2> | --12 <r> | --interleaved <i> | <s>} [<hit>]"
##  [3] ""                                                                                               
##  [4] "  <m1>    Comma-separated list of files containing upstream mates (or the"                      
##  [5] "          sequences themselves, if -c is set) paired with mates in <m2>"                        
##  [6] "  <m2>    Comma-separated list of files containing downstream mates (or the"                    
##  [7] "          sequences themselves if -c is set) paired with mates in <m1>"                         
##  [8] "  <r>     Comma-separated list of files containing Crossbow-style reads.  Can be"               
##  [9] "          a mixture of paired and unpaired.  Specify \"-\" for stdin."                          
## [10] "  <i>     Files with interleaved paired-end FASTQ reads."                                       
## [11] "  <s>     Comma-separated list of files containing unpaired reads, or the"                      
## [12] "          sequences themselves, if -c is set.  Specify \"-\" for stdin."                        
## [13] "  <hit>   File to write hits to (default: stdout)"                                              
## [14] "Input:"                                                                                         
## [15] "  -q                 query input files are FASTQ .fq/.fastq (default)"                          
## [16] "  -f                 query input files are (multi-)FASTA .fa/.mfa"                              
## [17] "  -r                 query input files are raw one-sequence-per-line"                           
## [18] "  -c                 query sequences given on cmd line (as <mates>, <singles>)"                 
## [19] "  -C                 reads and index are in colorspace"                                         
## [20] "  -Q/--quals <file>  QV file(s) corresponding to CSFASTA inputs; use with -f -C"                
## [21] "  --Q1/--Q2 <file>   same as -Q, but for mate files 1 and 2 respectively"                       
## [22] "  -s/--skip <int>    skip the first <int> reads/pairs in the input"                             
## [23] "  -u/--qupto <int>   stop after first <int> reads/pairs (excl. skipped reads)"                  
## [24] "  -5/--trim5 <int>   trim <int> bases from 5' (left) end of reads"                              
## [25] "  -3/--trim3 <int>   trim <int> bases from 3' (right) end of reads"                             
## [26] "  --phred33-quals    input quals are Phred+33 (default)"                                        
## [27] "  --phred64-quals    input quals are Phred+64 (same as --solexa1.3-quals)"                      
## [28] "  --solexa-quals     input quals are from GA Pipeline ver. < 1.3"                               
## [29] "  --solexa1.3-quals  input quals are from GA Pipeline ver. >= 1.3"                              
## [30] "  --integer-quals    qualities are given as space-separated integers (not ASCII)"               
## [31] "Alignment:"                                                                                     
## [32] "  -v <int>           report end-to-end hits w/ <=v mismatches; ignore qualities"                
## [33] "    or"                                                                                         
## [34] "  -n/--seedmms <int> max mismatches in seed (can be 0-3, default: -n 2)"                        
## [35] "  -e/--maqerr <int>  max sum of mismatch quals across alignment for -n (def: 70)"               
## [36] "  -l/--seedlen <int> seed length for -n (default: 28)"                                          
## [37] "  --nomaqround       disable Maq-like quality rounding for -n (nearest 10 <= 30)"               
## [38] "  -I/--minins <int>  minimum insert size for paired-end alignment (default: 0)"                 
## [39] "  -X/--maxins <int>  maximum insert size for paired-end alignment (default: 250)"               
## [40] "  --fr/--rf/--ff     -1, -2 mates align fw/rev, rev/fw, fw/fw (default: --fr)"                  
## [41] "  --nofw/--norc      do not align to forward/reverse-complement reference strand"               
## [42] "  --maxbts <int>     max # backtracks for -n 2/3 (default: 125, 800 for --best)"                
## [43] "  --pairtries <int>  max # attempts to find mate for anchor hit (default: 100)"                 
## [44] "  -y/--tryhard       try hard to find valid alignments, at the expense of speed"                
## [45] "  --chunkmbs <int>   max megabytes of RAM for best-first search frames (def: 64)"               
## [46] " --reads-per-batch   # of reads to read from input file at once (default: 16)"                  
## [47] "Reporting:"                                                                                     
## [48] "  -k <int>           report up to <int> good alignments per read (default: 1)"                  
## [49] "  -a/--all           report all alignments per read (much slower than low -k)"                  
## [50] "  -m <int>           suppress all alignments if > <int> exist (def: no limit)"                  
## [51] "  -M <int>           like -m, but reports 1 random hit (MAPQ=0); requires --best"               
## [52] "  --best             hits guaranteed best stratum; ties broken by quality"                      
## [53] "  --strata           hits in sub-optimal strata aren't reported (requires --best)"              
## [54] "Output:"                                                                                        
## [55] "  -t/--time          print wall-clock time taken by search phases"                              
## [56] "  -B/--offbase <int> leftmost ref offset = <int> in bowtie output (default: 0)"                 
## [57] "  --quiet            print nothing but the alignments"                                          
## [58] "  --refidx           refer to ref. seqs by 0-based index rather than name"                      
## [59] "  --al <fname>       write aligned reads/pairs to file(s) <fname>"                              
## [60] "  --un <fname>       write unaligned reads/pairs to file(s) <fname>"                            
## [61] "  --no-unal          suppress SAM records for unaligned reads"                                  
## [62] "  --max <fname>      write reads/pairs over -m limit to file(s) <fname>"                        
## [63] "  --suppress <cols>  suppresses given columns (comma-delim'ed) in default output"               
## [64] "  --fullref          write entire ref name (default: only up to 1st space)"                     
## [65] "Colorspace:"                                                                                    
## [66] "  --snpphred <int>   Phred penalty for SNP when decoding colorspace (def: 30)"                  
## [67] "     or"                                                                                        
## [68] "  --snpfrac <dec>    approx. fraction of SNP bases (e.g. 0.001); sets --snpphred"               
## [69] "  --col-cseq         print aligned colorspace seqs as colors, not decoded bases"                
## [70] "  --col-cqual        print original colorspace quals, not decoded quals"                        
## [71] "  --col-keepends     keep nucleotides at extreme ends of decoded alignment"                     
## [72] "SAM:"                                                                                           
## [73] "  -S/--sam           write hits in SAM format"                                                  
## [74] "  --mapq <int>       default mapping quality (MAPQ) to print for SAM alignments"                
## [75] "  --sam-nohead       supppress header lines (starting with @) for SAM output"                   
## [76] "  --sam-nosq         supppress @SQ header lines for SAM output"                                 
## [77] "  --sam-RG <text>    add <text> (usually \"lab=value\") to @RG line of SAM header"              
## [78] "Performance:"                                                                                   
## [79] "  -o/--offrate <int> override offrate of index; must be >= index's offrate"                     
## [80] "  -p/--threads <int> number of alignment threads to launch (default: 1)"                        
## [81] "  --mm               use memory-mapped I/O for index; many 'bowtie's can share"                 
## [82] "  --shmem            use shared mem for index; many 'bowtie's can share"                        
## [83] "Other:"                                                                                         
## [84] "  --seed <int>       seed for random number generator"                                          
## [85] "  --verbose          verbose output (for debugging)"                                            
## [86] "  --version          print version information and quit"                                        
## [87] "  -h/--help          print this usage message"

In the example below, readsFiles is the name of a file containing short reads to be aligned with bowtie, and samFiles specifies the name of the output file with the generated alignments.

readsFiles <- system.file(package="Rbowtie", "samples", "reads", "reads.fastq")
samFiles <- file.path(tempdir(), "alignments.sam")

bowtie(sequences=readsFiles, 
       index=file.path(indexDir, "index"), 
       outfile=samFiles, sam=TRUE,
       best=TRUE, force=TRUE)
strtrim(readLines(samFiles), 65)
## [1] "@HD\tVN:1.0\tSO:unsorted"                                                  
## [2] "@SQ\tSN:chr1\tLN:100000"                                                   
## [3] "@SQ\tSN:chr2\tLN:100000"                                                   
## [4] "@SQ\tSN:chr3\tLN:100000"                                                   
## [5] "@PG\tID:Bowtie\tVN:1.2.2\tCL:\"/tmp/RtmpX4sXpe/Rinstdde73ecfa1e/Rbowt"     
## [6] "HWUSI-EAS1513_0012:6:48:5769:946#0/1\t0\tchr1\t819\t255\t101M\t*\t0\t0\tTG"
## [7] "HWUSI-EAS1513_0012:6:48:6908:952#0/1\t0\tchr2\t1133\t255\t101M\t*\t0\t0\tA"
## [8] "HWUSI-EAS1513_0012:6:48:8070:953#0/1\t0\tchr1\t7543\t255\t101M\t*\t0\t0\tG"

3.3 Create spliced alignment with SpliceMap

While bowtie only generates ungapped alignments, the SpliceMap function can be used to generate spliced alignments. SpliceMap is itself using bowtie. To use it, it is necessary to create an index of the reference sequence as described in ??. SpliceMap parameters are specified in the form of a named list, which follows closely the configure file format of the original SpliceMap program(Au et al. 2010). Be aware that SpliceMap can only be used for reads that are at least 50bp long.

readsFiles <- system.file(package="Rbowtie", "samples", "reads", "reads.fastq")
refDir <- system.file(package="Rbowtie", "samples", "refs", "chr1.fa")
indexDir <- file.path(tempdir(), "refsIndex")
samFiles <- file.path(tempdir(), "splicedAlignments.sam")

cfg <- list(genome_dir=refDir,
            reads_list1=readsFiles,
            read_format="FASTQ",
            quality_format="phred-33",
            outfile=samFiles,
            temp_path=tempdir(),
            max_intron=400000,
            min_intron=20000,
            max_multi_hit=10,
            seed_mismatch=1,
            read_mismatch=2,
            num_chromosome_together=2,
            bowtie_base_dir=file.path(indexDir, "index"),
            num_threads=4,
            try_hard="yes",
            selectSingleHit=TRUE)
res <- SpliceMap(cfg)
## [SpliceMap] splitting reads into 25mers...done
## [SpliceMap] aligning 25mers...done
## [SpliceMap] sorting 25mer-alignments...done
## [SpliceMap] indexing 25mer-alignments...done
## [SpliceMap] finding spliced alignments for 1 chromosomes using 2 parallel processes...done
## [SpliceMap] extracting unmapped reads (single read mode)...done
## [SpliceMap] combining spliced alignments...done
## [SpliceMap] reordering final output...done
res
## [1] "/tmp/Rtmp31IA1d/splicedAlignments.sam"
strtrim(readLines(samFiles), 65)
## [1] "@HD\tVN:1.0\tSO:coordinate"                                                
## [2] "@SQ\tSN:chr1\tLN:100000"                                                   
## [3] "@PG\tID:SpliceMap\tVN:3.3.5.2 (55)"                                        
## [4] "HWUSI-EAS1513_0012:6:48:5769:946#0\t0\tchr1\t819\t255\t101M\t*\t0\t0\tTGGA"
## [5] "HWUSI-EAS1513_0012:6:48:6908:952#0\t4\t*\t0\t0\t*\t*\t0\t0\tAACATAGTGAAGAA"
## [6] "HWUSI-EAS1513_0012:6:48:8070:953#0\t0\tchr1\t7543\t255\t101M\t*\t0\t0\tGTC"
## [7] "HWUSI-EAS1513_0012:6:48:9942:949#0\t4\t*\t0\t0\t*\t*\t0\t0\tCGGTTCCTGTATCC"

4 Session information

The output in this vignette was produced under:

sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.3 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.10-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.10-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] Rbowtie_1.26.0   BiocStyle_2.14.0
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.2         bookdown_0.14      digest_0.6.22     
##  [4] magrittr_1.5       evaluate_0.14      rlang_0.4.1       
##  [7] stringi_1.4.3      rmarkdown_1.16     tools_3.6.1       
## [10] stringr_1.4.0      parallel_3.6.1     xfun_0.10         
## [13] yaml_2.2.0         compiler_3.6.1     BiocManager_1.30.9
## [16] htmltools_0.4.0    knitr_1.25

References

Au, K.F., H. Jiang, L. Lin, Y. Xing, and W.H. Wong. 2010. “Detection of Splice Junctions from Paired-End Rna-Seq Data by Splicemap.” Nucleic Acids Research 38 (14):4570–8.

Gaidatzis, D., A. Lerch, F. Hahne, and M.B. Stadler. 2015. “QuasR: Quantify and Annotate Short Reads in R.” Bioinformatics 31 (7):1130–2. https://doi.org/10.1093/bioinformatics/btu781.

Hahne, F., A. Lerch, and M.B. Stadler. 2012. “Rbowtie: A R Wrapper for Bowtie and SpliceMap Short Read Aligners.”

Langmead, B., C. Trapnell, M. Pop, and S.L. Salzberg. 2009. “Ultrafast and Memory-Efficient Alignment of Short Dna Sequences to the Human Genome.” Genome Biology 10 (3):R25.