Query bam file and assemble reads

[dustar1986]

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

[macrowave]

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.

[dustar1986]

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?

[nilshomer]

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.

[macrowave]

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

[cbeck]

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.

[dustar1986]

This would be a great time to learn perl/python/a-scripting-language.

I can use python actually...

[dustar1986]

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]

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.

[quinlana]

If you install bedtools, you could use the "mergeBed" program to merge overlapping reads on the same strand as follows:

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.

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

[cgjkjk]

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