This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Use tabix with a list of regions

Hi all,

I need to use tabix for a project and the following is the syntax that is used, for e.g.

tabix -fhftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20101123/interim_phase1_release/ALL.chr12.phase1.projectConsensus.genotypes.vcf.gz 12:2345295-2345295 > genotypes.vcf

The question I have is if i have a lot of regions in the same chromosome, is there a way that I can supply a file with the regions as input, what would the syntax be? In other words can I just replace the chromosome number and coordinate with a text file name with each region on a new line?

Thanks
Ashwin

tabix
• 22,754 views • link
updated by mg mg • written by win win

6 answers

Put your regions in a text file (one per line) and stream to tabix:

xargs -a regions.txt -I {} tabix -fh  ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20101123/interim_phase1_release/ALL.chr12.phase1.projectConsensus.genotypes.vcf.gz {} >> genotypes.vcf
• 89 views • link
written by JC JC

i just realized another issue, the vcf files i have are annotated with rs#, whereas the above approach will search using chromosomal coordinates, i wonder is there was a way to search by variant name using the above approach.

• 0 views • link
written by win win

you can convert the rs# to coordinates in dbSNP

• 0 views • link
written by JC JC

hi, i just tried to run in this command in a command window in ubuntu 12 and i get a message

[kftpconnectfile] 350 Restarting at 0. Send STORE or RETRIEVE to initiate transfer [main] fail to open the data file.

the command i typed was

xargs -a Chr12.txt -I {} tabix -fh ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20101123/interim_phase1_release/ALL.chr12.phase1.projectConsensus.genotypes.vcf.gz {} >> genotypes.vcf

any ideas?

• 0 views • link
written by win win

yes it does, remove that damn < > (I removed them from the response)

• 0 views • link
written by JC JC

hi jc, all of your input has certainly helped, but i had another question. When i run the aboe command i get a header line for each of the variant that was found. is it not possible to have just one header line accross the entire VCF file?

• 0 views • link
written by win win

you can remove the header removing the -h flag, but that removes all the headers, to obtain just one header, first generate the output file with a null base:

# generating the header
tabix -fh  ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20101123/interim_phase1_release/ALL.chr12.phase1.projectConsensus.genotypes.vcf.gz chr1:1 > genotypes.vcf
# then you can call all your positions
xargs -a regions.txt -I {} tabix -f  ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20101123/interim_phase1_release/ALL.chr12.phase1.projectConsensus.genotypes.vcf.gz {} >> genotypes.vcf
• 0 views • link
written by JC JC

Tabix has a -B option which allows you to input a .bed file of regions of interest:

tabix -fhB ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20101123/interim_phase1_release/ALL.chr12.phase1.projectConsensus.genotypes.vcf.gz regions.bed > genotypes.vcf
• 0 views • link
updated by Ram Ram • written by Leonor Palmeira Leonor Palmeira

I'm trying to get tabix to read in a .bed file as you say, so far without success. What should the format of the .bed lines be? I have tried multiple formats, including this specification at UCSC as well as the format suggested in an answer to similar Biostars question. Also, what about the comment by Rm on that answer, claiming that the -B option has been dropped from tabix? Is that out of date?

• 0 views • link
updated by Ram Ram • written by Greg P Greg P

ask a new question for your data please + show us a few lines of your bed.

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

I apologize, my only problem is that I was not careful in reading what people write. I thought that the .bed filename was part of the -B option. Of course, that's not how it works. The -B option just tells tabix to look for a bed file instead of 'region1'.

• 0 views • link
updated by Ram Ram • written by Greg P Greg P

@Greg P: I think Rm is talking about the outdated manual available here, I haven't seen that the -B option has been dropped (wouldn't make much sense to me either).

• 0 views • link
updated by Ram Ram • written by Leonor Palmeira Leonor Palmeira

-B has been dropped, and since version 1.2, tabix has -R option instead. See my answer for more info.

• 0 views • link
written by mg mg

Yet another way would be a BASH script as:

while read line;do
    tabix -fh ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20101123/interim_phase1_release/ALL.chr12.phase1.projectConsensus.genotypes.vcf.gz $line >> genotypes.vcf
done < coord_file

but I like the simplicity of the use of xargs in JC's answer otherwise.

• 14 views • link
written by cl.parallel • cl.parallel

thanks! I see you're in Seattle, where? I'm at ISB

• 0 views • link
written by JC JC

well no so far, maybe we can chat sometime

• 0 views • link
written by JC JC

Hi,

Thanks for this. I suppose this will generate as many files as the coord_file has right? Please let me know if I'm wrong. I'm trying to use tabix to extract information from vcf files for the coordinates I want, based on a text file with all my regions. I need to have as many files as my region file has. e.g if I have 100 regions, I will have 100 files as well. So in the example above, I will have 100 genotypes.$line.vcf.

Thanks

• 0 views • link
updated by Ram Ram • written by SGMS SGMS

I would rather curl the vcf and pipe the output with bedtools intersect

something like (I cannot check my command I'm away from my lab):

curl -s "ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/release/20101123/interim_phase1_release/ALL.chr12.phase1.projectConsensus.genotypes.vcf.gz"  | 
gunzip -c |
intersectBed -a stdin  -u -b list.bed
• 0 views • link
updated by Ram Ram • written by Pierre Lindenbaum Pierre Lindenbaum

I think it's better to use tabix. Curl forces you to download the whole file locally; while tabix can take advantage of the tabix indexes (which are already in the 1000genomes ftp), and download only the regions of interest.

• 0 views • link
written by Giovanni M Dall'Olio Giovanni M Dall'Olio

I don't think there is any magic on the server side. Tabix download the data and he knows how many bits he should 'skip' and when he should 'stop' scanning. But if all your regions are at the bottom (e.g chrY), it will have to fetch everything each time.

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Giovanni, I agree. But once tabix has downloaded the index, it must download the data itself from the server until it get the right offset.

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

This is a python script that I wrote as a wrapper to call tabix on a list of gene coordinates:

However, the xargs solution is more elegant. The only advantage of my code is that it doesn't download files if they are already present in the directory, so you can interrupt the download and continue it later, or add genes to the list after downloading them all.

• 0 views • link
written by Giovanni M Dall'Olio Giovanni M Dall'Olio

Since version 1.2, tabix uses -R option to accept file with regions of interest. From tabix's commandline help message:

-R, --regions FILE restrict to regions listed in the file

And tabix's doc has more info on how this file should look like:

-R, --regions FILE.      Restrict to regions listed in the FILE. The FILE can be BED file 
      (requires .bed, .bed.gz, .bed.bgz file name extension) or a TAB-delimited file with 
      CHROM, POS, and, optionally, POS_TO columns, where positions are 1-based and inclusive. 
      When this option is in use, the input file may not be sorted.

Note that, unlike typical bed file (where START is zero-based), both POS and POS_TO are one-based.

• 0 views • link
written by mg mg

Log in to answer this question.

More posts like this