I am looking for a tool/script/pony to correct the REF column in a vcf file whenever that nucleotide doesn't match the reference genome, as supplied in fasta format.
It sounds like a common task but I could not find something. I did find bcftools +fixref but that only works for SNPs. My vcf files are from structural variants.
Cheers,
W
Some code like this can help, but be aware that changing the reference column affects other fields in your VCF, such as the GT value.
#!/usr/bin/perl
use strict;
use warnings;
$ARGV[2] or die "use fixRefVCF.pl <GENOME_FASTA> <VCF> <FIX_VCF>\n";
my $genome = shift @ARGV;
my $orig_vcf = shift @ARGV;
my $new_vcf = shift @ARGV;
my %seq = ();
my ($chr, $pos, $ref, $obs);
warn "loading genome sequence from $genome\n";
open (my $fh, "<", $genome) or die "cannot read $genome\n";
while (<$fh>) {
chomp;
if (/>(\w+)/) {
$chr = $1;
} else {
$seq{$chr} .= $_;
}
}
close $fh;
open (my $vh, "<", $orig_vcf) or die "cannot read $orig_vcf\n";
open (my $oh, ">", $new_vcf) or die "cannot write $new_vcf\n";
while (<$vh>) {
if (/^#/) {
print $oh $_;
} else {
chomp;
my @var = split (/\t/, $_);
$chr = $var[0];
$pos = $var[1]; # VCF is 1-based, Perl string is 0-based
$ref = $var[3];
$obs = substr($seq{$chr}, $pos - 1, length $ref);
if ($ref ne $obs) {
warn "editing $chr:$pos $ref->$obs\n";
$var[3] = $obs;
}
print $oh join "\t", @var;
print $oh "\n";
}
}
close $vh;
close $oh;
Hello WouterDeCoster ,
there is another bcftools plugin that should do this:
$ bcftools +fill-from-fasta input.vcf -- -c REF -f genome.fa > output.vcf
fin swimmer
This is the python solution I came up with, using cyvcf (minimally requiring v0.10.2) and pyfaidx. It can use compressed vcf and fasta files.
Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
version 2.3.6
tagging: Pierre Lindenbaum
what does it currently look like in the REF/ALT column ? (why do you need to fill the REF column (thinking of a very large deletion...) in your downstream analysis ? )
My main issue is that EVA complains that the nucleotides are incorrect. It seems their software only uses POS to check the REF nucleotide, and not END.
Just some lines from the vcf:
Essentially a single nucleotide would be sufficient - and in one file 98% is correct. The other file has
Nfor every position so some more changes need to be made.I'm already working on a python script (now posted as answer) but I thought someone else should have done it already ;p
For your narrow purpose, where all your insertions are symbolic rather than explicitly spelled out, the existing answers should be a satisfactory workaround to "bcftools +fixref"'s insistence on SNP-only data.
However, it's worth noting that there is a very good reason that bcftools insists on only SNPs. It's actually impossible to fix ordinary short indels: the only difference between a 1-base deletion and an insertion of the same base at the same position is the REF/ALT order.