SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Dindel crash whilst parsing bam-sorted file michalkovac Bioinformatics 8 02-16-2012 07:56 PM
badly sorted BAM Filippo Bioinformatics 3 12-29-2011 01:39 PM
how to check whether a bam fille is sorted using picard in java jay2008 Bioinformatics 0 05-23-2011 04:14 PM
Sorted bam wangzkai Bioinformatics 3 05-07-2010 02:37 AM

Reply
 
Thread Tools
Old 09-28-2011, 09:18 PM   #1
dustar1986
Junior Member
 
Location: Sydney

Join Date: Sep 2011
Posts: 8
Default How to get all contig boundaries from a sorted bam file

Hi everyone,

Is there a smart way to assemble reads in the deep sequencing data?
After aligning reads to the genome using bowtie and samtools, I got the bam file. I want to extend the reads to get contigs if these reads have the same ref chromosome and are on the DNA strand.
Can I get the start position, end position, ref strand and chromosome number of all these contigs?

I wrote the following code but it ignore the strand (+/-) the reads come from. Reads from - strand might connected to a + strand one. That's not good.
Is that possible to extract ref strand info from samtools mpileup?

My code:

samtools mpileup input.sorted.bam |\
cut -d ' ' -f 1,2 |\
awk -F ' ' '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

And I really want to achieve this in the shell. Any suggestion?

Thanks a lot.

Dadi
dustar1986 is offline   Reply With Quote
Old 09-29-2011, 01:39 AM   #2
maubp
Peter (Biopython etc)
 
Location: Dundee, Scotland, UK

Join Date: Jul 2009
Posts: 1,543
Default

Ignore me: Thought you were working from the SAM/BAM file directly:
For the strand you'll have to parse the FLAG field, and restricting yourself to shell script makes this much more difficult than it need be. Is using a scripting language really not an option for you?
maubp is offline   Reply With Quote
Old 09-29-2011, 05:17 PM   #3
dustar1986
Junior Member
 
Location: Sydney

Join Date: Sep 2011
Posts: 8
Default

Quote:
Originally Posted by maubp View Post
Ignore me: Thought you were working from the SAM/BAM file directly:
For the strand you'll have to parse the FLAG field, and restricting yourself to shell script makes this much more difficult than it need be. Is using a scripting language really not an option for you?
No No. Using a scripting language, especially Python, is fine for me. Any suggestion on this please?
I just thought samtools might be better to deal with sequencing data coz is more "professional" than a script...
dustar1986 is offline   Reply With Quote
Old 09-30-2011, 01:31 AM   #4
maubp
Peter (Biopython etc)
 
Location: Dundee, Scotland, UK

Join Date: Jul 2009
Posts: 1,543
Default

You can use the samtools C API from a script, e.g. for Python there is PySam,
http://code.google.com/p/pysam/
maubp is offline   Reply With Quote
Reply

Tags
contigs, samtools

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 12:50 AM.


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