I have a BAM file and the resulting UMI text file generated from the codes at http://mccb.umassmed.edu/GUIDE-seq. However, the function getUniqueCleavageEvents() runs into an error, detailed below. The sample files on the website work with the function.

uniqueCleavageEvents <- getUniqueCleavageEvents(bamfile, umifile, n.cores.max = 1)

Error: subscript is a logical vector with out-of-bounds TRUE values
Traceback:

1. getUniqueCleavageEvents(alignment.inputfile = bamfile, umi.inputfile = umifile, 
 .     n.cores.max = 1)
2. importBAMAlignments(alignment.inputfile, min.mapping.quality, 
 .     keep.chrM, keep.R1only, keep.R2only, min.R1.mapped, min.R2.mapped, 
 .     apply.both.min.mapped, concordant.strand, same.chromosome, 
 .     max.paired.distance, distance.inter.chrom)
3. as(gal, "GAlignmentPairs")
4. asMethod(object)
5. GAlignmentPairs(ga[first], ga[last], isProperPair = isProperPair, 
 .     names = names(ga)[first])
6. is(first, "GAlignments")
7. ga[first]
8. ga[first]
9. subset_along_ROWS(x, i, drop = drop)
10. extractROWS(x, i)
11. extractROWS(x, i)
12. normalizeSingleBracketSubscript(i, x, as.NSBS = TRUE)
13. NSBS(i, x, exact = exact, strict.upper.bound = !allow.append, 
  .     allow.NAs = allow.NAs)
14. NSBS(i, x, exact = exact, strict.upper.bound = !allow.append, 
  .     allow.NAs = allow.NAs)
15. .subscript_error("subscript is a logical vector with out-of-bounds ", 
  .     "TRUE values")
16. stop(wmsg(...), call. = FALSE)

Hi Tiffany,

Would you be willing to try to narrow this down to a single chromosome of a single sample? I suspect there is just one or a few records causing this. Once you have identified the records/chromosome in a specific sample, it would be great if you could share the subset of data and I will try to reproduce the error and dig into the code.

In addition, could you please also share your R session info with sessionInfo()?

BTW, have you tried to visualize your bam files via IGV just to make sure that the alignment file is proper?

Thanks!

Best regards, Julie

Hi Tiffany,

If you do not feel comfortable sharing your BAM and UMI files, you can try to run the following script and let me know what you get.

bamdata <- GUIDEseq:::importBAMAlignments(file = "yourBAMfile")

Thanks!

Best,

Julie