@@ -1,6 +1,6 @@

 Package: crisprDesign

 Title: Comprehensive design of CRISPR gRNAs for nucleases and base editors

-Version: 1.1.6

+Version: 1.1.7

 Authors@R: c(

     person("Jean-Philippe", "Fortin", email = "[email protected]", role = c("aut", "cre")),

     person("Luke", "Hoberecht", email = "[email protected]", role = c("aut"))

@@ -133,6 +133,9 @@ exportMethods(targetOrigin)

 exportMethods(tssAnnotation)

 importClassesFrom(GenomeInfoDb,Seqinfo)

 importClassesFrom(GenomicRanges,GRanges)

+importClassesFrom(GenomicRanges,GRangesList)

+importClassesFrom(IRanges,DataFrameList)

+importClassesFrom(IRanges,IRanges)

 importFrom(AnnotationDbi,select)

 importFrom(BSgenome,getBSgenome)

 importFrom(BSgenome,getSeq)

@@ -353,9 +353,12 @@ setMethod("cutSites", "GuideSet",

     pamSites <- mcols(object)[["pam_site"]]

     nuc <- metadata(object)[["CrisprNuclease"]]

     strand <- as.character(strand(object))

-    out <- getCutSiteFromPamSite(pam_site=pamSites,

-                                 strand=strand,

-                                 nuclease=nuc)

+    ambiguousStrand <- strand == "*"

+    out <- rep(NA, length(object))

+    out[!ambiguousStrand] <- getCutSiteFromPamSite(

+        pam_site=pamSites[!ambiguousStrand],

+        strand=strand[!ambiguousStrand],

+        nuclease=nuc)

     names(out) <- names(object)

     return(out)

 })

