This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Split BAM into 1 BAM per contig

Hi,

I have mapped reads against a reference with bwa and got a big bam file now.

I would like to retrieve a separate bam file for every contig.

I do not mean "I want 1 bam for every chromosome". A chromosome is not a contig since my coverage is 0 in many parts of the genome (aDNA reads).

In this situation e.g. (I hope he doesn't mess up when I post this):

--------------------------------------
  ---   ---               ---   ---
    ----                    -----

I would like to get 2 bam files.

How do I do that?

sam contig
• 4,253 views • link
updated by Ram Ram • written by Marvin Marvin

1 answer

Use genomecov (http://bedtools.readthedocs.org/en/latest/content/tools/genomecov.html) to get a bed with >0 coverage. Then, loop over the bed to extract the reads.

awk '{printf("%s:%s-%s\n",$1,$2,$3);}' cov.bed | while read F; do samtools view -o "${F}.bam" -b your.bam "${F}"; done
• 0 views • link
updated by Ram Ram • written by Pierre Lindenbaum Pierre Lindenbaum

Thank you! Don't know how to use that bedtools genomecov but I will find out.

• 0 views • link
updated by Ram Ram • written by Marvin Marvin

I'm struggling, how do I achieve this?

• 0 views • link
updated by Ram Ram • written by Marvin Marvin

I just can't figure it out. The problem is that he opens a new bed line as soon as the coverage within a contig changes. But what I need is 1 line per contig. How is that possible with a single genomecov command? I know how to do it with a pipe but from your answer it sounds like it could be done immediately?

• 0 views • link
written by Marvin Marvin

Oh my god I finally made it.

Thank you for helping me alot by pointing me towards the right direction (bedtools is the right choice [didn't know that program before], then awk with samtools as posted by you). However that genomecov suggestion gave me struggles because it didn't get me what I wanted ... all I need is:

bedtools merge -i aln_sorted.bam

This gives me 1 line per contig. Afterwards I can apply the awk-samtools-pipe that you have posted :)

Thanks again!

• 0 views • link
written by Marvin Marvin

Can you share the code you used and the files you used as input? I would also like to obtain individual bam files per contig. I tried this process posted here but I am getting zero output files after running the awk command (without genomecov).

I instead used bedtools merge -i aln_sorted.bam > mycontigs.txt then used the awk command and instead of cov.bed I used mycontigs.txt

Do I need a bed file at all for this?

• 0 views • link
written by DNAngel DNAngel

Log in to answer this question.

More posts like this