How to sort a VCF file lexicographically by chromosome number?

Question

[Parth Patel]

Hello,

I have a large 155mb vcf file that I am trying to sort lexicographically by chromosome number (in the order 1,10, 11, 12 .. etc ) and then by position. I am trying to use vcf-sort from the vcf-tools perl package.

I have been running the command :

perl vcf-sort myfile.vcf > output.vcf

however, the output file loses all chromosome, position, and ref information for some reason and sorts variants by ID. According to the help page, this is suppose to sort the variants in lexicographic order.

I tried the command:

perl -c vcf-sort myfile.vcf > output.vcf

and the output file is sorted in chromosomally chronological order (1,2,3,4..) and all of the necessary information is retained.

Any idea on how to obtain a lexicographical sort by chromosome number (and then position)?

• Parth

[lavibharath]

Could you please post the command again...I am having some issues with the following command. Thanks.

grep -v '^#' in.vcf | LC_ALL=C sort -t '\t' -k1,1 -k2,2n >> out.vcf

error as follows:

sort: multi-character tab `\\t'

[Pierre Lindenbaum]

http://stackoverflow.com/questions/10627989/how-do-i-insert-a-tab-character-in-iterm

[Shaun Jackman]

Use $'\t'

grep '^#' in.vcf > out.vcf && grep -v '^#' in.vcf | LC_ALL=C sort -t $'\t' -k1,1 -k2,2n >> out.vcf

[lavibharath]

Oops... I got it. Thanks.

[lavibharath]

But the random chromosomes are really giving some problems...

Answer 1

[Pierre Lindenbaum]

grep '^#' in.vcf > out.vcf && grep -v '^#' in.vcf | LC_ALL=C sort -t $'\t' -k1,1 -k2,2n >> out.vcf

[dsmarcoantonio]

Running on Debian 8 is slight different

grep '^#' in.vcf > out.vcf && grep -v '^#' in.vcf |  sort -k1,1 -k2,2n >> out.vcf

Great solution Pierre

[Parth Patel]

Hi Pierre,

Thanks for the quick response. I tried your grep method and it yielded the same results as vcf-sort. The variants are sorted lexicographically by id, but the chromosome, position, and ref information is gone:

##fileformat=VCFv4.1
#CHROM POS   ID                REF   ALT   QUAL  FILTER INFO
             rs10110847_dbSNP        C    100   PASS
             rs10130831_dbSNP        G    100   PASS
             rs10132024_dbSNP        T    100   PASS
             rs10133030_dbSNP        A    100   PASS
             rs10135051_dbSNP        T    100   PASS
             rs10136467_dbSNP        C    100   PASS
             rs10147247_dbSNP        T    100   PASS
             rs10149947_dbSNP        T    100   PASS
             rs10152678_dbSNP        A    100   PASS
             rs10188587_dbSNP        T    100   PASS
             rs10277675_dbSNP        A    100   PASS

Is vcf-sort supposed to be doing this? I am using grep off of cygwin, so I am beginning to wonder if there is something wrong with my grep installation..

Parth

[Pierre Lindenbaum]

If

awk '/^\t/ {print}' in.vcf

produces any ouput then your VCF is corrupted.

[Parth Patel]

That was really helpful -- a formatting issue was in fact giving me a corrupt VCF file. Thanks!

[lavibharath]

Hi Pierre,

Thanks. Sorry, I am bit lost. Where to give tab (control+v+tab) as I still get the following error if cut and paste the aforesaid command. Thanks.

sort: multi-character tab `\\t'

[Pierre Lindenbaum]

You don't really to write the word '\t' but under linux: single-quote CTRL-V TAB single-quote

[Maximilian Haeussler]

Pierre, in bash, you can use $'\t' which is the same thing as "Ctrl-V tab" but can be copy/pasted from stackoverflow.

Answer 2

[ATpoint]

A solution without any grep:

cat in.vcf | awk '$1 ~ /^#/ {print $0;next} {print $0 | "LC_ALL=C sort -k1,1 -k2,2n"}' > sorted.vcf

[ibelcarri]

Is there any reason for using cat first and then awk instead of using directly awk and send the file as input?

[Pierre Lindenbaum]

no, but people often do this when they're exploring the data, inserting some tools in the pipeline. Imagine it could be gunzip -c in.vcf.gz instead of cat in.vcf