Thanks Valerie and Michael,
I’ll implement your suggestions and let you know how it goes.
As for the info section, last time I tried to modify the vcf object directly, that might have been quite different from setting again the whole info section.
As for the overlaps, following your indication I could also use overlapsAny(), which seems to work with two vcf objects as well.
thanks again,
Francesco


On 8 Jan 2014, at 19:09, Valerie Obenchain <vobencha@fhcrc.org<mailto:vobencha@fhcrc.org>> wrote:

Hi Francesco,

On 01/08/2014 05:59 AM, Francesco Lescai wrote:
Hi there,
following this suggestion I’ve converted the results of a number of
analyses into VCF files.
I believe however I never exploited the full potential of the VCF object
and I’m not sure which would be the most efficient way for a targeted
query, i.e. check if a variant from a VCF object is present in any of
the others, like I would to with an SQL database.

Ideally, I’d like to add a metadata field to one of my objects, that
could be reported in the INFO field like any other annotation when I
then write the VCF on a file.
To do this I’d have to add the field both in the info(header(vcf)) and
in the info(vcf), with/after the result of my query.

Two questions therefore:
1) how do you add an INFO field in a way that doesn’t generate error
from within R

> fl <- system.file("extdata", "ex2.vcf", package="VariantAnnotation")
> vcf <- readVcf(fl, "hg19")

Original info variables.

> vcf <- readVcf(fl, "hg19")
> names(info(vcf))
[1] "NS" "DP" "AF" "AA" "DB" "H2"
> dim(info(vcf))
[1] 5 6

Add new variables with the info<- setter.

> newInfo <- info(vcf)
> newInfo$foo <- 1:5
> info(vcf) <- newInfo
> names(info(vcf))
[1] "NS"  "DP"  "AF"  "AA"  "DB"  "H2"  "foo"
> info(vcf)$foo
[1] 1 2 3 4 5



2) what’s the most computationally efficient way to query for
presence/absence of a variant (I thought about a match on the names, a
countOverlaps on the GRanges etc, but that doesn’t exploit the tabix
index, and still seems to me a bit clumsy).

You mention the tabix index which is only applicable to the file on disk and not the VCF class in R. I'm not sure if you are interested in querying the file on disk or class in R so below I've addressed both.

Query vcf file on disk:

To match variant position in the file using a ScanVcfParam with GRanges is best if you want the VCF class as the return object.

Alternatively you could use filterVcf() where the output is a reduced vcf file on disk subset by the filter criteria. filterVcf() can subset by position, character string or specific field values and has options of 'prefilter' and 'filter'. A 'prefilter' searches unparsed lines from the file for a specific character string. 'filters' operate on the data once it has been parsed into the VCF object (this allows 'filters' to be defined in VCF-class syntax).

Query VCF class in R:

To match on position I'd use subsetByOverlaps() on the VCF GRanges if you want a GRanges back. If you only want a count of the hits (i.e., location if present) then use countOverlaps(). This is essentially what you're doing now but if chromosome and/or strand are important then this is the best way to go. Others free to chime in.

You could gain a little performance if the data are large and chromosome and strand are not important. Using countOverlaps() on the ranges() extracted from the GRanges will save a some time on big data.

To match on variant name I would match on the colnames(vcf) as you are doing.

Are any of these operations taking prohibitively long? Please let me know if they are.


Valerie

Begin forwarded message:

From: Michael Lawrence <lawrence.michael@gene.com<mailto:lawrence.michael@gene.com>>
Subject: Re: [BioC] VCF as custom variant annotation database
Date: 8 January 2014 16:56:09 CET
To: Francesco Lescai <francesco.lescai@hum-gen.au.dk<mailto:francesco.lescai@hum-gen.au.dk>>




On Wed, Jan 8, 2014 at 5:59 AM, Francesco Lescai <francesco.lescai@hum-gen.au.dk<mailto:francesco.lescai@hum-gen.au.dk>> wrote:
Hi there,
following this suggestion I’ve converted the results of a number of analyses into VCF files.
I believe however I never exploited the full potential of the VCF object and I’m not sure which would be the most efficient way for a targeted query, i.e. check if a variant from a VCF object is present in any of the others, like I would to with an SQL database.

