Biostar Beta. Not for public use.
Extract parts of reads from BAM file that overlap a specific region of genome
2
Entering edit mode
23 months ago
Chadi Saad • 60
France

Assuming that I have a BAM file with aligned reads, how can I extract the part of all reads that cover completely a specific region of a genome. (I tried with samtools view file.bam chr1:13974-22442 , but it reports all reads that intersect with this region, also it reports all the read sequence not only the intersection part...)

For example : (needed region to be extract between pipes)

=============|=====|========================= < Reference
    ---------|-----|--------------------------
-------------|-----|-----------
                         --------------------------------------------
           --|-----|--------------------------

Thank you

ADD COMMENTlink
0
Entering edit mode

why do you want only a part of the reads?

ADD REPLYlink
0
Entering edit mode

Hello Chadi Saad, I have the same request...I want to extract/trim the reads from bam file that exactly match with specific regions of a bed file...Did you get it?

Chiara

ADD REPLYlink
2
Entering edit mode
17 months ago
France/Nantes/Institut du Thorax - INSE…

using samjs :

samtools view  -bu -F 4  file.bam chr1:13974-22442 |\
java -jar dist/samjs.jar  -e 'record.alignmentStart <= 13974 && record.alignmentEnd >= 22442'

see also : Count reads within region

ADD COMMENTlink
0
Entering edit mode

thanks , worked perfectly. But it outputs all the reads sequence, while i want just the part of read that match in the required region. Can i do it in a simple way ?

ADD REPLYlink
1
Entering edit mode
ADD REPLYlink
0
Entering edit mode

It does what i want, but it have some bugs i think. Sometimes the read is clipped too much (see attachement screenshot, in purple the original read, in red the clipped one, it stops at 14776 instead of 14780). Sometimes reads are shifted by one position (see screenshot 2)...

enter image description here enter image description here

ADD REPLYlink
0
Entering edit mode

nice to know. please , send an issue to https://github.com/lindenb/jvarkit/issues w a minimal SAM file and a BED file.

ADD REPLYlink
0
Entering edit mode

I am also facing the same issue. I want the part of reads that match a specific region. How did you do it ?? Please help

ADD REPLYlink

Login before adding your answer.

Similar Posts
Loading Similar Posts
Powered by the version 2.3.1