0
$\begingroup$

I have a small sample of records from a bam file. I obtained them by intersecting with a bed file not enforcing strand-specificity. So I get them all back upon performing the same intersection again:

$ bedtools intersect -a sample.bam -b features.bed -wa -u -f 1.0 | samtools view
MN00481:17:000H2FH5N:1:11104:12480:15984    16  III 4499483 42  28M *   0   0   CAGGAACACGCACATCGTCGATTTTCAG    FFFFFAFFFFFFFFFFFFFFFFFFFFFF    AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:28 YT:Z:UU
MN00481:17:000H2FH5N:1:11103:9925:11522 16  IV  2388448 42  30M *   0   0   CACTGGAACTGGTCCGATCTTCAAGCCGGT  FFFFFFFFFFA/FFFFAFFFFFF/FFAFAF  AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:30 YT:Z:UU
MN00481:17:000H2FH5N:1:22109:11939:5243 16  IV  4749336 42  29M *   0   0   ATGTGGATGCGGTCAGCTTGGGCACGATG   FFFFFFFFFFFF/FFFF6FFFFFFFF/FF   AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:29 YT:Z:UU
MN00481:17:000H2FH5N:1:22104:17385:19431    16  V   4443064 42  29M *   0   0   CAAAATGTGAACAGGAACGTTTTCATGAG   F/FFFFFFFFFFFFFFFFFFFFFFFFFFF   AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:29 YT:Z:UU
MN00481:17:000H2FH5N:1:23105:22058:12968    0   X   981524  42  29M *   0   0   CCATCGATTCTCATTTTAATGCTGGTGCT   FFFFFFFFFFFFFFFFFFFFFFFFFFFFF   AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:29 YT:Z:UU
MN00481:17:000H2FH5N:1:12104:7867:7500  16  X   4214381 42  29M *   0   0   GCTGCAGGAATGACGGCTCTAGCAAACGT   FFFFFFFFFFFFFFFFFFFFFFFFFFFFF   AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:29 YT:Z:UU

Now I would like to enforce strand-specificity, but in either case, none of the records pass the filter:

Same strand (-s):

$ bedtools intersect -a sample.bam -b features.bed -wa -u -f 1.0 -s | samtools view
$ # nothing here

Opposite strand (-S):

$ bedtools intersect -a sample.bam -b features.bed -wa -u -f 1.0 -S | samtools view
$ # nothing here

What is the correct way of doing this?

Edit

Bedtools version:

$ bedtools --version
bedtools v2.27.1

My bed file has strand information:

$ tail features.bed 
X   17617826    17619458    "WBGene00015867"    +   .
X   17622777    17640743    "WBGene00019211"    -   .
X   17650177    17651594    "WBGene00044771"    +   .
X   17662457    17665357    "WBGene00003744"    +   .
X   17671028    17680790    "WBGene00020732"    +   .
X   17680910    17682294    "WBGene00044772"    -   .
X   17695231    17701744    "WBGene00017625"    +   .
X   17710485    17711674    "WBGene00019189"    -   .
X   17714455    17714832    "WBGene00007069"    +   .
X   17714967    17718717    "WBGene00007068"    +   .
$\endgroup$
2
$\begingroup$

Your bed file switch the score and strand columns. Strand should be column 6, score column 5. For example take this bed file with the 5th and 6th columns switched - zero results for strand specific intersection b/c strand is actually '.'

$ cat wrong.bed
chr1    30000000        31000000        noname  +       .
$ bedtools intersect -a mm10_500k_pe.bam -b wrong.bed -wa -u -f 1.0 \
    | samtools view | wc -l
166
$ bedtools intersect -a mm10_500k_pe.bam -b wrong.bed -wa -u -f 1.0 -s \
   | samtools view | wc -l                              
0

After fixing the bed file:

$ cat right.bed
chr1    30000000        31000000        noname  .       +
$ bedtools intersect -a mm10_500k_pe.bam -b right.bed -wa -u -f 1.0 \
      | samtools view | wc -l
166
$ bedtools intersect -a mm10_500k_pe.bam -b right.bed -wa -u -f 1.0 -s \
      | samtools view | wc -l
83
$\endgroup$
1
$\begingroup$

Your score and strand columns are swapped. If you correct that then you'll get stranded results.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.