SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Pileup / extract information from BAM/SAM files liu_xt005 Bioinformatics 4 01-19-2015 09:14 AM
Per base sequence coverage from sam/bam file? ewilbanks Bioinformatics 7 06-06-2012 01:03 PM
Going from transcriptome to genome coordinates with a bam file pbluescript Bioinformatics 5 05-23-2012 09:37 AM
Extract perfectly mapped reads from SAM/BAM file Graham Etherington Bioinformatics 2 07-21-2011 07:27 AM
Sam flags for bwa-aligned paired end reads with identical + / - strand coordinates spark Bioinformatics 0 03-09-2011 04:00 AM

Reply
 
Thread Tools
Old 08-14-2012, 06:21 AM   #1
pirates.genome
Member
 
Location: Earth

Join Date: Mar 2012
Posts: 11
Question Extract aligned sequence coordinates from SAM or BAM file

Hi,
I have a Illumina sequence reads which I mapped to HG19 using BWA. Resulting BAM file is indexed and sorted. Now what I want is to identify the region of alignment. Let me give an example, on chromosome 1 starting from base 1 to 10000, there are many reads aligned and then there is no any read aligned from 10000 to 20000 and then again from 20000 to 30000 there are many reads aligned and so on...

I want an output like
chr1 1 10000
chr1 20000 30000
etc...etc...

I searched for tools doing similar functions but didn't find anything.
Any help will be really appreciated.

Thanks.
pirates.genome is offline   Reply With Quote
Old 08-15-2012, 03:44 AM   #2
colindaven
Senior Member
 
Location: Germany

Join Date: Oct 2008
Posts: 401
Default

Have a look at a variation of:

samtools view -L bedfile.bed x.bam
colindaven is offline   Reply With Quote
Old 08-16-2012, 05:38 AM   #3
pirates.genome
Member
 
Location: Earth

Join Date: Mar 2012
Posts: 11
Default

Thnx colindaven,

I looked into samtools view -L option but that is not what I am looking for.

Let me explain the problem,
I have targeted DNA sequencing data from Illumina. After aligning to reference using BWA I got BAM files. From BAM files I want to extract a continuous stretch of reference region to which multiple reads are aligned. Assume that reads are aligned in following way:
read1 aligned to chr1:1-50
read2 chr1:10-60
read3 chr1:30-80
read4 chr1:150-200
read5 chr1:200-250
read6 chr1:250-300
read7 chr1:350-400
and so on....

Now I want continuous stretch of reference to which any read is aligned. Output should be:
chr1:1-80
chr1:150-300
chr1:350-400
and so on...

This output tells us that there are several number of reads aligned to chr1 position 1 to 80 in continuation of each other and then no any read aligned to position 80 to 150 and so on...

Is there an easy way to do so..??..

Thanks.
pirates.genome is offline   Reply With Quote
Old 08-16-2012, 06:54 AM   #4
Richard Finney
Senior Member
 
Location: bethesda

Join Date: Feb 2009
Posts: 699
Default

SIMPLE VERSION ... using samtools depth and awk ...
#ploc = previous location, pchr=previous chromsome
#Set SAMT and BAMF to samtools and your bamfile.
export SAMT=/h1/finneyr/samtools-0.1.18/samtools
export BAMF=98023.bam
$SAMT depth $BAMF | awk '{ if ($2!=(ploc+1)){if (ploc!=0){printf("%s %d-%d\n",$1,s,ploc);}s=$2} ploc=$2; }


More complicated version that handles chr1->chr2 transition and flushes at the end. (I think, this is not completely debugged and and is just a one-off. Not all corner conditions maybe addressed.)

$SAMT depth $BAMF | \
awk '
BEGIN{firsttime=1;}
{
if (pchr!=$1) { if (firsttime==1) { firsttime = 0;} else { printf("%s %d-%d\n",pchr,s,ploc);}s=$2}
else { if ($2!=(ploc+1)){if (ploc!=0){printf("%s %d-%d\n",$1,s,ploc);}s=$2} }
ploc=$2; pchr=$1
}
END{ printf("%s %d-%d\n",pchr,s,ploc);}
'

CAVEAT: samtools depth doesn't do cigar 'N' quite right so this won't work for RNA in the best way.

Last edited by Richard Finney; 08-16-2012 at 06:58 AM.
Richard Finney is offline   Reply With Quote
Old 08-18-2012, 11:07 AM   #5
pirates.genome
Member
 
Location: Earth

Join Date: Mar 2012
Posts: 11
Smile

Quote:
Originally Posted by Richard Finney View Post
SIMPLE VERSION ... using samtools depth and awk ...
#ploc = previous location, pchr=previous chromsome
#Set SAMT and BAMF to samtools and your bamfile.
export SAMT=/h1/finneyr/samtools-0.1.18/samtools
export BAMF=98023.bam
$SAMT depth $BAMF | awk '{ if ($2!=(ploc+1)){if (ploc!=0){printf("%s %d-%d\n",$1,s,ploc);}s=$2} ploc=$2; }


More complicated version that handles chr1->chr2 transition and flushes at the end. (I think, this is not completely debugged and and is just a one-off. Not all corner conditions maybe addressed.)

$SAMT depth $BAMF | \
awk '
BEGIN{firsttime=1;}
{
if (pchr!=$1) { if (firsttime==1) { firsttime = 0;} else { printf("%s %d-%d\n",pchr,s,ploc);}s=$2}
else { if ($2!=(ploc+1)){if (ploc!=0){printf("%s %d-%d\n",$1,s,ploc);}s=$2} }
ploc=$2; pchr=$1
}
END{ printf("%s %d-%d\n",pchr,s,ploc);}
'

CAVEAT: samtools depth doesn't do cigar 'N' quite right so this won't work for RNA in the best way.
Hi Richard Finney,
Thank you for very helpful code. First code worked good for me.
pirates.genome is offline   Reply With Quote
Old 08-20-2012, 08:06 AM   #6
ishmael
Member
 
Location: NY, US

Join Date: Jul 2008
Posts: 17
Default

I think BEDtools should also fit your demand.
1. convert the bam into bed by bam2bed
2. mergeBed.
ishmael is offline   Reply With Quote
Reply

Tags
aligned region, bam, sam

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 01:13 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2018, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO