Lifting over VCF's from one reference to another

delangel

liftOverVCF.pl

Contents 1 Introduction 2 Obtaining the Script 3 Example 4 Usage 5 Chain files

Introduction

This script converts a VCF file from one reference build to another. It runs 3 modules within our toolkit that are necessary for lifting over a VCF.
1. LiftoverVariants walker
2. sortByRef.pl to sort the lifted-over file
3. Filter out records whose ref field no longer matches the new reference

Obtaining the Script

The liftOverVCF.pl script is available in our public source repository under the 'perl' directory. Instructions for pulling down our source are available here.

Example

./liftOverVCF.pl -vcf calls.b36.vcf \
  -chain b36ToHg19.broad.over.chain \
  -out calls.hg19.vcf \
  -gatk /humgen/gsa-scr1/ebanks/Sting_dev
  -newRef /seq/references/Homo_sapiens_assembly19/v0/Homo_sapiens_assembly19
  -oldRef /humgen/1kg/reference/human_b36_both
  -tmp /broad/shptmp [defaults to /tmp]

Usage

Running the script with no arguments will show the usage:

Usage: liftOverVCF.pl
    -vcf        <input vcf>
    -gatk       <path to gatk trunk>
    -chain      <chain file>
    -newRef     <path to new reference prefix; we will need newRef.dict, .fasta, and .fasta.fai>
    -oldRef     <path to old reference prefix; we will need oldRef.fasta>
    -out        <output vcf>
    -tmp            <temp file location; defaults to /tmp>

• The 'tmp' argument is optional. It specifies the location to write the temporary file from step 1 of the process.

Chain files

Chain files from b36/hg18 to hg19 are located here within the Broad:

   /humgen/gsa-hpprojects/GATK/data/Liftover_Chain_Files/

External users can get them off our ftp site:

   location: ftp.broadinstitute.org
   username: gsapubftp-anonymous
   path:     Liftover_Chain_Files

Xiaomao

Where can I get the liftOverVCF.pl script? I am not able to find the "perl" dir ftp.broadinstitute.org/distribution. Thanks.

ebanks

It's available in our source code through github (instructions for downloading source are available on this forum).

Xiaomao

Got it, thanks!

lorania

Is the liftOverVCF.pl doing the right thing with the arguments for LiftoverVariants? When I'm using it according to the help, it fails complaining "Input files variant and reference have incompatible contigs: Found contigs with the same name but different lengths".

I think what happens is that it passes to LiftoverVariants old version of reference genome instead of a new one (in LiftoverVariants documentation refrence genome is not listed?) and then it clashes with the dict file. Below is what gets invoked by my command:

  ~/bin/GenomeAnalysisTK-2.1-11-g13c0244/liftOverVCF.pl  -vcf refIndels20111102_4str.vcf -gatk ~/bin/GenomeAnalysisTK-2.1-11-g13c0244/ -chain ../mm9/mm9ToMm10.over.chain -newRef ../mm10/mm10 -oldRef ../mm9/mm9 -out refIndels20111102_4str_mm10.vcf
(...)
INFO  13:33:53,122 HelpFormatter - Program Args: -T LiftoverVariants -R ../mm9/mm9.fa -V:variant refIndels20111102_4str.vcf -o /tmp/0.329585949265496.unsorted.vcf -chain ../mm9/mm9ToMm10.over.chain -dict ../mm10/mm10.dict  

Am I right or is something wrong with my input files?

thanks,

Ania

ebanks

Something is wrong with your input files.

philliprichmond

location: ftp.broadinstitute.org username: gsapubftp-anonymous path: Liftover_Chain_Files

This path no longer exists. Is there an updated location for the Chain files? I'm looking for mm9 to mm10.

Thanks, Phil

Geraldine_VdAuwera

Hi Phil,

I just checked and it's there. If you're accessing the ftp through a shell connection, try using a gui client instead. People have reported having problems with shell access; the ftp setup is not under our control so I can only recommend working around the problem.

Geraldine_VdAuwera

You can also find these files in our public repository (see Downloads page for link).

hoosier060

Hi I was wondering what went wrong with my command to lift a VCF file: The script stops with a message:

---

Re-sorting the vcf... Unknown option: tmp

The sorting step failed. Please correct the necessary errors before retrying.

---

So I did sorting using vcfsorter.pl, then ran again, but still getting the same error message. could I get some advise on how to fix this issue please? Thanks, -Charles

Geraldine_VdAuwera

Can you tell me what version you are using and what is your command line?

hoosier060

GATK version is (version 2.2-8-gec077cd) and the command line as follows:

