[Bioc-devel] VRanges-class positive strandness and locateVariants() strandawareness
Michael,
regarding our email exchange three weeks ago, I found a couple of places
in VariantAnnotation that IMO need to be updated to avoid enforcing
strand on VRanges.
these places occur when constructing and validating VRanges objects,
concretely at:
1. file R/methods-VRanges-class.R at the VRanges class constructor:
VRanges <-
function(seqnames = Rle(), ranges = IRanges(),
ref = character(), alt = NA_character_,
totalDepth = NA_integer_, refDepth = NA_integer_,
altDepth = NA_integer_, ..., sampleNames = NA_character_,
softFilterMatrix = FilterMatrix(matrix(nrow = length(gr),
ncol = 0L),
FilterRules()),
hardFilters = FilterRules())
{
gr <- GRanges(seqnames, ranges,
strand = .rleRecycleVector("*", length(ranges)), ...)
[...]
that precludes setting the strand at construction time:
library(VariantAnnotation)
VRanges(seqnames="chr1", ranges=IRanges(1, 5), ref="T", alt="C", strand="-")
Error in GRanges(seqnames, ranges, strand = .rleRecycleVector("*",
length(ranges)), :
formal argument "strand" matched by multiple actual arguments
2. R/AllClasses.R at the VRanges class validity function .valid.VRanges():
.valid.VRanges.strand <- function(object) {
if (any(object at strand == "-"))
paste("'strand' must always be '+' or '*'")
}
[...]
.valid.VRanges <- function(object)
{
c(.valid.VRanges.ref(object),
.valid.VRanges.alt(object),
.valid.VRanges.strand(object),
.valid.VRanges.depth(object))
}
that prompts an error when variants annotated on the negative strand are
detected:
library(VariantAnnotation)
example(VRanges)
strand(vr) <- "-"
c(vr)
Error in validObject(.Object) :
invalid class ?VRanges? object: 'strand' must always be '+' or '*'
cheers,
robert.
On 05/22/2015 09:49 PM, Michael Lawrence wrote:
This changed recently. VariantAnnotation in devel no longer enforces a
strand on VRanges, or at least it allows the "*" case.
On Fri, May 22, 2015 at 11:33 AM, Robert Castelo <robert.castelo at upf.edu
<mailto:robert.castelo at upf.edu>> wrote:
Hi,
I have encountered myself in a strange situation when using the
function locateVariants() from VariantAnnotation with an input
VRanges object. The problem is that some of the expected coding
annotations are not showing up when using locateVariants() with
default parameters.
After investigating this situation I think I found the reason, which
does not look like a bug but I would like that you give me some
clarification about the logic behind using locateVariants() with
VRanges objects.
The documentation of the VRanges-class says that in this class of
objects "The strand is always constrained to be positive (+).". I
guess there may be a good reason for this but I could not find it in
the documentation or googling about it.
This means that when you coerce a CollapsedVCF object (obtained, for
example, from a VCF file via readVcf()) to a VRanges object, even
though variants in the VCF may have no strand, they get a positive
strand in the VRanges object.
The problem arises then, when you call locateVariants() with this
VRanges object, because features on the negative strand are never
going to overlap with the variants since, by default, the argument
ignore.strand=FALSE.
Let me illustrate this with a toy example. Consider the SNP
rs1129038
(http://www.ncbi.nlm.nih.gov/projects/SNP/snp_ref.cgi?rs=1129038) at
chr15:28356859 with allels A/G. It is located on the 3' UTR of the
gene HERC2 coded on the negative strand of the human reference
genome. Let's build a toy VRanges object having this variant:
library(VariantAnnotation)
vr <- VRanges(seqnames="chr15",
ranges=IRanges(28356859, 28356859),
ref="A", alt="G",
refDepth=5, altDepth=7,
totalDepth=12, sampleNames="A")
strand(vr)
factor-Rle of length 1 with 1 run
Lengths: 1
Values : +
Levels(3): + - *
Let's build now its CollapsedVCF counterpart by using the
corresponding coercion method and set the strand to "*":
vcf <- asVCF(vr)
strand(vcf) <- "*"
Now run locateVariants() on both objects with UCSC annotations:
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
locateVariants(vcf, txdb, region=AllVariants())
GRanges object with 2 ranges and 9 metadata columns:
seqnames ranges strand | LOCATION LOCSTART
LOCEND QUERYID TXID CDSID
<Rle> <IRanges> <Rle> | <factor> <integer> <integer> <integer>
<character> <IntegerList>
[1] chr15 [28356859, 28356859] * | threeUTR 50 50
1 55386
[2] chr15 [28356859, 28356859] * | threeUTR 50 50
1 55387
GENEID PRECEDEID FOLLOWID
<character> <CharacterList> <CharacterList>
[1] 8924
[2] 8924
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
locateVariants(vr, txdb, region=AllVariants())
GRanges object with 1 range and 9 metadata columns:
seqnames ranges strand | LOCATION LOCSTART
LOCEND QUERYID TXID CDSID
<Rle> <IRanges> <Rle> | <factor> <integer> <integer> <integer>
<integer> <IntegerList>
[1] chr15 [28356859, 28356859] + | intergenic <NA> <NA>
1 <NA>
GENEID PRECEDEID FOLLOWID
<character> <CharacterList> <CharacterList>
[1] <NA> 100132565,100289656,100616223,... 2567
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
Note that while we get the 3' UTR annotation from the strandless VCF
object we do not get it from the VRanges object with the positive
strand. To make my point clear: this positive strand shows up when
you coerce a strandless VCF object to a VRanges one, because
positive strandness seems to be the convention for VRanges objects:
as(vcf, VRanges)
VRanges object with 1 range and 1 metadata column:
seqnames ranges strand ref
alt totalDepth refDepth altDepth
<Rle> <IRanges> <Rle> <character> <characterOrRle> <integerOrRle>
<integerOrRle> <integerOrRle>
[1] chr15 [28356859, 28356859] + A
G 12 5 7
sampleNames softFilterMatrix | QUAL
<factorOrRle> <matrix> | <numeric>
[1] A | <NA>
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
hardFilters: NULL
Of course, if I run locateVariants() with the argument
ignore.strand=TRUE, then I get the expected annotation:
locateVariants(vr, txdb, region=AllVariants(), ignore.strand=TRUE)
GRanges object with 2 ranges and 9 metadata columns:
seqnames ranges strand | LOCATION LOCSTART
LOCEND QUERYID TXID CDSID
<Rle> <IRanges> <Rle> | <factor> <integer> <integer> <integer>
<character> <IntegerList>
[1] chr15 [28356859, 28356859] + | threeUTR 677
677 1 55386
[2] chr15 [28356859, 28356859] + | threeUTR 677
677 1 55387
GENEID PRECEDEID FOLLOWID
<character> <CharacterList> <CharacterList>
[1] 8924
[2] 8924
-------
seqinfo: 1 sequence from an unspecified genome; no seqlengths
So, my question is, given that VRanges objects are enforced to have
a positive strand, would not be better to have ignore.strand=TRUE as
default in locateVariants?
Alternatively, I would recommend that locateVariants() issues a
warning, or maybe an error, when the input object is VRanges and
ignore.strand=FALSE.
Finally, out of curiosity, why a VRanges object enforces the
positive strand in all its genomic ranges? Would not be better just
taking the strand of the CollapsedVCF object when coercing the
CollapsedVCF object to VRanges?
thanks!!
robert.
_______________________________________________
Bioc-devel at r-project.org <mailto:Bioc-devel at r-project.org> mailing list
https://stat.ethz.ch/mailman/listinfo/bioc-devel
Robert Castelo, PhD Associate Professor Dept. of Experimental and Health Sciences Universitat Pompeu Fabra (UPF) Barcelona Biomedical Research Park (PRBB) Dr Aiguader 88 E-08003 Barcelona, Spain telf: +34.933.160.514 fax: +34.933.160.550