Read genotype information with VariantAnnotation
1
1
Entering edit mode
4.6 years ago
jeni ▴ 90

Hi,

is there any way to read genotype information (FORMAT field) GT:AD:AF using VariantAnnotation? Moreover, as I have used a somatic caller, I have 2 format fields for each variant (one from normal and one from tumour sample); so, how can I preserve sample name as a column name when I transform this vcf to a dataframe?

As an example:

#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Normal_sample   Tumour_sample
1   52052   .   C   T   .   PASS    CONTQ=93;DP=56 GT:AD:AF 0/0:36,0:0.026  0/1:15,4:0.238

I am reading the vcf with this command:

as.data.frame(cbind(vcf@rowRanges@seqnames,vcf@rowRanges@ranges@start, vcf@fixed, info(vcf)))

So I get the following dataframe:

vcf.rowRanges.ranges.start  REF ALT QUAL    FILTER  CONTQ   DP  
52052    C    T  pass 93   56

I would like to get something like this:

vcf.rowRanges.ranges.start REF ALT QUAL FILTER CONTQ DP NGT NAD NAF TGT TAD TAF
52052    C    T  pass 93   56  0/0  36,0 0.026  0/1 15,4:0 0.238
SNP R vcf Bioconductor VariantAnnotation • 1.4k views
ADD COMMENT
1
Entering edit mode

Please use the formatting bar (especially the code option) to present your post better. You can use backticks for inline code (`text` becomes text), or select a chunk of text and use the highlighted button to format it as a code block. I've done it for you this time.
code_formatting

ADD REPLY
0
Entering edit mode
4.6 years ago

This is pretty well illustrated in the manual... Most standard entries can be accessed using the info() or geno() function call. I don't think you are going to be able to simply cbind all this information together without checking the format of each of the different columns you want to append. For instance AD provides two values - so you would likely have to make separate data.frames before merging.

> your.vcf <- readVcf('your.vcf')
> geno(your.vcf)
List of length 5
names(5): GT AD DP GQ PL
> head(geno(your.vcf)$GT)[[1]]
[1] "0/1"
> geno(your.vcf)$AD[[1]]
[1] 10  2
> info(your.vcf)$AF[[1]]
[1] 0.5
ADD COMMENT

Login before adding your answer.

Traffic: 1964 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6