[BioC] writeVcf bug

Valerie Obenchain vobencha at fhcrc.org
Fri Apr 19 00:03:13 CEST 2013


Hi Richard,

Thanks for reporting this. I've gone back and completely rewritten the 
code that parses/packages the geno data for writeVcf(). I didn't do a 
good job with this the first time around and it's now much cleaner and 
streamlined. Patched versions are 1.6.3 and 1.7.8 and should be 
available via biocLite() tomorrow (Friday) after ~9am PST.

You mentioned you're 'still' having problems with writeVcf(). Are there 
other problems that haven't been addressed (fixed) yet?

Valerie

On 04/18/2013 03:55 AM, Richard Pearson wrote:
> Hi
>
> I've recently upgraded to R 3.0.0, but am still having a (now slightly
> different) problem with writeVcf. When there is only one element in
> geno(vcf), each variant gets a separate row for each sample. For
> example, in the following code, out1.vcf has 5 variants, each with
> genotypes for 3 samples, but out2.vcf has 15 variants, each with
> genotypes for just one sample.
>
>      library(VariantAnnotation)
>      fl <- system.file("extdata", "ex2.vcf", package="VariantAnnotation")
>      in1 <- readVcf(fl, "hg19")
>      writeVcf(in1, "out1.vcf")
>      in2 <- readVcf(fl, "hg19", ScanVcfParam(geno="GT"))
>      writeVcf(in2, "out2.vcf")
>
>  > sessionInfo()
> R version 3.0.0 (2013-04-03)
> Platform: x86_64-unknown-linux-gnu (64-bit)
>
> locale:
>   [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C LC_TIME=en_GB.UTF-8
> LC_COLLATE=en_GB.UTF-8 LC_MONETARY=en_GB.UTF-8 LC_MESSAGES=en_GB.UTF-8
> LC_PAPER=C                 LC_NAME=C
>   [9] LC_ADDRESS=C               LC_TELEPHONE=C
> LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
>
> attached base packages:
> [1] parallel  stats     graphics  grDevices utils     datasets methods base
>
> other attached packages:
> [1] VariantAnnotation_1.6.2 Rsamtools_1.12.1 Biostrings_2.28.0
> GenomicRanges_1.12.2 IRanges_1.18.0          BiocGenerics_0.6.0
> devtools_1.2            BiocInstaller_1.10.0
>
> loaded via a namespace (and not attached):
>   [1] AnnotationDbi_1.22.1   Biobase_2.20.0 biomaRt_2.16.0
> bitops_1.0-5           BSgenome_1.28.0 DBI_0.2-5 digest_0.6.3
> evaluate_0.4.3 GenomicFeatures_1.12.0 httr_0.2
> [11] memoise_0.1            RCurl_1.95-4.1 RSQLite_0.11.3
> rtracklayer_1.20.1     stats4_3.0.0 stringr_0.6.2 tools_3.0.0
> whisker_0.1 XML_3.96-1.1           zlibbioc_1.6.0
>
> Cheers
>
> Richard
>
>
> On 29/11/2012 17:36, Richard Pearson wrote:
>> Hi Michael
>>
>> My devel sessionInfo is below. It seems biocLite has given me version
>> 1.5.17of the package, but I see that .makeVcfGenolooks different in the
>> latest from SVN (1.5.19) - so I'm sure you're right and I just need to
>> wait for things to filter through. Apologies for the noise.
>>
>> Every time I upgrade VariantAnnotation it is an even better package than
>> before - many thanks for all the hard work!
>>
>> Richard
>>
>>
>>   > sessionInfo()
>> R Under development (unstable) (2012-11-27 r61172)
>> Platform: x86_64-unknown-linux-gnu (64-bit)
>>
>> locale:
>>    [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C LC_TIME=en_GB.UTF-8
>> LC_COLLATE=en_GB.UTF-8 LC_MONETARY=en_GB.UTF-8
>> LC_MESSAGES=en_GB.UTF-8 LC_PAPER=C                 LC_NAME=C
>>    [9] LC_ADDRESS=C               LC_TELEPHONE=C
>> LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
>>
>> attached base packages:
>> [1] parallel  stats     graphics  grDevices utils datasets  methods
>> base
>>
>> other attached packages:
>> [1] VariantAnnotation_1.5.17 Rsamtools_1.11.11 Biostrings_2.27.7
>> GenomicRanges_1.11.15 IRanges_1.17.15          BiocGenerics_0.5.1
>> devtools_0.8             BiocInstaller_1.9.4
>>
>> loaded via a namespace (and not attached):
>>    [1] AnnotationDbi_1.21.7   Biobase_2.19.1 biomaRt_2.15.0
>> bitops_1.0-5 BSgenome_1.27.1        DBI_0.2-5 digest_0.6.0
>> evaluate_0.4.2 GenomicFeatures_1.11.5 httr_0.2
>> [11] memoise_0.1            plyr_1.7.1 RCurl_1.95-3
>> RSQLite_0.11.2 rtracklayer_1.19.6     stats4_2.16.0
>> stringr_0.6.1          tools_2.16.0 whisker_0.1            XML_3.95-0.1
>> [21] zlibbioc_1.5.0
>>
>>
>> On 29/11/2012 17:22, Michael Lawrence wrote:
>>> I am pretty sure this bug has been fixed in devel. Would you mind
>>> pasting your devel session info?
>>>
>>> Thanks,
>>> Michael
>>>
>>>
>>> On Thu, Nov 29, 2012 at 8:44 AM, Richard Pearson
>>> <rpearson at well.ox.ac.uk <mailto:rpearson at well.ox.ac.uk>> wrote:
>>>
>>>      Hi
>>>
>>>      I've noticed a bug with writeVcf when trying to write out
>>>      relatively simple vcf files, e.g. a file containing only GT in the
>>>      FORMAT fields. The following is a (hopefully) reproducible example
>>>      based on examples in ?writeVcf
>>>
>>>      library(VariantAnnotation)
>>>      fl <- system.file("extdata", "ex2.vcf",
>>> package="VariantAnnotation")
>>>      out1.vcf <- tempfile()
>>>      out2.vcf <- tempfile()
>>>      in1 <- readVcf(fl, "hg19")
>>>      writeVcf(in1, out1.vcf)
>>>      in2 <- readVcf(fl, "hg19", ScanVcfParam(geno="GT"))
>>>      writeVcf(in2, out2.vcf)
>>>
>>>      This last line gives me:
>>>
>>>      Error in row(genoMat) :
>>>        a matrix-like object is required as argument to 'row'
>>>
>>>      Looking at the following code snippet from .makeVcfGeno, it seems
>>>      genoMat is first getting created as a matrix, but then changed to
>>>      something that is not a matrix by "genoMat <- genoMatFlat". The
>>>      error then comes from the "row(genoMat)" on the last line shown. I
>>>      don't really understand what is going on in this code so can't
>>>      offer a patch I'm afraid.
>>>
>>>          genoMat <- matrix(unlist(as.list(geno), use.names = FALSE,
>>>              recursive = FALSE), nsub * nrec, length(geno))
>>>          genoMatFlat <- as.character(unlist(genoMat))
>>>          genoMatFlat[is.na <http://is.na>(genoMatFlat)] <- "."
>>>          if (is.list(genoMat)) {
>>>              genoMatList <- relist(genoMatFlat,
>>> PartitioningByEnd(genoMat))
>>>              genoMatFlat <- .pasteCollapse(genoMatList, ",")
>>>              genoMat <- matrix(genoMatFlat, nrow(genoMat),
>>> ncol(genoMat))
>>>          }
>>>          else genoMat <- genoMatFlat
>>>          formatMatPerSub <- matrix(rep(t(formatMat), nsub), nsub *
>>>              nrec, length(geno), byrow = TRUE)
>>>          keep <- !is.na <http://is.na>(formatMatPerSub)
>>>          genoListBySub <- seqsplit(genoMat[keep], row(genoMat)[keep])
>>>
>>>      FWIW, as a temporary fix I made an extra, slightly more complex
>>>      FORMAT field using the following:
>>>
>>>      geno(in2)[["DUMMY"]] <-
>>>      matrix(sapply(as.vector(geno(in2)[["GT"]]), function(x) list(c(x,
>>>      x))), ncol=ncol(geno(in2)[["GT"]]),
>>>      dimnames=dimnames(geno(in2)[["GT"]]))
>>>      writeVcf(in2, out2.vcf)
>>>
>>>      I've checked and it seems this bug is also present in latest
>>>      devel. Any chance of a fix?
>>>
>>>      Thanks
>>>
>>>      Richard
>>>
>>>      > sessionInfo()
>>>      R version 2.15.1 (2012-06-22)
>>>      Platform: x86_64-unknown-linux-gnu (64-bit)
>>>
>>>      locale:
>>>       [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C LC_TIME=en_GB.UTF-8
>>>           LC_COLLATE=en_GB.UTF-8 LC_MONETARY=en_GB.UTF-8
>>>       LC_MESSAGES=en_GB.UTF-8 LC_PAPER=C                 LC_NAME=C
>>>       [9] LC_ADDRESS=C               LC_TELEPHONE=C
>>>      LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C
>>>
>>>      attached base packages:
>>>      [1] stats     graphics  grDevices utils     datasets  methods base
>>>
>>>      other attached packages:
>>>      [1] VariantAnnotation_1.4.5 Rsamtools_1.10.2 Biostrings_2.26.2
>>>        GenomicRanges_1.10.5 IRanges_1.16.4          BiocGenerics_0.4.0
>>>      devtools_0.8  BiocInstaller_1.8.3
>>>
>>>      loaded via a namespace (and not attached):
>>>       [1] AnnotationDbi_1.20.3   Biobase_2.18.0 biomaRt_2.14.0
>>>      bitops_1.0-5 BSgenome_1.26.1        DBI_0.2-5 digest_0.6.0
>>>        evaluate_0.4.2 GenomicFeatures_1.10.1 httr_0.2
>>>      [11] memoise_0.1            parallel_2.15.1 plyr_1.7.1
>>>      RCurl_1.95-3 RSQLite_0.11.2         rtracklayer_1.18.1
>>>      stats4_2.15.1          stringr_0.6.1 tools_2.15.1 whisker_0.1
>>>      [21] XML_3.95-0.1           zlibbioc_1.4.0
>>>
>>>      _______________________________________________
>>>      Bioconductor mailing list
>>>      Bioconductor at r-project.org <mailto:Bioconductor at 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]]
>>
>> _______________________________________________
>> Bioconductor mailing list
>> Bioconductor at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/bioconductor
>> Search the archives:
>> http://news.gmane.org/gmane.science.biology.informatics.conductor
>>
>>
>
> _______________________________________________
> Bioconductor mailing list
> Bioconductor at r-project.org
> https://stat.ethz.ch/mailman/listinfo/bioconductor
> Search the archives:
> http://news.gmane.org/gmane.science.biology.informatics.conductor



More information about the Bioconductor mailing list