1 Installation

crisprBwa can be installed from Bioconductor using the following commands in a fresh R session:

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("crisprBwa")

2 Overview of crisprBwa

crisprBwa provides two main functions to align short DNA sequences to a reference genome using the short read aligner BWA-backtrack (Li and Durbin 2009) and return the alignments as R objects: runBwa and runCrisprBwa``. It utilizes the Bioconductor packageRbwato access the BWA program in a platform-independent manner. This means that users do not need to install BWA prior to usingcrisprBwa`.

The latter function (runCrisprBwa) is specifically designed to map and annotate CRISPR guide RNA (gRNA) spacer sequences using CRISPR nuclease objects and CRISPR genomic arithmetics defined in the Bioconductor crisprBase package. This enables a fast and accurate on-target and off-target search of gRNA spacer sequences for virtually any type of CRISPR nucleases.

3 Building a bwa index

To use runBwa or runCrisprBwa, users need to first build a BWA genome index. For a given genome, this step has to be done only once. The Rbwa package conveniently provides the function bwa_build_index to build a BWA index from any custom genome from a FASTA file.

As an example, we build a BWA index for a small portion of the human chromosome 12 (chr12.fa file provided in the crisprBwa package) and save the index file as myIndex to a temporary directory:

library(Rbwa)
fasta <- system.file(package="crisprBwa", "example/chr12.fa")
outdir <- tempdir()
index <- file.path(outdir, "chr12")
Rbwa::bwa_build_index(fasta,
                      index_prefix=index)

4 Alignment using runCrisprBwa

As an example, we align 5 spacer sequences (of length 20bp) to the custom genome built above, allowing a maximum of 3 mismatches between the spacer and protospacer sequences.

We specify that the search is for the wildtype Cas9 (SpCas9) nuclease by providing the CrisprNuclease object SpCas9 available through the crisprBase package. The argument canonical=FALSE specifies that non-canonical PAM sequences are also considered (NAG and NGA for SpCas9). The function getAvailableCrisprNucleases in crisprBase returns a character vector of available crisprNuclease objects found in crisprBase.

We also need to provide a BSgenome object corresponding to the reference genome used for alignment to extract protospacer and PAM sequences of the target sequences.

library(crisprBwa)
library(BSgenome.Hsapiens.UCSC.hg38)
## Loading required package: BSgenome
## Loading required package: BiocGenerics
## 
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
## 
##     Filter, Find, Map, Position, Reduce, anyDuplicated, append,
##     as.data.frame, basename, cbind, colnames, dirname, do.call,
##     duplicated, eval, evalq, get, grep, grepl, intersect, is.unsorted,
##     lapply, mapply, match, mget, order, paste, pmax, pmax.int, pmin,
##     pmin.int, rank, rbind, rownames, sapply, setdiff, sort, table,
##     tapply, union, unique, unsplit, which.max, which.min
## Loading required package: S4Vectors
## Loading required package: stats4
## 
## Attaching package: 'S4Vectors'
## The following objects are masked from 'package:base':
## 
##     I, expand.grid, unname
## Loading required package: IRanges
## Loading required package: GenomeInfoDb
## Loading required package: GenomicRanges
## Loading required package: Biostrings
## Loading required package: XVector
## 
## Attaching package: 'Biostrings'
## The following object is masked from 'package:base':
## 
##     strsplit
## Loading required package: rtracklayer
data(SpCas9, package="crisprBase")
crisprNuclease <- SpCas9
bsgenome <- BSgenome.Hsapiens.UCSC.hg38
spacers <- c("AGCTGTCCGTGGGGGTCCGC",
             "CCCCTGCTGCTGTGCCAGGC",
             "ACGAACTGTAAAAGGCTTGG",
             "ACGAACTGTAACAGGCTTGG",
             "AAGGCCCTCAGAGTAATTAC")
runCrisprBwa(spacers,
             bsgenome=bsgenome,
             crisprNuclease=crisprNuclease,
             n_mismatches=3,
             canonical=FALSE,
             bwa_index=index)
## [runCrisprBwa] Using BSgenome.Hsapiens.UCSC.hg38 
## [runCrisprBwa] Searching for SpCas9 protospacers
##                 spacer          protospacer pam   chr pam_site strand
## 1 AAGGCCCTCAGAGTAATTAC AAGGCCCTCAGAGTAATTAC AGA chr12   170636      +
## 2 ACGAACTGTAAAAGGCTTGG ACGAACTGTAAAAGGCTTGG AGG chr12   170815      -
## 3 ACGAACTGTAACAGGCTTGG ACGAACTGTAAAAGGCTTGG AGG chr12   170815      -
## 4 AGCTGTCCGTGGGGGTCCGC AGCTGTCCGTGGGGGTCCGC AGG chr12   170585      +
## 5 CCCCTGCTGCTGTGCCAGGC CCCCTGCTGCTGTGCCAGGC CGG chr12   170609      +
##   n_mismatches canonical
## 1            0     FALSE
## 2            0      TRUE
## 3            1      TRUE
## 4            0      TRUE
## 5            0      TRUE

5 Applications beyond CRISPR

The function runBwa is similar to runCrisprBwa, but does not impose constraints on PAM sequences. It can be used to search for any short read sequence in a genome.

5.1 Example using RNAi (siRNA design)

Seed-related off-targets caused by mismatch tolerance outside of the seed region is a well-studied and characterized problem observed in RNA interference (RNAi) experiments. runBWa can be used to map shRNA/siRNA seed sequences to reference genomes to predict putative off-targets:

seeds <- c("GTAAGCGGAGTGT", "AACGGGGAGATTG")
runBwa(seeds,
       n_mismatches=2,
       bwa_index=index)
##           query   chr    pos strand n_mismatches
## 1 AACGGGGAGATTG chr12  68337      -            2
## 2 AACGGGGAGATTG chr12   1666      -            2
## 3 AACGGGGAGATTG chr12 123863      +            2
## 4 AACGGGGAGATTG chr12 151731      -            2
## 5 AACGGGGAGATTG chr12 110901      +            2
## 6 GTAAGCGGAGTGT chr12 101550      -            2

6 Session info

sessionInfo()
## R version 4.2.0 RC (2022-04-19 r82224)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.15-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.15-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB              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] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] BSgenome.Hsapiens.UCSC.hg38_1.4.4 BSgenome_1.64.0                  
##  [3] rtracklayer_1.56.0                Biostrings_2.64.0                
##  [5] XVector_0.36.0                    GenomicRanges_1.48.0             
##  [7] GenomeInfoDb_1.32.0               IRanges_2.30.0                   
##  [9] S4Vectors_0.34.0                  BiocGenerics_0.42.0              
## [11] crisprBwa_1.0.0                   Rbwa_1.0.0                       
## [13] BiocStyle_2.24.0                 
## 
## loaded via a namespace (and not attached):
##  [1] lattice_0.20-45             Rsamtools_2.12.0           
##  [3] digest_0.6.29               utf8_1.2.2                 
##  [5] R6_2.5.1                    evaluate_0.15              
##  [7] pillar_1.7.0                zlibbioc_1.42.0            
##  [9] rlang_1.0.2                 jquerylib_0.1.4            
## [11] Matrix_1.4-1                rmarkdown_2.14             
## [13] BiocParallel_1.30.0         readr_2.1.2                
## [15] stringr_1.4.0               bit_4.0.4                  
## [17] RCurl_1.98-1.6              DelayedArray_0.22.0        
## [19] compiler_4.2.0              xfun_0.30                  
## [21] pkgconfig_2.0.3             htmltools_0.5.2            
## [23] tidyselect_1.1.2            SummarizedExperiment_1.26.0
## [25] tibble_3.1.6                GenomeInfoDbData_1.2.8     
## [27] bookdown_0.26               matrixStats_0.62.0         
## [29] XML_3.99-0.9                fansi_1.0.3                
## [31] crayon_1.5.1                tzdb_0.3.0                 
## [33] GenomicAlignments_1.32.0    bitops_1.0-7               
## [35] grid_4.2.0                  jsonlite_1.8.0             
## [37] lifecycle_1.0.1             magrittr_2.0.3             
## [39] vroom_1.5.7                 cli_3.3.0                  
## [41] stringi_1.7.6               crisprBase_1.0.0           
## [43] bslib_0.3.1                 ellipsis_0.3.2             
## [45] vctrs_0.4.1                 rjson_0.2.21               
## [47] restfulr_0.0.13             tools_4.2.0                
## [49] bit64_4.0.5                 Biobase_2.56.0             
## [51] glue_1.6.2                  purrr_0.3.4                
## [53] hms_1.1.1                   MatrixGenerics_1.8.0       
## [55] parallel_4.2.0              fastmap_1.1.0              
## [57] yaml_2.3.5                  BiocManager_1.30.17        
## [59] knitr_1.39                  sass_0.4.1                 
## [61] BiocIO_1.6.0

References

Li, Heng, and Richard Durbin. 2009. “Fast and Accurate Short Read Alignment with Burrows–Wheeler Transform.” Bioinformatics 25 (14): 1754–60.