1
$\begingroup$

From the bedtools intersect man (for version 2.30.0):

-a BAM/BED/GFF/VCF file “A”. Each feature in A is compared to B in search of overlaps. 
   Use “stdin” if passing A with a UNIX pipe.

The stdin part, however, does not apply to version 2.27.1. The help on this version gives no advice for piping.

bam mem ref fq1.fastq fq2.fastq | samtools view -b | bedtools intersect -a stdin -b blacklist.bed -v > align_black.sam
Error: unable to open file or unable to determine types for file stdin

- Please ensure that your file is TAB delimited (e.g., cat -t FILE).
- Also ensure that your file has integer chromosome coordinates in the
  expected columns (e.g., cols 2 and 3 for BED).

I also tried:

(map to bam) | bedtools intersect -a "stdin" -b blacklist.bed -v > align_black.sam
(map to bam) | bedtools intersect -a -b blacklist.bed -v > align_black.sam
(map to bam) | bedtools intersect -b blacklist.bed -v > align_black.sam

But these don't work either, either for the same or a different reason. Here is a site that tells me to use the above syntax: https://bedtools.readthedocs.io/en/latest/content/example-usage.html

What is the solution to this if I want to pipe everything together for my users?

$\endgroup$
3
$\begingroup$

With v2.27.1 (or any version I ever used) instead of stdin use -

(...) | bedtools intersect -a - -b blacklist.bed -v > out.bam

Be aware that the output of this is a BAM, not SAM file. It would be easier though to simply make a list of regions you want to keep, e.g. the complement of the backlist file with the entire genome (bedtools complement), and then use samtools view -L option to only keep overlapping reads. That saves you the run through bedtools.

$\endgroup$
2
  • $\begingroup$ Thanks, I wonder why this isn't in the docs. I had my own workaround before you commented, but I'll accept this as the answer. Is there a reason you suggest not using bedtools? $\endgroup$
    – Jeff
    Apr 19 at 13:59
  • $\begingroup$ @Jeff As samtools can do the same and you anyway run it you save one additional tool/pipe/core. That is the only reason. Not sure how the performance compares. $\endgroup$
    – ATpoint
    Apr 20 at 7:34
1
$\begingroup$

Here's a way to skip conversion to BAM and pipe between tools, which performs the same set operations:

$ bam mem ref fq1.fastq fq2.fastq | sam2bed - | bedops -e 1 - blacklist.bed > blacklisted_reads.bed

With this toolkit, the hyphen (-) follows a Unix convention of serving as a placeholder for stdin.

$\endgroup$
1
$\begingroup$

While I have chosen ATpoint as the official answer, another possible syntax to use is

bedtools intersect -a stdin -b blacklist.bed -v < `(map to bam)` > align_black.bam

The main issue with this syntax is that it obfuscates the direction of flow for the pipeline, piping step 2 after step 1 is better in the official answer.

$\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.