Filtering bam file with bedtools intersect using strand specificity

Question

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"    +   .

Answer 1

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

Answer 2

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