1
$\begingroup$

I'm having some problems in exporting the annotations of some variants from a VCF file, after using the "VariantAnnotation" package from bioconductor. Below the script I've used:

#!/usr/bin/env Rscript

library(VariantAnnotation)
library(BSgenome.Scerevisiae.UCSC.sacCer3)
library(TxDb.Scerevisiae.UCSC.sacCer3.sgdGene)

args = commandArgs(trailingOnly=TRUE)

vcf <- readVcf(args[1], "sacCer3")

sc <- BSgenome.Scerevisiae.UCSC.sacCer3
txdb <- TxDb.Scerevisiae.UCSC.sacCer3.sgdGene

coding <- data.frame(predictCoding(vcf, txdb, seqSource=sc))

The coding dataframe looks like this:

  seqnames start end width strand paramRangeID REF     ALT QUAL FILTER
1     chrI   455 455     1      +         <NA>   C       T   40   PASS
2     chrI   455 455     1      +         <NA>   C           40   PASS
3     chrI   455 455     1      +         <NA>   C       G   40   PASS
4     chrI   455 455     1      +         <NA>   C       A   40   PASS
5     chrI   456 456     1      +         <NA>   T       A   40   PASS
6     chrI   456 456     1      +         <NA>   C CCTGCTC   40   PASS
  varAllele CDSLOC.start CDSLOC.end CDSLOC.width PROTEINLOC QUERYID TXID CDSID
1         T          121        121            1         41       1    1     1
2                    121        121            1         41       2    1     1
3         G          121        121            1         41       3    1     1
4         A          121        121            1         41       4    1     1
5         A          122        122            1         41       5    1     1
6   CCTGCTC          122        122            1         41       6    1     1
   GENEID    CONSEQUENCE REFCODON  VARCODON REFAA VARAA
1 YAL069W  nonsynonymous      CTC       TTC     L     F
2 YAL069W not translated      CTC        TC            
3 YAL069W  nonsynonymous      CTC       GTC     L     V
4 YAL069W  nonsynonymous      CTC       ATC     L     I
5 YAL069W  nonsynonymous      CTC       CAC     L     H
6 YAL069W  nonsynonymous      CTC CCCTGCTCC     L   PCS

But if I try to save it to file, for instance using the following command:

write.table(coding, 'test.tsv', sep='\t', row.names=F)

The output table looks like this:

"seqnames"      "start" "end"   "width" "strand"        "paramRangeID"  "REF"   "ALT"   "QUAL"  "FILTER"        "varAllele"     "CDSLOC.start"  "CDSLOC.end"    "CDSLOC.width"  "PROTEINLOC"    "QUERYID"       "TXID"  "CDSID" "GENEID"        "CONSEQUENCE"   "REFCODON"      "VARCODON"      "REFAA" "VARAA"             
"chrI"  455     455     1       "+"     NA      "C"     <S4 object of class "DNAStringSet">     40      "PASS"  "T"     121     121     1       41      1       "1"     1       "YAL069W"       "nonsynonymous" "CTC"   "TTC"   "L"     "F"                            
"chrI"  455     455     1       "+"     NA      "C"     <S4 object of class "DNAStringSet">     40      "PASS"  ""      121     121     1       41      2       "1"     1       "YAL069W"       "not translated" "CTC"   "TC"    ""      ""

Notice the <S4 object of class "DNAStringSet"> bit. Is there any way to force R to give me the same representation as head? I couldn't find anything in Biostrings' documentation.

Improve this question
$\endgroup$
1
  • $\begingroup$ You might be able to convert the DNAStringSet to character and then store it. But it might be a function available to store a GRange to a file I am not aware $\endgroup$
    – llrs
    Jan 22, 2018 at 12:04

1 Answer 1

Reset to default
2
$\begingroup$

Is there any way to force R to give me the same representation as head?

head just cuts off everything after the first six entries. The magic, as it were, happens in the default print method for the specific type. But the representation generated by it is only useful for ad-hoc display, not for storage in a structured data file.

The best solution is therefore to manually convert the column types. Since it’s only one column:

coding$ALT = sapply(coding$ALT, as.character)

I’m a bit puzzled that the REF and ALT column have different types, anyway: that’s probably a bug in the {VariantAnnotation} package that should be reported.

Furthermore, you probably want to use quote = FALSE for the write.table. In general, write.table has terrible defaults. readr::write_tsv is almost always a better choice.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Unfortunately, if I perform as.character(coding$ALT) (not REF, that gets written fine) I get a "<S4 object of class \"DNAStringSet\">" $\endgroup$
    – mgalardini
    Jan 22, 2018 at 12:58
  • $\begingroup$ I removed the cast to data.frame around the predictCoding function and now I get this message when trying to cast to character: Error in as.vector(x, mode = "character") : no method for coercing this S4 class to a vector $\endgroup$
    – mgalardini
    Jan 22, 2018 at 13:07
  • $\begingroup$ @mgalardini That’s bizarre, it works for me. But maybe try loading the {Biostrings} package explicitly — it sounds like the method call is incorrectly dispatched (but it works for me even without attaching the package; as it should). $\endgroup$
    – Konrad Rudolph
    Jan 22, 2018 at 13:44
  • 1
    $\begingroup$ If this continues to not work you can also try to invoke as.character on individual DNAStrings, like so: sapply(coding$REF, as.character). However, if the former doesn’t work I fear the latter won’t either. $\endgroup$
    – Konrad Rudolph
    Jan 22, 2018 at 13:59
  • $\begingroup$ This last suggestion actually worked; I ignore while the others didn't... I'll update your answer, and thanks a lot! $\endgroup$
    – mgalardini
    Jan 22, 2018 at 14:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.