bedtools intersect zero output

Question

[tnminh89]

Hi everyone,

I am trying to see chromosome coordinates (.bed file) are in a gene in the same direction, opposite direction, or not overlapping a gene. I downloaded the genome file from biomart. Here is the command

for i in *_sorted_removed.bed; do bedtools intersect -wa -S -a $i -b biomart_mm39.bed > "${i%_sorted_removed.bed}"_opposite_gene.bed; bedtools intersect -wa -s -a $i -b biomart_mm39.bed > "${i%_sorted_removed.bed}"_same_gene.bed; bedtools intersect -v -a $i -b biomart_mm39.bed > ${i%_sorted_removed.bed}_no_overlap_gene.bed; done

Here is the example of the input files:

chr4    91687896        91694326        UID-13571       -       2643    741     6430    0.1152411
chr9    123953239       123959891       UID-1573        -       143     5194    6652    0.7808178000000001
chr10   121917511       121924041       UID-10126       -       142     4126    6530    0.631853
chr10   33797594        33803682        UID-10121       -       132     4690    6088    0.7703679

Here is the example of biomart file:

chr1    3143476 3144545 ENSUSG00000102693       +       4933401J01Rik   TEC     .
chr1    3172239 3172348 ENSUSG00000064842       +       Gm26206 snRNA   .
chr1    3276124 3741721 ENSUSG00000051951       -       Xkr4    protein_coding  .
chr1    3322980 3323459 ENSUSG00000102851       +       Gm18956 processed_pseudogene    .

While there was no error when I ran that command, the outputs for chromosome coordinates in the same-direction and opposite-direction gene were all 0. There were outputs for chromosomal coordinates not overlapping a gene. When I check these coordinates in IGV, there should have been outputs for both same-direction and opposite-direction files, and the outputs for no overlapping were wrong. For example:

chr1    5048503 5054639 UID-8777        -       10      1205    6136    0.1963820

This should have been in the output for the same-direction gene output. In IGV it is in the same-direction gene (Rgs20).

I have tried to change the format of the biomart genome file (changing columns, checking tabs, make sure the numbers of columns are the same, sorting biomart files) but the errors were still there. I tried this command and genome file with other samples and this still happened.

Anyone has any idea of what is going on here? Thank you so much!!

Answer 1

[tnminh89]

I figured out what the problem was. It turned out bedtools can only really work for my bed files if the strand information is on the 6th column. Once I reorganized my biomart bed file, it worked.

Answer 2

[lethalfang]

Did you miss a do before the first bedtools?

[tnminh89]

Yes I copied the command wrong, sorry! But the errors still remained