Question: Recode VCF file according to ancestral alleles (AA)
1
gravatar for abascalfederico
4.8 years ago by
abascalfederico1.1k
Spain
abascalfederico1.1k wrote:

I am working with 1000G vcf files. These files contain an AA field indicating which is the Ancestral Allele. I would like to use this information to recode the vcf file. When the REF allele is different from the AA,  REF and ALT alleles should be switched,  AF fields recalculated (derived AF=1-AF), and genotypes be properly recoded (1/1 --> 0/0, 0/1 --> 0/1...).

Are you aware of any tool doing this? If there were no such tool, I would probably try myself and share the script here.

Thanks!

Federico

vcf • 2.6k views
ADD COMMENTlink modified 4.8 years ago • written 4.8 years ago by abascalfederico1.1k
2
gravatar for abascalfederico
4.8 years ago by
abascalfederico1.1k
Spain
abascalfederico1.1k wrote:

Here is the solution I (dirty) coded - in case it is useful to others.

while(<>) {
    if(/^#/) {
        print $_;
        next;
    }
    chomp;
    my($chr,$site,$rs,$ref,$alt,$info) = (split(/\t/,$_))[0,1,2,3,4,7];
    my $ancestral_allele;
    $ancestral_allele = (split(/AA=/,$info))[1];
    $ancestral_allele = (split(/\|/,$ancestral_allele))[0];
    $ancestral_allele = uc($ancestral_allele);
    if($ancestral_allele !~ /[ACGT]/) {
        print STDERR "Ancestral allele ($ancestral_allele) is not ACGT: $_\n";
        $ancestral_allele = $ref;
    }

    if($ancestral_allele eq $ref) { #ok as it is
        print "$_\n";
        next;
    } else { #Switch alleles
        my @tmp = split(/;/,$info);
        foreach my $t ( @tmp ) {
            next if($t !~ /AF/);
            my($key,$val) = split(/=/,$t);
            my $new_AF = 1-$val;
            $_ =~ s/$key=$val/$key=$new_AF/;
        }
        @tmp = split(/\t/,$_);
        my $str = "";
        foreach my $t ( @tmp ) {
            if($t =~ /^[01][\|\/][01]/) {
                my($a,$b) = (split(/\|/,$t))[0,1];
                if($a == 0)    {    $a = 1;    }
                elsif($a == 1) {    $a = 0;    }
                if($b == 0)    {    $b = 1;    }
                elsif($b == 1) {    $b = 0;    }
                $t = "$a|$b";
            }
            $str .= "$t\t";
        }
        $str =~ s/\t$//;
        print "$str\n";
    }
}
ADD COMMENTlink modified 10 months ago by RamRS30k • written 4.8 years ago by abascalfederico1.1k
I made a post here, could you also help me answer it. Thank you very much. Ancestral Allele for the Genus Gallus
ADD REPLYlink written 4.8 years ago by lawalakinyanju0
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: 2101 users visited in the last hour