How to get all contig boundaries from a sorted bam file

[dustar1986]

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

[maubp]

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?

[dustar1986]

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

[maubp]

You can use the samtools C API from a script, e.g. for Python there is PySam,

http://code.google.com/p/pysam/