plotting AF in vcf files
2
0
Entering edit mode
11 weeks ago
rheab1230 ▴ 30

Hello everyone, I want to plot AF distribution for all my vcf files. Is there any software to do that. I tried using vcfstats but its showing error in my vcf files.

AF vcf plotting • 799 views
ADD COMMENT
0
Entering edit mode

I got AF from vcf files using vcftools software. Now I am trying to plot these values using the command:

plot(density(var_qual$X.FREQ.))
Error in density.default(var_qual$X.FREQ.) : argument 'x' must be numeric

The thing is that in my case I have two frequency value being available at one position.

22 chr
51238249 pos
2 alleles
930  N_CHR
0.910753    0.0892473 freq
ADD REPLY
0
Entering edit mode

Hi, what is the output of str(var_qual)?

I would use bcftools query to output the AF values, and then import these to R where I would plot them via plot(), density(), and/or hist()

ADD REPLY
0
Entering edit mode

This is the output of str(var_qual)

data.frame':    1 obs. of  6 variables:
 $ X.cichlid_subset.: chr "cichlid_subset"
 $ X.CHROM.         : chr "CHROM"
 $ X.POS.           : chr "POS"
 $ X.N_ALLELES.     : chr "N_ALLELES"
 $ X.N_CHR.         : chr "N_CHR"
 $ X..FREQ..        : chr "{FREQ}"
ADD REPLY
2
Entering edit mode
11 weeks ago
sbstevenlee ▴ 240

Here is a Python API solution using the pyvcf.VcfFrame.plot_hist method I wrote.

Below is a simple example:

from fuc import common, pyvcf
common.load_dataset('pyvcf')
vcf_file = '~/fuc-data/pyvcf/normal-tumor.vcf'
vf = pyvcf.VcfFrame.from_file(vcf_file)
vf.plot_hist('DP')

enter image description here

We can draw multiple histograms with hue mapping:

annot_file = '~/fuc-data/pyvcf/normal-tumor-annot.tsv'
af = common.AnnFrame.from_file(annot_file, sample_col='Sample')
vf.plot_hist('DP', af=af, group_col='Tissue')

enter image description here

We can show AF instead of DP:

vf.plot_hist('AF')

enter image description here

ADD COMMENT
0
Entering edit mode

Thank you so much for this. I will definitely try this.

ADD REPLY
1
Entering edit mode

Hello, I am using the command to plot AF distribution for vcf files. But I am getting some error. The command:

#!/bin/python
from fuc import common, pyvcf
common.load_dataset('pyvcf')
vcf_file = 'GEUVADIS.chr22.genotype.vcf'
vf = pyvcf.VcfFrame.from_file(vcf_file)
vf.plot_hist('AF')

the error:

file.python:5: DtypeWarning: Columns (5) have mixed types.Specify dtype option on import or set low_memory=False.
vf = pyvcf.VcfFrame.from_file(vcf_file)
Traceback (most recent call last):
File "file.python", line 6, in
vf.plot_hist('AF')
File "/home/anaconda3/envs/var/lib/python3.7/site-packages/fuc/api/pyvcf.py", line 1991, in plot_hist
df = self.extract(k, as_nan=True, func=d[k])
File "/home/anaconda3/envs/var/lib/python3.7/site-packages/fuc/api/pyvcf.py", line 3912, in extract
df = self.df.apply(one_row, axis=1)
File "/home/anaconda3/envs/var/lib/python3.7/site-packages/pandas/core/frame.py", line 8736, in apply
return op.apply()
File "/home/anaconda3/envs/var/lib/python3.7/site-packages/pandas/core/apply.py", line 688, in apply
return self.apply_standard()
File "/home/anaconda3/envs/var/lib/python3.7/site-packages/pandas/core/apply.py", line 805, in apply_standard
results, res_index = self.apply_series_generator()
File "/home/anaconda3/envs/var/lib/python3.7/site-packages/pandas/core/apply.py", line 821, in apply_series_generator
results[i] = self.f(v)
File "/home//anaconda3/envs/var/lib/python3.7/site-packages/fuc/api/pyvcf.py", line 3901, in one_row
i = r.FORMAT.split(':').index(k)
ValueError: 'AF' is not in list
ADD REPLY
1
Entering edit mode

I responded to your GitHub issue here: https://github.com/sbslee/fuc/issues/39

ADD REPLY
0
Entering edit mode

Yes, I saw that. Thank you.

ADD REPLY
1
Entering edit mode
11 weeks ago

I want to plot AF distribution for all my vcf files. Is there any software to do that.

 bcftools stats --af-bins 0.1,0.2,0.3,0.5,1.0 input.vcf

followed by

plot-vcfstats

http://samtools.github.io/bcftools/bcftools.html#plot-vcfstats

ADD COMMENT
0
Entering edit mode

Okay, I will do that. Thank You.

ADD REPLY

Login before adding your answer.

Traffic: 2639 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