SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Write the subset of reads from BAM file into new SAM/BAM file, using R tools. Old Pioneer Bioinformatics 0 01-27-2016 04:41 AM
Extract multiple fasta sequences from a fasta file based on sequenes entomology Bioinformatics 38 12-19-2015 06:28 PM
Drag-and-drop alignment of sequences to form SAM file Alex8 Bioinformatics 0 05-05-2011 12:10 AM
BAM/SAM to a gapped multiple sequence alignment query Bioinformatics 6 04-06-2011 05:07 AM
BWA: specifying SAM/BAM file header fields before read alignment? nora Bioinformatics 3 12-04-2010 09:11 PM

Reply
 
Thread Tools
Old 02-09-2016, 05:10 AM   #1
sequence_hard
Junior Member
 
Location: Sweden

Join Date: Feb 2016
Posts: 5
Default Separating a .sam or .bam file that is based on alignment to multiple sequences

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!
sequence_hard is offline   Reply With Quote
Old 02-09-2016, 05:20 AM   #2
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,937
Default

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/
GenoMax is offline   Reply With Quote
Old 02-09-2016, 05:22 AM   #3
WhatsOEver
Senior Member
 
Location: Germany

Join Date: Apr 2012
Posts: 215
Default

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.
WhatsOEver is offline   Reply With Quote
Old 02-09-2016, 05:33 AM   #4
sequence_hard
Junior Member
 
Location: Sweden

Join Date: Feb 2016
Posts: 5
Default

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

@GenoMax: Thanks, that looks usefull!

Last edited by sequence_hard; 02-09-2016 at 05:41 AM.
sequence_hard is offline   Reply With Quote
Old 02-09-2016, 05:48 AM   #5
WhatsOEver
Senior Member
 
Location: Germany

Join Date: Apr 2012
Posts: 215
Default

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
WhatsOEver is offline   Reply With Quote
Old 02-09-2016, 05:54 AM   #6
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,937
Default

Samtools depth may work in a pinch. You could also use Qualimap to get visual maps of coverage.
GenoMax is offline   Reply With Quote
Old 02-09-2016, 05:54 AM   #7
sequence_hard
Junior Member
 
Location: Sweden

Join Date: Feb 2016
Posts: 5
Default

@WhatsOEver: Awesome, thanks! I am already using bedtools genomecov but I did not know it had that function. Nice!
sequence_hard is offline   Reply With Quote
Old 02-09-2016, 06:05 AM   #8
WhatsOEver
Senior Member
 
Location: Germany

Join Date: Apr 2012
Posts: 215
Default

Quote:
Originally Posted by sequence_hard View Post
@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.
WhatsOEver is offline   Reply With Quote
Reply

Tags
mapping, multi fasta, sam bam

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 04:52 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2019, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO