SEQanswers

Go Back   SEQanswers > Sequencing Technologies/Companies > SOLiD



Similar Threads
Thread Thread Starter Forum Replies Last Post
extracting unmapped reads from BAM using Samtools? Lspoor Bioinformatics 17 08-25-2013 12:22 AM
Extracting (Slicing) a BAM file ashkot Bioinformatics 4 12-13-2011 11:02 AM
Extracting reads from BAM files alpesh Bioinformatics 1 10-12-2011 02:59 PM
Unique Reads from BAM/SAM Gen2007 Bioinformatics 6 01-25-2011 08:29 AM
Extracting Reference Sequence from a bam File andy11 Bioinformatics 6 12-13-2010 02:06 PM

Reply
 
Thread Tools
Old 05-27-2010, 08:01 AM   #1
JohnK
Senior Member
 
Location: Los Angeles, China.

Join Date: Feb 2010
Posts: 106
Default Extracting unique reads from a .ma or .bam file?

Hi all,

I wanted to know the most efficient way to extract uniquely mapped reads from BioScope's resequencing mapping pipeline output in either the .ma or .bam format (those are the files I have to work with)? By unique, I mean reads that mapped to one location alone. Thanks!

J
JohnK is offline   Reply With Quote
Old 05-27-2010, 08:22 AM   #2
drio
Senior Member
 
Location: 4117'49"N / 24'42"E

Join Date: Oct 2008
Posts: 323
Default

.ma way:

They keep all the coordinates for all the mapping positions available. You can use that
to drop reads that map to multiple locations.
Actually, they should have another file (.ma.unique.*) that does that for you already.

.bam way:

Look to see if they have any tag (NH for example) to indicate how many alignments where done for that read.

You can post some records here we can take a look.
__________________
-drd
drio is offline   Reply With Quote
Old 05-27-2010, 10:04 AM   #3
JohnK
Senior Member
 
Location: Los Angeles, China.

Join Date: Feb 2010
Posts: 106
Default

Quote:
Originally Posted by drio View Post
.ma way:

They keep all the coordinates for all the mapping positions available. You can use that
to drop reads that map to multiple locations.
Actually, they should have another file (.ma.unique.*) that does that for you already.

.bam way:

Look to see if they have any tag (NH for example) to indicate how many alignments where done for that read.

You can post some records here we can take a look.

I was looking for the .unique.ma, but I didn't find it or any options with BioScope after mapping my resequencing sample. I understand that Corona produced that automatically. From what I understand, unique reads in bioscope are a matter of score > 5.



>2_29_53_F3,9_-110926286.224.2.15):q7
T33333203331202321313113130023213300030321303003300
>2_29_134_F3
T23331313312033202331333300333110330110312303003300
>2_29_155_F3
T11322113002233330002013020033212300020300303003302
>2_29_160_F3
T11331323300033211333223330333013320210301303003300
>2_29_200_F3,17_-90327533.241.5.0):q21,8_-43471518.229.3.0):q0
T01030013331133333333303333033303333330333303303300
>2_29_217_F3
T13330023311233323333300302133123220020302303002300



******

Is there a way to pick out a uniquely mapped read from a .bam file with SAM/Bam tools?
JohnK is offline   Reply With Quote
Old 05-27-2010, 11:31 AM   #4
drio
Senior Member
 
Location: 4117'49"N / 24'42"E

Join Date: Oct 2008
Posts: 323
Default

Quote:
Originally Posted by JohnK View Post
Is there a way to pick out a uniquely mapped read from a .bam file with SAM/Bam tools?
No within the mandatory fields. That's why I was asking if you can post here
some SAM entries from Bioscope so we can take a look.
__________________
-drd
drio is offline   Reply With Quote
Old 05-27-2010, 02:21 PM   #5
lh3
Senior Member
 
Location: Boston

Join Date: Feb 2008
Posts: 693
Default

http://sourceforge.net/apps/mediawik...from_SAM.2FBAM.
lh3 is offline   Reply With Quote
Old 05-28-2010, 06:03 AM   #6
drio
Senior Member
 
Location: 4117'49"N / 24'42"E

Join Date: Oct 2008
Posts: 323
Default

He wants to get the unique mapped reads.
But yes I agree with this (from the FAQ):

Quote:
We prefer to say an alignment is `reliable' rather than `unique' as `uniqueness' is not well defined in general cases. You can get reliable alignments by setting a threshold on mapping quality
__________________
-drd
drio is offline   Reply With Quote
Old 05-28-2010, 08:49 AM   #7
JohnK
Senior Member
 
Location: Los Angeles, China.

Join Date: Feb 2010
Posts: 106
Default

Quote:
Originally Posted by drio View Post
He wants to get the unique mapped reads.
But yes I agree with this (from the FAQ):
I had some trouble converting between bam back to sam. It appears to be incorrect. Here's my usage:

samtools view in.bam > name.sam

