Efficient way of vcf files handling
1
0
Entering edit mode
7.9 years ago
Earendil ▴ 50

Hi all

I have a .tgz file in which there are folders, one for each individual. Each folder contains vcf files, one for each chromosome. It's in this fashion: .tgz has ./individual1/chr1.vcf ./individual1/chr2.vcf etc, ./individual2/chr1.vcf ./individual2/chr2.vcf etc.

I want to create vcf files, one for each chromosome, containing individuals of a population using bcftools merge, so that later I can use the vcf files for an Fst estimation among populations using vcftools.

Having no previous experience for this kind of work, the only thing I can think is to extract the .tgz file, compress each vcf with bgzip and index it with tabix, run bcftools merge for each chromosome and population and finally run the Fst test on vcftools.

But is there a better way to do this? Extracting the .tgz file will chew up a lot of storage space.

vcf bcftools genome vcftools • 2.8k views
3
Entering edit mode
7.9 years ago

It's a not-widely-appreciated feature of tar that you can just untar the files you want; so that, for instance, you could do this by chromosome if disk space is the concern:

#!/bin/bash
mkdir -p merged
for chrom in X Y $( seq 1 22 ) do tar xzvf vcfs.tgz individual{1,2,3,4,5}/chr${chrom}.vcf
for individual in individual*
do
cd $individual bgzip chr${chrom}.vcf
tabix -p vcf chr${chrom}.vcf.gz cd .. done bcftools merge -m all --force-samples -O z individual*/chr${chrom}.vcf.gz  > merged/chr${chrom}.vcf.gz rm -f individual*/chr${chrom}.*
done

0
Entering edit mode

I don't see the "tar's under-appreciated feature" part. Really good use of a bash loop and the bash 3.0 {} syntax, but all bash commands can "just work on the files one wants".

1
Entering edit mode

It's been my experience lately that a number of people - including many who are quite proficient in the linux shell environment - don't realize you can just pull out the individual files you want out of a tar file, rather than untarring the whole thing. That may indeed be already familiar to you - I assumed everyone knew that myself until recently.

0
Entering edit mode

Apologies - I see what you're highlighting now. I'd misread OP's question. It is indeed an underused feature, because we rarely have need for that today. In fact, it was only when I had to pick a script from a large project archive on my HPC that I discovered this nifty little feature :-)

0
Entering edit mode

Thank you very much for your answer! I am just a starter on the linux shell environment so I need to explore various ways to do things. I have been experimenting with scripts based on your answer and the storage problem is solved.

To try to optimize it a bit more, would there be a way to make the process faster? With brace and parameter expansion, the tar xzvf command runs anew each time for every different individual/chromosome combination, so it needs to sweep the whole tar file each time. This makes things slow; is there any possible way to avoid it? I could leave bcftools merge aside for the time being if that is a problem.

Thanks again