The most recently developed genome editing system, CRISPR-Cas9, has greater inherent flexibility than prior programmable nuclease platforms because sequence-specific recognition resides primarily within the associated sgRNA, which permits a simple alteration of its recognition sequence. The short Protospacer Adjacent Motif (PAM), which is recognized by Cas9, is the chief constraint on the target site design density. Because of its simplicity and efficacy, this technology is revolutionizing biological studies and holds tremendous promise for therapeutic applications (Ledford 2015; Cox et al. 2015).
However, imperfect cleavage specificity of CRISPR/Cas9 nuclease within the genome is a cause for concern for its therapeutic application. S. pyogenes Cas9 (SpyCas9)-based nucleases can cleave an imperfect heteroduplex formed between the guide sequence and a DNA sequence containing a functional PAM where the number, position and type of base mismatches can impact the level of activity (Hsu et al. 2013; Mali et al. 2013; Fu et al. 2013). This degree of promiscuity is problematic for therapeutic applications, since the generation of DNA breaks at unintended (off-target) sites has the potential to alter gene expression and function through direct mutagenesis or the generation of genomic rearrangements. Experimentally defining the number and activity of off-target sites for a given Cas9/sgRNA complex genome-wide is critical to assess and improve nuclease precision.
A new suite of genome-wide off-target detection methods have recently been described that can identify sites with low cleavage activity within a population of nuclease-treated cells. One of the most sensitive and straightforward methods to employ is GUIDE-seq (Tsai et al. 2015). This method relies on erroneous NHEJ-mediated DNA repair to capture co-introduced blunt-ended double stranded oligonucleotides (dsODNs) at Cas9-induced breakpoints within the genome. The GUIDE-seq dsODNs display high insertion frequency (up to 50% of the measured indel rate (Tsai et al. 2015)) at Cas9-induced DSBs, thereby tagging these loci for selective amplification and subsequent deep sequencing. The method is quite sensitive as off-target sites with >0.1% indel frequency can be detected, and the frequency of dsODN insertion appears to be correlated with the frequency of Cas9-induced lesions at each site (Tsai et al. 2015). This method has been used successfully to evaluate the precision of Cas9 and its variants (tru-sgRNAs (Tsai et al. 2015) or PAM variants (Kleinstiver et al. 2015)). Given its favorable properties, GUIDE-seq could become a standard in the nuclease field for off-target analysis.
While the GUIDE-seq method is straightforward to employ, to date no bioinformatic tools have been released to the community to support the analysis of this data. We developed the GUIDEseq package (Zhu et al. 2017) to facilitate the analysis of GUIDE-seq datasets, including retaining one read per unique molecular identifier (UMI), filtering reads lacking the integration oligo sequence (dsODNs), identifying peak locations (cleavage sites) and heights, merging cleavage sites from the plus strand and those from the minus strand, and performing target and off-target search of the input gRNA. This analysis leverages our ChIPpeakAnno package (Zhu et al. 2010) for merging cleavage sites from the plus strand and minus strand, and the CRISPRseek package (Zhu et al. 2014) for defining the homology of any identified off-target site to the guide sequence and Cas9 PAM specificity.
Here is the workflow of GUIDE-seq data analysis with human sequence. First load GUIDEseq and BSgenome.Hsapiens.UCSC.hg19.
To find the BSgenome of other species, please refer to
available.genomes in the BSgenome
package. For example, BSgenome.Hsapiens.UCSC.hg19
for hg19, BSgenome.Mmusculus.UCSC.mm10
for mm10, BSgenome.Celegans.UCSC.ce6
for ce6, BSgenome.Rnorvegicus.UCSC.rn5
for rn5, BSgenome.Drerio.UCSC.danRer7
for Zv9, and BSgenome.Dmelanogaster.UCSC.dm3
for dm3.
Then specify the alignment file path as
alignment.inputfile, and a umi file path as
umi.inputfile containing the unique molecular identifier
for each sequence.
library(GUIDEseq)
umifile <- system.file("extdata", "UMI-HEK293_site4_chr13.txt",
package = "GUIDEseq")
bedfile <- system.file("extdata", "bowtie2.HEK293_site4_chr13.sort.bed",
package = "GUIDEseq")
bamfile <- system.file("extdata", "bowtie2.HEK293_site4_chr13.sort.bam",
package = "GUIDEseq")The alignment.inputfile is an alignment file in BED
format containing CIGAR information. It contains chromosome, start, end,
readID, mapping quality, strand and CIGAR information as a tab-delimited
file. Here is an example line:
chr13 27629253 27629403 HWI-M01326:156:1:113:4572:6938/1 44 + 150MWhen fastq files are available, scripts for binning reads, removing adaptors and mapping to the genome are available at http://mccb.umassmed.edu/GUIDE-seq/. Otherwise, a one-line, 6-argument pipeline GS-Preprocess at https://github.com/umasstr/GS-Preprocess (Rodriguez2021?) can be used to generate the needed input files for GUIDEseq. The input file for GS-Preprocess is the standard raw data output in BCL file format.
The umi.inputfile is a tab-delimited file containing at
least two columns, read IDs and corresponding unique molecular
identifiers (UMI). A script for creating the umi.inputfile
is available at http://mccb.umassmed.edu/GUIDE-seq/getUmi.pl. An example
input file is at http://mccb.umassmed.edu/GUIDE-seq/testGetUmi/. Please
make sure to use R1 reads as input to getUmi.pl.
PCR amplification often leads to biased representation of the
starting sequence population. To track the sequence tags present in the
initial sequence library, unique molecular identifiers (UMI) are added
to the 5 prime of each sequence in the starting library. The function
getUniqueCleavageEvents uses the UMI sequence in the
umi.inputfile (optionally containing the UMI plus the first
few sequence from the R1 reads) to obtain the starting sequence library.
It also filters out reads that do not contain the integration oligo
sequence, are too short, or are not in the right paired
configuration.
For detailed parameter settings for the function
getUniqueCleavageEvents, please type
help(getUniqueCleavageEvents).
uniqueCleavageEvents <- getUniqueCleavageEvents(bamfile, umifile,
n.cores.max = 1)
uniqueCleavageEvents$cleavage.gr
#> GRanges object with 3841 ranges and 2 metadata columns:
#> seqnames ranges strand | total umi
#> <Rle> <IRanges> <Rle> | <numeric> <character>
#> [1] chr13 69979989 + | 1 CGTATTTT
#> [2] chr13 39262929 + | 1 CTTATAGA
#> [3] chr13 94623776 + | 1 TCATCTTG
#> [4] chr13 28089785 + | 1 TTACCTCC
#> [5] chr13 27629410 + | 1 GCTGCTAT
#> ... ... ... ... . ... ...
#> [3837] chr13 51894561 - | 1 GTTAAACT
#> [3838] chr13 93491789 - | 1 ATTATTCA
#> [3839] chr13 27629409 - | 1 GGAGTTCT
#> [3840] chr13 98262510 - | 1 ACAACAGG
#> [3841] chr13 37190177 - | 1 AATAAATT
#> -------
#> seqinfo: 25 sequences from an unspecified genome; no seqlengthsCalling the function getPeaks with the results from
getUniqueCleavageEvents outputs summarized cleavage events
for each moving window with at least min.reads of cleavage
events.
By default, window size is set to 20, step is set to 20, and
min.reads is set to 2. For detailed parameter settings
using the function getPeaks, please type
help(getPeaks).
peaks <- getPeaks(uniqueCleavageEvents$cleavage.gr, min.reads = 80)
#> computing coverage for plus strand ...
#> computing coverage for minus strand ...
#> call peaks ...
#> finding local max for chromosome: chr13
peaks.gr <- peaks$peaks
peaks.gr
#> GRanges object with 4 ranges and 5 metadata columns:
#> seqnames ranges strand | count bg p.value
#> <Rle> <IRanges> <Rle> | <integer> <numeric> <numeric>
#> [1] chr13 27629413-27629420 + | 146 0.584 1.49081e-291
#> [2] chr13 39262922-39262939 + | 156 0.624 3.20197e-311
#> [3] chr13 27629400-27629416 - | 103 0.412 5.74504e-207
#> [4] chr13 39262918-39262920 - | 157 0.636 2.53405e-312
#> SNratio adjusted.p.value
#> <numeric> <numeric>
#> [1] 250.000 1.49081e-291
#> [2] 250.000 3.20197e-311
#> [3] 250.000 5.74504e-207
#> [4] 246.855 2.53405e-312
#> -------
#> seqinfo: 25 sequences from an unspecified genome; no seqlengthsCalling the function mergePlusMinusPeaks with the output
from getPeaks merges peaks from the plus strand and minus
strand with a specific orientation and within a certain distance
apart.
By default, plus.strand.start.gt.minus.strand.end is set
to TRUE and distance.threshold is set to 40,
i.e., twice the window size. For detailed parameter settings using the
function mergePlusMinusPeaks, please type
help(mergePlusMinusPeaks).
mergedPeaks <- mergePlusMinusPeaks(peaks.gr = peaks.gr,
output.bedfile = "mergedPeaks.bed")
mergedPeaks$mergedPeaks.gr
#> GRanges object with 2 ranges and 2 metadata columns:
#> seqnames ranges
#> <Rle> <IRanges>
#> chr13+:27629413:27629420:chr13-:27629400:27629416 chr13 27629400-27629420
#> chr13+:39262922:39262939:chr13-:39262918:39262920 chr13 39262918-39262939
#> strand | count
#> <Rle> | <numeric>
#> chr13+:27629413:27629420:chr13-:27629400:27629416 + | 249
#> chr13+:39262922:39262939:chr13-:39262918:39262920 + | 313
#> bg
#> <numeric>
#> chr13+:27629413:27629420:chr13-:27629400:27629416 0.996
#> chr13+:39262922:39262939:chr13-:39262918:39262920 1.260
#> -------
#> seqinfo: 25 sequences from an unspecified genome; no seqlengths
head(mergedPeaks$mergedPeaks.bed)
#> seqnames minStart maxEnd names
#> 1 chr13 27629400 27629420 chr13+:27629413:27629420:chr13-:27629400:27629416
#> 2 chr13 39262918 39262939 chr13+:39262922:39262939:chr13-:39262918:39262920
#> totalCount strand
#> 1 249 +
#> 2 313 +Calling the function offTargetAnalysisOfPeakRegions with
the input gRNA, peaks and genome of interest annotates the identified
cleavage sites with sequence homology to the input gRNA. For detailed
parameter settings using the function
offTargetAnalysisOfPeakRegions, please type
help(offTargetAnalysisOfPeakRegions).
library(BSgenome.Hsapiens.UCSC.hg19)
peaks <- system.file("extdata", "T2plus100OffTargets.bed",
package = "CRISPRseek")
gRNAs <- system.file("extdata", "T2.fa",
package = "CRISPRseek")
outputDir <- getwd()
offTargets <- offTargetAnalysisOfPeakRegions(gRNA = gRNAs, peaks = peaks,
format = c("fasta", "bed"),
peaks.withHeader = TRUE, BSgenomeName = Hsapiens,
upstream = 50, downstream = 50, PAM.size = 3, gRNA.size = 20,
PAM = "NGG", PAM.pattern = "(NAG|NGG|NGA)$", max.mismatch = 2,
outputDir = outputDir,
orderOfftargetsBy = "predicted_cleavage_score",
allowed.mismatch.PAM = 2, overwrite = TRUE
)
#> search for gRNAs for input file1...
#> [1] "Scoring ..."
#> finish off-target search in sequence 2
#> finish off-target search in sequence 1
#> finish feature vector building
#> finish score calculation
#> [1] "Done!"The function GUIDEseqAnalysis is a wrapper function that
uses the UMI sequence (or the UMI plus the first few bases of each
sequence from the R1 reads) to estimate the starting sequence library,
piles up reads with a user-defined window and step size, identifies the
cleavage sites, merges cleavage sites from the plus strand and minus
strand, and performs off-target analysis of the extended regions around
the identified cleavage sites. For detailed parameter settings using the
function GUIDEseqAnalysis, please type
help(GUIDEseqAnalysis).
gRNA.file <- system.file("extdata", "gRNA.fa", package = "GUIDEseq")
system.time(guideSeqRes <- GUIDEseqAnalysis(
alignment.inputfile = bamfile,
umi.inputfile = umifile, gRNA.file = gRNA.file,
orderOfftargetsBy = "peak_score",
descending = TRUE, n.cores.max = 1,
BSgenomeName = Hsapiens, min.reads = 1))
#> search for gRNAs for input file1...
#> [1] "Scoring ..."
#> finish off-target search in sequence 2
#> finish off-target search in sequence 1
#> finish feature vector building
#> finish score calculation
#> [1] "Done!"
#> Extract PAM sequence and n.PAM.mismatch.
#> Done with offtarget search!
#> Add gene and exon information to offTargets ....
#> Order offtargets.
#> Save offtargets.
#> user system elapsed
#> 5.540 0.053 5.594
names(guideSeqRes)
#> [1] "offTargets" "merged.peaks" "peaks" "uniqueCleavages"
#> [5] "read.summary" "sequence.depth"#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [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
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] BSgenome.Hsapiens.UCSC.hg19_1.4.3 BSgenome_1.81.0
#> [3] rtracklayer_1.73.0 BiocIO_1.23.3
#> [5] Biostrings_2.81.5 XVector_0.53.0
#> [7] GUIDEseq_1.43.1 GenomicRanges_1.65.1
#> [9] Seqinfo_1.3.0 IRanges_2.47.2
#> [11] S4Vectors_0.51.5 BiocGenerics_0.59.10
#> [13] generics_0.1.4 BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] RColorBrewer_1.1-3 sys_3.4.3
#> [3] jsonlite_2.0.0 magrittr_2.0.5
#> [5] GenomicFeatures_1.65.0 farver_2.1.2
#> [7] rmarkdown_2.31 vctrs_0.7.3
#> [9] multtest_2.69.0 memoise_2.0.1
#> [11] Rsamtools_2.29.0 RCurl_1.98-1.19
#> [13] base64enc_0.1-6 htmltools_0.5.9
#> [15] S4Arrays_1.13.0 BiocBaseUtils_1.15.1
#> [17] progress_1.2.3 lambda.r_1.2.4
#> [19] curl_7.1.0 Rhdf5lib_2.1.0
#> [21] SparseArray_1.13.2 rhdf5_2.57.1
#> [23] sass_0.4.10 bslib_0.11.0
#> [25] httr2_1.2.3 futile.options_1.0.1
#> [27] cachem_1.1.0 rio_1.3.0
#> [29] buildtools_1.0.0 GenomicAlignments_1.49.1
#> [31] ChIPpeakAnno_3.47.2 whisker_0.4.1
#> [33] lifecycle_1.0.5 mltools_0.3.5
#> [35] pkgconfig_2.0.3 Matrix_1.7-5
#> [37] R6_2.6.1 fastmap_1.2.0
#> [39] MatrixGenerics_1.25.0 digest_0.6.39
#> [41] patchwork_1.3.2 AnnotationDbi_1.75.0
#> [43] regioneR_1.45.0 RSQLite_3.53.3
#> [45] filelock_1.0.3 tfruns_1.5.4
#> [47] httr_1.4.8 abind_1.4-8
#> [49] compiler_4.6.1 withr_3.0.3
#> [51] bit64_4.8.2 S7_0.2.2
#> [53] BiocParallel_1.47.0 DBI_1.3.0
#> [55] tensorflow_2.20.0 biomaRt_2.69.0
#> [57] MASS_7.3-65 rappdirs_0.3.4
#> [59] DelayedArray_0.39.3 rjson_0.2.23
#> [61] gtools_3.9.5 tools_4.6.1
#> [63] otel_0.2.0 zip_3.0.0
#> [65] glue_1.8.1 VennDiagram_1.8.2
#> [67] restfulr_0.0.17 InteractionSet_1.41.0
#> [69] rhdf5filters_1.25.0 grid_4.6.1
#> [71] keras_2.16.1 ade4_1.7-24
#> [73] seqinr_4.2-44 gtable_0.3.6
#> [75] tidyr_1.3.2 ensembldb_2.37.3
#> [77] data.table_1.18.4 hms_1.1.4
#> [79] pillar_1.11.1 stringr_1.6.0
#> [81] limma_3.69.2 splines_4.6.1
#> [83] dplyr_1.2.1 BiocFileCache_3.3.0
#> [85] lattice_0.22-9 survival_3.8-6
#> [87] bit_4.6.0 universalmotif_1.31.46
#> [89] tidyselect_1.2.1 RBGL_1.89.0
#> [91] maketools_1.3.2 knitr_1.51
#> [93] ProtGenerics_1.45.0 SummarizedExperiment_1.43.0
#> [95] futile.logger_1.4.9 xfun_0.59
#> [97] Biobase_2.73.1 statmod_1.5.2
#> [99] matrixStats_1.5.0 stringi_1.8.7
#> [101] UCSC.utils_1.9.0 lazyeval_0.2.3
#> [103] yaml_2.3.12 evaluate_1.0.5
#> [105] codetools_0.2-20 cigarillo_1.3.0
#> [107] tibble_3.3.1 hash_2.2.6.4
#> [109] BiocManager_1.30.27 graph_1.91.0
#> [111] cli_3.6.6 reticulate_1.46.0
#> [113] jquerylib_0.1.4 Rcpp_1.1.2
#> [115] GenomeInfoDb_1.49.1 CRISPRseek_1.53.0
#> [117] zeallot_0.2.0 dbplyr_2.6.0
#> [119] png_0.1-9 XML_3.99-0.23
#> [121] parallel_4.6.1 ggplot2_4.0.3
#> [123] blob_1.3.0 prettyunits_1.2.0
#> [125] AnnotationFilter_1.37.0 bitops_1.0-9
#> [127] pwalign_1.9.1 scales_1.4.0
#> [129] purrr_1.2.2 openxlsx_4.2.8.1
#> [131] crayon_1.5.3 rlang_1.3.0
#> [133] KEGGREST_1.53.5 formatR_1.14