Answers.Body: <p>Your score and strand columns are swapped. If you correct that then you'll get stranded results.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-04-13T12:55:29.990Z
Answers.LastActivityDate: 2018-04-13T12:55:29.990Z
Answers.UserDisplayName: Devon Ryan
Answers.Body: <p>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 '.'</p>
  
  <pre><code>$ 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
  </code></pre>
  
  <p>After fixing the bed file:</p>
  
  <pre><code>$ 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
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-04-13T13:02:51.327Z
Answers.LastActivityDate: 2018-04-13T13:09:11.667Z
Answers.LastEditDate: 2018-04-13T13:09:11.667Z
Answers.UserDisplayName: Wolfgang Resch
Body: <p>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:</p>
  
  <pre><code>$ 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
  </code></pre>
  
  <p>Now I would like to enforce strand-specificity, but in either case, none of the records pass the filter:</p>
  
  <p>Same strand (<code>-s</code>):</p>
  
  <pre><code>$ bedtools intersect -a sample.bam -b features.bed -wa -u -f 1.0 -s | samtools view
  $ # nothing here
  </code></pre>
  
  <p>Opposite strand (<code>-S</code>):</p>
  
  <pre><code>$ bedtools intersect -a sample.bam -b features.bed -wa -u -f 1.0 -S | samtools view
  $ # nothing here
  </code></pre>
  
  <p><strong>What is the correct way of doing this?</strong></p>
  
  <h3>Edit</h3>
  
  <p>Bedtools version:</p>
  
  <pre><code>$ bedtools --version
  bedtools v2.27.1
  </code></pre>
  
  <p>My bed file has strand information:</p>
  
  <pre><code>$ 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"    +   .
  </code></pre>
  
Tags: <sam><bedtools>
Title: Filtering bam file with bedtools intersect using strand specificity
ThreadText: # Filtering bam file with bedtools intersect using strand specificity
  
  (Asked by: username_0 on 2018-04-12)
  
  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"    +   .
  
  ```
  
  **Tags:** sam, bedtools
  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_2 on 2018-04-13. Score: 2 | Accepted answer
  
  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 by username_1 on 2018-04-13. Score: 1
  
  Your score and strand columns are swapped. If you correct that then you'll get stranded results.
  
  
  ---
  
  