Ideally, I’d like to add a metadata field to one of my objects, that could be reported in the INFO field like any other annotation when I then write the VCF on a file.
To do this I’d have to add the field both in the info(header(vcf)) and in the info(vcf), with/after the result of my query.

Two questions therefore:
1) how do you add an INFO field in a way that doesn’t generate error from within R


You might want to look into the VRanges class. It's just a GRanges for variants. Add an mcol like with GRanges. By default, the mcols end up in the 'geno' field, but you can change this during coercion to VCF:

vcf <- asVCF(vr, info="myColumn")
writeVcf(vcf, "my.vcf")


2) what’s the most computationally efficient way to query for presence/absence of a variant (I thought about a match on the names, a countOverlaps on the GRanges etc, but that doesn’t exploit the tabix index, and still seems to me a bit clumsy).


Directly from a file? No explicit support for that, but it's an interesting idea. Right now, say you have a VRanges subset to some variants of interest, called "variants", then:

vr <- readVcfAsVRanges("my.vcf", genome=genome(variants)[1], param=variants)
variants %in% vr

I'll make a method for %in%,VRanges,TabixFile to make this easier.

thanks for any advice on how to best use this class,
best
Francesco





thanks for any advice on how to best use this class,
best
Francesco



On 3 Oct 2013, at 15:53, Vincent Carey <stvjc@channing.harvard.edu<mailto:stvjc@channing.harvard.edu>
<mailto:stvjc@channing.harvard.edu>> wrote:




On Thu, Oct 3, 2013 at 7:18 AM, Francesco Lescai
<francesco.lescai@hum-gen.au.dk<mailto:francesco.lescai@hum-gen.au.dk>
<mailto:francesco.lescai@hum-gen.au.dk>> wrote:

   Hi guys,
   I would like to store lists of variants (SNPs and INDELs) in a
   convenient format to be used frequently with other bioconductor
   packages, and I was thinking to store them as Annotation databases.
   In this way I could store and query my own list of variants with
   annotations coming from my work and/or previous experiments.


   Most of the examples/vignettes I could find like SQLForge are
   however gene-centric.
   Can I create a custom one, maybe using SNPlocs as template?


Have you considered tabix-indexed VCF?  This can be stored in a
package and Rsamtools/VariantAnnotation facilities can be used for
targeted querying.  The SNPlocs containers are likely to be redesigned
in the near future.

   thanks for any suggestions,
   Francesco



           [[alternative HTML version deleted]]

   _______________________________________________
   Bioconductor mailing list
   Bioconductor@r-project.org<mailto:Bioconductor@r-project.org> <mailto:Bioconductor@r-project.org>
   https://stat.ethz.ch/mailman/listinfo/bioconductor
   Search the archives:
   http://news.gmane.org/gmane.science.biology.informatics.conductor







--
Valerie Obenchain

Program in Computational Biology
Division of Public Health Sciences
Fred Hutchinson Cancer Research Center
1100 Fairview Ave. N, M1-B155
P.O. Box 19024
Seattle, WA 98109-1024

E-mail: vobencha@fhcrc.org<mailto:vobencha@fhcrc.org>
Phone:  (206) 667-3158
Fax:    (206) 667-1319

________________________________

[X]

Francesco Lescai, PhD, EDBT
Associate Professor

Department of Biomedicine, Human Genetics
Building 1243 - Bioinformatics
Wilhelm Meyers Alle 4
8000 Aarhus C, Denmark


Tel: +45 871 68496
francesco.lescai@hum-gen.au.dk<mailto:francesco.lescai@hum-gen.au.dk>

iSEQ
Centre for iSequencing
iPSYCH
Initiative for Integrative Psychiatric Research

www.ipsych.dk<http://www.ipsych.dk>




	[[alternative HTML version deleted]]

