Question: Extract parts of reads from BAM file that overlap a specific region of genome
2
gravatar for Chadi Saad
2.4 years ago by
Chadi Saad60
France
Chadi Saad60 wrote:

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

sequencing bam alignment • 1.8k views
ADD COMMENTlink modified 2.4 years ago by Pierre Lindenbaum124k • written 2.4 years ago by Chadi Saad60

why do you want only a part of the reads?

ADD REPLYlink written 2.4 years ago by guillaume.rbt730

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 written 2.4 years ago by chiaraunivpm0
2
gravatar for Pierre Lindenbaum
2.4 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum124k wrote:

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 modified 2.4 years ago • written 2.4 years ago by Pierre Lindenbaum124k

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 written 2.4 years ago by Chadi Saad60
1

see this post: Limiting variant calls to amplicon target regions?

ADD REPLYlink written 2.4 years ago by Pierre Lindenbaum124k

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 modified 2.4 years ago • written 2.4 years ago by Chadi Saad60

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 written 2.4 years ago by Pierre Lindenbaum124k

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 written 8 months ago by msumaira360
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1692 users visited in the last hour