intersectBed only intersects the first 9 chromosomes and X Y chromoses

[Angel Ba]

Hi all,

I am dealing with a normal-size VCF file and a big BED file (4GB). Both are sorted as indicated in the documentation. The thing is that when I intersect them, it works fine for the first 9 chromosomes and sex chromosomes, but doesn't interesect the rest (i.e. chr10 to chr22)

Has anyone experienced the same? Am I doing something wrong?

PS: the commands I already tried to use are the following (and all procude the same oputput... when it finished with chr9, it stops for a while till it arrives to chrX and chrY)

zcat myVCFfile.vcf.gz | intersectBed -header -sorted -a stdin -b mybigBEDfile.sorted.bed

or

zcat myVCFfile.vcf.gz | bedtools intersect -header -a stdin -b mybigBEDfile.sorted.bed -sorted

[Aaron Quinlan]

Hi Angel,

Could you provide a sense of how the chromosomes are ordered in the two files with the following two commands:

zcat myVCFfile.vcf.gz | grep -v ^# | cut -f 1 | uniq -c

cut -f 1 mybigBEDfile.sorted.bed | uniq -c

This will help me to understand what might be happening.  Also, what version of bedtools are you using?

- Aaron

quinlanlab.org

--
You received this message because you are subscribed to the Google Groups "bedtools-discuss" group.
To unsubscribe from this group and stop receiving emails from it, send an email to bedtools-discu...@googlegroups.com.
For more options, visit https://groups.google.com/groups/opt_out.

[Angel Ba]

Sure,

when I apply "zcat myVCFfile.vcf.gz | grep -v ^# | cut -f 1 | uniq -c" to my VCF file it appears:

   4048 chr1

   4911 chr2

   3309 chr3

   3746 chr4

   3455 chr5

   3603 chr6

   3033 chr7

   3101 chr8

   2139 chr9

   3050 chr10

   2296 chr11

   2262 chr12

   2037 chr13

   1727 chr14

   1535 chr15

   1598 chr16

   1612 chr17

   1807 chr18

   1029 chr19

   1320 chr20

   1001 chr21

    958 chr22

    813 chrX

     12 chrY

And when I apply "cut -f 1 mybigBEDfile.sorted.bed | uniq -c" to mybigBEDfile.bed it appears:

14020064 chr1

15690104 chr2

13089776 chr3

12802487 chr4

11751924 chr5

11179073 chr6

9524406 chr7

9475744 chr8

7049615 chr9

8153483 chr10

8270912 chr11

8350866 chr12

6664851 chr13

5678683 chr14

4825831 chr15

4347627 chr16

4276806 chr17

5116545 chr18

2487346 chr19

3727116 chr20

2204808 chr21

1841339 chr22

    249 chrM

8300227 chrX

 510526 chrY

Version:

Tool:    bedtools intersect (aka intersectBed)

Version: v2.17.0

Summary: Report overlaps between two feature files.

bedtools v2.17.0

El diumenge 16 de febrer de 2014 13:59:42 UTC+1, Aaron Quinlan va escriure:

[Aaron Quinlan]

Hi Angel,

Version 2.17 only supported lexicographical chromosome sorting (i.e. chr1, chr10, chr11 … chr2, chr3, etc.). Your files are "version sorted".  Staring with version 2.18 and continuing with the current version (2.19), alternative sorting criteria are supported, but you need to provide a "genome file" that specifies the expected chromosome ordering in your data files.  You could make such a file from the data files themselves:

$ zcat myVCFfile.vcf.gz | grep -v ^# | cut -f 1 | awk '{print $0"\t1"}' > my.genome

Then, using >= 2.18, you would do the following:

$ bedtools intersect -header -a myVCFfile.vcf.gz -b mybigBEDfile.sorted.bed -sorted -g my.genome

I hope this helps!

- Aaron