davetang / bedr Goto Github PK

Hi bedr authors,

I am trying to use your package for selecting tabix-indexed regions and do some additional stuff. However, when I tried to apply the "tabix" function both on my own files and also your provided test files, I got the following error: "tabix: invalid option -- 'B' ".

To me, it seems that you are specifying an option --B for argument "params" in your code tabix.R although tabix does not have this option. I would appreciate any help or hint how to solve this. Thanks!

I am using tabix-0.2.6 version.

Here is the whole input / output of your provided example:

if (check.binary("tabix")) {
query.regions <- c("1:1000-100000", "1:1000000-1100000")
cosmic.vcf.example <- system.file(
"extdata/CosmicCodingMuts_v66_20130725_ex.vcf.gz",
package = "bedr"
)
cosmic.query <- tabix(query.regions, cosmic.vcf.example, check.chr = FALSE)
}

• Checking path for tabix... PASS
  /usr/bin/tabix
  TABIX-QUERY
• Processing input (1):
  CONVERT TO BED
• Checking input type... PASS
  Input is in index format
  VALIDATE REGIONS
• Check if index is a string... PASS
• Check index pattern... PASS
• Check for missing values... PASS
• Check for larger start position... PASS.
• Check if zero based... PASS
• Checking sort order... PASS
• Checking for overlapping 'contiguous' regions... PASS

bash -c 'tabix -B /data/home/zane/R/x86_64-pc-linux-gnu-library/3.4/bedr/extdata/CosmicCodingMuts_v66_20130725_ex.vcf.gz /tmp/Rtmpbthdnq/_2dd44b53421a.bed '
tabix: invalid option -- 'B'

Version: 1.1
Usage: tabix [OPTIONS] [FILE] [REGION [...]]
Options:
-0, --zero-based coordinates are zero-based
-b, --begin INT column number for region start [4]
-c, --comment CHAR skip comment lines starting with CHAR [null]
-e, --end INT column number for region end (if no end, set INT to -b) [5]
-f, --force overwrite existing index without asking
-h, --print-header print also the header lines
-H, --only-header print only the header lines
-l, --list-chroms list chromosome names
-m, --min-shift INT set the minimal interval size to 1<<INT; 0 for the old tabix index [0]
-p, --preset STR gff, bed, sam, vcf, bcf, bam
-r, --reheader FILE replace the header with the content of FILE
-s, --sequence INT column number for sequence names (suppressed by -p) [1]
-S, --skip-lines INT skip first INT lines [0]

Head of region...

Version: 1.1
Usage: tabix [OPTIONS] [FILE] [REGION [...]]
Options:
-0, --zero-based coordinates are zero-based
-b, --begin INT column number for region start [4]
-c, --comment CHAR skip comment lines starting with CHAR [null]
-e, --end INT column number for region end (if no end, set INT to -b) [5]
-f, --force overwrite existing index without asking
-h, --print-header print also the header lines
-H, --only-header print only the header lines
-l, --list-chroms list chromosome names
-m, --min-shift INT set the minimal interval size to 1<<INT; 0 for the old tabix index [0]
-p, --preset STR gff, bed, sam, vcf, bcf, bam
-r, --reheader FILE replace the header with the content of FILE
-s, --sequence INT column number for sequence names (suppressed by -p) [1]
-S, --skip-lines INT skip first INT lines [0]

ERROR: Looks like tabix had a problem
Error in tabix(query.regions, cosmic.vcf.example, check.chr = FALSE) :
In addition: Warning message:
running command 'bash -c 'tabix -B /data/home/zane/R/x86_64-pc-linux-gnu-library/3.4/bedr/extdata/CosmicCodingMuts_v66_20130725_ex.vcf.gz /tmp/Rtmpbthdnq/_2dd44b53421a.bed '' had status 1

sessionInfo()
R version 3.4.1 (2017-06-30)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 8 (jessie)

Matrix products: default
BLAS: /usr/lib/atlas-base/atlas/libblas.so.3.0
LAPACK: /usr/lib/atlas-base/atlas/liblapack.so.3.0

locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=C LC_COLLATE=C
[5] LC_MONETARY=C LC_MESSAGES=C LC_PAPER=C LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=C LC_IDENTIFICATION=C

attached base packages:
[1] stats graphics grDevices utils datasets methods base

other attached packages:
[1] bedr_1.0.3 RMySQL_0.10.11 DBI_0.7

loaded via a namespace (and not attached):
[1] crayon_1.3.2 R.methodsS3_1.7.1 grid_3.4.1 R6_2.2.2
[5] futile.options_1.0.0 magrittr_1.5 VennDiagram_1.6.17 data.table_1.10.4
[9] testthat_1.0.2 R.oo_1.21.0 R.utils_2.5.0 futile.logger_1.4.3
[13] lambda.r_1.1.9 tools_3.4.1 parallel_3.4.1 yaml_2.1.14
[17] compiler_3.4.1

I'm trying to use bedr to walk through the bedtools tutorial to understand how bedtools code translates into bedr and to compare their output.

However I'm having an issue with the main bedtools wrapper function bedr. I'm looking to run a basic intersection of two files cpg.bed and exons.bed provided on the link above (too large to upload here). The basic unix command for this would be:

bedtools intersect -a cpg.bed -b exons.bed | head -5

To emulate this line of code I tried to run the following:

bedr(engine = "bedtools", 
     params = "| head -5",
     input = list(a = read.table("cpg.bed"), b = read.table("exons.bed")),
     method = intersect,
     tmpDir = NULL,
     deleteTmpDir = TRUE,
     outputDir = NULL,
     outputFile = NULL,
     check.chr = TRUE,
     check.zero.based = TRUE,
     check.valid = TRUE,
     check.sort = TRUE,
     check.merge = TRUE,
     verbose = TRUE
)

This produces the following error message:

Error in method == "sort" : comparison (1) is possible only for atomic and list

Help messages on the web suggest this error may be due to R and R Studio versions being out of sync, but I have since updated to the latest versions of each and still get the same error.

It's not clear to me from the help files and vignette provided with bedr how to input files correctly, particularly if using large files i.e. do I need to use the read.table function as shown above? I also tried using R objects to input the files:

cpg <- read.table("cpg.bed")
exons <- read.table("exons.bed")

Then altered the input parameter in the bedr command to:

input = list(a = cpg, b = exons)

But I got the same error.

These files are located in my working directory. I also tried setting PARAMS to NULL, and omitting all superfluous arguments one by one, but the same error is always generated. I'm a bit stumped.

My understanding of what this command does may be flawed, but from my reading of the help files, the first 4 arguments of the bedr function effectively constructs the line of unix code I've written above, and this is run via unix, eventually generating a new R object/or file that contains the intersect information.