SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
How to replace select reads in a bam file? Heisman Bioinformatics 8 01-02-2012 02:49 PM
Getting Reads from Bam file empyrean Bioinformatics 4 10-12-2011 12:57 PM
Obtaining reference identical reads from a BAM file Sakti Bioinformatics 2 05-17-2011 10:40 AM
Query bam files, awk question NearyJL78 Bioinformatics 5 05-06-2011 10:58 AM
Extracting unique reads from a .ma or .bam file? JohnK SOLiD 14 06-04-2010 12:32 AM

Reply
 
Thread Tools
Old 09-27-2011, 07:21 PM   #1
dustar1986
Junior Member
 
Location: Sydney

Join Date: Sep 2011
Posts: 8
Default Query bam file and assemble reads

Hi,

I am quite a newbie in Bioinformatic and maybe asking something stupid.

I am trying to generate a four-column file for assembling overlapped sequencing reads into longer contigs from a sorted bam file. The file needs to contain the following information in this format:

chr#: start position of alignment - stop position of alignment: strand (+/-)

Currently I am trying to awk the bam files to obtain the information, my code is here:

samtools mpileup input.sorted.bam |\
cut -d '\t' -f 1,2 |\
awk -F '\t' 'BEGIN {chr="";start=-1;end=-1} {if(chr!=$1 || int($2)!=end+1) { if(chr!="") {printf("%s:%d-%d\n",chr,start,end);} chr=$1;start=int($2);end=int($2);} else { end=end+1;}} END {if(chr!="") {printf("%s:%d-%d\n",chr,start,end); } }'>output.txt

It can work except with strand info. Actually, it ignores the strand info. If a read from + strand overlapped with a read from - strand, it will form a contig and that's not what I want. I want to assemble contigs in the same strand.

How can I improve my code to take in strand info and make the assembly according to strand?

Please help. Thank you very much.

Dadi Gao
dustar1986 is offline   Reply With Quote
Old 09-28-2011, 09:19 AM   #2
macrowave
Member
 
Location: New York

Join Date: May 2010
Posts: 13
Default

I think Bio:B::Sam is the way to go, if you know a bit perl programming. You can use its methods to access sam/bam for any information, including map positions and strands.
macrowave is offline   Reply With Quote
Old 09-28-2011, 04:17 PM   #3
dustar1986
Junior Member
 
Location: Sydney

Join Date: Sep 2011
Posts: 8
Default

Quote:
Originally Posted by macrowave View Post
I think Bio:B::Sam is the way to go, if you know a bit perl programming. You can use its methods to access sam/bam for any information, including map positions and strands.
Unfortunately I know nothing about perl
Is there another way to do it directly from linux shell?
dustar1986 is offline   Reply With Quote
Old 09-29-2011, 04:32 AM   #4
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

Quote:
Originally Posted by dustar1986 View Post
Unfortunately I know nothing about perl
Is there another way to do it directly from linux shell?
This would be a great time to learn perl/python/a-scripting-language.
nilshomer is offline   Reply With Quote
Old 09-29-2011, 09:34 AM   #5
macrowave
Member
 
Location: New York

Join Date: May 2010
Posts: 13
Default

Maybe you can dig something out from the sam bitwise flags, I'm sure the strand information would be encoded there. For details, you can see the sam specifications here:
http://samtools.sourceforge.net/SAM1.pdf. And there is an easy tool to explain the flags here: http://picard.sourceforge.net/explain-flags.html
macrowave is offline   Reply With Quote
Old 09-29-2011, 09:44 AM   #6
cbeck
Junior Member
 
Location: Montreal

Join Date: Sep 2011
Posts: 6
Default

You might also try the bedtools suite - bamtobed and then mergebed. Not sure if that will deal with the strands,you might be able to filter the two strands at the bam output level.
cbeck is offline   Reply With Quote
Old 09-29-2011, 04:10 PM   #7
dustar1986
Junior Member
 
Location: Sydney

Join Date: Sep 2011
Posts: 8
Default

Quote:
Originally Posted by nilshomer View Post
This would be a great time to learn perl/python/a-scripting-language.
I can use python actually...
dustar1986 is offline   Reply With Quote
Old 09-29-2011, 04:12 PM   #8
dustar1986
Junior Member
 
Location: Sydney

Join Date: Sep 2011
Posts: 8
Default

Quote:
Originally Posted by macrowave View Post
Maybe you can dig something out from the sam bitwise flags, I'm sure the strand information would be encoded there. For details, you can see the sam specifications here:
http://samtools.sourceforge.net/SAM1.pdf. And there is an easy tool to explain the flags here: http://picard.sourceforge.net/explain-flags.html
Thanks. But the flag only appear in samtools -view not in samtools -mpileup.
I want to use the latter one coz it's easier for assembly.
dustar1986 is offline   Reply With Quote
Old 09-29-2011, 04:12 PM   #9
dustar1986
Junior Member
 
Location: Sydney

Join Date: Sep 2011
Posts: 8
Default

Quote:
Originally Posted by cbeck View Post
You might also try the bedtools suite - bamtobed and then mergebed. Not sure if that will deal with the strands,you might be able to filter the two strands at the bam output level.
Thanks, I'll have a look at that.
dustar1986 is offline   Reply With Quote
Old 09-29-2011, 04:53 PM   #10
quinlana
Senior Member
 
Location: Charlottesville

Join Date: Sep 2008
Posts: 119
Default

If you install bedtools, you could use the "mergeBed" program to merge overlapping reads on the same strand as follows:
Code:
bamToBed -i aln.bam | mergeBed -i stdin -s > merged.intervals.bed
Alternatively, you could use the "genomeCoverageBed" tool to create strand-specific BEDGRAPH files, one for each strand. You could then combine them into a single file.
Code:
bamToBed -i aln.bam | \
          awk '$6=="+" | \
          genomeCoverageBed -i stdin -g chrom.sizes -bg \
          > merged.intervals.fwd.bedg

bamToBed -i aln.bam | \
          awk '$6=="-" | \
          genomeCoverageBed -i stdin -g chrom.sizes -bg \
          > merged.intervals.rev.bedg

cat merged.intervals.fwd.bedg merged.intervals.rev.bedg > merged.intervals.both.bedg
quinlana is offline   Reply With Quote
Old 09-29-2011, 07:48 PM   #11
cgjkjk
Member
 
Location: Roskilde

Join Date: Oct 2008
Posts: 14
Default

BAM file is compressed, you need do "samtools view -u -bh aln.bam > my.bam" first.
cgjkjk is offline   Reply With Quote
Reply

Tags
awk question, query bam files

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 07:36 AM.


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