perl liftOverVCF.pl -vcf my.vcf \ -chain Liftover_Chain_Files/b37tohg19.chain \ -out mynew.vcf \ -gatk /usr/local/GATK \ -newRef hg19_canonical \ -oldRef human_g1k_v37 \ -tmp temporary_folder

Thanks, -Charles

ebanks

My guess is that you don't have the sortByRef.pl program in $gatk/public/perl/

hoosier060

liftOverVCF.pl and sortByRef.pl is under the same directory as gatk directory I used.

Geraldine_VdAuwera

Can you try without specifying a tmp folder?

hoosier060

without specifying a tmp folder still gives the same error...

Geraldine_VdAuwera

This is a weird error. I'll put this on my to-do bugfix list, but I won't be able to get to it until Tuesday at the earliest, sorry.

danielecook

Is there a liftover chain file called 'b36ToHg19.broad.over.chain' as you see in the example? I have looked in the 'Liftover Chain Files' and dont see it.

Geraldine_VdAuwera

Some of the liftover chain files may only be accessible internally at the Broad, sorry.

danielecook

Okay so how can I go about lifting over from b36 to hg19 in that case?

Geraldine_VdAuwera

I'll see if we can make the other liftover files publicly available. In the meantime I'm afraid I can't offer you any solutions, sorry.

PeteHaitch

Are the .chain files different to those provided by the UCSC genome browser, e.g. http://hgdownload.cse.ucsc.edu/goldenPath/mm9/liftOver/mm9ToMm10.over.chain.gz? If so, are you able to make the mm9 to mm10 files publicly available? Thanks

Geraldine_VdAuwera

I believe they are the same format, so the UCSC ones should work with our tool if that's what you're asking.

PeteHaitch

Thanks, Geraldine. That's indeed what I was asking. It might be worth noting that in the documentation.

nikmal

I would like to liftover from hg19 to GRCh37, but there is no chain file for this conversion on the FTP server. Is it possible to get access to such a file?

If I have understood this correctly, the difference is that hg19 uses 0-based coordinates and GRCh37 uses 1-based coordinates, as well as 'chr1' instead of '1' for the contig names. If I were to just write a script that removes all the 'chr' from the contig names and adds one to the coordinates, would that do the job correctly?

ebanks

Sure, I just added it to the ftp site (and we'll add it to the codebase on github too).

Incidentally, hg19 does NOT use 0-based coordinates. Also, your renaming proposal wouldn't work with some of the contigs (e.g. chrM -> MT).

nikmal

@ebanks said: Sure, I just added it to the ftp site (and we'll add it to the codebase on github too).

Incidentally, hg19 does NOT use 0-based coordinates. Also, your renaming proposal wouldn't work with some of the contigs (e.g. chrM -> MT).

Oh right, thank you for enlightening me. :) And thanks a lot for providing the file!

nikmal

Hm, I get a strange error when I try to do the liftover:

Bad input: there is a problem with the chain file you are using: chain line does not start with 'chain' in chain file hg19tob37.chain at line 1

I checked the file and line 1 does indeed start with 'chain'. I compared with the b37tohg19 chain file and the only difference is the spacing; single space in hg19tob37 and tab space in b37tohg19. I replaced single space to tab, since b37tohg19 worked for me before, but it didn't help.

Full command:

./liftOverVCF.pl -vcf callsToLiftOver.vcf -chain hg19tob37.chain -out calls_b37.vcf -gatk ~/software/GATK_2.5-2 -newRef GRCh37-lite -oldRef ucsc.hg19 -tmp ~/tmp/

What is wrong?

EDIT:

I figured out what was wrong: the file was stored with some strange encoding, when I read the first line of 'hg19tob37.chain' in Python I got this: '\xff\xfec\x00h\x00a\x00i\x00n\x00 \x001\x00 \x00c\x00h\x00r\x001\x00 \x002\x004\x009\x002\x005\x000\x006\x002\x001\x00 \x00+\x00 \x000\x00 \x002\x004\x009\x002\x005\x000\x006\x002\x001\x00 \x001\x00 \x002\x004\x009\x002\x005\x000\x006\x002\x001\x00 \x00+\x00 \x000\x00 \x002\x004\x009\x002\x005\x000\x006\x002\x001\x00 \x001\x00\n'

And when I opened the file in GEdit (Ubuntu 12.04) it looked normal. I simply copy-pasted the contents into a new text-file and saved it in UTF-8 encoding so now the liftover works. I hope I didn't mess anything up by doing this, but so far everything seems to have worked out well.

valentin

Thanks for the diagnosis. I have updated the file in the ftp site with a version that should work.