This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Get sequence from aligned BAM file according to interval in BED file

Hello,

I have a BED file which define an interval of 3'UTR. The example of my BED is like this:

7   127591299   127591705   ENSG00000004059 0   +
12  8940364 8941817 ENSG00000003056 0   -
11  64315966    64316738    ENSG00000173153 0   +
12  2803258 2805423 ENSG00000004478 0   +

The order is : Chromosome Start End Gene Name Length Strand

I need to extract sequence in aligned BAM with corresponding interval per gene name. Is there any tools to do that? Thank you.

Update :

I use samtool view to see specific sequence in the given interval. So, according to my BED file from the first line, I use this command :

samtools view 33_realigned.bam 7:127591429-127591705

And it gives subset of BAM (in SAM format) corresponds to that interval but not merged. Is it possible to get the merged sequence? Thank you

alignment sequence
• 7,147 views • link
updated by george.ry george.ry • written by bharata1803 bharata1803

It looks like you are trying to make consensus sequence out of a bam file.

• 0 views • link
written by GouthamAtla GouthamAtla

2 answers

Samtools itself will do that for you (the new versions, anyway):

samtools view -L your.bed -b your.bam > subset.bam
• 0 views • link
written by george.ry george.ry

Make your BED a tab delimited file and use bedtools. For this task, you would use

intersectBed -abam aligned.bam -b ToExtract.bed > subset.bam
• 0 views • link
written by venu venu

Thank you for your suggestion. I tried it and it seems the result is different from what I expect. The bam result converted to fastq doesn't give me what I need. So, I want an output which will consist of the sequence corresponds to the BED. For example for BED with line: 

7   127591299   127591705   ENSG00000004059 0   +

I want to get string of ACGT between those region. So, probably the output is like this :

7   127591299   127591705   ENSG00000004059 0   + ACTGACAGT..

Is it possible?

The result I want is like if I use

bedtools getfasta -fi fasta.in -bed bed.in -fo csv.out -tab -s -fullHeader

I don't have the fasta.in and all I have is bam file.

• 0 views • link
written by bharata1803 bharata1803

Log in to answer this question.

More posts like this