Hi, 

I am working on RNA-seq data. I have gotten a VCF file using exactSNP function in Rsubread package, Then I want to use predictCoding function in VariantAnnotation to compute amino acid coding changes, there is an error. Thanks in advance for your help. here is the scripts and the error:

The scripts:

library(BSgenome.Btaurus.UCSC.bosTau6)
library(VariantAnnotation)
library(GenomicFeatures)
vcf <- readVcf(file="/usr/home/qgg/flz/RNA-seq/SNP/SNP.VCF","bosTau6")##SNP.VCF gotten #from   exactSNP function

bosTau6 <- makeTranscriptDbFromUCSC(genome = "bosTau6", tablename = "ensGene")

saveDb(bosTau6, file="/usr/home/qgg/flz/RNA-seq/bosTau6.sqlite")

txdb <- loadDb("/usr/home/qgg/flz/RNA-seq/bosTau6.sqlite")

vcf <- keepSeqlevels(vcf,c(seqlevels(vcf)[1:30]))
txdb <- keepSeqlevels(txdb,c(seqlevels(txdb)[1:30]))
seqlevels(vcf)<-gsub("Chr","chr",seqlevels(vcf))
aa <- predictCoding(vcf,txdb,Btaurus)

The error: 

Error in .Call2("DNAStringSet_translate", x, skip_code, dna_codes[codon_alphabet],  : 
  in 'x[[591]]': not a base at pos 2
Calls: predictCoding ... .predictCodingGRangesList -> translate -> translate -> .Call2 -> .Call

Hi,

Looks like a problem with the 'alt' variable.  predictCoding() calls translate() which expects the alternate alleles to be valid DNA or RNA bases.  See ?translate. The error message says you have an invalid (i.e., non-translatable) value in 'alt' at position 591. You should be able to see the offending character by subsetting the DNAStringSetList,

alt(vcf)[[591]]

fyi, these functions are also useful for checking / inspecting DNA strings:

hasOnlyBaseLetters(unlist(alt(vcf)))

colSums(alphabetFrequency(unlist(alt(vcf))))

Once you've identified the non-base entries you can remove them from the VCF object and re-run predictCoding(). For example,

invalid <- c(591, 600, 802)

sub_vcf <- vcf[!invalid]

Valerie