Hi-

I'm trying to figure out the best method for adding custom annotations to
multi-sample VCF files using the VariantAnnotation package.

I have an already functionally-annotated multi-sample VCF file that I would like to
match up with another VCF file containing different annotations (not multi-sample). Where the chr, start, end, ref, and alt fields are identical, I would like to add the annotation
from the custom vcf file to the existing annotations on the multi-file.

So far I've been grappling with pulling in each vcf and coercing them into a data frame (which is really no fun with DNAStringSetList objects), and attempting to manually match up the fields. However, the ALT field can contain a comma-separated list in both fields, which makes things trickier, and I just feel that there has to be a better way. I'd really appreciate any help or feedback; would findoverlap and Vranges work better? 

Here's what I've tried so far: 

Data frame matching: 

##############
## customvcf ##
##############

humie_ref = "hg19"

##read in vcf with custom annotations
custom <- readVcf(custom_db, humie_ref)

##rename chromosomes to match with vcf files
custom <- renameSeqlevels(custom, c("1"="chr1"))

##make data frame of needed fields
cust.df = cbind(chr= seqlevels(custom), as.data.frame(ranges(custom), stringsToFactors=FALSE), 
               ref=as.character(values(rowData(custom))[["REF"]]), alt=values(rowData(custom))[["ALT"]], info(custom))
##convert alts col to character list for matching
cust.df$alt.value = as.character(cust.df$alt.value)

#############
## multi-vcf ##
#############

##read in everything but sample data for speediness
vcf_param = ScanVcfParam(samples=NA)

##read in multi-sample vcf
vcf <- readVcf(multi_vcf, humie_ref, param=vcf_param)

##data frame vcf info 
vcf.df = cbind(chr= seqlevels(vcf), as.data.frame(ranges(vcf), stringsToFactors=FALSE), fixed(vcf), info(vcf))

vcf.df = cbind(chr= seqlevels(vcf), as.data.frame(ranges(vcf), stringsToFactors=FALSE), 
               ref=as.character(values(rowData(vcf))[["REF"]]), alt=values(rowData(vcf))[["ALT"]], info(vcf))

##convert alt alleles to character for matching
vcf.df$alt.value = as.character(vcf.df$alt.value)

###############
## match-vcf ##
###############
## Here's where I get a little stuck because no matter what I try, I end up with things like "AA|T" matching up to "G", etc. 

##subset for all fields except alt
match1 = data.frame(merge(vcf.df, kav.df, by=c("chr", "start", "end", "ref")))

##keep matches with same alt allele
matches = unique(match1[unlist(sapply(match1$ALT.x, grep, match1$ALT.y, fixed=TRUE)),])

Please let me know if you'd also like me to include example output. 

Thanks,
Summer

> sessionInfo()
R version 3.1.2 (2014-10-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] BiocInstaller_1.16.1     data.table_1.9.4         VariantAnnotation_1.12.9 Rsamtools_1.18.2        
 [5] Biostrings_2.34.1        XVector_0.6.0            GenomicRanges_1.18.4     GenomeInfoDb_1.2.4      
 [9] IRanges_2.0.1            S4Vectors_0.4.0          BiocGenerics_0.12.1     

loaded via a namespace (and not attached):
 [1] AnnotationDbi_1.28.1    base64enc_0.1-2         BatchJobs_1.5           BBmisc_1.9             
 [5] Biobase_2.26.0          BiocParallel_1.0.3      biomaRt_2.22.0          bitops_1.0-6           
 [9] brew_1.0-6              BSgenome_1.34.1         checkmate_1.5.1         chron_2.3-45           
[13] codetools_0.2-10        DBI_0.3.1               digest_0.6.8            fail_1.2               
[17] foreach_1.4.2           GenomicAlignments_1.2.1 GenomicFeatures_1.18.3  iterators_1.0.7        
[21] packrat_0.4.3           plyr_1.8.1              Rcpp_0.11.3             RCurl_1.95-4.5         
[25] reshape2_1.4            RSQLite_1.0.0           rtracklayer_1.26.2      sendmailR_1.2-1        
[29] stringr_0.6.2           tools_3.1.2             XML_3.98-1.1            zlibbioc_1.12.0  

Hi Summer,

I think you want one ALT per row in your data frame so the vcf data needs to be expanded. Take a look at ?expand on a VCF object or try reading the data in with ?readVcfAsVRanges. 

Once you've created your merged data frame (match1) compare the ALT fields directly:

match1$ALT.x == match1$ALT.y

This will give you a logical vector that you can subset on.

As a side note, DNAStringSet objects can be compared directly and you don't need to coerce to character for comparison. Here's an example of subsetting a VCF object on matching alleles:

Read in vcf files:

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

Change ALT in vcf2 for sake of example:

> alt(vcf2)
DNAStringSetList of length 5
[[1]] A
[[2]] A
[[3]] G T
[[4]] 
[[5]] G GTCT

> alt(vcf2)[[4]] <- "TT"

> alt(vcf2)
DNAStringSetList of length 5
[[1]] A
[[2]] A
[[3]] G T
[[4]] TT
[[5]] G GTCT

Compares every ALT in the list:

> alt(vcf1) == alt(vcf2)
LogicalList of length 5
[[1]] TRUE
[[2]] TRUE
[[3]] TRUE TRUE
[[4]] FALSE
[[5]] TRUE TRUE

Return TRUE only if all match:

> idx <- all(alt(vcf1) == alt(vcf2))
> idx
[1]  TRUE  TRUE  TRUE FALSE  TRUE

Subset the VCF:

> dim(vcf1)
[1] 5 3
> dim(vcf1[idx, ])
[1] 4 3

Valerie