Split BAM into 1 BAM per contig

Question: Split BAM into 1 BAM per contig

5.6 years ago by

Marvin • 170

Australia

Marvin • 170 wrote:

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 • 2.3k views

ADD COMMENT • link •

modified 5.6 years ago by Pierre Lindenbaum ♦ 133k • written 5.6 years ago by Marvin • 170

5.6 years ago by

Pierre Lindenbaum ♦ 133k

France/Nantes/Institut du Thorax - INSERM UMR1087

Pierre Lindenbaum ♦ 133k wrote:

use genomecov ( http://bedtools.readthedocs.org/en/latest/content/tools/genomecov.html ) 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

ADD COMMENT • link written 5.6 years ago by Pierre Lindenbaum ♦ 133k

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

ADD REPLY • link modified 5.6 years ago • written 5.6 years ago by Marvin • 170

I'm struggling, how do I achieve this ... !?

ADD REPLY • link written 5.6 years ago by Marvin • 170

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?

ADD REPLY • link written 5.6 years ago by Marvin • 170

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!

ADD REPLY • link written 5.6 years ago by Marvin • 170

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?

ADD REPLY • link written 15 months ago by DNAngel • 90

Similar posts • Search »