Probably not right usage here. Also, I'm re-formatting the matobam.ini file so it calcs 'unique' (or more reliable) reads based on the QV > 5. If anybody knows if this threshold is optimal, your help would be greatly appreciated? I'll then use bed tools to convert to bed and get the coverage off the unique. Hopefully, this will work and I can get a rough estimate of concordance between the unique reads from Corona versus the unique reads from BioScope.
JohnK is offline   Reply With Quote
Old 05-28-2010, 11:21 AM   #8
drio
Senior Member
 
Location: 4117'49"N / 24'42"E

Join Date: Oct 2008
Posts: 323
Default

Your samtools cmd is correct. What error/problem are you getting?
__________________
-drd
drio is offline   Reply With Quote
Old 05-28-2010, 01:13 PM   #9
JohnK
Senior Member
 
Location: Los Angeles, China.

Join Date: Feb 2010
Posts: 106
Default

Quote:
Originally Posted by drio View Post
Your samtools cmd is correct. What error/problem are you getting?
17_1177_738 16 chr1 2999999 7 25H25M * 0 0 NNGAATTCTTTTCTATGATTTAGTT """IIIIIIIIIIIIIIIIIIIII! RG:Z:20100521202751742 NH:i:1 CS:Z:T30123003213322000220302023010232222031211331023323 CQ:Z:BB?@;A@??B@=4>B<6<<A<?=?A=96=?>289>;7/=?04<0B3+(=> CM:i:2
774_1259_2007 16 chr1 2999999 7 25H25M * 0 0 NNGAATTCTTTTCTATGATTTAGTT """IIIIIIIIIIIIIIIIIIIII! RG:Z:20100521202751742 NH:i:1 CS:Z:T30123003213322000220302023010232222031211331023123 CQ:Z:B:AABA7BBBA?@?B=;A>B;<B:A?7=.@?2<;A:=<<?;<?5>:15?@ CM:i:2
JohnK is offline   Reply With Quote
Old 05-28-2010, 05:04 PM   #10
drio
Senior Member
 
Location: 4117'49"N / 24'42"E

Join Date: Oct 2008
Posts: 323
Default

So what is the problem? What were you expecting as a output?
__________________
-drd
drio is offline   Reply With Quote
Old 05-28-2010, 06:40 PM   #11
JohnK
Senior Member
 
Location: Los Angeles, China.

Join Date: Feb 2010
Posts: 106
Default

Quote:
Originally Posted by drio View Post
So what is the problem? What were you expecting as a output?
Nah- I realized it's correct now. It's just my first time really messing with and having to look at the SAM format meticulously. So this is what I'm expecting as output, but I just want to grab the unique reads. I'll keep reading over the manuals and what not and see if I can figure this out.
JohnK is offline   Reply With Quote
Old 05-29-2010, 12:40 PM   #12
JohnK
Senior Member
 
Location: Los Angeles, China.

Join Date: Feb 2010
Posts: 106
Default

Quote:
Originally Posted by drio View Post
So what is the problem? What were you expecting as a output?
Just a follow up, if anyone cares to know. You can extract unique (aka more 'reliable' reads) via the matobam.output.parameter='unique' setting. Then with BedTools convert from bam to bed and so forth... I appreciate the help, all!
JohnK is offline   Reply With Quote
Old 06-03-2010, 02:18 AM   #13
bair
Member
 
Location: London

Join Date: Jan 2010
Posts: 65
Default

My bam file includes 2x35bp mate-pair and 2x101 pair-end, any efficient way to split it into MP.bam and PE.bam files?
bair is offline   Reply With Quote
Old 06-03-2010, 11:46 AM   #14
drio
Senior Member
 
Location: 4117'49"N / 24'42"E

Join Date: Oct 2008
Posts: 323
Default

Quote:
Originally Posted by bair View Post
My bam file includes 2x35bp mate-pair and 2x101 pair-end, any efficient way to split it into MP.bam and PE.bam files?
If you are confortable with C or java, use samtools or picard to iterate over your bam and dump the data to PE or MP bam based on the read length.
That's the most efficient way.

There are other ports of samtools for lisp and python (maybe more) if you
want a more higher level language.

Another option would be to stream the bam data to sam with samtools and pipe it to your split script. After that you can convert back to BAM. It would be slower but it'd work just fine.

Something like:

Code:
#!/bin/bash
#
samtools view ./i.bam  | ruby -a -ne 'puts $_ if $F[9].size == 35' > 35.sam
samtools view ./i.bam  | ruby -a -ne 'puts $_ if $F[9].size == 120' > 120.sam 
samtools view -bS ./35.sam > 35.bam
samtools view -bS ./120.sam > 120.bam
__________________
-drd
drio is offline   Reply With Quote
Old 06-04-2010, 12:32 AM   #15
bair
Member
 
Location: London

Join Date: Jan 2010
Posts: 65
Default

Thank you Drio,

Your suggestions and script code are very helpful, like to use the simple code.

Thanks again.
bair is offline   Reply With Quote
Reply

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 08:09 PM.


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