Question: vcftools and multiple alternative allelles
0
gravatar for oddmundn
3.4 years ago by
oddmundn0
oddmundn0 wrote:

I am trying to extract only variants which as an ALT allele frequency above a certain value, using vcftools and the paramrter --non-ref-af-any. Several of my variants have several alternative alleles, but when using the --non-ref-af-any it doesn't filter at all. Every variant is included, despite having AF=0 for the only ALT alllele. When using the --non-ref-af option, only single ALT alleles with AF > x is extracted.

Is this a bug in vcftools? I am using VCF files from Ion Torrent Suite.

 

Thanks for any help!

myposts • 1.9k views
ADD COMMENTlink modified 3.4 years ago by abascalfederico1.1k • written 3.4 years ago by oddmundn0
0
gravatar for abascalfederico
3.4 years ago by
abascalfederico1.1k
Spain
abascalfederico1.1k wrote:

I think it is not a bug. The documentation says that the site is not filtered if "any" of the variants meet the specified criteria. I had a similar problem and after struggling with different programs finally wrote my own perl script to remove AF=0 variants. It works with multi-allelic variants - only the affected allele will be removed. But it only works for the specific case of alleles with AF=0 (after Beagle refinement in my case). Here it is in case you want to try it.

#!/usr/bin/perl -w

use strict;
$|=1;

while(<>) {
    if(/^#/) {
        print;
        next;
    }
    chomp;
    my @tmp = split(/\t/,$_);
    my %allele_counts;
    my($chr,$pos,$id,$ref,$alt,$qual,$filter,$info,$format) = (@tmp)[0..8];
    my $num_alt;
    my @alts;
    $alts[0] = $ref;
    if($alt =~ /,/) {
        my @tmp2 = split(/,/,$alt);
        $num_alt = scalar(@tmp2)+1;
        push(@alts,@tmp2);
    } else {
        $num_alt = 2;
        $alts[1] = $alt;
    }
    for(my $i=9; $i<=$#tmp; $i++) {
        my $genotype = (split(/:/,$tmp[$i]))[0];
        my($a,$b) = (split(/\|/,$genotype))[0,1];
        $allele_counts{$a}++;
        $allele_counts{$b}++;
    }
    if(scalar(keys(%allele_counts)) == 1) {
        print STDERR "REMOVED, AF=0: $_\n";
        next;
    }
    if(!exists($allele_counts{'0'})) {
        print STDERR "REMOVED, REF ALLELE HAS AF=0: $_\n";
        next;
    }
    if(scalar(keys(%allele_counts)) != $num_alt) {
        print STDERR "REMOVED: There are multiple alleles but at least one has to be removed: $_\n";
        for(my $i=1; $i<=$#alts; $i++) {
            if(!exists($allele_counts{$i})) {
                if($i==1) {
                    $tmp[4] =~ s/$alts[$i],//;
                    for(my $j=9; $j<=$#tmp; $j++) {
                        my $recode = $i+1; #al siguiente alelo tenemos que bajarle la cuenta
                        $tmp[$j] =~ s/\|$recode/\|1/;
                        $tmp[$j] =~ s/$recode\|/1\|/;
                    }
                } else {
                    $tmp[4] =~ s/,$alts[$i]//;
                    print STDERR "\tRemoving allele $i=$alts[$i]\n";
                }
            }
        }
    }
    print "@tmp\n";
}    


__END__
ADD COMMENTlink modified 3.4 years ago • written 3.4 years ago by abascalfederico1.1k

 

Thank you very much!

Unfortunately, your filter removed all entries in the VCF file, even though some of them have AF > 0.  Used this way: scriptname filename.vcf > outfile.vcf

The VCF file is genereated by Ion Torrent Suite software (Life Technologies/Thermo).

 

Oddmund

 

ADD REPLYlink written 3.4 years ago by oddmundn0

Could you send me the first 1,000 or 10,000 lines of your vcf to test the script? You can email me to fa8 - at - sanger.ac.uk

 

ADD REPLYlink written 3.4 years ago by abascalfederico1.1k

Earlier today I wrote my own python script to filter NOCALLS and AF=0 from my files, which works fine on my VCF files. That means; I don't need help anymore. Please tell me if you would like  the "head" of my Ion Torrent VCF file anyway?

ADD REPLYlink written 3.4 years ago by oddmundn0

Glad to hear that. I don't need it, but thanks

 

ADD REPLYlink written 3.4 years ago by abascalfederico1.1k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1765 users visited in the last hour