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
•
written
by
abascalfederico
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";
}
}
Log in to answer this question.
