Peter,
Could you please post your future question at the official support site at https://support.bioconductor.org/ for others to benefit/contribute? Thanks!
If you could try to download the dev version of CRISPRseek and see if you encounter the same problem, that would be great! Thanks!
http://www.bioconductor.org/packages/devel/bioc/html/CRISPRseek.html
Best,
Julie
On 10/27/14 1:02 PM, "Peter Waltman" <pwaltman@mail.rockefeller.edu> wrote:
Hi Lihua -
Thanks for putting together this package! During testing of it, I've found that I'm getting an error thrown in the filterOffTargets function.
Specifically, when using the fasta file that I've attached, and calling the offTargetAnalysis function below, I get the following error message:
call to offTargetAnalysis:
unpaired.no_re.results <- offTargetAnalysis( "../test.fa", format="fasta", findgRNAs=TRUE, exportAllgRNAs=c("all", "fasta", "genbank", "no"), findgRNAsWithREcutOnly=FALSE, minREpatternSize=6, overlap.gRNA.positions = c(17, 18), findPairedgRNAOnly=FALSE, min.gap=0, max.gap=20, gRNA.name.prefix="gRNA", PAM.size=3, gRNA.size=20, PAM="NGG", BSgenomeName = Hsapiens, chromToSearch="all", max.mismatch = 3, PAM.pattern = "N[A|G]G$", gRNA.pattern = "", min.score = 0.5, topN = 100, topN.OfftargetTotalScore = 10, annotateExon = TRUE, txdb = TxDb.Hsapiens.UCSC.hg19.knownGene, outputDir="./unpaired.grnas.no_re.max_3", fetchSequence = TRUE, upstream = 200, downstream = 200, orgAnn = org.Hs.egSYMBOL)
output:
...
>>> DONE searching
>>> Finding all hits in sequence chrUn_gl000241 ...
>>> DONE searching
>>> Finding all hits in sequence chrUn_gl000242 ...
>>> DONE searching
>>> Finding all hits in sequence chrUn_gl000243 ...
>>> DONE searching
>>> Finding all hits in sequence chrUn_gl000244 ...
>>> DONE searching
>>> Finding all hits in sequence chrUn_gl000245 ...
>>> DONE searching
>>> Finding all hits in sequence chrUn_gl000246 ...
>>> DONE searching
>>> Finding all hits in sequence chrUn_gl000247 ...
>>> DONE searching
>>> Finding all hits in sequence chrUn_gl000248 ...
>>> DONE searching
>>> Finding all hits in sequence chrUn_gl000249 ...
>>> DONE searching
Building feature vectors for scoring ...
Calculating scores ...
Annotating, filtering and generating reports ...
[1] "55870" "79136" "79136" "7920" "79136" "79136" "7920" "79136" "79136" "7920" "79136" "79136" "7920" "79136"
[15] "79136" "7920" "79136" "79136" "7920" "79136" "79136" "7920"
Error in this.score$symbol[query.ind] = overlapGenes.symbol :
NAs are not allowed in subscripted assignments
>
> traceback()
2: filterOffTarget(scores = scores, outputDir = outputDir, BSgenomeName = BSgenomeName,
fetchSequence = fetchSequence, txdb = txdb, orgAnn = orgAnn,
min.score = min.score, topN = topN, topN.OfftargetTotalScore = topN.OfftargetTotalScore,
upstream = upstream, downstream = downstream, annotateExon = annotateExon,
baseBeforegRNA = baseBeforegRNA, baseAfterPAM = baseAfterPAM,
featureWeightMatrixFile = featureWeightMatrixFile)
1: offTargetAnalysis("../test.fa", format = "fasta", findgRNAs = TRUE,
exportAllgRNAs = c("all", "fasta", "genbank", "no"), findgRNAsWithREcutOnly = FALSE,
minREpatternSize = 6, overlap.gRNA.positions = c(17, 18),
findPairedgRNAOnly = FALSE, min.gap = 0, max.gap = 20, gRNA.name.prefix = "gRNA",
PAM.size = 3, gRNA.size = 20, PAM = "NGG", BSgenomeName = Hsapiens,
chromToSearch = "all", max.mismatch = 3, PAM.pattern = "N[A|G]G$",
gRNA.pattern = "", min.score = 0.5, topN = 100, topN.OfftargetTotalScore = 10,
annotateExon = TRUE, txdb = TxDb.Hsapiens.UCSC.hg19.knownGene,
outputDir = "./unpaired.grnas.no_re.max_3", fetchSequence = TRUE,
upstream = 200, downstream = 200, orgAnn = org.Hs.egSYMBOL)
I've traced the error to the following line (including my debug code):
tryres <- try( this.score$symbol[query.ind] <- overlapGenes.symbol )
if (class(tryres)=="try-error" ) browser()
Although I'm not sure if I have the expertise to be able to determine why query.ind would be set to NA in that case.
Thanks,
Peter Waltman
Peter,
Could you please send me the input sequence file so that I can try to replicate the error? Also could you please let me know the version of CRISPRseek you used? Thanks!
FYI, the most recent version is 1.7.10. http://www.bioconductor.org/packages/devel/bioc/html/CRISPRseek.html
Best regards,
Julie