Separating a .sam or .bam file that is based on alignment to multiple sequences

[sequence_hard]

Hey people!
I have the following setup. I have a .sam file based on the alignment of E.coli genome reads to a multifasta file containing the sequences of several plasmids. Now I want to separate the .sam file into several sam files, each containing the header and mapping information for one plasmid, so that I can later calculate the % coverage per plasmid.

I would do it such that I create a dictionary in python, append the header of a sequence as key and the corresponding sequence (identifiable because the key should also be present in the line of the alignment/mapping block). Then I write that to a file, convert to .bam and so on.

Now my questions: Is this an appropriate way to do what I want to do or would I lose information? Is there a better way to do this when the file is in .bam format?

Thanks for your help!

[GenoMax]

Unless you want to write your own code you can use "bamtools split" (or other programs) as noted in this thread: https://www.biostars.org/p/46327/

[WhatsOEver]

Which program did you for aligning your reads? If you have used a multi-fasta as reference, you should normally have the sequence headers of this multi-fasta as header in your sam file. Likewise, you should have the sequence ids as "reference id" in the sam line. Then, its quite straightforward to use the sam/bam file as it is and compute the coverage of the individual reference seqs.

[sequence_hard]

@WhatsOEver: I used bowtie2. But how can I calculate the coverage of each reference sequence using only my file?

@GenoMax: Thanks, that looks usefull!

[WhatsOEver]

Either you use a different mapper which has this capability (eg bbmap does this with the scafstats parameter) or (what I prefer) you use a program like bedtools genomecov (eg bedtools genomecov -ibam yourAlignedData.bam -g yourMultiFasta.fasta). The function of bedtools is nicely explained here: http://bedtools.readthedocs.org/en/l...genomecov.html

[GenoMax]

Samtools depth may work in a pinch. You could also use Qualimap to get visual maps of coverage.

[sequence_hard]

@WhatsOEver: Awesome, thanks! I am already using bedtools genomecov but I did not know it had that function. Nice!

[WhatsOEver]

@WhatsOEver: Awesome, thanks! I am already using bedtools genomecov but I did not know it had that function. Nice!

You are welcome
As to my experience, there is pretty much nothing you can't do with bedtools when it comes to sequence coverage.