[Bioc-devel] VRanges-class positive strandness and locateVariants() strandawareness

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.