2
$\begingroup$

In my gzipped VCF file I have many variants, and they report allele frequency such as:

##INFO=<ID=CAF,Number=.,Type=String,Description="An ordered, comma delimited list of allele frequencies based on 1000Genomes, starting with the reference allele followed by alternate alleles as ordered in the ALT column. Where a 1000Genomes alternate allele is not in the dbSNPs alternate allele set, the allele is added to the ALT column. The minor allele is the second largest value in the list, and was previuosly reported in VCF as the GMAF. This is the GMAF reported on the RefSNP and EntrezSNP pages and VariationReporter">

and

CAF=0.5747,0.4253

I would like to filter out the variants >= 1% allelle frequency. How do I do that?

$\endgroup$
5
  • $\begingroup$ Did you try any program or regex? How big is your VCF file ? (perhaps you could post some lines with what you want and what you don't want) to test it $\endgroup$
    – llrs
    Jun 1 '18 at 7:05
  • $\begingroup$ @Llopis I want all lines that the allele frequency of the alternative allele is >= 0.01. In my example, I'd keep it because 0.4253 > 0.01. $\endgroup$
    – SmallChess
    Jun 1 '18 at 7:06
  • $\begingroup$ @Llopis My VCF files are big and gzipped. Direct regression expression won't work. $\endgroup$
    – SmallChess
    Jun 1 '18 at 7:07
  • 1
    $\begingroup$ I agree with @Llopis, you should post a few lines of your vcf file to give a workable example. You should probably use AWK. $\endgroup$
    – benn
    Jun 1 '18 at 7:54
  • 1
    $\begingroup$ How do you want to deal with multi-allelic variants (e.g. A C,T) which will have >2 CAF values? Do you want all of them to be <1% or any? And as others have said, we can't really do this without a few example lines to work with. Ideally, give us a few lines including a multiallelic one and show us your desired output. $\endgroup$
    – terdon
    Jun 1 '18 at 8:36
4
$\begingroup$

using vcffilterjdk:

you want to keep the variant if there is no INFO/CAF.

Otherwise, get All values of CAF as a String list (because it is declared as a string in the VCF header...), stream those data, skip the first value for REF , convert each string to double, you want to keep the variant if any allele is rare (< 0.01)

java -jar dist/vcffilterjdk.jar -e 'return !variant.hasAttribute("CAF") || variant.variant.getAttributeAsStringList("CAF","1.0").stream().skip(1L).map(A->Double.parseDouble(A)).anyMatch(A->A<0.01);' input.vcf
$\endgroup$
1
  • $\begingroup$ I have not encountered this tool before, but in my experience, if Pierre suggests something it will work as advertised. I was going to suggest a tool that I’d read about but never got around until last week: bioawk. It was blazingly fast. $\endgroup$
    – mdperry
    Sep 22 '18 at 21:55
0
$\begingroup$

This could be done with bcftools view. To include variants where the first element in CAF is > x:

bcftools view -i 'INFO/CAF[0] > 0.1' ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/common_all_20180418.vcf.gz

However, it appears that CAF is encoded as a list of strings in this input so the > operator cannot be used. In addition, you may want more control on the filtering. In this case, I would suggest using a dedicated library in your favourite language for dealing with VCF files. In python, you could use pysam along these lines:

python -c "import pysam
invcf= pysam.VariantFile('ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/VCF/common_all_20180418.vcf.gz')

# Write to stdout
outvcf= pysam.VariantFile('-', mode= 'w', header= invcf.header)

for line in invcf:
    if float(line.info['CAF'][0]) > 0.1:
        outvcf.write(line)
outvcf.close()
invcf.close()
"

I would stay away from parsing VCF files with regex, awk, or custom text processing scripts as even if the format is quite well defined they can be quite tricky.

$\endgroup$
1
  • $\begingroup$ The first value is for the reference allele. Presumably, the OP is interested in filtering by the variant's frequency instead. So it would have to be based on elements 1 through the end of the array. $\endgroup$
    – terdon
    Jun 2 '18 at 13:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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