Extract parts of reads from BAM file that overlap a specific region of genome

Question

[Chadi Saad]

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

[guillaume.rbt]

why do you want only a part of the reads?

[chiaraunivpm]

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

Answer 1

[Pierre Lindenbaum]

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

[Chadi Saad]

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 ?

[Pierre Lindenbaum]

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

[Chadi Saad]

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)...

[Pierre Lindenbaum]

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

[msumaira36]

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

Answer 2

[Prakki Rama]

@Chadi Saad @chiaraunivpm @msumaira36

This should do, what you want.

 library(GenomicAlignments) # Load the package

 stack <- stackStringsFromBam("test_4_seqs_sorted.bam", 
                              param=GRanges("Reference_barcodes:54-78")) # Extract the stacked reads from the position

 write.table(stack, Aligned_Reads_at_Region_Of_Interest.txt, quote = FALSE, col.names = F, row.names = F) # save into a text file

Source: Extract from a .bam file the reads that map to specific area in the reference genome