Easy way to split VCF file by chromosome

Question

[ConvolutedGenome]

Hi,

Im trying to submit a job on the TOPMED/Michigan imputation server, but it returns an error saying that I need to split my VCF by chromosome

Is there an easy way to do this? Will bcftools help?

Answer 1

[Pierre Lindenbaum]

bcftools index -s in.vcf.gz | cut -f 1 | while read C; do bcftools view -O z -o split.${C}.vcf.gz in.vcf.gz "${C}" ; done

[ConvolutedGenome]

omg thanks!!

[rheab1230]

Hello,

I tried to do this

bcftools index -s /scratch/Databases/GTEx/downloaded/dbGAP/GTEx_Analysis_2017-06-05_v8_WholeExomeSeq_979Indiv_VEP_annot.vcf.gz | cut -f 1 | while read C; do bcftools view -O z -o split.${C}.vcf.gz in.vcf.gz "${C}" ; done

but it throws an error:

[E::hts_open_format] Failed to open file in.vcf.gz
Failed to open in.vcf.gz: No such file or directory

[Pierre Lindenbaum]

of course it doesn't work. You should have a look at my command line and _understand_ how it works.

[Alex]

You should _explain_ your command

[rheab1230]

#!/bin/bash
source ~/miniconda3/bin/activate
conda activate train
cd /scratch/downloaded/dbGAP/
for C in {1 .. 22}
do
  bcftools index -s in.vcf.gz | cut -f 1 | while read C
  do
    bcftools view -O z -o /home/data_download/split.${C}.vcf.gz in.vcf.gz "${C}"
  done

I tried to do this but its throwing error:

line 9: syntax error: unexpected end of file

[nitin]

This will Work!

[Pierre Lindenbaum]

why a screenshot when you can just copy and paste the text ?

[Alex]

There are nicer ways of saying that

[Pierre Lindenbaum]

sorry for my French

[Pierre Lindenbaum]

This is just wrong. Why do you use two loops when all the chromosomes are available in bcftools index -s

[BroNobo]

I also tried this, but I am getting the following error:

[E::idx_find_and_load] Could not retrieve index file for 'in.vcf.gz'
Could not load index for VCF: in.vcf.gz

I saw elsewhere bcftools might be bugged (https://github.com/samtools/bcftools/issues/881)? Or am I doing this wrong? The workarounds seem a bit too involved for me to understand.

The command do bcftools view -O z -o split.${C}.vcf.gz in.vcf.gz "${C}" works great. I think it's just something with bcftools index.

For more context, I had some Illumina genotyping array data (.bed, .bim, .fam, .map, .ped) that I converted to .vcf using Plink (and then bgzip). Maybe that's why I'm losing the index somewhere?

Thanks for the help.

[Joana]

You need to run bcftools index on your vcf file before running the suggested command.

[Ram]

Or the more popular tabix -p vcf <vcf_file>.

Answer 2

[LTDavid]

See the explanation below the code

vcf_in=/scratch/Databases/GTEx/downloaded/dbGAP/GTEx_Analysis_2017-06-05_v8_WholeExomeSeq_979Indiv_VEP_annot.vcf.gz

vcf_out_stem=/scratch/Databases/GTEx/downloaded/dbGAP/by_chrom

for i in {1..22}
do
bcftools view ${vcf_in} --regions ${i} -o ${vcf_out_stem}_${i}.vcf.gz -Oz
done

1. First, define vcf_in using your path to your input file. Here, I assumed that your input file was /scratch/Databases/GTEx/downloaded/dbGAP/GTEx_Analysis_2017-06-05_v8_WholeExomeSeq_979Indiv_VEP_annot.vcf.gz.
2. Deine vcf_out_stem using the path of the new file. vcf_out_stem should only include the portion of the path and file name leading up to the chromosome number, assuming you want the chromosome number in the path name. Here, I assumed that it would have the same path as the input vcf.
3. You don't need to make changes to the rest of the code. The for line means that it's going to repeat the code 22 times. The range {1..22} corresponds to the chromosomes you want to subset. This is going to loop through; each time, it's going to take your input file and subset the chromosome using --regions ${i}. If your naming convention for your chromosome is different, you may need to change it to that. For example, if your chromosome name in your input put is chr1, then you may need to use --regions chr${i}. And then -o ${vcf_out_stem}_${i}.vcf.gz is going to give you the new file name by using whatever you defined for vcf_out_stem and adding on the chromosome number using i.

[gebk]

LTDavid:

My data have chromosomes 1..22, X, Y, and MT.

How do I modify the syntax for i in {1..22}

Thanks

[Ram]

To stick to the for, understand that {1..5} (for example) simply expands to 1 2 3 4 5, which means you can always add more space-separated values like {1..22} X Y MT.

for i in {1..22} X Y MT
do
  bcftools view ${vcf_in} --regions ${i} -o ${vcf_out_stem}_${i}.vcf.gz -Oz
done

[gebk]

Thank you! I will let you know.

[gebk]

The code works good. It created vcf.gz files by chromosome. I'm using Michigan Imputation Server. I uploaded all the files and run the imputation.

During the input validation phase, I got the following error message.

No valid chromosomes found!

Any idea why it gives the error message?

Thanks, again

[Pierre Lindenbaum]

it's unrelated to the original question.

[gebk]

Yes, it is but the next step after my earlier question. Shall I remove it and repost it on its own?

[gebk]

Hi Ram, could you please comment on this too? I want to try convert 1..22 to chr1..chr22 and run the imputation. How do I modify the above syntax? I really appreciate your help.

[Ram]

If you want to use chr1...chr22, change the command text like so:

for i in {1..22} X Y MT
do
  bcftools view ${vcf_in} --regions chr${i} -o ${vcf_out_stem}_chr${i}.vcf.gz -Oz
done

Just be aware that you'll probably need to use M and not MT as MT becomes chrM.

[gebk]

Thank you very much, Ram.

[gebk]

Ram, can you look at the syntax again? Something is not right. All the separate files have the same size which are probably close to empty. Thanks.

[Ram]

Run the bcftools command without the loop just for chr22 and see the results. If things look OK, add an echo before the bcftools command in the loop and check that the echo'd commands look fine. This is a great opportunity for you to get debugging.

[Pierre Lindenbaum]

(seq 1 22  && echo X Y MT) | while read i; do ....

[Ram]

Please do not add answers unless you're answering the top level question. Use Add Comment or Add Reply as appropriate. I've moved your post to a comment this time.