[BUG] --ref-genome flag does not work

[wenjingk]

./GPSeq-RadiCal-0.0.9/gpseq-radical.R "$metadata" /home/${run_id}/test \
>        --chromosome-wide -b 1e6:1e5,5e5:5e4,1e5:1e4,5e4:5e4,25e3:25e3 \
>        --normalize-by lib \
>        --ref-genome mm10 --chrom-tag 19:X,Y \
>        --threads 6 --export-level 1 \
>        --site-domain universe --site-bed "$cutsite_path" \
>        --mask-bed "$mask_path"
2022-03-23 14:04:50 INFO::This log will be stored at '/home/WW319/test/gpseq-radical.log'.
2022-03-23 14:04:50 INFO::Created output folder '/home/WW319/test'.
2022-03-23 14:04:50 INFO::Exported input parameters to '/home/WW319/test/gpseq-radical.opts.rds'.
2022-03-23 14:04:50 INFO::Parsing metadata from '/home/B319_centrality_GPSeq_meta.tsv'.
2022-03-23 14:04:50 INFO::Storing metadata.
2022-03-23 14:04:50 INFO::Opening UCSC browser session...
2022-03-23 14:05:00 INFO::Querying UCSC for 'mm10' chromosome info...
Loading required package: BiocGenerics

Attaching package: ‘BiocGenerics’

The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs

The following objects are masked from ‘package:base’:

    anyDuplicated, append, as.data.frame, basename, cbind, colnames,
    dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
    grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
    rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which.max, which.min

Loading required package: S4Vectors
Loading required package: stats4

Attaching package: ‘S4Vectors’

The following objects are masked from ‘package:base’:

    expand.grid, I, unname

Error in h(simpleError(msg, call)) :
  error in evaluating the argument 'object' in selecting a method for function 'getTable': Table 'chromInfo' is unavailable
Calls: <Anonymous> ... tableName<- -> tableName<- -> .local -> normArgTable
Execution halted

[ggirelli]

Hi @wenjingk, thanks for reporting this.

If you could share the input to your command I can see if I am able to reproduce this.

From a first look, this seems related to the getTable function of the rtracklayer package. Which version of the package do you have installed?

In the meantime, you can always provide a chromosome size table yourself using the --cinfo-path option instead. As explained in the help page, the script queries UCSC only when you do not provide a chromosome size table. You can retrieve your chromosome size table of interest from UCSC TableBrowser. I would recommend including only the chromosomes you are interested into in the table.

For more details check out the related code:

GPSeq-RadiCal/gpseq-radical.R

Lines 788 to 818 in 4ec245f

	# Read chromosome info bed -----------------------------------------------------
	if (is.na(args$bin_bed)) {
	cinfo = NULL
	if (is.na(args$cinfo_path)) {
	logging::loginfo("Opening UCSC browser session...")
	ucsc = rtracklayer::browserSession("UCSC")
	rtracklayer::genome(ucsc) = args$ref_genome
	logging::loginfo(sprintf(
	"Querying UCSC for '%s' chromosome info...",
	rtracklayer::genome(ucsc)))
	cinfo = data.table::data.table(rtracklayer::getTable(
	rtracklayer::ucscTableQuery(ucsc,
	table="chromInfo")))
	cinfo = cinfo[, .(start=1, end=size), by=chrom]
	} else {
	assert(file.exists(args$cinfo_path),
	sprintf("Cannot find chromosome info bed file '%s'.",
	args$cinfo_path))
	logging::loginfo(sprintf(
	"Reading chromosome info from '%s'.", args$cinfo_path))
	cinfo = data.table::as.data.table(
	rtracklayer::import.bed(args$cinfo_path))
	data.table::setnames(cinfo, "seqnames", "chrom")
	cinfo[, c("width", "strand") := NULL]
	}
	assert(!is.null(cinfo), "Failed to build or retrieve chromosome info.")
	assert(nrow(cinfo[start == 1]) == nrow(cinfo),
	"Chromosome start should be 1.")
	}

[ggirelli]

Hi, I am closing this because of lack of feedback. Just drop a comment below if the issue is not solved yet 😄