crisprDesign is the core package of the crisprVerse ecosystem, and
plays the role of a one-stop shop for designing and annotating CRISPR
guide RNA (gRNA) sequences. This includes the characterization of
on-targets and off-targets using different aligners, on- and off-target
scoring, gene context annotation, SNP annotation, sequence feature
characterization, repeat annotation, and many more.
The software was developed to be as applicable and generalizable as
possible.
It currently support five types of CRISPR modalities (modes of perturbations): CRISPR knockout (CRISPRko), CRISPR activation (CRISPRa), CRISPR interference (CRISPRi), CRISPR base editing (CRISPRbe), and CRISPR knockdown (CRISPRkd) (see Kampmann (2018) for a review of CRISPR modalities).
It utilizes the crisprBase package to enable gRNA design
for any CRISPR nuclease and base editor via the
CrisprNuclease and BaseEditor classes,
respectively. Nucleases that are commonly used in the field are
provided, including DNA-targeting nucleases (e.g. SpCas9, AsCas12a) and
RNA-targeting nucleases (e.g. CasRx (RfxCas13d)).
crisprDesign is fully developed to work with the genome
of any organism, and can also be used to design gRNAs targeting custom
DNA sequences.
Finally, more specialized gRNA design functionalities are also available, including design for optical pooled screening (OPS), paired gRNA design, and gRNA filtering and ranking functionalities.
This vignette is meant to be an overview of the main features included in the package, using toy examples for the sake of time (the vignette has to compile within a few minutes, as required by Bioconductor). For detailed and comprehensive tutorials, please visit our crisprVerse tutorials page.
crisprDesign can be installed from from the Bioconductor
devel branch using the following commands in a fresh R session:
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install(version="devel")
BiocManager::install("crisprDesign")Users interested in contributing to crisprDesign might
want to look at the following CRISPR-related package dependencies:
bowtieBWAYou can contribute to the package by submitting pull requests to our GitHub repo.
CRISPR nucleases are examples of RNA-guided endonucleases. They require two binding components for cleavage. First, the nuclease needs to recognize a constant nucleotide motif in the target DNA called the protospacer adjacent motif (PAM) sequence. Second, the gRNA, which guides the nuclease to the target sequence, needs to bind to a complementary sequence adjacent to the PAM sequence, called the protospacer sequence. The latter can be thought of as a variable binding motif that can be specified by designing corresponding gRNA sequences.
The spacer sequence is used in the gRNA construct to guide the CRISPR nuclease to the target protospacer sequence in the host genome.
For DNA-targeting nucleases, the nucleotide sequence of the spacer and protospacer are identical. For RNA-targeting nucleases, they are the reverse complement of each other.
While a gRNA spacer sequence may not always uniquely target the host genome (i.e. it may map to multiple protospacers in the host genome), we can, for a given reference genome, uniquely identify a protospacer sequence with a combination of 3 attributes:
chr: chromosome namestrand: forward (+) or reverse (-)pam_site: genomic coordinate of the first nucleotide of
the nuclease-specific PAM sequence (e.g. for SpCas9, the “N” in the NGG
PAM sequence; for AsCas12a, the first “T” of the TTTV PAM sequence)For CRISPRko, we use an additional genomic coordinate, called
cut_site, to represent where the double-stranded break
(DSB) occurs. For SpCas9, the cut site (blunt-ended dsDNA break) is
located 4nt upstream of the pam_site (PAM-proximal editing). For
AsCas12a, the 5nt 5’ overhang dsDNA break will cause a cut 19nt after
the PAM sequence on the targeted strand, and 23nt after the PAM sequence
on the opposite strand (PAM-distal editing).
We will illustrate the main functionalities of
crisprDesign by performing a common task: designing gRNAs
to knock out a coding gene. In our example, we will design gRNAs for the
wildtype SpCas9 nuclease, with spacers having a length of 20nt.
The crisprBase package provides functionalities to
create objects that store information about CRISPR nucleases, and
functions to interact with those objects (see the
crisprBase vignette). It also provides commonly-used CRISPR
nucleases. Let’s look at the SpCas9 nuclease object:
## Class: CrisprNuclease
## Name: SpCas9
## Target type: DNA
## Metadata: list of length 1
## PAMs: NGG, NAG, NGA
## Weights: 1, 0.2593, 0.0694
## Spacer length: 20
## PAM side: 3prime
## Distance from PAM: 0
## Prototype protospacers: 5'--SSSSSSSSSSSSSSSSSSSS[NGG]--3', 5'--SSSSSSSSSSSSSSSSSSSS[NAG]--3', 5'--SSSSSSSSSSSSSSSSSSSS[NGA]--3'The three motifs (NGG, NAG and NGA) represent the recognized PAM sequences by SpCas9, and the weights indicate a recognition score. The canonical PAM sequence NGG is fully recognized (weight of 1), while the two non-canonical PAM sequences NAG and NGA are much less tolerated.
The spacer sequence is located on the 5-prime end with respect to the
PAM sequence, and the default spacer sequence length is 20 nucleotides.
If necessary, we can change the spacer length using the function
crisprBase::spacerLength. Let’s see what the protospacer
construct looks like by using prototypeSequence:
## [1] "5'--SSSSSSSSSSSSSSSSSSSS[NGG]--3'"As an example, we will design gRNAs that knockout the human gene IQSEC3 by finding all protospacer sequences located in the coding region (CDS) of IQSEC3.
To do so, we need to create a GRanges object that
defines the genomic coordinates of the CDS of IQSEC3 in a reference
genome.
The toy dataset grListExample object in
crisprDesign contains gene coordinates in hg38 for exons of
all human IQSEC3 isoforms, and was obtained by converting an Ensembl
TxDb object into a GRangesList object using
the TxDb2GRangesList convenience function in
crisprDesign.
The queryTxObject function allows us to query such
objects for a specific gene and feature. Here, we obtain a
GRanges object containing the CDS coordinates of
IQSEC3:
gr <- queryTxObject(txObject=grListExample,
featureType="cds",
queryColumn="gene_symbol",
queryValue="IQSEC3")We will only consider the first exon to speed up design:
findSpacers is the main function to obtain a list of all
possible spacer sequences targeting protospacers located in the target
DNA sequence(s). If a GRanges object is provided as input,
a BSgenome object (object containing sequences of a
reference genome) will need to be provided as well:
library(BSgenome.Hsapiens.UCSC.hg38)
bsgenome <- BSgenome.Hsapiens.UCSC.hg38
guideSet <- findSpacers(gr,
bsgenome=bsgenome,
crisprNuclease=SpCas9)
guideSet## GuideSet object with 123 ranges and 5 metadata columns:
## seqnames ranges strand | protospacer pam
## <Rle> <IRanges> <Rle> | <DNAStringSet> <DNAStringSet>
## spacer_1 chr12 66893 - | CGCGCACCGGATTCTCCAGC AGG
## spacer_2 chr12 66896 + | GGGCGGCATGGAGAGCCTGC TGG
## spacer_3 chr12 66905 + | GGAGAGCCTGCTGGAGAATC CGG
## spacer_4 chr12 66906 - | AGGTAGAGCACGGCGCGCAC CGG
## spacer_5 chr12 66916 - | GAGCTCCTTGAGGTAGAGCA CGG
## ... ... ... ... . ... ...
## spacer_119 chr12 67407 + | CACAAATCCCCCTCCGCCCT CGG
## spacer_120 chr12 67412 + | ATCCCCCTCCGCCCTCGGCA AGG
## spacer_121 chr12 67413 + | TCCCCCTCCGCCCTCGGCAA GGG
## spacer_122 chr12 67421 - | CTCACTCAGGTCTCCTGCTC AGG
## spacer_123 chr12 67426 + | TCGGCAAGGGCGTCCTGAGC AGG
## pam_site cut_site region
## <numeric> <numeric> <character>
## spacer_1 66893 66896 region_1
## spacer_2 66896 66893 region_1
## spacer_3 66905 66902 region_1
## spacer_4 66906 66909 region_1
## spacer_5 66916 66919 region_1
## ... ... ... ...
## spacer_119 67407 67404 region_1
## spacer_120 67412 67409 region_1
## spacer_121 67413 67410 region_1
## spacer_122 67421 67424 region_1
## spacer_123 67426 67423 region_1
## -------
## seqinfo: 711 sequences (1 circular) from hg38 genome
## crisprNuclease: SpCas9This returns a GuideSet object that stores genomic
coordinates for all spacer sequences found in the regions provided by
gr. The GuideSet object is an extension of a
GenomicRanges object that stores additional information
about gRNAs.
For the subsequent sections, we will only work with a random subset of 20 spacer sequences:
Several accessor functions are provided to extract information about the spacer sequences:
## DNAStringSet object of length 20:
## width seq names
## [1] 20 CCGAGTTGCTGCGCTGCTGC spacer_107
## [2] 20 GCTCTGCTGGTTCTGCACGA spacer_9
## [3] 20 CGGCCGCCGCGTCAGCACCA spacer_74
## [4] 20 GCCCTTGCCGAGGGCGGAGG spacer_112
## [5] 20 GGCCCCGCTGGGGCTGCTCC spacer_76
## ... ... ...
## [16] 20 TCCCCCTCCGCCCTCGGCAA spacer_121
## [17] 20 CGGCAGCGGGGCCGATGACG spacer_34
## [18] 20 GACGAGCCCGGGCGGAGGCT spacer_24
## [19] 20 CTCGTCGATACGCTCTCGCT spacer_13
## [20] 20 CAGTCGCCCCACAAGCATCT spacer_95## DNAStringSet object of length 20:
## width seq names
## [1] 20 CCGAGTTGCTGCGCTGCTGC spacer_107
## [2] 20 GCTCTGCTGGTTCTGCACGA spacer_9
## [3] 20 CGGCCGCCGCGTCAGCACCA spacer_74
## [4] 20 GCCCTTGCCGAGGGCGGAGG spacer_112
## [5] 20 GGCCCCGCTGGGGCTGCTCC spacer_76
## ... ... ...
## [16] 20 TCCCCCTCCGCCCTCGGCAA spacer_121
## [17] 20 CGGCAGCGGGGCCGATGACG spacer_34
## [18] 20 GACGAGCCCGGGCGGAGGCT spacer_24
## [19] 20 CTCGTCGATACGCTCTCGCT spacer_13
## [20] 20 CAGTCGCCCCACAAGCATCT spacer_95## DNAStringSet object of length 20:
## width seq names
## [1] 3 CGG spacer_107
## [2] 3 TGG spacer_9
## [3] 3 CGG spacer_74
## [4] 3 GGG spacer_112
## [5] 3 AGG spacer_76
## ... ... ...
## [16] 3 GGG spacer_121
## [17] 3 GGG spacer_34
## [18] 3 GGG spacer_24
## [19] 3 GGG spacer_13
## [20] 3 GGG spacer_95## spacer_107 spacer_9 spacer_74 spacer_112 spacer_76 spacer_55
## 67371 66943 67233 67396 67244 67153## spacer_107 spacer_9 spacer_74 spacer_112 spacer_76 spacer_55
## 67368 66946 67230 67399 67247 67156The genomic locations stored in the IRanges represent the PAM site locations in the reference genome.
There are specific spacer sequence features, independent of the genomic context of the protospacer sequence, that can reduce or even eliminate gRNA activity:
Use the function addSequenceFeatures to adds these
spacer sequence characteristics to the GuideSet object:
## GuideSet object with 6 ranges and 12 metadata columns:
## seqnames ranges strand | protospacer pam
## <Rle> <IRanges> <Rle> | <DNAStringSet> <DNAStringSet>
## spacer_107 chr12 67371 + | CCGAGTTGCTGCGCTGCTGC CGG
## spacer_9 chr12 66943 - | GCTCTGCTGGTTCTGCACGA TGG
## spacer_74 chr12 67233 + | CGGCCGCCGCGTCAGCACCA CGG
## spacer_112 chr12 67396 - | GCCCTTGCCGAGGGCGGAGG GGG
## spacer_76 chr12 67244 - | GGCCCCGCTGGGGCTGCTCC AGG
## spacer_55 chr12 67153 - | CTGGTCCTGGAGAGGTTCCT GGG
## pam_site cut_site region percentGC polyA polyC
## <numeric> <numeric> <character> <numeric> <logical> <logical>
## spacer_107 67371 67368 region_1 70 FALSE FALSE
## spacer_9 66943 66946 region_1 60 FALSE FALSE
## spacer_74 67233 67230 region_1 80 FALSE FALSE
## spacer_112 67396 67399 region_1 80 FALSE FALSE
## spacer_76 67244 67247 region_1 85 FALSE TRUE
## spacer_55 67153 67156 region_1 60 FALSE FALSE
## polyG polyT startingGGGGG NNGG
## <logical> <logical> <logical> <character>
## spacer_107 FALSE FALSE FALSE CCGG
## spacer_9 FALSE FALSE FALSE ATGG
## spacer_74 FALSE FALSE FALSE ACGG
## spacer_112 FALSE FALSE FALSE GGGG
## spacer_76 TRUE FALSE FALSE CAGG
## spacer_55 FALSE FALSE FALSE TGGG
## -------
## seqinfo: 711 sequences (1 circular) from hg38 genome
## crisprNuclease: SpCas9In order to select gRNAs that are most specific to our target of interest, it is important to avoid gRNAs that target additional loci in the genome with either perfect sequence complementarity (multiple on-targets), or imperfect complementarity through tolerated mismatches (off-targets).
For instance, both the SpCas9 and AsCas12a nucleases can be tolerant to mismatches between the gRNA spacer sequence (RNA) and the protospacer sequence (DNA), thereby making it critical to characterize off-targets to minimize the introduction of double-stranded breaks (DSBs) beyond our intended target.
The addSpacerAlignments function appends a list of
putative on- and off-targets to a GuideSet object using one
of three methods. The first method uses the fast aligner bowtie (Langmead et al. 2009) via the
crisprBowtie package to map spacer sequences to a specified
reference genome. This can be done by specifying
aligner="bowtie" in addSpacerAlignments.
The second method uses the fast aligner BWA via the crisprBwa
package to map spacer sequences to a specified reference genome. This
can be done by specifying aligner="bwa" in
addSpacerAlignments. Note that this is not available for
Windows machines.
The third method uses the package Biostrings to search
for similar sequences in a set of DNA coordinates sequences, usually
provided through a BSGenome object. This can be done by
specifying aligner="biostrings" in
addSpacerAlignments. This is extremely slow, but can be
useful when searching for off-targets in custom short DNA sequences.
We can control the alignment parameters and output using several
function arguments. n_mismatches sets the maximum number of
permitted gRNA:DNA mismatches (up to 3 mismatches).
n_max_alignments specifies the maximum number of alignments
for a given gRNA spacer sequence (1000 by default). The
n_max_alignments parameter may be overruled by setting
all_Possible_alignments=TRUE, which returns all possible
alignments. canonical=TRUE filters out protospacer
sequences that do not have a canonical PAM sequence.
Finally, the txObject argument in
addSpacerAlignmentsused allows users to provide a
TxDb object, or a TxDb object converted in a
GRangesList using the TxDb2GRangesList
function, to annotate genomic alignments with a gene model annotation.
This is useful to understand whether or not off-targets are located in
the CDS of another gene, for instance.
For the sake of time here, we will search only for on- and off-targets located in the beginning of human chr12 where IQSEC3 is located. We note note that users should always perform a genome-wide search as shown in the [CRISPRko design tutorial](https://github.com/crisprVerse/Tutorials/tree/master/Design_CRISPRko_Cas9].
We will use the bowtie method, with a maximum of 2 mismatches. First,
we need to build a bowtie index sequence using the fasta file provided
in crisprDesign. We use the RBowtie package to
build the index:
library(Rbowtie)
fasta <- system.file(package="crisprDesign", "fasta/chr12.fa")
outdir <- tempdir()
Rbowtie::bowtie_build(fasta,
outdir=outdir,
force=TRUE,
prefix="chr12")
bowtie_index <- file.path(outdir, "chr12")For genome-wide off-target search, users will need to create a bowtie index on the whole genome. This is explained in this tutorial.
Finally, we also need to specify a BSgenome object
storing DNA sequences of the human reference genome:
We are now ready to search for on- and off-targets:
guideSet <- addSpacerAlignments(guideSet,
txObject=grListExample,
aligner_index=bowtie_index,
bsgenome=bsgenome,
n_mismatches=2)## Loading required namespace: crisprBwaLet’s look at what was added to the GuideSet:
## GuideSet object with 20 ranges and 19 metadata columns:
## seqnames ranges strand | protospacer pam
## <Rle> <IRanges> <Rle> | <DNAStringSet> <DNAStringSet>
## spacer_107 chr12 67371 + | CCGAGTTGCTGCGCTGCTGC CGG
## spacer_9 chr12 66943 - | GCTCTGCTGGTTCTGCACGA TGG
## spacer_74 chr12 67233 + | CGGCCGCCGCGTCAGCACCA CGG
## spacer_112 chr12 67396 - | GCCCTTGCCGAGGGCGGAGG GGG
## spacer_76 chr12 67244 - | GGCCCCGCTGGGGCTGCTCC AGG
## ... ... ... ... . ... ...
## spacer_121 chr12 67413 + | TCCCCCTCCGCCCTCGGCAA GGG
## spacer_34 chr12 67093 - | CGGCAGCGGGGCCGATGACG GGG
## spacer_24 chr12 67069 - | GACGAGCCCGGGCGGAGGCT GGG
## spacer_13 chr12 66976 - | CTCGTCGATACGCTCTCGCT GGG
## spacer_95 chr12 67308 + | CAGTCGCCCCACAAGCATCT GGG
## pam_site cut_site region percentGC polyA polyC
## <numeric> <numeric> <character> <numeric> <logical> <logical>
## spacer_107 67371 67368 region_1 70 FALSE FALSE
## spacer_9 66943 66946 region_1 60 FALSE FALSE
## spacer_74 67233 67230 region_1 80 FALSE FALSE
## spacer_112 67396 67399 region_1 80 FALSE FALSE
## spacer_76 67244 67247 region_1 85 FALSE TRUE
## ... ... ... ... ... ... ...
## spacer_121 67413 67410 region_1 75 FALSE TRUE
## spacer_34 67093 67096 region_1 80 FALSE FALSE
## spacer_24 67069 67072 region_1 80 FALSE FALSE
## spacer_13 66976 66979 region_1 60 FALSE FALSE
## spacer_95 67308 67305 region_1 60 FALSE TRUE
## polyG polyT startingGGGGG NNGG n0 n1
## <logical> <logical> <logical> <character> <numeric> <numeric>
## spacer_107 FALSE FALSE FALSE CCGG 1 0
## spacer_9 FALSE FALSE FALSE ATGG 1 0
## spacer_74 FALSE FALSE FALSE ACGG 1 0
## spacer_112 FALSE FALSE FALSE GGGG 1 0
## spacer_76 TRUE FALSE FALSE CAGG 1 0
## ... ... ... ... ... ... ...
## spacer_121 FALSE FALSE FALSE AGGG 1 0
## spacer_34 TRUE FALSE FALSE GGGG 1 0
## spacer_24 FALSE FALSE FALSE TGGG 1 0
## spacer_13 FALSE FALSE FALSE TGGG 1 0
## spacer_95 FALSE FALSE FALSE TGGG 1 0
## n2 n0_c n1_c n2_c alignments
## <numeric> <numeric> <numeric> <numeric> <GRangesList>
## spacer_107 0 1 0 0 chr12:67371:+
## spacer_9 0 1 0 0 chr12:66943:-
## spacer_74 0 1 0 0 chr12:67233:+
## spacer_112 0 1 0 0 chr12:67396:-
## spacer_76 0 1 0 0 chr12:67244:-
## ... ... ... ... ... ...
## spacer_121 0 1 0 0 chr12:67413:+
## spacer_34 0 1 0 0 chr12:67093:-
## spacer_24 0 1 0 0 chr12:67069:-
## spacer_13 0 1 0 0 chr12:66976:-
## spacer_95 0 1 0 0 chr12:67308:+
## -------
## seqinfo: 711 sequences (1 circular) from hg38 genome
## crisprNuclease: SpCas9A few columns were added to the GuideSet object to
summarize the number of on- and off-targets for each spacer sequence,
taking into account genomic context:
To look at the individual on- and off-targets and their context, use
the alignments function to retrieve a table of all genomic
alignments stored in the GuideSet object:
## GRanges object with 20 ranges and 14 metadata columns:
## seqnames ranges strand | spacer
## <Rle> <IRanges> <Rle> | <DNAStringSet>
## spacer_107 chr12 67371 + | CCGAGTTGCTGCGCTGCTGC
## spacer_9 chr12 66943 - | GCTCTGCTGGTTCTGCACGA
## spacer_74 chr12 67233 + | CGGCCGCCGCGTCAGCACCA
## spacer_112 chr12 67396 - | GCCCTTGCCGAGGGCGGAGG
## spacer_76 chr12 67244 - | GGCCCCGCTGGGGCTGCTCC
## ... ... ... ... . ...
## spacer_121 chr12 67413 + | TCCCCCTCCGCCCTCGGCAA
## spacer_34 chr12 67093 - | CGGCAGCGGGGCCGATGACG
## spacer_24 chr12 67069 - | GACGAGCCCGGGCGGAGGCT
## spacer_13 chr12 66976 - | CTCGTCGATACGCTCTCGCT
## spacer_95 chr12 67308 + | CAGTCGCCCCACAAGCATCT
## protospacer pam pam_site n_mismatches
## <DNAStringSet> <DNAStringSet> <numeric> <integer>
## spacer_107 CCGAGTTGCTGCGCTGCTGC CGG 67371 0
## spacer_9 GCTCTGCTGGTTCTGCACGA TGG 66943 0
## spacer_74 CGGCCGCCGCGTCAGCACCA CGG 67233 0
## spacer_112 GCCCTTGCCGAGGGCGGAGG GGG 67396 0
## spacer_76 GGCCCCGCTGGGGCTGCTCC AGG 67244 0
## ... ... ... ... ...
## spacer_121 TCCCCCTCCGCCCTCGGCAA GGG 67413 0
## spacer_34 CGGCAGCGGGGCCGATGACG GGG 67093 0
## spacer_24 GACGAGCCCGGGCGGAGGCT GGG 67069 0
## spacer_13 CTCGTCGATACGCTCTCGCT GGG 66976 0
## spacer_95 CAGTCGCCCCACAAGCATCT GGG 67308 0
## canonical cut_site cds fiveUTRs threeUTRs
## <logical> <numeric> <character> <character> <character>
## spacer_107 TRUE 67368 IQSEC3 <NA> <NA>
## spacer_9 TRUE 66946 IQSEC3 <NA> <NA>
## spacer_74 TRUE 67230 IQSEC3 <NA> <NA>
## spacer_112 TRUE 67399 IQSEC3 <NA> <NA>
## spacer_76 TRUE 67247 IQSEC3 <NA> <NA>
## ... ... ... ... ... ...
## spacer_121 TRUE 67410 IQSEC3 <NA> <NA>
## spacer_34 TRUE 67096 IQSEC3 <NA> <NA>
## spacer_24 TRUE 67072 IQSEC3 <NA> <NA>
## spacer_13 TRUE 66979 IQSEC3 <NA> <NA>
## spacer_95 TRUE 67305 IQSEC3 <NA> <NA>
## exons introns intergenic intergenic_distance
## <character> <character> <character> <integer>
## spacer_107 IQSEC3 <NA> <NA> <NA>
## spacer_9 IQSEC3 <NA> <NA> <NA>
## spacer_74 IQSEC3 <NA> <NA> <NA>
## spacer_112 IQSEC3 <NA> <NA> <NA>
## spacer_76 IQSEC3 <NA> <NA> <NA>
## ... ... ... ... ...
## spacer_121 IQSEC3 <NA> <NA> <NA>
## spacer_34 IQSEC3 <NA> <NA> <NA>
## spacer_24 IQSEC3 <NA> <NA> <NA>
## spacer_13 IQSEC3 <NA> <NA> <NA>
## spacer_95 IQSEC3 <NA> <NA> <NA>
## -------
## seqinfo: 711 sequences (1 circular) from hg38 genomeThe functions onTargets and offTargets will
return on-target alignments (no mismatch) and off-target alignment (with
at least one mismatch), respectively. See
?addSpacerAlignments for more details about the different
options and the following outputs:
gRNAs that align to hundreds of different locations are highly
unspecific and undesirable. This can also cause
addSpacerAlignments to be slow. To mitigate this, we
provide addSpacerAlignmentsIterative, an iterative version
of addSpacerAlignments that curtails alignment searches for
gRNAs having more hits than the user-defined threshold (see
?addSpacerAlignmentsIterative).
To remove protospacer sequences located in repeats or low-complexity
DNA sequences (regions identified by RepeatMasker), which are usually
not of interest due to their low specificity, we provide the convenience
function removeRepeats:
After retrieving a list of putative off-targets and on-targets for a
given spacer sequence, we can use addOffTargetScores to
predict the likelihood of the nuclease to cut at the off-targets based
on mismatch tolerance. Currently, only off-target scoring for the SpCas9
nuclease are available (MIT and CFD algorithms):
## GuideSet object with 17 ranges and 22 metadata columns:
## seqnames ranges strand | protospacer pam
## <Rle> <IRanges> <Rle> | <DNAStringSet> <DNAStringSet>
## spacer_107 chr12 67371 + | CCGAGTTGCTGCGCTGCTGC CGG
## spacer_9 chr12 66943 - | GCTCTGCTGGTTCTGCACGA TGG
## spacer_74 chr12 67233 + | CGGCCGCCGCGTCAGCACCA CGG
## spacer_112 chr12 67396 - | GCCCTTGCCGAGGGCGGAGG GGG
## spacer_76 chr12 67244 - | GGCCCCGCTGGGGCTGCTCC AGG
## ... ... ... ... . ... ...
## spacer_71 chr12 67218 - | TGTCCGTGGTGCTGACGCGG CGG
## spacer_121 chr12 67413 + | TCCCCCTCCGCCCTCGGCAA GGG
## spacer_24 chr12 67069 - | GACGAGCCCGGGCGGAGGCT GGG
## spacer_13 chr12 66976 - | CTCGTCGATACGCTCTCGCT GGG
## spacer_95 chr12 67308 + | CAGTCGCCCCACAAGCATCT GGG
## pam_site cut_site region percentGC polyA polyC
## <numeric> <numeric> <character> <numeric> <logical> <logical>
## spacer_107 67371 67368 region_1 70 FALSE FALSE
## spacer_9 66943 66946 region_1 60 FALSE FALSE
## spacer_74 67233 67230 region_1 80 FALSE FALSE
## spacer_112 67396 67399 region_1 80 FALSE FALSE
## spacer_76 67244 67247 region_1 85 FALSE TRUE
## ... ... ... ... ... ... ...
## spacer_71 67218 67221 region_1 70 FALSE FALSE
## spacer_121 67413 67410 region_1 75 FALSE TRUE
## spacer_24 67069 67072 region_1 80 FALSE FALSE
## spacer_13 66976 66979 region_1 60 FALSE FALSE
## spacer_95 67308 67305 region_1 60 FALSE TRUE
## polyG polyT startingGGGGG NNGG n0 n1
## <logical> <logical> <logical> <character> <numeric> <numeric>
## spacer_107 FALSE FALSE FALSE CCGG 1 0
## spacer_9 FALSE FALSE FALSE ATGG 1 0
## spacer_74 FALSE FALSE FALSE ACGG 1 0
## spacer_112 FALSE FALSE FALSE GGGG 1 0
## spacer_76 TRUE FALSE FALSE CAGG 1 0
## ... ... ... ... ... ... ...
## spacer_71 FALSE FALSE FALSE GCGG 1 0
## spacer_121 FALSE FALSE FALSE AGGG 1 0
## spacer_24 FALSE FALSE FALSE TGGG 1 0
## spacer_13 FALSE FALSE FALSE TGGG 1 0
## spacer_95 FALSE FALSE FALSE TGGG 1 0
## n2 n0_c n1_c n2_c alignments inRepeats
## <numeric> <numeric> <numeric> <numeric> <GRangesList> <logical>
## spacer_107 0 1 0 0 chr12:67371:+ FALSE
## spacer_9 0 1 0 0 chr12:66943:- FALSE
## spacer_74 0 1 0 0 chr12:67233:+ FALSE
## spacer_112 0 1 0 0 chr12:67396:- FALSE
## spacer_76 0 1 0 0 chr12:67244:- FALSE
## ... ... ... ... ... ... ...
## spacer_71 0 1 0 0 chr12:67218:- FALSE
## spacer_121 0 1 0 0 chr12:67413:+ FALSE
## spacer_24 0 1 0 0 chr12:67069:- FALSE
## spacer_13 0 1 0 0 chr12:66976:- FALSE
## spacer_95 0 1 0 0 chr12:67308:+ FALSE
## score_cfd score_mit
## <numeric> <numeric>
## spacer_107 1 1
## spacer_9 1 1
## spacer_74 1 1
## spacer_112 1 1
## spacer_76 1 1
## ... ... ...
## spacer_71 1 1
## spacer_121 1 1
## spacer_24 1 1
## spacer_13 1 1
## spacer_95 1 1
## -------
## seqinfo: 711 sequences (1 circular) from hg38 genome
## crisprNuclease: SpCas9Note that this will only work after calling
addSpacerAlignments, as it requires a list of off-targets
for each gRNA entry. The returned GuideSet object has now
the additional columns score_mit and score_cfd
representing the gRNA-level aggregated off-target specificity scores.
The off-target table also contains a cutting likelihood score for each
gRNA and off-target pair:
## GRanges object with 6 ranges and 16 metadata columns:
## seqnames ranges strand | spacer
## <Rle> <IRanges> <Rle> | <DNAStringSet>
## spacer_107 chr12 67371 + | CCGAGTTGCTGCGCTGCTGC
## spacer_9 chr12 66943 - | GCTCTGCTGGTTCTGCACGA
## spacer_74 chr12 67233 + | CGGCCGCCGCGTCAGCACCA
## spacer_112 chr12 67396 - | GCCCTTGCCGAGGGCGGAGG
## spacer_76 chr12 67244 - | GGCCCCGCTGGGGCTGCTCC
## spacer_55 chr12 67153 - | CTGGTCCTGGAGAGGTTCCT
## protospacer pam pam_site n_mismatches
## <DNAStringSet> <DNAStringSet> <numeric> <integer>
## spacer_107 CCGAGTTGCTGCGCTGCTGC CGG 67371 0
## spacer_9 GCTCTGCTGGTTCTGCACGA TGG 66943 0
## spacer_74 CGGCCGCCGCGTCAGCACCA CGG 67233 0
## spacer_112 GCCCTTGCCGAGGGCGGAGG GGG 67396 0
## spacer_76 GGCCCCGCTGGGGCTGCTCC AGG 67244 0
## spacer_55 CTGGTCCTGGAGAGGTTCCT GGG 67153 0
## canonical cut_site cds fiveUTRs threeUTRs
## <logical> <numeric> <character> <character> <character>
## spacer_107 TRUE 67368 IQSEC3 <NA> <NA>
## spacer_9 TRUE 66946 IQSEC3 <NA> <NA>
## spacer_74 TRUE 67230 IQSEC3 <NA> <NA>
## spacer_112 TRUE 67399 IQSEC3 <NA> <NA>
## spacer_76 TRUE 67247 IQSEC3 <NA> <NA>
## spacer_55 TRUE 67156 IQSEC3 <NA> <NA>
## exons introns intergenic intergenic_distance score_cfd
## <character> <character> <character> <integer> <numeric>
## spacer_107 IQSEC3 <NA> <NA> <NA> 1
## spacer_9 IQSEC3 <NA> <NA> <NA> 1
## spacer_74 IQSEC3 <NA> <NA> <NA> 1
## spacer_112 IQSEC3 <NA> <NA> <NA> 1
## spacer_76 IQSEC3 <NA> <NA> <NA> 1
## spacer_55 IQSEC3 <NA> <NA> <NA> 1
## score_mit
## <numeric>
## spacer_107 1
## spacer_9 1
## spacer_74 1
## spacer_112 1
## spacer_76 1
## spacer_55 1
## -------
## seqinfo: 711 sequences (1 circular) from hg38 genomeaddOnTargetScores adds scores from all on-target
efficiency algorithms available in the R package
crisprScore and appends them to the GuideSet.
By default, scores for all available methods for a given nuclease will
be computed. Here, for the sake of time, let’s add only the CRISPRater
score:
## GuideSet object with 6 ranges and 23 metadata columns:
## seqnames ranges strand | protospacer pam
## <Rle> <IRanges> <Rle> | <DNAStringSet> <DNAStringSet>
## spacer_107 chr12 67371 + | CCGAGTTGCTGCGCTGCTGC CGG
## spacer_9 chr12 66943 - | GCTCTGCTGGTTCTGCACGA TGG
## spacer_74 chr12 67233 + | CGGCCGCCGCGTCAGCACCA CGG
## spacer_112 chr12 67396 - | GCCCTTGCCGAGGGCGGAGG GGG
## spacer_76 chr12 67244 - | GGCCCCGCTGGGGCTGCTCC AGG
## spacer_55 chr12 67153 - | CTGGTCCTGGAGAGGTTCCT GGG
## pam_site cut_site region percentGC polyA polyC
## <numeric> <numeric> <character> <numeric> <logical> <logical>
## spacer_107 67371 67368 region_1 70 FALSE FALSE
## spacer_9 66943 66946 region_1 60 FALSE FALSE
## spacer_74 67233 67230 region_1 80 FALSE FALSE
## spacer_112 67396 67399 region_1 80 FALSE FALSE
## spacer_76 67244 67247 region_1 85 FALSE TRUE
## spacer_55 67153 67156 region_1 60 FALSE FALSE
## polyG polyT startingGGGGG NNGG n0 n1
## <logical> <logical> <logical> <character> <numeric> <numeric>
## spacer_107 FALSE FALSE FALSE CCGG 1 0
## spacer_9 FALSE FALSE FALSE ATGG 1 0
## spacer_74 FALSE FALSE FALSE ACGG 1 0
## spacer_112 FALSE FALSE FALSE GGGG 1 0
## spacer_76 TRUE FALSE FALSE CAGG 1 0
## spacer_55 FALSE FALSE FALSE TGGG 1 0
## n2 n0_c n1_c n2_c alignments inRepeats
## <numeric> <numeric> <numeric> <numeric> <GRangesList> <logical>
## spacer_107 0 1 0 0 chr12:67371:+ FALSE
## spacer_9 0 1 0 0 chr12:66943:- FALSE
## spacer_74 0 1 0 0 chr12:67233:+ FALSE
## spacer_112 0 1 0 0 chr12:67396:- FALSE
## spacer_76 0 1 0 0 chr12:67244:- FALSE
## spacer_55 0 1 0 0 chr12:67153:- FALSE
## score_cfd score_mit score_crisprater
## <numeric> <numeric> <numeric>
## spacer_107 1 1 0.782780
## spacer_9 1 1 0.834319
## spacer_74 1 1 0.764870
## spacer_112 1 1 0.795745
## spacer_76 1 1 0.755493
## spacer_55 1 1 0.711902
## -------
## seqinfo: 711 sequences (1 circular) from hg38 genome
## crisprNuclease: SpCas9See the crisprScore vignette for a full description of
the different scores.
Restriction enzymes are usually involved in the gRNA library
synthesis process. Removing gRNAs that contain specific restriction
sites is often necessary. We provide the function
addRestrictionEnzymes to indicate whether or not gRNAs
contain restriction sites for a user-defined set of enzymes:
When no enzymes are specified, the function adds annotation for the
following default enzymes: EcoRI, KpnI, BsmBI, BsaI, BbsI, PacI, ISceI
and MluI. The function also has two additional arguments,
flanking5 and flanking3, to specify nucleotide
sequences flanking the spacer sequence (5’ and 3’, respectively) in the
lentiviral cassette that will be used for gRNA delivery. The function
will effectively search for restriction sites in the full sequence
[flanking5][spacer][flanking3].
The enzymeAnnotation function can be used to retrieve
the added annotation:
## DataFrame with 6 rows and 7 columns
## EcoRI KpnI BsmBI BsaI BbsI PacI
## <logical> <logical> <logical> <logical> <logical> <logical>
## spacer_107 FALSE FALSE FALSE FALSE FALSE FALSE
## spacer_9 FALSE FALSE FALSE FALSE FALSE FALSE
## spacer_74 FALSE FALSE FALSE FALSE FALSE FALSE
## spacer_112 FALSE FALSE FALSE FALSE FALSE FALSE
## spacer_76 FALSE FALSE FALSE FALSE FALSE FALSE
## spacer_55 FALSE FALSE FALSE FALSE FALSE FALSE
## MluI
## <logical>
## spacer_107 FALSE
## spacer_9 FALSE
## spacer_74 FALSE
## spacer_112 FALSE
## spacer_76 FALSE
## spacer_55 FALSEThe function addGeneAnnotation adds transcript- and
gene-level contextual information to gRNAs from a TxDb-like
object:
The gene annotation can be retrieved using the function
geneAnnotation:
## DataFrame with 17 rows and 24 columns
## chr anchor_site strand gene_symbol gene_id
## <factor> <integer> <factor> <character> <character>
## spacer_107 chr12 67368 + IQSEC3 ENSG00000120645
## spacer_9 chr12 66946 - IQSEC3 ENSG00000120645
## spacer_74 chr12 67230 + IQSEC3 ENSG00000120645
## spacer_112 chr12 67399 - IQSEC3 ENSG00000120645
## spacer_76 chr12 67247 - IQSEC3 ENSG00000120645
## ... ... ... ... ... ...
## spacer_71 chr12 67221 - IQSEC3 ENSG00000120645
## spacer_121 chr12 67410 + IQSEC3 ENSG00000120645
## spacer_24 chr12 67072 - IQSEC3 ENSG00000120645
## spacer_13 chr12 66979 - IQSEC3 ENSG00000120645
## spacer_95 chr12 67305 + IQSEC3 ENSG00000120645
## tx_id protein_id exon_id cut_cds
## <character> <character> <character> <logical>
## spacer_107 ENST00000538872 ENSP00000437554 ENSE00002310174 TRUE
## spacer_9 ENST00000538872 ENSP00000437554 ENSE00002310174 TRUE
## spacer_74 ENST00000538872 ENSP00000437554 ENSE00002310174 TRUE
## spacer_112 ENST00000538872 ENSP00000437554 ENSE00002310174 TRUE
## spacer_76 ENST00000538872 ENSP00000437554 ENSE00002310174 TRUE
## ... ... ... ... ...
## spacer_71 ENST00000538872 ENSP00000437554 ENSE00002310174 TRUE
## spacer_121 ENST00000538872 ENSP00000437554 ENSE00002310174 TRUE
## spacer_24 ENST00000538872 ENSP00000437554 ENSE00002310174 TRUE
## spacer_13 ENST00000538872 ENSP00000437554 ENSE00002310174 TRUE
## spacer_95 ENST00000538872 ENSP00000437554 ENSE00002310174 TRUE
## cut_fiveUTRs cut_threeUTRs cut_introns percentCDS aminoAcidIndex
## <logical> <logical> <logical> <numeric> <numeric>
## spacer_107 FALSE FALSE FALSE 13.7 162
## spacer_9 FALSE FALSE FALSE 1.8 22
## spacer_74 FALSE FALSE FALSE 9.8 116
## spacer_112 FALSE FALSE FALSE 14.6 173
## spacer_76 FALSE FALSE FALSE 10.3 122
## ... ... ... ... ... ...
## spacer_71 FALSE FALSE FALSE 9.6 113
## spacer_121 FALSE FALSE FALSE 14.9 176
## spacer_24 FALSE FALSE FALSE 5.4 64
## spacer_13 FALSE FALSE FALSE 2.7 33
## spacer_95 FALSE FALSE FALSE 11.9 141
## downstreamATG percentTx nIsoforms totalIsoforms percentIsoforms
## <numeric> <numeric> <integer> <numeric> <numeric>
## spacer_107 1 8.5 1 2 50
## spacer_9 0 2.5 1 2 50
## spacer_74 0 6.5 1 2 50
## spacer_112 1 8.9 1 2 50
## spacer_76 0 6.8 1 2 50
## ... ... ... ... ... ...
## spacer_71 0 6.4 1 2 50
## spacer_121 1 9.1 1 2 50
## spacer_24 0 4.3 1 2 50
## spacer_13 0 3.0 1 2 50
## spacer_95 1 7.6 1 2 50
## isCommonExon nCodingIsoforms totalCodingIsoforms
## <logical> <integer> <numeric>
## spacer_107 FALSE 1 2
## spacer_9 FALSE 1 2
## spacer_74 FALSE 1 2
## spacer_112 FALSE 1 2
## spacer_76 FALSE 1 2
## ... ... ... ...
## spacer_71 FALSE 1 2
## spacer_121 FALSE 1 2
## spacer_24 FALSE 1 2
## spacer_13 FALSE 1 2
## spacer_95 FALSE 1 2
## percentCodingIsoforms isCommonCodingExon
## <numeric> <logical>
## spacer_107 50 FALSE
## spacer_9 50 FALSE
## spacer_74 50 FALSE
## spacer_112 50 FALSE
## spacer_76 50 FALSE
## ... ... ...
## spacer_71 50 FALSE
## spacer_121 50 FALSE
## spacer_24 50 FALSE
## spacer_13 50 FALSE
## spacer_95 50 FALSEIt contains a lot of information that contextualizes the genomic location of the protospacer sequences.
The ID columns (tx_id, gene_id,
protein_id, exon_id) give Ensembl IDs. The
exon_rank gives the order of the exon for the transcript,
for example “2” indicates it is the second exon (from the 5’ end) in the
mature transcript.
The columns cut_cds, cut_fiveUTRs,
cut_threeUTRs and cut_introns indicate whether
the guide sequence overlaps with CDS, 5’ UTR, 3’ UTR, or an intron,
respectively.
percentCDS gives the location of the
cut_site within the transcript as a percent from the 5’ end
to the 3’ end. aminoAcidIndex gives the number of the
specific amino acid in the protein where the cut is predicted to occur.
downstreamATG shows how many in-frame ATGs are downstream
of the cut_site (and upstream from the defined percent
transcript cutoff, met_cutoff), indicating a potential
alternative translation initiation site that may preserve protein
function.
For more information about the other columns, type
?addGeneAnnotation.
Similarly, one might want to know which protospacer sequences are located within promoter regions of known genes:
data(tssObjectExample, package="crisprDesign")
guideSet <- addTssAnnotation(guideSet,
tssObject=tssObjectExample)
tssAnnotation(guideSet)## DataFrame with 10 rows and 11 columns
## chr anchor_site strand tx_id gene_id
## <factor> <integer> <factor> <character> <character>
## spacer_9 chr12 66946 - ENST00000538872 ENSG00000120645
## spacer_74 chr12 67230 + ENST00000538872 ENSG00000120645
## spacer_76 chr12 67247 - ENST00000538872 ENSG00000120645
## spacer_55 chr12 67156 - ENST00000538872 ENSG00000120645
## spacer_72 chr12 67224 - ENST00000538872 ENSG00000120645
## spacer_54 chr12 67145 + ENST00000538872 ENSG00000120645
## spacer_15 chr12 66995 + ENST00000538872 ENSG00000120645
## spacer_71 chr12 67221 - ENST00000538872 ENSG00000120645
## spacer_24 chr12 67072 - ENST00000538872 ENSG00000120645
## spacer_13 chr12 66979 - ENST00000538872 ENSG00000120645
## gene_symbol promoter tss_id tss_strand tss_pos dist_to_tss
## <character> <character> <character> <character> <integer> <numeric>
## spacer_9 IQSEC3 P1 IQSEC3_P1 + 66767 179
## spacer_74 IQSEC3 P1 IQSEC3_P1 + 66767 463
## spacer_76 IQSEC3 P1 IQSEC3_P1 + 66767 480
## spacer_55 IQSEC3 P1 IQSEC3_P1 + 66767 389
## spacer_72 IQSEC3 P1 IQSEC3_P1 + 66767 457
## spacer_54 IQSEC3 P1 IQSEC3_P1 + 66767 378
## spacer_15 IQSEC3 P1 IQSEC3_P1 + 66767 228
## spacer_71 IQSEC3 P1 IQSEC3_P1 + 66767 454
## spacer_24 IQSEC3 P1 IQSEC3_P1 + 66767 305
## spacer_13 IQSEC3 P1 IQSEC3_P1 + 66767 212For more information, type ?addTssAnnotation.
Common single-nucleotide polymorphisms (SNPs) can change the
on-target and off-target properties of gRNAs by altering the binding.
The function addSNPAnnotation annotates gRNAs with respect
to a reference database of SNPs (stored in a VCF file), specified by the
vcf argument.
VCF files for common SNPs (dbSNPs) can be downloaded from NCBI on the dbSNP website. We include in this package an example VCF file for common SNPs located in the proximity of human gene IQSEC3. This was obtained using the dbSNP151 RefSNP database obtained by subsetting around IQSEC.
vcf <- system.file("extdata",
file="common_snps_dbsnp151_example.vcf.gz",
package="crisprDesign")
guideSet <- addSNPAnnotation(guideSet, vcf=vcf)
snps(guideSet)## DataFrame with 0 rows and 9 columnsThe rs_site_rel gives the relative position of the SNP
with respect to the pam_site. allele_ref and
allele_minor report the nucleotide of the reference and
minor alleles, respectively. MAF_1000G and
MAF_TOPMED report the minor allele frequency (MAF) in the
1000Genomes and TOPMED populations.
Once gRNAs are fully annotated, it is easy to filter out any unwanted
gRNAs since GuideSet objects can be subsetted like regular
vectors in R.
As an example, suppose that we only want to keep gRNAs that have percent GC between 20% and 80% and that do not contain a polyT stretch. This can be achieved using the following lines:
guideSet <- guideSet[guideSet$percentGC>=20]
guideSet <- guideSet[guideSet$percentGC<=80]
guideSet <- guideSet[!guideSet$polyT]Similarly, it is easy to rank gRNAs based on a set of criteria using
the regular order function.
For instance, let’s sort gRNAs by the CRISPRater on-target score:
# Creating an ordering index based on the CRISPRater score:
# Using the negative values to make sure higher scores are ranked first:
o <- order(-guideSet$score_crisprater)
# Ordering the GuideSet:
guideSet <- guideSet[o]
head(guideSet)## GuideSet object with 6 ranges and 28 metadata columns:
## seqnames ranges strand | protospacer pam
## <Rle> <IRanges> <Rle> | <DNAStringSet> <DNAStringSet>
## spacer_9 chr12 66943 - | GCTCTGCTGGTTCTGCACGA TGG
## spacer_112 chr12 67396 - | GCCCTTGCCGAGGGCGGAGG GGG
## spacer_107 chr12 67371 + | CCGAGTTGCTGCGCTGCTGC CGG
## spacer_74 chr12 67233 + | CGGCCGCCGCGTCAGCACCA CGG
## spacer_76 chr12 67244 - | GGCCCCGCTGGGGCTGCTCC AGG
## spacer_121 chr12 67413 + | TCCCCCTCCGCCCTCGGCAA GGG
## pam_site cut_site region percentGC polyA polyC
## <numeric> <numeric> <character> <numeric> <logical> <logical>
## spacer_9 66943 66946 region_1 60 FALSE FALSE
## spacer_112 67396 67399 region_1 80 FALSE FALSE
## spacer_107 67371 67368 region_1 70 FALSE FALSE
## spacer_74 67233 67230 region_1 80 FALSE FALSE
## spacer_76 67244 67247 region_1 85 FALSE TRUE
## spacer_121 67413 67410 region_1 75 FALSE TRUE
## polyG polyT startingGGGGG NNGG n0 n1
## <logical> <logical> <logical> <character> <numeric> <numeric>
## spacer_9 FALSE FALSE FALSE ATGG 1 0
## spacer_112 FALSE FALSE FALSE GGGG 1 0
## spacer_107 FALSE FALSE FALSE CCGG 1 0
## spacer_74 FALSE FALSE FALSE ACGG 1 0
## spacer_76 TRUE FALSE FALSE CAGG 1 0
## spacer_121 FALSE FALSE FALSE AGGG 1 0
## n2 n0_c n1_c n2_c alignments inRepeats
## <numeric> <numeric> <numeric> <numeric> <GRangesList> <logical>
## spacer_9 0 1 0 0 chr12:66943:- FALSE
## spacer_112 0 1 0 0 chr12:67396:- FALSE
## spacer_107 0 1 0 0 chr12:67371:+ FALSE
## spacer_74 0 1 0 0 chr12:67233:+ FALSE
## spacer_76 0 1 0 0 chr12:67244:- FALSE
## spacer_121 0 1 0 0 chr12:67413:+ FALSE
## score_cfd score_mit score_crisprater enzymeAnnotation
## <numeric> <numeric> <numeric> <SplitDataFrameList>
## spacer_9 1 1 0.834319 FALSE:FALSE:FALSE:...
## spacer_112 1 1 0.795745 FALSE:FALSE:FALSE:...
## spacer_107 1 1 0.782780 FALSE:FALSE:FALSE:...
## spacer_74 1 1 0.764870 FALSE:FALSE:FALSE:...
## spacer_76 1 1 0.755493 FALSE:FALSE:FALSE:...
## spacer_121 1 1 0.741315 FALSE:FALSE:FALSE:...
## geneAnnotation tssAnnotation hasSNP
## <SplitDataFrameList> <SplitDataFrameList> <logical>
## spacer_9 chr12:66946:-:... chr12:66946:-:... FALSE
## spacer_112 chr12:67399:-:... :...,... FALSE
## spacer_107 chr12:67368:+:... :...,... FALSE
## spacer_74 chr12:67230:+:... chr12:67230:+:... FALSE
## spacer_76 chr12:67247:-:... chr12:67247:-:... FALSE
## spacer_121 chr12:67410:+:... :...,... FALSE
## snps
## <SplitDataFrameList>
## spacer_9 :...,...
## spacer_112 :...,...
## spacer_107 :...,...
## spacer_74 :...,...
## spacer_76 :...,...
## spacer_121 :...,...
## -------
## seqinfo: 711 sequences (1 circular) from hg38 genome
## crisprNuclease: SpCas9One can also sort gRNAs using several annotation columns. For instance, let’s sort gRNAs using the CRISPRrater score, but also by prioritizing first gRNAs that have no 1-mismatch off-targets:
o <- order(guideSet$n1, -guideSet$score_crisprater)
# Ordering the GuideSet:
guideSet <- guideSet[o]
head(guideSet)## GuideSet object with 6 ranges and 28 metadata columns:
## seqnames ranges strand | protospacer pam
## <Rle> <IRanges> <Rle> | <DNAStringSet> <DNAStringSet>
## spacer_9 chr12 66943 - | GCTCTGCTGGTTCTGCACGA TGG
## spacer_112 chr12 67396 - | GCCCTTGCCGAGGGCGGAGG GGG
## spacer_107 chr12 67371 + | CCGAGTTGCTGCGCTGCTGC CGG
## spacer_74 chr12 67233 + | CGGCCGCCGCGTCAGCACCA CGG
## spacer_76 chr12 67244 - | GGCCCCGCTGGGGCTGCTCC AGG
## spacer_121 chr12 67413 + | TCCCCCTCCGCCCTCGGCAA GGG
## pam_site cut_site region percentGC polyA polyC
## <numeric> <numeric> <character> <numeric> <logical> <logical>
## spacer_9 66943 66946 region_1 60 FALSE FALSE
## spacer_112 67396 67399 region_1 80 FALSE FALSE
## spacer_107 67371 67368 region_1 70 FALSE FALSE
## spacer_74 67233 67230 region_1 80 FALSE FALSE
## spacer_76 67244 67247 region_1 85 FALSE TRUE
## spacer_121 67413 67410 region_1 75 FALSE TRUE
## polyG polyT startingGGGGG NNGG n0 n1
## <logical> <logical> <logical> <character> <numeric> <numeric>
## spacer_9 FALSE FALSE FALSE ATGG 1 0
## spacer_112 FALSE FALSE FALSE GGGG 1 0
## spacer_107 FALSE FALSE FALSE CCGG 1 0
## spacer_74 FALSE FALSE FALSE ACGG 1 0
## spacer_76 TRUE FALSE FALSE CAGG 1 0
## spacer_121 FALSE FALSE FALSE AGGG 1 0
## n2 n0_c n1_c n2_c alignments inRepeats
## <numeric> <numeric> <numeric> <numeric> <GRangesList> <logical>
## spacer_9 0 1 0 0 chr12:66943:- FALSE
## spacer_112 0 1 0 0 chr12:67396:- FALSE
## spacer_107 0 1 0 0 chr12:67371:+ FALSE
## spacer_74 0 1 0 0 chr12:67233:+ FALSE
## spacer_76 0 1 0 0 chr12:67244:- FALSE
## spacer_121 0 1 0 0 chr12:67413:+ FALSE
## score_cfd score_mit score_crisprater enzymeAnnotation
## <numeric> <numeric> <numeric> <SplitDataFrameList>
## spacer_9 1 1 0.834319 FALSE:FALSE:FALSE:...
## spacer_112 1 1 0.795745 FALSE:FALSE:FALSE:...
## spacer_107 1 1 0.782780 FALSE:FALSE:FALSE:...
## spacer_74 1 1 0.764870 FALSE:FALSE:FALSE:...
## spacer_76 1 1 0.755493 FALSE:FALSE:FALSE:...
## spacer_121 1 1 0.741315 FALSE:FALSE:FALSE:...
## geneAnnotation tssAnnotation hasSNP
## <SplitDataFrameList> <SplitDataFrameList> <logical>
## spacer_9 chr12:66946:-:... chr12:66946:-:... FALSE
## spacer_112 chr12:67399:-:... :...,... FALSE
## spacer_107 chr12:67368:+:... :...,... FALSE
## spacer_74 chr12:67230:+:... chr12:67230:+:... FALSE
## spacer_76 chr12:67247:-:... chr12:67247:-:... FALSE
## spacer_121 chr12:67410:+:... :...,... FALSE
## snps
## <SplitDataFrameList>
## spacer_9 :...,...
## spacer_112 :...,...
## spacer_107 :...,...
## spacer_74 :...,...
## spacer_76 :...,...
## spacer_121 :...,...
## -------
## seqinfo: 711 sequences (1 circular) from hg38 genome
## crisprNuclease: SpCas9The rankSpacers function is a convenience function that
implements our recommended rankings for the SpCas9, enAsCas12a and CasRx
nucleases. For a detailed description of our recommended rankings, see
the documentation of rankSpacers by typing
?rankSpacers.
If an Ensembl transcript ID is provided, the ranking function will also take into account the position of the gRNA within the target CDS of the transcript ID in the ranking procedure. Our recommendation is to specify the Ensembl canonical transcript as the representative transcript for the gene. In our example, ENST00000538872 is the canonical transcript for IQSEC3:
For CRISPRa and CRISPRi applications, the CRISPR nuclease is
engineered to lose its endonuclease activity, therefore should not
introduce double-stranded breaks (DSBs). We will use the dead SpCas9
(dSpCas9) nuclease as an example here. Note that users don’t have to
distinguish between dSpCas9 and SpCas9 when specifying the nuclease in
crisprDesign and crisprBase as they do not
differ in terms of the characteristics stored in the
CrisprNuclease object.
CRISPRi: Fusing dSpCas9 with a Krüppel-associated box (KRAB) domain has been shown to be effective at repressing transcription in mammalian cells (Gilbert et al. 2013). The dSpCas9-KRAB fused protein is a commonly-used construct to conduct CRISPR inhibition (CRISPRi) experiments. To achieve optimal inhibition, gRNAs are usually designed targeting the region directly downstream of the gene transcription starting site (TSS).
CRISPRa: dSpCas9 can also be used to activate gene expression by coupling the dead nuclease with activation factors. The technology is termed CRISPR activation (CRISPRa), and several CRISPRa systems have been developed (see Kampmann (2018) for a review). For optimal activation, gRNAs are usually designed to target the region directly upstream of the gene TSS.
crisprDesign provides functionalities to be able to take
into account design rules that are specific to CRISPRa and CRISPRi
applications. The queryTss function allows to specify
genomic coordinates of promoter regions. The
addTssAnnotation annotates gRNAs for known TSSs, and
includes a column named dist_to_tss that indicates the
distance between the TSS position and the PAM site of the gRNA. For
CRISPRi, we recommend targeting the 25-75bp region downstream of the TSS
for optimal inhibition. For CRISPRa, we recommend targeting the region
75-150bp upstream of the TSS for optimal activation; see (Sanson et al. 2018) for more information.
For more information, please see the following two tutorials:
We illustrate the CRISPR base editing (CRISPRbe) functionalities of
crisprDesign by designing and characterizing gRNAs
targeting IQSEC3 using the cytidine base editor BE4max (Koblan et al. 2018).
We first load the BE4max BaseEditor object available in
crisprBase:
## Class: BaseEditor
## CRISPR Nuclease name: SpCas9
## Target type: DNA
## Metadata: list of length 2
## PAMs: NGG, NAG, NGA
## Weights: 1, 0.2593, 0.0694
## Spacer length: 20
## PAM side: 3prime
## Distance from PAM: 0
## Prototype protospacers: 5'--SSSSSSSSSSSSSSSSSSSS[NGG]--3', 5'--SSSSSSSSSSSSSSSSSSSS[NAG]--3', 5'--SSSSSSSSSSSSSSSSSSSS[NGA]--3'
## Base editor name: BE4max
## Editing strand: original
## Maximum editing weight: C2T at position -15The editing probabilities of the base editor BE4max are stored in a
matrix where rows correspond to the different nucleotide substitutions,
and columns correspond to the genomic coordinate relative to the PAM
site. The editingWeights function from
crisprBase allows to retrieve those probabilities. One can
see that C to T editing is optimal around 15 nucleotides upstream of the
PAM site for the BE4max base editor:
## -36 -35 -34 -33 -32 -31 -30 -29 -28 -27 -26 -25 -24
## 0.007 0.007 0.008 0.018 0.010 0.020 0.014 0.012 0.023 0.013 0.024 0.022 0.034
## -23 -22 -21 -20 -19 -18 -17 -16 -15 -14 -13 -12 -11
## 0.022 0.021 0.035 0.058 0.162 0.318 0.632 0.903 1.000 0.870 0.620 0.314 0.163
## -10 -9 -8 -7 -6 -5 -4 -3 -2 -1
## 0.100 0.056 0.033 0.019 0.018 0.024 0.017 0.005 0.002 0.001We obtain a GuideSet object using the first exon of the
IQSEC3 gene and retain only the first 2 gRNAs for the sake of time:
gr <- queryTxObject(txObject=grListExample,
featureType="cds",
queryColumn="gene_symbol",
queryValue="IQSEC3")
gs <- findSpacers(gr[1],
bsgenome=bsgenome,
crisprNuclease=BE4max)
gs <- gs[1:2]The function addEditedAlleles finds, characterizes, and
scores predicted edited alleles for each gRNA, for a chosen transcript.
It requires a transcript-specific annotation that can be obtained using
the function getTxInfoDataFrame. Here, we will perform the
analysis using the main isoform of IQSEC3 (transcript id
ENST00000538872).
We first get the transcript table for ENST00000538872,
txid <- "ENST00000538872"
txTable <- getTxInfoDataFrame(tx_id=txid,
txObject=grListExample,
bsgenome=bsgenome)
head(txTable)## DataFrame with 6 rows and 10 columns
## chr pos nuc aa aa_number exon pos_plot
## <character> <numeric> <character> <character> <integer> <integer> <integer>
## 1 chr12 66737 G NA NA NA 1
## 2 chr12 66738 G NA NA NA 2
## 3 chr12 66739 C NA NA NA 3
## 4 chr12 66740 A NA NA NA 4
## 5 chr12 66741 G NA NA NA 5
## 6 chr12 66742 T NA NA NA 6
## pos_mrna pos_cds region
## <integer> <integer> <character>
## 1 NA NA Upstream
## 2 NA NA Upstream
## 3 NA NA Upstream
## 4 NA NA Upstream
## 5 NA NA Upstream
## 6 NA NA Upstreamand then add the edited alleles annotation to the
GuideSet:
editingWindow <- c(-20,-8)
gs <- addEditedAlleles(gs,
baseEditor=BE4max,
txTable=txTable,
editingWindow=editingWindow)## [addEditedAlleles] Obtaining edited alleles at each gRNA target site.## [addEditedAlleles] Adding functional consequences to alleles.The editingWindow argument specifies the window of
editing that we are interested in. When not provided, it uses the
default window provided in the BaseEditor object. Note that
providing large windows can exponentially increase computing time as the
number of possible alleles grows exponentially.Let’s retrieve the edited
alleles for the first gRNA:
It is a DataFrame object that contains useful metadata
information:
## list()The wildtypeAllele reports the unedited nucleotide
sequence of the region specified by the editing window (with respect to
the gRNA PAM site). It is always reported from the 5’ to 3’ direction on
the strand corresponding to the gRNA strand. The start and
end specify the corresponding coordinates on the
transcript.
Let’s look at the edited alleles:
## DNAStringSet object of length 6:
## width seq
## [1] 13 CGCGTATTGGATT
## [2] 13 CGCGTATCGGATT
## [3] 13 CGTGTATTGGATT
## [4] 13 CGTGTATCGGATT
## [5] 13 CGCGTACTGGATT
## [6] 13 CGCGCATTGGATTThe DataFrame is ordered so that the top predicted
alleles (based on the score column) are shown first. The
score represents the likelihood of the edited allele to
occur relative to all possible edited alleles, and is calculated using
the editing weights stored in the BE4max object. The
seq column represents the edited nucleotide sequences.
Similar to the wildtypeAllele above, they are always
reported from the 5’ to 3’ direction on the strand corresponding to the
gRNA strand. The variant column indicates the functional
consequence of the editing event (silent, nonsense or missense
mutation). In case an edited allele leads to multiple editing events,
the most detrimental mutation (nonsense over missense, missense over
silent) is reported. The aa column reports the result
edited amino acid sequence.
Note that several gRNA-level aggregate scores have also been added to
the GuideSet object when calling
addEditedAlleles:
## GuideSet object with 2 ranges and 11 metadata columns:
## seqnames ranges strand | protospacer pam
## <Rle> <IRanges> <Rle> | <DNAStringSet> <DNAStringSet>
## spacer_1 chr12 66893 - | CGCGCACCGGATTCTCCAGC AGG
## spacer_2 chr12 66896 + | GGGCGGCATGGAGAGCCTGC TGG
## pam_site cut_site region
## <numeric> <numeric> <character>
## spacer_1 66893 66896 region_1
## spacer_2 66896 66893 region_1
## editedAlleles
## <list>
## spacer_1 CGCGTATTGGATT:0.247151:missense:...,CGCGTATCGGATT:0.161844:missense:...,CGTGTATTGGATT:0.105779:missense:...,...
## spacer_2 GGGTGGTATGGAG:0.4644396:silent:...,GGGCGGTATGGAG:0.2976235:silent:...,GGGTGGCATGGAG:0.0699329:silent:...,...
## score_missense score_nonsense score_silent maxVariant
## <numeric> <numeric> <numeric> <character>
## spacer_1 0.9020188 0 0.0745221 missense
## spacer_2 0.0036734 0 0.9514897 silent
## maxVariantScore
## <numeric>
## spacer_1 0.902019
## spacer_2 0.951490
## -------
## seqinfo: 711 sequences (1 circular) from hg38 genome
## crisprNuclease: SpCas9The score_missense, score_nonsense and
score_silent columns represent aggregated scores for each
of the mutation type. They were obtained by summing adding up all scores
for a given mutation type across the set of edited alleles for a given
gRNA. The maxVariant column indicates the most likely to
occur mutation type for a given gRNA, and is based on the maximum
aggregated score, which is stored in maxVariantScore. For
instance, for spacer_1, the higher score is the
score_missense, and therefore maxVariant is
set to missense.
For more information, please see the following tutorial:
It is also possible to design gRNAs for RNA-targeting nucleases using
crisprDesign. In contrast to DNA-targeting nucleases, the
target spacer is composed of mRNA sequences instead of DNA genomic
sequences.
We illustrate the functionalities of crisprDesign for
RNA-targeting nucleases by designing gRNAs targeting IQSEC3 using the
CasRx (RfxCas13d) nuclease (Konermann et al.
2018).
We first load the CasRx CrisprNuclease object from
crisprBase:
## Class: CrisprNuclease
## Name: CasRx
## Target type: RNA
## Metadata: list of length 2
## PFS: N
## Weights: 1
## Spacer length: 23
## PFS side: 3prime
## Distance from PFS: 0
## Prototype protospacers: 5'--SSSSSSSSSSSSSSSSSSSSSSS[N]--3'The PFS sequence (the equivalent of a PAM sequence for RNA-targeting
nucleases) for CasRx is N, meaning that there is no
specific PFS sequences preferred by CasRx.
We will now design CasRx gRNAs for the transcript ENST00000538872 of IQSEC3.
Let’s first extract all mRNA sequences for IQSEC3:
txid <- c("ENST00000538872","ENST00000382841")
mrnas <- getMrnaSequences(txid=txid,
bsgenome=bsgenome,
txObject=grListExample)
mrnas## DNAStringSet object of length 2:
## width seq names
## [1] 2701 AAGCCCCTCCCCTTCTCTGGGCC...AAAGTTACTGCTAGCATGGGTAA ENST00000382841
## [2] 7087 AGGCTGGGCCGGTGGGAGAGGGA...TTATATTGAAAGATGTCACTTGA ENST00000538872We can use the usual function findSpacers to design
gRNAs, and we only consider a random subset of 100 gRNAs for the sake of
time:
## GuideSet object with 6 ranges and 5 metadata columns:
## seqnames ranges strand | protospacer
## <Rle> <IRanges> <Rle> | <DNAStringSet>
## spacer_1000 region_1 1023 + | TTGACCTAAAGAATAAACAGATT
## spacer_1001 region_1 1024 + | TGACCTAAAGAATAAACAGATTG
## spacer_1002 region_1 1025 + | GACCTAAAGAATAAACAGATTGA
## spacer_1003 region_1 1026 + | ACCTAAAGAATAAACAGATTGAA
## spacer_1004 region_1 1027 + | CCTAAAGAATAAACAGATTGAAA
## spacer_1005 region_1 1028 + | CTAAAGAATAAACAGATTGAAAT
## pam pam_site cut_site region
## <DNAStringSet> <numeric> <numeric> <character>
## spacer_1000 G 1023 NA region_1
## spacer_1001 A 1024 NA region_1
## spacer_1002 A 1025 NA region_1
## spacer_1003 A 1026 NA region_1
## spacer_1004 T 1027 NA region_1
## spacer_1005 G 1028 NA region_1
## -------
## seqinfo: 1 sequence from custom genome
## crisprNuclease: CasRxNote that all protospacer sequences are located on the original strand of the mRNA sequence. For RNA-targeting nucleases, the spacer and protospacer sequences are the reverse complement of each other:
## DNAStringSet object of length 6:
## width seq names
## [1] 23 AATCTGTTTATTCTTTAGGTCAA spacer_1000
## [2] 23 CAATCTGTTTATTCTTTAGGTCA spacer_1001
## [3] 23 TCAATCTGTTTATTCTTTAGGTC spacer_1002
## [4] 23 TTCAATCTGTTTATTCTTTAGGT spacer_1003
## [5] 23 TTTCAATCTGTTTATTCTTTAGG spacer_1004
## [6] 23 ATTTCAATCTGTTTATTCTTTAG spacer_1005## DNAStringSet object of length 6:
## width seq names
## [1] 23 TTGACCTAAAGAATAAACAGATT spacer_1000
## [2] 23 TGACCTAAAGAATAAACAGATTG spacer_1001
## [3] 23 GACCTAAAGAATAAACAGATTGA spacer_1002
## [4] 23 ACCTAAAGAATAAACAGATTGAA spacer_1003
## [5] 23 CCTAAAGAATAAACAGATTGAAA spacer_1004
## [6] 23 CTAAAGAATAAACAGATTGAAAT spacer_1005The addSpacerAlignments can be used to perform an
off-target search across all mRNA sequences using the argument
custom_seq. Here, for the sake of time, we only perform an
off-target search to the 2 isoforms of IQSEC3 specified by the
mRNAs object:
gs <- addSpacerAlignments(gs,
aligner="biostrings",
txObject=grListExample,
n_mismatches=1,
custom_seq=mrnas)
tail(gs)## GuideSet object with 6 ranges and 10 metadata columns:
## seqnames ranges strand | protospacer
## <Rle> <IRanges> <Rle> | <DNAStringSet>
## spacer_1095 region_1 1118 + | CGCCAATACCAGCTCAGCAAGAA
## spacer_1096 region_1 1119 + | GCCAATACCAGCTCAGCAAGAAC
## spacer_1097 region_1 1120 + | CCAATACCAGCTCAGCAAGAACT
## spacer_1098 region_1 1121 + | CAATACCAGCTCAGCAAGAACTT
## spacer_1099 region_1 1122 + | AATACCAGCTCAGCAAGAACTTC
## spacer_1100 region_1 1123 + | ATACCAGCTCAGCAAGAACTTCG
## pam pam_site cut_site region n0_tx
## <DNAStringSet> <numeric> <numeric> <character> <numeric>
## spacer_1095 C 1118 NA region_1 2
## spacer_1096 T 1119 NA region_1 2
## spacer_1097 T 1120 NA region_1 2
## spacer_1098 C 1121 NA region_1 2
## spacer_1099 G 1122 NA region_1 2
## spacer_1100 A 1123 NA region_1 2
## n1_tx n0_gene n1_gene
## <numeric> <numeric> <numeric>
## spacer_1095 0 1 0
## spacer_1096 0 1 0
## spacer_1097 0 1 0
## spacer_1098 0 1 0
## spacer_1099 0 1 0
## spacer_1100 0 1 0
## alignments
## <GRangesList>
## spacer_1095 ENST00000382841:505:+,ENST00000538872:1118:+
## spacer_1096 ENST00000382841:506:+,ENST00000538872:1119:+
## spacer_1097 ENST00000382841:507:+,ENST00000538872:1120:+
## spacer_1098 ENST00000382841:508:+,ENST00000538872:1121:+
## spacer_1099 ENST00000382841:509:+,ENST00000538872:1122:+
## spacer_1100 ENST00000382841:510:+,ENST00000538872:1123:+
## -------
## seqinfo: 1 sequence from custom genome
## crisprNuclease: CasRxThe columns n0_gene and n0_tx report the
number of on-targets at the gene- and transcript-level, respectively.
For instance, spacer_1095 maps to the two isoforms of
IQSEC3 has n0_tx is equal to 2:
## GRanges object with 2 ranges and 9 metadata columns:
## seqnames ranges strand | spacer
## <Rle> <IRanges> <Rle> | <character>
## spacer_1095 ENST00000382841 505 + | TTCTTGCTGAGCTGGTATTG..
## spacer_1095 ENST00000538872 1118 + | TTCTTGCTGAGCTGGTATTG..
## protospacer pam pam_site n_mismatches
## <DNAStringSet> <DNAStringSet> <numeric> <numeric>
## spacer_1095 CGCCAATACCAGCTCAGCAAGAA C 505 0
## spacer_1095 CGCCAATACCAGCTCAGCAAGAA C 1118 0
## canonical cut_site gene_id gene_symbol
## <logical> <numeric> <character> <character>
## spacer_1095 TRUE NA ENSG00000120645 IQSEC3
## spacer_1095 TRUE NA ENSG00000120645 IQSEC3
## -------
## seqinfo: 2 sequences from custom genomeNote that one can also use the bowtie aligner to perform
an off-target search to a set of mRNA sequences. This requires building
a transcriptome bowtie index first instead of building a genome index.
See the crisprBowtie vignette for more detail.
For more information, please see the following tutorial:
Optical pooled screening (OPS) combines image-based sequencing (in situ sequencing) of gRNAs and optical phenotyping on the same physical wells (Feldman et al. 2019). In such experiments, gRNA spacer sequences are partially sequenced from the 5 prime end. From a gRNA design perspective, additional gRNA design constraints are needed to ensure sufficient dissimilarity of the truncated spacer sequences. The length of the truncated sequences, which corresponds to the number of sequencing cycles, is fixed and chosen by the experimentalist.
To illustrate the functionalities of crisprDesign for
designing OPS libraries, we use the guideSetExample. We
will design an OPS library with 8 cycles.
We add the 8nt OPS barcodes to the GuideSet using the
addOpsBarcodes function:
data(guideSetExample, package="crisprDesign")
guideSetExample <- addOpsBarcodes(guideSetExample,
n_cycles=n_cycles)
head(guideSetExample$opsBarcode)## DNAStringSet object of length 6:
## width seq names
## [1] 8 CGCGCACC spacer_1
## [2] 8 GGGCGGCA spacer_2
## [3] 8 GGAGAGCC spacer_3
## [4] 8 AGGTAGAG spacer_4
## [5] 8 GAGCTCCT spacer_5
## [6] 8 CGATGGCC spacer_6The function getBarcodeDistanceMatrix calculates the
nucleotide distance between a set of query barcodes and a set of target
barcodes. The type of distance (hamming or levenshtein) can be specified
using the dist_method argument. The Hamming distance
(default) only considers substitutions when calculating distances, while
the Levenshtein distance allows insertions and deletions.
When the argument binnarize is set to
FALSE, the return object is a matrix of pairwise distances
between query and target barcodes:
barcodes <- guideSetExample$opsBarcode
dist <- getBarcodeDistanceMatrix(barcodes[1:5],
barcodes[6:10],
binnarize=FALSE)
print(dist)## 5 x 5 sparse Matrix of class "dgCMatrix"
## CGATGGCC GCGCGCCG GCTCTACC GCTCTGCT GGGTGTGG
## CGCGCACC 4 7 5 7 7
## GGGCGGCA 4 3 5 4 4
## GGAGAGCC 3 6 5 5 6
## AGGTAGAG 5 6 8 7 4
## GAGCTCCT 7 3 4 3 6When binnarize is set to TRUE (default),
the matrix of distances is binnarized so that 1 indicates similar
barcodes, and 0 indicates dissimilar barcodes. The
min_dist_edit argument specifies the minimal distance
between two barcodes to be considered dissimilar:
dist <- getBarcodeDistanceMatrix(barcodes[1:5],
barcodes[6:10],
binnarize=TRUE,
min_dist_edit=4)
print(dist)## 5 x 5 sparse Matrix of class "dtCMatrix"
## CGATGGCC GCGCGCCG GCTCTACC GCTCTGCT GGGTGTGG
## CGCGCACC . . . . .
## GGGCGGCA . 1 . . .
## GGAGAGCC 1 . . . .
## AGGTAGAG . . . . .
## GAGCTCCT . 1 . 1 .The designOpsLibrary allows users to perform a complete
end-to-end library design; see ?designOpsLibrary for
documentation.
For more information, please see the following tutorial:
The findSpacerPairs function in
crisprDesign enables the design of pairs of gRNAs and works
similar to findSpacers. As an example, we will design
candidate pairs of gRNAs that target a small locus located on chr12 in
the human genome:
library(GenomicRanges)
library(BSgenome.Hsapiens.UCSC.hg38)
library(crisprBase)
bsgenome <- BSgenome.Hsapiens.UCSC.hg38We first specify the genomic locus:
and find all pairs using the function
findSpacerPairs:
The first and second arguments of the function specify the which
genomic region the first and second gRNA should target, respectively. In
our case, we are targeting the same region with both gRNAs. The other
arguments of the function are similar to the findSpacers
function described below.
The output object is a PairedGuideSet, which can be
thought of a list of two GuideSet:
## PairedGuideSet object with 2626 pairs and 4 metadata columns:
## first second | pamOrientation pamDistance
## <GuideSet> <GuideSet> | <character> <numeric>
## [1] chr12:22224025:- chr12:22224033:+ | out 8
## [2] chr12:22224025:- chr12:22224055:- | rev 30
## [3] chr12:22224033:+ chr12:22224055:- | in 22
## [4] chr12:22224025:- chr12:22224056:- | rev 31
## [5] chr12:22224033:+ chr12:22224056:- | in 23
## ... ... ... . ... ...
## [2622] chr12:22224937:- chr12:22224994:+ | out 57
## [2623] chr12:22224938:- chr12:22224994:+ | out 56
## [2624] chr12:22224944:- chr12:22224994:+ | out 50
## [2625] chr12:22224950:+ chr12:22224994:+ | fwd 44
## [2626] chr12:22224958:- chr12:22224994:+ | out 36
## spacerDistance cutLength
## <integer> <numeric>
## [1] -32 2
## [2] 11 30
## [3] 24 28
## [4] 12 31
## [5] 25 29
## ... ... ...
## [2622] 17 51
## [2623] 16 50
## [2624] 10 44
## [2625] 25 44
## [2626] -4 30The first and second GuideSet store information about
gRNAs at position 1 and position 2, respectively. They can be accessed
using the first and second functions:
## GuideSet object with 2626 ranges and 5 metadata columns:
## seqnames ranges strand | protospacer pam
## <Rle> <IRanges> <Rle> | <DNAStringSet> <DNAStringSet>
## spacer_1 chr12 22224025 - | ATTAGTACAACCTTTCTTTT AGG
## spacer_1 chr12 22224025 - | ATTAGTACAACCTTTCTTTT AGG
## spacer_2 chr12 22224033 + | CTTTTGTTTTCCTAAAAGAA AGG
## spacer_1 chr12 22224025 - | ATTAGTACAACCTTTCTTTT AGG
## spacer_2 chr12 22224033 + | CTTTTGTTTTCCTAAAAGAA AGG
## ... ... ... ... . ... ...
## spacer_68 chr12 22224937 - | GGCTGCCAGTCATTGGATCA GGG
## spacer_69 chr12 22224938 - | AGGCTGCCAGTCATTGGATC AGG
## spacer_70 chr12 22224944 - | TTTATAAGGCTGCCAGTCAT TGG
## spacer_71 chr12 22224950 + | GTGAGCCCTGATCCAATGAC TGG
## spacer_72 chr12 22224958 - | CACTGTTTTTTCTTTTTATA AGG
## pam_site cut_site region
## <numeric> <numeric> <character>
## spacer_1 22224025 22224028 region_1
## spacer_1 22224025 22224028 region_1
## spacer_2 22224033 22224030 region_1
## spacer_1 22224025 22224028 region_1
## spacer_2 22224033 22224030 region_1
## ... ... ... ...
## spacer_68 22224937 22224940 region_1
## spacer_69 22224938 22224941 region_1
## spacer_70 22224944 22224947 region_1
## spacer_71 22224950 22224947 region_1
## spacer_72 22224958 22224961 region_1
## -------
## seqinfo: 711 sequences (1 circular) from hg38 genome
## crisprNuclease: SpCas9## GuideSet object with 2626 ranges and 5 metadata columns:
## seqnames ranges strand | protospacer pam
## <Rle> <IRanges> <Rle> | <DNAStringSet> <DNAStringSet>
## spacer_2 chr12 22224033 + | CTTTTGTTTTCCTAAAAGAA AGG
## spacer_3 chr12 22224055 - | TATTCTCATGCACTGCTAGT GGG
## spacer_3 chr12 22224055 - | TATTCTCATGCACTGCTAGT GGG
## spacer_4 chr12 22224056 - | ATATTCTCATGCACTGCTAG TGG
## spacer_4 chr12 22224056 - | ATATTCTCATGCACTGCTAG TGG
## ... ... ... ... . ... ...
## spacer_73 chr12 22224994 + | CAGTGACATAGATCATACAT AGG
## spacer_73 chr12 22224994 + | CAGTGACATAGATCATACAT AGG
## spacer_73 chr12 22224994 + | CAGTGACATAGATCATACAT AGG
## spacer_73 chr12 22224994 + | CAGTGACATAGATCATACAT AGG
## spacer_73 chr12 22224994 + | CAGTGACATAGATCATACAT AGG
## pam_site cut_site region
## <numeric> <numeric> <character>
## spacer_2 22224033 22224030 region_1
## spacer_3 22224055 22224058 region_1
## spacer_3 22224055 22224058 region_1
## spacer_4 22224056 22224059 region_1
## spacer_4 22224056 22224059 region_1
## ... ... ... ...
## spacer_73 22224994 22224991 region_1
## spacer_73 22224994 22224991 region_1
## spacer_73 22224994 22224991 region_1
## spacer_73 22224994 22224991 region_1
## spacer_73 22224994 22224991 region_1
## -------
## seqinfo: 711 sequences (1 circular) from hg38 genome
## crisprNuclease: SpCas9The pamOrientation function returns the PAM orientation
of the pairs:
## [1] "out" "rev" "in" "rev" "in" "rev"and takes 4 different values: in (for PAM-in
configuration) out (for PAM-out configuration),
fwd (both gRNAs target the forward strand) and
rev (both gRNAs target the reverse strand).
The function pamDistance returns the distance between
the PAM sites of the two gRNAs. The function cutLength
returns the distance between the cut sites of the two gRNAs. The
function spacerDistance returns the distance between the
two spacer sequences of the gRNAs.
For more information, please see the following tutorial:
crisprDesign also allows gRNA design for DNA sequences
without genomic context (such as a synthesized DNA construct). See
?findSpacers for more information, and here’s an
example:
seqs <- c(seq1="AGGCGGAGGCCCGACCCGGGCGCGGGGCGGCGC",
seq2="AGGCGGAGGCCCGACCCGGGCGCGGGAAAAAAAGGC")
gs <- findSpacers(seqs)
head(gs)## GuideSet object with 6 ranges and 5 metadata columns:
## seqnames ranges strand | protospacer pam
## <Rle> <IRanges> <Rle> | <DNAStringSet> <DNAStringSet>
## spacer_1 seq1 12 - | CGCCGCCCCGCGCCCGGGTC GGG
## spacer_2 seq1 13 - | GCGCCGCCCCGCGCCCGGGT CGG
## spacer_3 seq1 23 + | GCGGAGGCCCGACCCGGGCG CGG
## spacer_4 seq1 24 + | CGGAGGCCCGACCCGGGCGC GGG
## spacer_5 seq1 25 + | GGAGGCCCGACCCGGGCGCG GGG
## spacer_6 seq1 28 + | GGCCCGACCCGGGCGCGGGG CGG
## pam_site cut_site region
## <numeric> <numeric> <character>
## spacer_1 12 15 seq1
## spacer_2 13 16 seq1
## spacer_3 23 20 seq1
## spacer_4 24 21 seq1
## spacer_5 25 22 seq1
## spacer_6 28 25 seq1
## -------
## seqinfo: 2 sequences from custom genome
## crisprNuclease: SpCas9One can also search for off-targets in a custom sequence as follows:
ontarget <- "AAGACCCGGGCGCGGGGCGGGGG"
offtarget <- "TTGACCCGGGCGCGGGGCGGGGG"
gs <- findSpacers(ontarget)
gs <- addSpacerAlignments(gs,
aligner="biostrings",
n_mismatches=2,
custom_seq=offtarget)
head(alignments(gs))## GRanges object with 1 range and 7 metadata columns:
## seqnames ranges strand | spacer
## <Rle> <IRanges> <Rle> | <DNAStringSet>
## spacer_1 custom_seq1 21 + | AAGACCCGGGCGCGGGGCGG
## protospacer pam pam_site n_mismatches canonical
## <DNAStringSet> <DNAStringSet> <numeric> <numeric> <logical>
## spacer_1 TTGACCCGGGCGCGGGGCGG GGG 21 2 TRUE
## cut_site
## <numeric>
## spacer_1 18
## -------
## seqinfo: 1 sequence from custom genomeFor more information, please see the following tutorial:
It is common to include non-targeting controls (NTCs) in screening experiments. Such controls are usually designed to be random sequences that do not align anywhere in the target genome, and therefore can serve as non-cutting negative controls. Given that they cannot be represented by a genomic range specific to a chromosome, they require a special treatment.
The function addNtcs can be used to add a named vector
of ntcs to a targeting GuideSet object as follows:
gs <- guideSetExample
ntcs <- c(ntc_1="AGTGCTGTGTGTGTGATGCT",
ntc_2="GGGTGCCTTTTTACTCGATG")
gs <- addNtcs(gs, ntcs)
print(gs)## GuideSet object with 1253 ranges and 6 metadata columns:
## seqnames ranges strand | protospacer pam
## <Rle> <IRanges> <Rle> | <DNAStringSet> <DNAStringSet>
## spacer_1 chr12 66893 - | CGCGCACCGGATTCTCCAGC AGG
## spacer_2 chr12 66896 + | GGGCGGCATGGAGAGCCTGC TGG
## spacer_3 chr12 66905 + | GGAGAGCCTGCTGGAGAATC CGG
## spacer_4 chr12 66906 - | AGGTAGAGCACGGCGCGCAC CGG
## spacer_5 chr12 66916 - | GAGCTCCTTGAGGTAGAGCA CGG
## ... ... ... ... . ... ...
## spacer_1249 chr12 175019 + | ACCGCTACTCCAGTGGCTCA AGG
## spacer_1250 chr12 175026 + | CTCCAGTGGCTCAAGGAGCC TGG
## spacer_1251 chr12 175026 - | GGTGGGGCAGAGTCTACACC AGG
## ntc_1 ntc_1 0 * | AGTGCTGTGTGTGTGATGCT NA
## ntc_2 ntc_2 0 * | GGGTGCCTTTTTACTCGATG NA
## pam_site cut_site region opsBarcode
## <numeric> <numeric> <character> <list>
## spacer_1 66893 66896 region_1 C,G,C,...
## spacer_2 66896 66893 region_1 G,G,G,...
## spacer_3 66905 66902 region_1 G,G,A,...
## spacer_4 66906 66909 region_1 A,G,G,...
## spacer_5 66916 66919 region_1 G,A,G,...
## ... ... ... ... ...
## spacer_1249 175019 175016 region_14 A,C,C,...
## spacer_1250 175026 175023 region_14 C,T,C,...
## spacer_1251 175026 175029 region_14 G,G,T,...
## ntc_1 0 NA <NA>
## ntc_2 0 NA <NA>
## -------
## seqinfo: 642 sequences (3 circular) from 2 genomes (hg38, ntc)
## crisprNuclease: SpCas9One can see that additional sequence names (ntc_1 and ntc_2) were
added to the set of human chromosomes to allow NTCs to be represented in
the GuideSet object, with strand set to *, and
pam_site set to 0.
Other design functions can be used as usual on such
GuideSet objects; let’s use
addSequenceFeatures as an example:
## GuideSet object with 1253 ranges and 13 metadata columns:
## seqnames ranges strand | protospacer pam
## <Rle> <IRanges> <Rle> | <DNAStringSet> <DNAStringSet>
## spacer_1 chr12 66893 - | CGCGCACCGGATTCTCCAGC AGG
## spacer_2 chr12 66896 + | GGGCGGCATGGAGAGCCTGC TGG
## spacer_3 chr12 66905 + | GGAGAGCCTGCTGGAGAATC CGG
## spacer_4 chr12 66906 - | AGGTAGAGCACGGCGCGCAC CGG
## spacer_5 chr12 66916 - | GAGCTCCTTGAGGTAGAGCA CGG
## ... ... ... ... . ... ...
## spacer_1249 chr12 175019 + | ACCGCTACTCCAGTGGCTCA AGG
## spacer_1250 chr12 175026 + | CTCCAGTGGCTCAAGGAGCC TGG
## spacer_1251 chr12 175026 - | GGTGGGGCAGAGTCTACACC AGG
## ntc_1 ntc_1 0 * | AGTGCTGTGTGTGTGATGCT NA
## ntc_2 ntc_2 0 * | GGGTGCCTTTTTACTCGATG NA
## pam_site cut_site region opsBarcode percentGC polyA
## <numeric> <numeric> <character> <list> <numeric> <logical>
## spacer_1 66893 66896 region_1 C,G,C,... 70 FALSE
## spacer_2 66896 66893 region_1 G,G,G,... 75 FALSE
## spacer_3 66905 66902 region_1 G,G,A,... 60 FALSE
## spacer_4 66906 66909 region_1 A,G,G,... 70 FALSE
## spacer_5 66916 66919 region_1 G,A,G,... 55 FALSE
## ... ... ... ... ... ... ...
## spacer_1249 175019 175016 region_14 A,C,C,... 60 FALSE
## spacer_1250 175026 175023 region_14 C,T,C,... 65 FALSE
## spacer_1251 175026 175029 region_14 G,G,T,... 65 FALSE
## ntc_1 0 NA <NA> 50 FALSE
## ntc_2 0 NA <NA> 50 FALSE
## polyC polyG polyT startingGGGGG NNGG
## <logical> <logical> <logical> <logical> <character>
## spacer_1 FALSE FALSE FALSE FALSE CAGG
## spacer_2 FALSE FALSE FALSE FALSE CTGG
## spacer_3 FALSE FALSE FALSE FALSE CCGG
## spacer_4 FALSE FALSE FALSE FALSE CCGG
## spacer_5 FALSE FALSE FALSE FALSE ACGG
## ... ... ... ... ... ...
## spacer_1249 FALSE FALSE FALSE FALSE AAGG
## spacer_1250 FALSE FALSE FALSE FALSE CTGG
## spacer_1251 FALSE TRUE FALSE FALSE CAGG
## ntc_1 FALSE FALSE FALSE FALSE <NA>
## ntc_2 FALSE FALSE TRUE FALSE <NA>
## -------
## seqinfo: 642 sequences (3 circular) from 2 genomes (hg38, ntc)
## crisprNuclease: SpCas9## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 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.26.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=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
##
## 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] Rbowtie_1.45.0 BSgenome.Hsapiens.UCSC.hg38_1.4.5
## [3] BSgenome_1.73.1 rtracklayer_1.65.0
## [5] BiocIO_1.17.0 Biostrings_2.75.0
## [7] XVector_0.45.0 GenomicRanges_1.57.2
## [9] GenomeInfoDb_1.41.2 IRanges_2.39.2
## [11] S4Vectors_0.43.2 BiocGenerics_0.53.0
## [13] crisprDesign_1.9.0 crisprBase_1.9.1
## [15] BiocStyle_2.35.0
##
## loaded via a namespace (and not attached):
## [1] DBI_1.2.3 bitops_1.0-9
## [3] httr2_1.0.5 biomaRt_2.63.0
## [5] rlang_1.1.4 magrittr_2.0.3
## [7] matrixStats_1.4.1 compiler_4.4.1
## [9] RSQLite_2.3.7 GenomicFeatures_1.57.1
## [11] dir.expiry_1.13.1 txdbmaker_1.1.2
## [13] png_0.1-8 vctrs_0.6.5
## [15] stringr_1.5.1 pkgconfig_2.0.3
## [17] crayon_1.5.3 fastmap_1.2.0
## [19] dbplyr_2.5.0 utf8_1.2.4
## [21] Rsamtools_2.21.2 rmarkdown_2.28
## [23] tzdb_0.4.0 UCSC.utils_1.1.0
## [25] bit_4.5.0 xfun_0.48
## [27] randomForest_4.7-1.2 zlibbioc_1.51.2
## [29] cachem_1.1.0 jsonlite_1.8.9
## [31] progress_1.2.3 blob_1.2.4
## [33] DelayedArray_0.31.14 BiocParallel_1.39.0
## [35] prettyunits_1.2.0 parallel_4.4.1
## [37] R6_2.5.1 VariantAnnotation_1.51.2
## [39] bslib_0.8.0 stringi_1.8.4
## [41] reticulate_1.39.0 jquerylib_0.1.4
## [43] Rcpp_1.0.13 SummarizedExperiment_1.35.5
## [45] knitr_1.48 readr_2.1.5
## [47] Matrix_1.7-1 tidyselect_1.2.1
## [49] abind_1.4-8 yaml_2.3.10
## [51] codetools_0.2-20 curl_5.2.3
## [53] lattice_0.22-6 tibble_3.2.1
## [55] withr_3.0.2 Biobase_2.67.0
## [57] basilisk.utils_1.19.0 KEGGREST_1.45.1
## [59] evaluate_1.0.1 crisprScoreData_1.9.0
## [61] BiocFileCache_2.15.0 xml2_1.3.6
## [63] ExperimentHub_2.13.1 pillar_1.9.0
## [65] BiocManager_1.30.25 filelock_1.0.3
## [67] MatrixGenerics_1.17.1 crisprScore_1.9.3
## [69] generics_0.1.3 vroom_1.6.5
## [71] RCurl_1.98-1.16 BiocVersion_3.21.1
## [73] hms_1.1.3 glue_1.8.0
## [75] maketools_1.3.1 tools_4.4.1
## [77] AnnotationHub_3.15.0 sys_3.4.3
## [79] crisprBwa_1.9.0 GenomicAlignments_1.41.0
## [81] buildtools_1.0.0 XML_3.99-0.17
## [83] grid_4.4.1 AnnotationDbi_1.69.0
## [85] GenomeInfoDbData_1.2.13 basilisk_1.19.0
## [87] restfulr_0.0.15 cli_3.6.3
## [89] rappdirs_0.3.3 fansi_1.0.6
## [91] S4Arrays_1.5.11 dplyr_1.1.4
## [93] Rbwa_1.9.0 crisprBowtie_1.9.0
## [95] sass_0.4.9 digest_0.6.37
## [97] SparseArray_1.5.45 rjson_0.2.23
## [99] memoise_2.0.1 htmltools_0.5.8.1
## [101] lifecycle_1.0.4 httr_1.4.7
## [103] bit64_4.5.2