A VCF file describes how to translate the reference genome into (one or more) alternative haplotypes. I need a VCF file that translates one of these haplotypes back into the reference. In other words, I need to "invert" the VCF file I have.
If only SNPs were present, one would simply swap reference and alternative haplotype columns. Indels make it necessary to update the position. Of course, larger structural variants would also change the coordinates (sometimes dramatically), but for now I only care about SNPs and indels.
It does not seem too difficult to code, but would be quite time consuming (for me) to do it right.
Is there a tool for it?
•
written
by
Michael Kosicki •
Ok, I think I've found the solution myself:
## get mouse chr19
curl ftp://ftp.ensembl.org/pub/release-91/fasta/mus_musculus/dna/Mus_musculus.GRCm38.dna.chromosome.19.fa.gz > chr19.fa.gz
zcat chr19.fa.gz > chr19.fa
## get BALBc SNPs
curl ftp://ftp-mouse.sanger.ac.uk/current_snps/strain_specific_vcfs/BALB_cJ.mgp.v5.indels.dbSNP142.normed.vcf.gz > dir.raw.vcf.gz
zcat dir.raw.vcf.gz > dir.raw.vcf
## save header
grep '^#' dir.raw.vcf > dir.vcf
## collapse to first haplotype (phased or not)
grep -v '^#' dir.raw.vcf | sed 's/,[GTCA]*//' > dir.hap.vcf
## clean up (ie turn CCCA>CA to CCA>A etc)
cat dir.hap.vcf | awk 'BEGIN{OFS="\t"} \
{ if(length($4)>length($5)) {min=length($5)} else { min=length($4) } } \
{ print $1,$2+min-1,$3,substr($4,min),substr($5,min),$6,$7,$8,$9,$10 }' \
> dir.hap.clean.vcf
## resort
sort -k1,1 -k2,2n dir.hap.clean.vcf > dir.hap.sort.vcf
## remove overlapping indels
grep -v '^#' dir.hap.sort.vcf | awk 'BEGIN{pos=0; chr=1} \
{ if($1!=chr) {chr=$1; pos=0} } \
{if($2>pos) { print $0; pos=$2-(length($5)-length($4)) } }' \
>> dir.vcf
## zip and index
bgzip -f dir.vcf && tabix -p vcf dir.vcf.gz
## convert ref to alt
bcftools consensus -f chr19.fa dir.vcf.gz > chr19.alt.fa
## invert vcf
zcat dir.vcf.gz | grep '^#' > alt.vcf && # get header
zcat dir.vcf.gz | grep -v '^#' | awk 'BEGIN{OFS="\t"; chr=1; i=0} \
{ if($1!=chr) {chr=$1; i=0} } \
{print $1,$2+i,$3,$5,$4,$6,$7,$8,$9,$10; i=i+length($5)-length($4)}' \
>> alt.vcf
## zip and index
bgzip -f alt.vcf && tabix -p vcf alt.vcf.gz
## convert alt to ref
bcftools consensus -f chr19.alt.fa alt.vcf.gz > chr19.dir2.fa
## check if orignal restored, no difference
cmp --silent chr19.dir2.fa chr19.fa && echo "Files are identical." || echo "Files are different."
This solutions deals with (i.e. removes) alternative haplotypes and overlapping indels (i.e. alternative haplotypes in practice). Of course, one might want to handle these differently, e.g. mask the offending (likely heterozygous) regions completely.
Looks good - if you are certain that it does what you need, then I can move it to an Answer, and then you can Accept it. Let me know.
Sure, go ahead, thanks.
Done - you should be able to Accept the answer now
Log in to answer this question.
I guess you can use FastaAlternateReferenceMaker from GATK, but other way around.
Could you elaborate? This tool just "merges" reference fasta and a vcf file, and does not produce a new vcf file, as far as I can tell.
Here is a solution but using coding, and tested on a 1000 Genomes BCF file that i have. You can see that it functions correctly based on the fact that I had previously set the ID field to unique values based on CHROM:POS:REF:ALT:
I need to also "invert" the indels, this script seems to ignore them and only "invert" the SNPs.
Yes, I wrote it to just flip the SNPs / SNVs. You had mentioned how, by flipping indels, the position would change. I'm actually looking at it right now and I believe that the position stays the same for indels. This code will do it:
...or, even shorter:
I don't get it. Consider ATC>ATGCT with indels 2 T TG and 3 C CT. "Inverse" indels are 2 TG T and 4 CT C - second one shifted due to the first one "adding" sequence.
It is for those cases that you should always normalise your VCF variants by left-aligning indels and splitting multi-allelic calls into multiple lines. This helps to drastically reduce ambiguity. Even the publicly-available 1000 Genomes data is not normalised.
This can easily be done with:
When you do this and you flip the REF and ALT alleles, you do not need to subsequently change the POS field. Please test it by loading the different variant positions in the UCSC Genome Browser from my example above. If you additionally wanted to record the end-position of the indel, which is not recorded as standard in VCF format, then that would of course change.
Edit: if you look at the example that I posted above, I normalised that 1000 Genomes BCF file - see entries 8 and 9 in the header:
Can you please confirm whether this has helped or not by up-voting and/or accepting (also for the post of grant.hovhannisyan)? We are volunteers here and give up free time to help others. We appreciate acknowledgements and, more importantly, it's good practice to show others with the same problem the solutions / tips that worked.
Thanks, Kevin
I think I might not have made myself clear.
Consider a fasta reference "ref1" and a vcf file "vcf1" such that replacing all the variants in "ref1" with ones in "vcf1" gives me an alternative reference "ref2", which may have more or less nucleotides than "ref1" depending on the types of indels in "vcf1". What I need is such a vcf file "vcf2", that when I replace all the variants in "ref2" with ones in "vcf2" I would get "ref1" again. In pseudo-mathematic form: ref1+vcf1=ref2, ref2+vcf2=ref1.
Another example (same as ATC>ATGCT, but more dramatic): let's say there's two variants in vcf1: a 1000nt insertion and a SNV on the next position. Ref2 would be 1000nt longer than ref1. If I want to convert ref2 to ref1 again, the SNV entry would have to have a position 1002, with respect to the insertion. No matter, which way one aligns the indels, that's what it's going to be.
And yes, I've checked, with normalization or without (bcftools norm does not modify my vcf file, because it's already normed) and it does not work. Applying vcf-convert to ref2 (obtained with vcf-convert on ref1 and vcf1) and vcf2 (obtained with your script, which does swap column 4 and 5 as you say) yields an error of "the ref does not agree at...", as expected.
If you have only SNPs, then just swap the REF and ALT column in initial vcf file. The same with indels, but you have to recalculate coordinates in fasta and vcf (here you have to code, don't know any tools for that)
Hi Michael, thanks for getting back. That does indeed seem like a complex and very specific problem! It's not impossible but I think that you may have to go the coding route.
A quick search through Google reveals possible leads:
https://bioinformatics.stackexchange.com/questions/2825/converting-a-vcf-into-a-fasta-given-a-reference-with-python-r
https://bioinformatics.stackexchange.com/questions/2223/how-to-manipulate-a-reference-fasta-or-bam-to-include-variants-from-a-vcf
*the final link is what Grant had suggested
Kevin