@@ -80,10 +80,18 @@ addEditedAlleles <- function(guideSet,

                 "each gRNA target site.")

     }

     alleles <- lapply(seq_along(guideSet), function(guide){

-        .getEditedAllelesPerGuide(gs=guideSet[guide],

-                                  baseEditor=baseEditor,

-                                  editingWindow=editingWindow,

-                                  nMaxAlleles=nMaxAlleles)

+        seqname <- as.character(GenomeInfoDb::seqnames(guideSet[guide]))

+        genome <- GenomeInfoDb::genome(guideSet[guide])

+        genome <- genome[seqname]

+        if (genome == "ntc"){

+            S4Vectors::DataFrame(seq=DNAStringSet(character(0)),

+                                 score=numeric(0))

+        } else {

+            .getEditedAllelesPerGuide(gs=guideSet[guide],

+                                      baseEditor=baseEditor,

+                                      editingWindow=editingWindow,

+                                      nMaxAlleles=nMaxAlleles)

+        }

     })

     if (addFunctionalConsequence){

         if (verbose){

@@ -175,7 +183,6 @@ addEditedAlleles <- function(guideSet,

                              editingWindow)

     nucChanges <- .getPossibleNucChanges(ws)

-

     # Getting gRNA information:

     pamSite <- pamSites(gs)

     strand <- as.character(strand(gs))

@@ -188,7 +195,7 @@ addEditedAlleles <- function(guideSet,

     pos <- seq(editingWindow[1],

                editingWindow[2])

     names(nucs) <- pos

-

+    

     # Getting scores for the edited nucleotides:

     nucsReduced <- nucs[nucs %in% names(nucChanges)]

@@ -423,6 +430,12 @@ addEditedAlleles <- function(guideSet,

 .addFunctionalConsequences <- function(editedAlleles,

                                        txTable

 ){

+    if (nrow(editedAlleles) == 0){

+        editedAlleles$variant <- character(0)

+        editedAlleles$aa <- character(0)

+        return(editedAlleles)

+    }

+    

     if (txTable$chr[[1]]!=metadata(editedAlleles)$chr){

         stop("editedAlleles are not on the same chromosome.")

     }

@@ -519,8 +532,8 @@ addEditedAlleles <- function(guideSet,

 ){

     ws <- editingWeights(baseEditor)

     ws <- .rescaleWeights(ws)

-    ws <- ws[, as.numeric(colnames(ws))>=editingWindow[1],drop=FALSE]

-    ws <- ws[, as.numeric(colnames(ws))<=editingWindow[2],drop=FALSE]

+    ws <- ws[, as.numeric(colnames(ws)) >= editingWindow[1], drop=FALSE]

+    ws <- ws[, as.numeric(colnames(ws)) <= editingWindow[2], drop=FALSE]

     ws <- crisprBase:::.getReducedEditingMatrix(ws)

     ws <- .addWildtypeWeights(ws)

     return(ws)

@@ -30,10 +30,13 @@ setMethod("addEditingSites",

     }

     pamSites <- pamSites(object)

     strand <- as.character(strand(object))

-    editingSite <- getEditingSiteFromPamSite(pam_site=pamSites,

-                                             strand=strand,

-                                             baseEditor=nuc,

-                                             substitution=substitution)

+    ambiguousStrand <- strand == "*"

+    editingSite <- rep(NA, length(object))

+    editingSite[!ambiguousStrand] <- getEditingSiteFromPamSite(

+        pam_site=pamSites[!ambiguousStrand],

+        strand=strand[!ambiguousStrand],

+        baseEditor=nuc,

+        substitution=substitution)

     mcols(object)$editing_site <- editingSite 

     return(object)

 })

@@ -185,6 +185,8 @@ setMethod("addNtcs", "NULL", function(object,

 #' @importFrom S4Vectors mcols mcols<-

 #' @importFrom methods is

+#' @importClassesFrom GenomicRanges GRangesList

+#' @importClassesFrom IRanges DataFrameList

 .mergeNtcGuideSet <- function(object,

                               ntcGuideSet

 ){

@@ -196,15 +198,15 @@ setMethod("addNtcs", "NULL", function(object,

         if (!is.atomic(dataColumn)){

             S4Vectors::mcols(ntcGuideSet)[[i]] <-

                 lapply(seq_along(ntcGuideSet), function(x){

-                    x <- S4Vectors::mcols(object)[[i]][0]

-                    if (methods::is(x, "GRangesList")){

-                        x <- unlist(x)

-                    }

-                    x

+                    unlist(S4Vectors::mcols(object)[[i]][0])

                 })

-            if (is(dataColumn, "GRangesList")){

+            ## set generic list to appropriate type

+            if (methods::is(dataColumn, "GRangesList")){

                 S4Vectors::mcols(ntcGuideSet)[[i]] <-

-                    GRangesList(S4Vectors::mcols(ntcGuideSet)[[i]])

+                    GenomicRanges::GRangesList(S4Vectors::mcols(ntcGuideSet)[[i]])

+            } else if (methods::is(dataColumn, "DFrameList")){

+                S4Vectors::mcols(ntcGuideSet)[[i]] <-

+                    IRanges::DataFrameList(S4Vectors::mcols(ntcGuideSet)[[i]])

             }

         }

     }

@@ -241,11 +241,11 @@ convertToProtospacerGRanges <- function(guideSet){

-#' @title Convert a GuideSet object into a GRanges storing the range of 

-#'     all gRNAs.

+#' @title Convert a GuideSet object into a GRanges containing the range of 

+#'     all targeting gRNAs.

 #' @description  Convert a GuideSet object into a GRanges object containing 

-#'    the minimum and maximum coordinates of all gRNAs. 

+#'    the minimum and maximum coordinates for all targeting gRNAs. 

 #' 

 #' @param guideSet A \linkS4class{GuideSet} object.

 #' @param anchor A character string specifying which gRNA-specific coordinate

@@ -260,18 +260,25 @@ convertToProtospacerGRanges <- function(guideSet){

 #' data(guideSetExample, package="crisprDesign")

 #' gr <- convertToMinMaxGRanges(guideSetExample)

 #' 

-#' @author Jean-Philippe Fortin

+#' @author Jean-Philippe Fortin, Luke Hoberecht

 #' 

-#' @importFrom GenomicRanges GRanges

-#' @importFrom IRanges IRanges

-#' @importFrom GenomeInfoDb seqlevels seqlevels<-

-#' @importFrom GenomeInfoDb seqinfo seqinfo<-

+#' @importClassesFrom GenomicRanges GRanges

+#' @importClassesFrom IRanges IRanges

+#' @importFrom GenomeInfoDb seqlevels seqlevels<- genome dropSeqlevels

+#' @importFrom GenomeInfoDb seqinfo seqinfo<- seqnames

 #' @export

 convertToMinMaxGRanges <- function(guideSet,

                                    anchor=c("cut_site", "pam_site")

 ){

     anchor <- match.arg(anchor)

-    grs <- split(guideSet, f=as.character(seqnames(guideSet)))

+    genomeSeqlevels <- GenomeInfoDb::genome(guideSet)

+    ntc_seqs <- names(genomeSeqlevels)[genomeSeqlevels == "ntc"]

+    if (length(ntc_seqs) > 0){

+        guideSet <- GenomeInfoDb::dropSeqlevels(guideSet,

+                                                ntc_seqs,

+                                                pruning.mode="coarse")

+    }

+    grs <- split(guideSet, f=as.character(GenomeInfoDb::seqnames(guideSet)))

     grs <- lapply(grs, function(gr){

         if (anchor=="cut_site"){

             start <- min(cutSites(gr), na.rm=TRUE)

@@ -281,11 +288,11 @@ convertToMinMaxGRanges <- function(guideSet,

             end   <- max(pamSites(gr), na.rm=TRUE)

         }

-        chr <- as.character(seqnames(gr))[1]

-        out <- GRanges(chr,

-                       IRanges(start=start,end=end))

-        seqlevels(out) <- seqlevels(gr)

-        seqinfo(out) <- seqinfo(gr)

+        chr <- as.character(GenomeInfoDb::seqnames(gr))[1]

+        out <- GenomicRanges::GRanges(chr,

+                                      IRanges::IRanges(start=start,end=end))

+        GenomeInfoDb::seqlevels(out) <- GenomeInfoDb::seqlevels(gr)

+        GenomeInfoDb::seqinfo(out) <- GenomeInfoDb::seqinfo(gr)

         out

     })

     gr <- Reduce(c,grs)

@@ -2,8 +2,8 @@

 % Please edit documentation in R/completeSpacers.R

 \name{convertToMinMaxGRanges}

 \alias{convertToMinMaxGRanges}

-\title{Convert a GuideSet object into a GRanges storing the range of 

-    all gRNAs.}

+\title{Convert a GuideSet object into a GRanges containing the range of 

+    all targeting gRNAs.}

 \usage{

 convertToMinMaxGRanges(guideSet, anchor = c("cut_site", "pam_site"))

 }

@@ -21,7 +21,7 @@ A GRanges object with start and end coordinates

 }

 \description{

 Convert a GuideSet object into a GRanges object containing 

-   the minimum and maximum coordinates of all gRNAs.

+   the minimum and maximum coordinates for all targeting gRNAs.

 }

 \examples{

 data(guideSetExample, package="crisprDesign")

@@ -29,5 +29,5 @@ gr <- convertToMinMaxGRanges(guideSetExample)

 }

 \author{

-Jean-Philippe Fortin

+Jean-Philippe Fortin, Luke Hoberecht

 }

@@ -98,7 +98,11 @@ test_that("added ntcs have NA or empty list annotations", {

     out_ntcs <- out[intersect(names(out), names(all_ntcs))]

     lapply(listcols, function(i){

         lapply(seq_along(out_ntcs), function(ii){

-            expect_equal(length(mcols(out_ntcs)[[i]][[ii]]), 0)

+            if (is(mcols(out_ntcs)[[i]], "GRangesList")){

+                expect_equal(length(mcols(out_ntcs)[[i]][[ii]]), 0)

+            } else {

+                expect_equal(nrow(mcols(out_ntcs)[[i]][[ii]]), 0)

+            }

         })

     })

     lapply(atomicCols, function(x){

@@ -129,9 +133,8 @@ test_that("ntcs names must be unique and distinct from object ids, seqnames", {

 test_that("crisprDesign computational functions handle ntcs in GuideSet", {

     out <- addNtcs(guideSetExample, all_ntcs)

-    # getPAMSequence(seqnames(out), pamSites(out), strand(out)) # error

-    # getSpacerSequence(seqnames(out), pamSites(out), strand(out)) # error

-    # convertToMinMaxGRanges(out) # issue with strand

+    expect_error(convertToMinMaxGRanges(out),

+                 regex=NA)

     expect_error(convertToProtospacerGRanges(out),

                  regexp=NA) # no error, but meaningless output for ntcs

 })

@@ -147,8 +150,63 @@ out <- addNtcs(head(guideSetExample), all_ntcs[1])

 out_full <- addNtcs(head(guideSetExampleFullAnnotation), all_ntcs[1])

+## ACCESSOR FUNCTIONS =========================================================

+

+test_that("atomic accessor functions handle ntcs in GuideSet gracefully", {

+    expect_error(spacers(out_full), regexp=NA)

+    expect_error(pams(out_full), regexp=NA)

+    expect_error(pamSites(out_full), regexp=NA)

+    expect_error(cutSites(out_full), regexp=NA)

+    expect_error(protospacers(out_full), regexp=NA)

+})

+

+test_that("snps accessor handles ntcs in GuideSet gracefully", {

+    expect_error(snps(out_full),

+                 regexp=NA)

+    expect_false(is.null(snps(out_full)))

+})

+

+test_that("alignments accessors handles ntcs in GuideSet gracefully", {

+    expect_error(alignments(out_full),

+                 regexp=NA)

+    expect_error(onTargets(out_full),

+                 regexp=NA)

+    expect_error(offTargets(out_full),

+                 regexp=NA)

+    expect_false(is.null(alignments(out_full)))

+    expect_false(is.null(onTargets(out_full)))

+    expect_false(is.null(offTargets(out_full)))

+})

+

+test_that("geneAnnotation accessor handles ntcs in GuideSet gracefully", {

+    expect_error(geneAnnotation(out_full),

+                 regexp=NA)

+    expect_false(is.null(geneAnnotation(out_full)))

+})

+

+test_that("tssAnnotation accessor handles ntcs in GuideSet gracefully", {

+    expect_error(tssAnnotation(out_full),

+                 regexp=NA)

+    expect_false(is.null(tssAnnotation(out_full)))

+})

+

+test_that("enzymeAnnotation accessor handles ntcs in GuideSet gracefully", {

+    expect_error(enzymeAnnotation(out_full),

+                 regexp=NA)

+    expect_false(is.null(enzymeAnnotation(out_full)))

+})

+

+## guideSetExampleFullAnnotation lacks editedAlleles annotation

+# test_that("editedAlleles accessor handles ntcs in GuideSet gracefully", {

+#     expect_error(editedAlleles(out_full),

+#                  regexp=NA)

+#     expect_false(is.null(editedAlleles(out_full)))

+# })

+

+## ANNOTATION FUNCTIONS =======================================================

+

 ## split for each crisprDesign function (so more helpful message if test fails)

 test_that("addCutSites handles ntcs in GuideSet gracefully", {

     expect_error(addCutSites(out),

@@ -167,20 +225,29 @@ test_that("addSNPAnnotation handles ntcs in GuideSet gracefully", {

     VCF_PATH <- system.file("extdata",

                             file="common_snps_dbsnp151_example.vcf.gz",

                             package="crisprDesign")

-    expect_error(addSNPAnnotation(out, vcf=VCF_PATH),

+    expect_error(res <- addSNPAnnotation(out, vcf=VCF_PATH),

+                 regexp=NA)

+    expect_error(snps(res),

                  regexp=NA)

+    expect_false(is.null(snps(res)))

 })

 test_that("addGeneAnnotation handles ntcs in GuideSet gracefully", {

-    expect_error(addGeneAnnotation(out, txObject=txdb_human),

+    expect_error(res <- addGeneAnnotation(out, txObject=txdb_human),

                  regexp=NA)

+    expect_error(geneAnnotation(res),

+                 regexp=NA)

+    expect_false(is.null(geneAnnotation(res)))

 })

 test_that("addTssAnnotation handles ntcs in GuideSet gracefully", {

-    expect_error(addTssAnnotation(out, tssObject=tss_human),

+    expect_error(res <- addTssAnnotation(out, tssObject=tss_human),

+                 regexp=NA)

+    expect_error(tssAnnotation(res),

                  regexp=NA)

+    expect_false(is.null(tssAnnotation(res)))

 })

@@ -197,8 +264,11 @@ test_that("addPamScores handles ntcs in GuideSet gracefully", {

 test_that("addRestrictionEnzymes handles ntcs in GuideSet gracefully", {

-    expect_error(addRestrictionEnzymes(out),

+    expect_error(res <- addRestrictionEnzymes(out),

                  regexp=NA)

+    expect_error(enzymeAnnotation(res),

+                 regexp=NA)

+    expect_false(is.null(enzymeAnnotation(res)))

 })

@@ -211,7 +281,7 @@ test_that("addSpacerAlignments/Iterative handles ntcs in GuideSet gracefully", {

                           force=TRUE,

                           prefix="tempIndex")

     index <- file.path(outdir, "tempIndex")

-    expect_error(addSpacerAlignments(

+    expect_error(res <- addSpacerAlignments(

         out,

         txObject=txdb_human,

         tssObject=tss_human,

@@ -225,12 +295,15 @@ test_that("addSpacerAlignments/Iterative handles ntcs in GuideSet gracefully", {

         aligner_index=index,

         bsgenome=BSgenome.Hsapiens.UCSC.hg38),

                  regexp=NA)

+    expect_error(alignments(res),

+                 regexp=NA)

+    expect_false(is.null(alignments(res)))

 })

 test_that("addOnTargetScores handles ntcs in GuideSet gracefully", {

     expect_error(addOnTargetScores(out, methods=c("deephf")),

-                 regexp=NA) # need to test all methods?

+                 regexp=NA) # need to test all methods

 })

@@ -246,20 +319,99 @@ test_that("GuideSet with ntcs can be converted to data.frame", {

 })

+test_that("addEditingSites handles ntcs in GuideSet gracefully", {

+    gs <- out

+    metadata(gs)$CrisprNuclease <- BE4max

+    expect_error(addEditingSites(gs, "C2T"),

+                 regexp=NA)

+})

+

+

+test_that("addExonTable handles ntcs in the GuideSet gracefully", {

+    outga <- addGeneAnnotation(out, txObject=txdb_human)

+    expect_error(addExonTable(outga,

+                              gene_id="ENSG00000120645",

+                              txObject=txdb_human),

+                 regexp=NA)

+})

+

+

+## uses local file

+# test_that("addConservationScores handles ntcs in the GuideSet gracefully", {

+#     conservationFile <- getConservationFiles("human")

+#     expect_error(addConservationScores(out,

+#                                        conservationFile=conservationFile),

+#                  regexp=NA)

+# })

+

+

+test_that("addDistanceToTss handles ntcs in the GuideSet gracefully", {

+    tss_id <- "ENSG00000120645_P1"

+    expect_error(addDistanceToTss(out_full, tss_id),

+                 regexp=NA)

+})

+

+

+test_that("addTxTable handles ntcs in the GuideSet gracefully", {

+    gene_id <- "ENSG00000120645"

+    expect_error(addTxTable(out_full, gene_id, txdb_human),

+                 regexp=NA)

+})

+

+

+## uses local files

+test_that("addCrispraiScores handles ntcs in the GuideSet gracefully", {

+    # gr <- queryTss(tss_human,

+    #                "gene_symbol",

+    #                "IQSEC3")

+    # gs <- findSpacers(gr,

+    #                   crisprNuclease=SpCas9,

+    #                   bsgenome=BSgenome.Hsapiens.UCSC.hg38)

+    # gs <- addNtcs(head(gs), all_ntcs[1])

+    # chromatinFiles <- "~/crisprIndices/chromatin/hg38"

+    # chromatinFiles <- file.path(chromatinFiles, list.files(chromatinFiles))

+    # names(chromatinFiles) <- c("dnase", "faire", "mnase")

+    # fastaFile <- "~/crisprIndices/genomes/hg38/hg38.fa.gz"

+    # addCrispraiScores(gs,

+    #                   gr=gr,

+    #                   tssObject=tss_human,

+    #                   chromatinFiles=chromatinFiles,

+    #                   fastaFile=fastaFile)

+})

+

+

+test_that("addEditedAlleles handles ntcs in the GuideSet gracefully", {

+    gs <- out

+    metadata(gs)$CrisprNuclease <- BE4max

+    txTable <- getTxInfoDataFrame(tx_id="ENST00000538872",

+                                  txObject=txdb_human,

+                                  bsgenome=BSgenome.Hsapiens.UCSC.hg38)

+    

+    expect_error(res <- addEditedAlleles(gs,

+                                         baseEditor=BE4max,

+                                         txTable=txTable,

+                                         editingWindow=c(-20, -8)),

+                 regexp=NA)

+    expect_error(editedAlleles(res),

+                 regexp=NA)

+    expect_false(is.null(editedAlleles(res)))

+})

+

+test_that("addIsoformAnnotation handles ntcs in the GuideSet gracefully", {

+    tx_id <- "ENST00000538872"

+    expect_error(addIsoformAnnotation(out_full,

+                                      tx_id="ENST00000538872"),

+                 regexp=NA)

+})

-## other crisprDesign functions that may need testing

-# addEditingSites

-# addExonTable_consensusIsoform

-# addConservationScores

-# addDistanceToTss

-# addExonTable

-# addTxTable

-# addCrispraiScores

-# addEditedAlleles

-# addExonTable_allIsoforms

-# addIsoformAnnotation

-# addPfamDomains

+test_that("addPfamDomains handles ntcs in the GuideSet gracefully", {

+    pfamTable <- preparePfamTable(txdb_human,

+                                  mart_dataset="hsapiens_gene_ensembl")

+    expect_error(addPfamDomains(out_full,

+                                pfamTable=pfamTable),

+                 regexp=NA)

+    

+})

-## check that function can also handle PairedGuideSet?