How to extract and cut sequence from a given region from a BAM file

Question

[nbargues]

Hi,

I have an aligned BAM file and I want to extract sequence from a given region and just that, not what's come after or before the specific region.

I try samtools with:

samtools view -b alignment_sort_JX888104.bam "JX888104.1:569-895" > filt_aln_grp1_JX888104.bam

but it give me all the read that comes in this interval even if it's one base.

After that I try bedtools:

bedtools intersect -wa -a alignment_sort_JX888104.bam -b grp1.bed > test_intersect_grp1.bam

with

grp1.bed -> JX888104.1 569 895

but it give me an error that my bed is not recognize.

Do you have another idea?

Thanks in advance

[Pierre Lindenbaum]

see Limiting variant calls to amplicon target regions? ; tool for bam soft clipping reads within bed file regions

[Pierre Lindenbaum]

but it give me an error that my bed is not recognize.

what is the output of

file  grp1.bed
tr "\t" "#" < grp1.bed

?

[nbargues]

ASCII text

JX888140.1 569 895

[Pierre Lindenbaum]

you're not using a tab as the column separator but spaces.

Output should have been

ASCII text
JX888140.1#569#895

[Pierre Lindenbaum]

furthermore, if your bam interval was "JX888104.1:569-895" , then your bed should have been

JX888104.1  568 895

see What Are The Advantages/Disadvantages Of One-Based Vs. Zero-Based Genome Coordinate Systems

[nbargues]

you're right, I thought that with atom tab was imply but apparently not, thanks. I fix that and I don't have the bedtools error anymore but when I look the data on IGV, I still have all my reads.

am I doing something wrong with my bedtools command ?

Cheers