# GUIDEseq User’s Guide

# Introduction

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.

# Workflow of GUIDE-seq data analysis

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 + 150M

When 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`.

## Step 1: Remove PCR bias and obtain unique cleavage events

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 seqlengths

## Step 2: Summarize cleavage events

Calling 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 seqlengths

## Step 3: Merge peaks from plus and minus strand

Calling 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      +

## Step 4: Off-target analysis of extended regions around the identified cleavage sites

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!"

## Run all steps in one workflow function

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"

# Session info

    #> 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

# References

Cox, D. B., R. J. Platt, and F. Zhang. 2015. “Therapeutic Genome Editing: Prospects and Challenges.” Journal Article. Nat Med 21 (2): 121–31. https://doi.org/10.1038/nm.3793.
Fu, Y., J. A. Foden, C. Khayter, et al. 2013. “High-Frequency Off-Target Mutagenesis Induced by CRISPR-Cas Nucleases in Human Cells.” Journal Article. Nat Biotechnol 31 (9): 822–26. https://doi.org/10.1038/nbt.2623.
Hsu, P. D., D. A. Scott, J. A. Weinstein, et al. 2013. “DNA Targeting Specificity of RNA-Guided Cas9 Nucleases.” Journal Article. Nat Biotechnol 31 (9): 827–32. https://doi.org/10.1038/nbt.2647.
Kleinstiver, B. P., M. S. Prew, S. Q. Tsai, et al. 2015. “Engineered CRISPR-Cas9 Nucleases with Altered PAM Specificities.” Journal Article. Nature 523 (7561): 481–85. https://doi.org/10.1038/nature14592.
Ledford, H. 2015. “CRISPR, the Disruptor.” Journal Article. Nature 522 (7554): 20–24. https://doi.org/10.1038/522020a.
Mali, P., J. Aach, P. B. Stranges, et al. 2013. “CAS9 Transcriptional Activators for Target Specificity Screening and Paired Nickases for Cooperative Genome Engineering.” Journal Article. Nat Biotechnol 31 (9): 833–38. https://doi.org/10.1038/nbt.2675.
Tsai, S. Q., Z. Zheng, N. T. Nguyen, et al. 2015. “GUIDE-Seq Enables Genome-Wide Profiling of Off-Target Cleavage by CRISPR-Cas Nucleases.” Journal Article. Nat Biotechnol 33 (2): 187–97. https://doi.org/10.1038/nbt.3117.
Zhu, L. J., C. Gazin, N. D. Lawson, et al. 2010. “ChIPpeakAnno: A Bioconductor Package to Annotate ChIP-Seq and ChIP-Chip Data.” Journal Article. BMC Bioinformatics 11: 237. https://doi.org/10.1186/1471-2105-11-237.
Zhu, L. J., B. R. Holmes, N. Aronin, and M. H. Brodsky. 2014. “CRISPRseek: A Bioconductor Package to Identify Target-Specific Guide RNAs for CRISPR-Cas9 Genome-Editing Systems.” Journal Article. PLoS One 9 (9): e108424. https://doi.org/10.1371/journal.pone.0108424.
Zhu, L. J., M. Lawrence, A. Gupta, et al. 2017. “GUIDEseq: A Bioconductor Package to Analyze GUIDE-Seq Datasets for CRISPR-Cas Nucleases.” Journal Article. BMC Genomics 18 (1): 379. https://doi.org/10.1186/s12864-017-3746-y.