Hi,
I found a bug in the GOTHiC package. The bug occurs when different seqlevels are present in the read-pairs - even though when making both seqlevels the same.
Please find below a short example of the error.
Best wishes,
Bettina
library(GOTHiC)
library(BSgenome.Dmelanogaster.UCSC.dm3)
# create GRanges object
> sampleName <- "Test"
> paired_reads_1 <- GRanges(seqnames = Rle(c("chr2L", "chr3R", "chrX")),
+ ranges = IRanges(c(2000, 2000, 2000), c(2100, 2100, 2100)),
+ strand = c("+", "-", "+"))
> seqlevels(paired_reads_1) <- seqlevels(BSgenome.Dmelanogaster.UCSC.dm3)
> paired_reads_2 <- GRanges(seqnames = Rle(c("chr2L", "chrX", "chrX")),
+ ranges = IRanges(c(3000, 3000, 3000), c(3100, 3100, 3100)),
+ strand = c("+", "-", "+"))
> seqlevels(paired_reads_2) <- seqlevels(BSgenome.Dmelanogaster.UCSC.dm3)
> filtered <- GRangesList("paired_reads_1" = paired_reads_1, "paired_reads_2" = paired_reads_2)
> outputfilename <- paste(sampleName, "paired_filtered", sep = "_")
> save(filtered, file = outputfilename)
# map reads to fragments
> mapped <- mapReadsToRestrictionSites(filtered,
+ sampleName=sampleName,
+ BSgenomeName="BSgenome.Dmelanogaster.UCSC.dm3",
+ genome=BSgenome.Dmelanogaster.UCSC.dm3,
+ restrictionSite="^GATC",
+ enzyme="DpnII",
+ parallel=FALSE,
+ cores=1)
# run GOTHiC
> binom <- GOTHiC(fileName1, fileName2,
+ sampleName=sampleName,
+ res=10000,
+ BSgenomeName="BSgenome.Dmelanogaster.UCSC.dm3",
+ genome=BSgenome.Dmelanogaster.UCSC.dm3,
+ restrictionSite="^GATC",
+ enzyme="DpnII",
+ cistrans='all',
+ filterdist=10000,
+ DUPLICATETHRESHOLD=1,
+ fileType='BAM', parallel=FALSE, cores=NULL)
Loading paired reads file ...
Loading mapped reads file ...
Computing binomial ...
Error in as.vector(ifelse(df_filtered$chr1 == df_filtered$chr2, abs(df_filtered$locus1 - :
error in evaluating the argument 'x' in selecting a method for function 'as.vector': Error in ifelse(df_filtered$chr1 == df_filtered$chr2, abs(df_filtered$locus1 - :
error in evaluating the argument 'test' in selecting a method for function 'ifelse': Error in Ops.factor(df_filtered$chr1, df_filtered$chr2) :
level sets of factors are different
> traceback()
3: as.vector(ifelse(df_filtered$chr1 == df_filtered$chr2, abs(df_filtered$locus1 -
df_filtered$locus2), Inf))
2: .binomialHiC(interactions, res, sampleName, BSgenomeName, genome,
restrictionSite, enzyme, parallel, cores, cistrans = cistrans,
filterdist = filterdist)
1: GOTHiC(fileName1, fileName2, sampleName = sampleName, res = 10000,
BSgenomeName = "BSgenome.Dmelanogaster.UCSC.dm3", genome = BSgenome.Dmelanogaster.UCSC.dm3,
restrictionSite = "^GATC", enzyme = "DpnII", cistrans = "all",
filterdist = 10000, DUPLICATETHRESHOLD = 1, fileType = "BAM",
parallel = FALSE, cores = NULL)
> sessionInfo()
R version 3.2.0 (2015-04-16)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 7 (wheezy)
locale:
[1] LC_CTYPE=en_GB.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_GB.UTF-8 LC_COLLATE=en_GB.UTF-8
[5] LC_MONETARY=en_GB.UTF-8 LC_MESSAGES=en_GB.UTF-8
[7] LC_PAPER=en_GB.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats4 parallel stats graphics grDevices utils datasets
[8] methods base
other attached packages:
[1] BSgenome.Dmelanogaster.UCSC.dm3_1.4.0 GOTHiC_1.4.0
[3] data.table_1.9.4 BSgenome_1.36.0
[5] rtracklayer_1.28.2 Biostrings_2.36.0
[7] XVector_0.8.0 GenomicRanges_1.20.3
[9] GenomeInfoDb_1.4.0 IRanges_2.2.1
[11] S4Vectors_0.6.0 BiocGenerics_0.14.0
loaded via a namespace (and not attached):
[1] Rcpp_0.11.5 RColorBrewer_1.1-2 futile.logger_1.4.1
[4] plyr_1.8.2 bitops_1.0-6 futile.options_1.0.0
[7] tools_3.2.0 zlibbioc_1.14.0 digest_0.6.8
[10] gtable_0.1.2 lattice_0.20-31 proto_0.3-10
[13] stringr_0.6.2 hwriter_1.3.2 grid_3.2.0
[16] Biobase_2.28.0 XML_3.98-1.1 BiocParallel_1.2.0
[19] latticeExtra_0.6-26 reshape2_1.4.1 ggplot2_1.0.1
[22] lambda.r_1.1.7 Rsamtools_1.20.1 scales_0.2.4
[25] GenomicAlignments_1.4.1 MASS_7.3-40 ShortRead_1.26.0
[28] colorspace_1.2-6 RCurl_1.95-4.6 munsell_0.4.2
[31] chron_2.3-45
