bedtools intersect -c flag output
0
0
Entering edit mode
9.2 years ago
bruce.moran ▴ 970

I have been using BEDtools intersect to determine on-target proportions of reads from BAM files in a BED. I used the -c flag which specifies that "For each entry in A, report the number of hits in B" which is exactly what I want. Output is in the form:

chr1    9999    10062   HISEQ:529:C71NKANXX:3:1305:12386:7195/2 0 ...

So in this case, the read does not overlap any feature in my BED file. I looked into this as colleagues were querying the method and found that the output from BEDtools stops at 60. Is this a local maximum to save time? There is nothing about this in documentation.

Also while I am asking, is there a better way to determine on-target reads and bases? Intersect seemed to be able to give me both (using -wao flag and increment on col 16 for bases). Picard HSmetrics is good but lacks on-target reads values.

intersect BEDtools • 3.2k views
ADD COMMENT
0
Entering edit mode

I do not think there is such limit:

printf 'chr1\t1\t100\n%.0s' {1..70} > A.bed
printf 'chr1\t20\t80\n%.0s' {1..1} > B.bed
bedtools intersect -a B.bed -b A.bed -c

Output:

chr1    1    100    70
ADD REPLY
0
Entering edit mode

I did test as you have and get no limit, using the -abam flag for BAM input may output MAPQ, limit of which is 60. Haven't had time to make a 'fake' BAM to tet if this is what is output, had hoped someone else came across this issue and could clarify.

Thanks for you input.

ADD REPLY
0
Entering edit mode

Bedtools intersect has -abam option ?

ADD REPLY
0
Entering edit mode

If you want to get the on target reads, then maybe

samtools view -h input.bam -L region.bed | samtools view -bSh - > output.bam

That should give you all the reads in the required region (e.g. on target reads)

ADD REPLY
0
Entering edit mode

I used:

bedtools intersect -a input.bam -b region.bed -v

and just divide by total reads to get off-target, which gives on-target obviously (half time of samtools method).

Thanks for suggestion though.

ADD REPLY

Login before adding your answer.

Traffic: 1665 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6