On Thu, Jul 19, 2012 at 8:59 PM, Timothy Duff <timothy.duff@ncf.edu> wrote:

> Hi. I am trying to determine, from a filtered set of Ilumina 450k probes,
> which of them occur with a specified frequency in a given population. I was
> referred to the Variant Annotation package by this mailing list. While the
> readVcf function seems to handle small loads nicely, looping over the
> regions of interest seems to cause allocation troubles. R tells me "Realloc
> could not re-allocate memory (0 bytes)" after about 4 iterations. Below is
> the relevent code, and below it the output of sessionInfo(). If anyone
> might sugget some diagnostic measures or an alternate way of doing this I
> would appreciate it. Thanks.
>
> ------
>
> library(VariantAnnotation)
> library(IlluminaHumanMethylation450kprobe)
>
> filename <- "
>
> ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20100804/supporting/EUR.2of4intersection_allele_freq.20100804.genotypes.vcf.gz
> "
>
> load("rsids.Rdata") # a data frame containing probe id, rs id, and
> chromosome
> data(IlluminaHumanMethylation450kprobe)
> colnames(rsids) <- c("Probe_ID", "RS_ID", "CHR")
> m <- merge(IlluminaHumanMethylation450kprobe, rsids, by="Probe_ID")
>
> filename <- "
>
> ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20100804/supporting/EUR.2of4intersection_allele_freq.20100804.genotypes.vcf.gz
> "
>
>
Briefly, at a minimum you should set up your param to be a GRanges with
multiple ranges, instead
of running readVcf one range at a time.  How many ranges, and how to best
use the yield capability, I cannot
comment on at the moment.  There is no reason to do this with record-wise
looping, but some kind of chunking
may be in order.  Have a look at scanVcf which may give you more control
over the consumption as the VCF is traversed.

for (i in 1:length(m$Probe_ID)) {
>   snprange <- readVcf(TabixFile(filename), "hg19",
> param=GRanges(as.character(m$CHR[i]), IRanges(as.integer(m$start[i]),
> as.integer(m$end[i]))))
>   freq <- elementMetadata(info(snprange))["EUR_R2"][1,1]
>   if (is.na(freq) == FALSE & freq < .99 & freq > .01) {
>     m$CpGs[i] <- 1
>     }
>   else {
>     m$CpGs[i] <- 0
>   }
> }
>
>
> ----
>
> > sessionInfo()
> R version 2.15.1 (2012-06-22)
> Platform: x86_64-pc-linux-gnu (64-bit)
>
> locale:
>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
>  [7] LC_PAPER=C                 LC_NAME=C
>  [9] LC_ADDRESS=C               LC_TELEPHONE=C
> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] stats     graphics  grDevices utils     datasets  methods   base
>
> other attached packages:
>  [1] IlluminaHumanMethylation450kprobe_2.0.6
>  [2] AnnotationDbi_1.18.1
>  [3] Biobase_2.16.0
>  [4] BiocInstaller_1.4.7
>  [5] VariantAnnotation_1.2.9
>  [6] Rsamtools_1.8.5
>  [7] Biostrings_2.24.1
>  [8] GenomicRanges_1.8.7
>  [9] IRanges_1.14.4
> [10] BiocGenerics_0.2.0
>
> loaded via a namespace (and not attached):
>  [1] biomaRt_2.2.0         bitops_1.0-4.1        BSgenome_1.24.0
>  [4] DBI_0.2-5             GenomicFeatures_1.8.2 grid_2.15.1
>  [7] lattice_0.20-6        Matrix_1.0-6          RCurl_1.4-3
> [10] RSQLite_0.9-2         rtracklayer_1.16.3    snpStats_1.6.0
> [13] splines_2.15.1        stats4_2.15.1         survival_2.36-14
> [16] tools_2.15.1          XML_3.9-4             zlibbioc_1.2.0
>
>
>
> --
> Tim
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor@r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
> http://news.gmane.org/gmane.science.biology.informatics.conductor
>

	[[alternative HTML version deleted]]

