[Bioc-devel] seqlevelsStyle() warning-message by message-warning, and NCBI/Ensembl seq styles

Robert Castelo robert.castelo at upf.edu
Mon Jan 25 11:39:19 CET 2016


hi,

i would like to ask if current line #142 of GenomeInfoDb/R/seqlevelsStyle.R:

         message("warning! Multiple seqlevels styles found.")

could be replaced by

         warning("Multiple seqlevels styles found.")

since, after all, the message is a warning.

the reason for my request is that a recent update on GenomeInfoDb added 
the 'Ensembl' sequence style and this executes the previous line when i 
try to figure out the sequence style of a BAM file produced from human 
sequence data and GATK. for instance:

path2bam <- file.path(system.file("extdata", 
package="VariantFiltering"), "NA12878.subset.bam")

hdr <- scanBamHeader(path2bam)
names(hdr[[1]]$targets)
  [1] "1"          "2"          "3"          "4"          "5"
  [6] "6"          "7"          "8"          "9"          "10"
[11] "11"         "12"         "13"         "14"         "15"
[16] "16"         "17"         "18"         "19"         "20"
[21] "21"         "22"         "X"          "Y"          "MT"
[26] "GL000207.1" "GL000226.1" "GL000229.1" "GL000231.1" "GL000210.1"
[31] "GL000239.1" "GL000235.1" "GL000201.1" "GL000247.1" "GL000245.1"
[36] "GL000197.1" "GL000203.1" "GL000246.1" "GL000249.1" "GL000196.1"
[41] "GL000248.1" "GL000244.1" "GL000238.1" "GL000202.1" "GL000234.1"
[46] "GL000232.1" "GL000206.1" "GL000240.1" "GL000236.1" "GL000241.1"
[51] "GL000243.1" "GL000242.1" "GL000230.1" "GL000237.1" "GL000233.1"
[56] "GL000204.1" "GL000198.1" "GL000208.1" "GL000191.1" "GL000227.1"
[61] "GL000228.1" "GL000214.1" "GL000221.1" "GL000209.1" "GL000218.1"
[66] "GL000220.1" "GL000213.1" "GL000211.1" "GL000199.1" "GL000217.1"
[71] "GL000216.1" "GL000215.1" "GL000205.1" "GL000219.1" "GL000224.1"
[76] "GL000223.1" "GL000195.1" "GL000212.1" "GL000222.1" "GL000200.1"
[81] "GL000193.1" "GL000194.1" "GL000225.1" "GL000192.1" "NC_007605"
[86] "hs37d5"

seqlevelsStyle(names(hdr[[1]]$targets))
warning! Multiple seqlevels styles found.
[1] "NCBI"    "Ensembl"

this is all fine in interactive mode so that the user is aware that 
he/she may have to take action on  this fact. however, at a specific 
place of my package VariantFiltering i'd like to take action without 
telling anything to the user. for that purpose i would like to use 
suppressWarnings() but does not work currently with the message() call. 
i could use suppressMessages() but would not prefer to.

on a side note, both NCBI and Ensembl sequence style seem to be identical:

genomeStyles()$Homo_sapiens
    circular  auto   sex NCBI  UCSC dbSNP Ensembl
1     FALSE  TRUE FALSE    1  chr1   ch1       1
2     FALSE  TRUE FALSE    2  chr2   ch2       2
3     FALSE  TRUE FALSE    3  chr3   ch3       3
4     FALSE  TRUE FALSE    4  chr4   ch4       4
5     FALSE  TRUE FALSE    5  chr5   ch5       5
6     FALSE  TRUE FALSE    6  chr6   ch6       6
7     FALSE  TRUE FALSE    7  chr7   ch7       7
8     FALSE  TRUE FALSE    8  chr8   ch8       8
9     FALSE  TRUE FALSE    9  chr9   ch9       9
10    FALSE  TRUE FALSE   10 chr10  ch10      10
11    FALSE  TRUE FALSE   11 chr11  ch11      11
12    FALSE  TRUE FALSE   12 chr12  ch12      12
13    FALSE  TRUE FALSE   13 chr13  ch13      13
14    FALSE  TRUE FALSE   14 chr14  ch14      14
15    FALSE  TRUE FALSE   15 chr15  ch15      15
16    FALSE  TRUE FALSE   16 chr16  ch16      16
17    FALSE  TRUE FALSE   17 chr17  ch17      17
18    FALSE  TRUE FALSE   18 chr18  ch18      18
19    FALSE  TRUE FALSE   19 chr19  ch19      19
20    FALSE  TRUE FALSE   20 chr20  ch20      20
21    FALSE  TRUE FALSE   21 chr21  ch21      21
22    FALSE  TRUE FALSE   22 chr22  ch22      22
23    FALSE FALSE  TRUE    X  chrX   chX       X
24    FALSE FALSE  TRUE    Y  chrY   chY       Y
25     TRUE FALSE FALSE   MT  chrM  chMT      MT

this leads to the previous use case with a BAM file, and also to this 
other one:

library(BSgenome.Hsapiens.NCBI.GRCh38)

seqlevelsStyle(BSgenome.Hsapiens.NCBI.GRCh38)
warning! Multiple seqlevels styles found.
[1] "NCBI"    "Ensembl"

note that now if you have a UCSC-style GRanges object like this one:

gr <- GRanges(c("chr1", "chr2"), IRanges(c(10, 20), c(30, 40)))
seqlevelsStyle(gr)
[1] "UCSC"

that you want to use with the BSgenome object, the following simple 
operation will not work anymore:

seqlevelsStyle(gr) <- seqlevelsStyle(BSgenome.Hsapiens.NCBI.GRCh38)
warning! Multiple seqlevels styles found.
Error in mapSeqlevels(x_seqlevels, value, drop = FALSE) :
   the supplied seqname style must be a single string

of course, this has an easy fix:

seqlevelsStyle(gr) <- seqlevelsStyle(BSgenome.Hsapiens.NCBI.GRCh38)[1]
warning! Multiple seqlevels styles found.

but in my view, we loose simplicity in the grammar of moving back and 
forth between sequence styles. one workaround to keep the former simpler 
syntax would be that 'seqlevelsStyle()<-' takes the first value and 
gives a warning instead of prompting an error.

in any case, just out of curiosity, what is the reason to have two 
identical sequence styles under different names?


thanks!

robert.



More information about the Bioc-devel mailing list