SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
bowtie SAM mapq field rgregor Bioinformatics 2 12-19-2012 05:04 PM
SAM Format - SEQ field '=' Bio.X2Y Bioinformatics 3 04-25-2012 05:26 AM
extracting certain scaffolds from SAM file lukas1848 Bioinformatics 1 01-18-2012 09:33 AM
snpEff and extracting EFF field kasthuri Bioinformatics 2 12-30-2011 01:32 PM
Is * really a valid value for a SAM FLAG field? derobins Bioinformatics 1 01-20-2011 10:06 AM

Reply
 
Thread Tools
Old 12-07-2009, 07:23 AM   #1
jdrum00
Member
 
Location: Houston, TX, USA

Join Date: Dec 2009
Posts: 16
Default Extracting one field from a SAM file

I just became a bioinformatician (is that what we're called?) a few weeks ago, and am surprised that I can't find anything in samtools that will output a subset of fields from a SAM/BAM file. It seems to be entirely geared toward filtering "horizontally", for subsets of alignment records. I know that's a major need, but it's not actually what I need at the moment!

What I'm looking for is a derivative file with, for example, just one long list of CIGAR strings. I don't care what reads or alignments they're from; I just want to do some statistics on them in the aggregate.

What I'll do in the meantime is just parse the SAM with python to make my derivative file. But surely samtools would be better at it than I am. Anyone know how to make it do what I want? Thanks in advance!

(I'm excited to see the new python wrapper for samtools, which I'll probably check out soon.)

ANSWER from below, thanks to nilshomer: Just operate directly on the SAM-formatted stream that comes out of samtools view. E.g., samtools view <bamname> | awk <whatever field you need> | <further processing>.

Last edited by jdrum00; 01-04-2010 at 09:50 AM. Reason: Added answer to original question, for the benefit of searchers.
jdrum00 is offline   Reply With Quote
Old 12-07-2009, 08:31 AM   #2
krobison
Senior Member
 
Location: Boston area

Join Date: Nov 2007
Posts: 747
Default

Doing that sort of slice is trivial from a SAM file with standard UNIX tools - cut, awk, perl, etc. For example "cut -f1" would pull out the first field. Count over (index origin 1) to CIGAR and put that value for the -f argument. Similarly, "perl -n -e '@a=split(/\t/); print $a[0],"\n"; (note that Perl is index origin 0, as most languages -- cut is wierd).

The binary BAM file would be a bit trickier; if you are more C-fluent than I you could use the decompress routines to decompress and then parse things out.

With the Perl (and probably Python; I'm just haven't dug into it) wrapper you can certainly iterate over all records & only grab the field(s) you want. Probably a significant overhead penalty, but it will get the job done.
krobison is offline   Reply With Quote
Old 12-07-2009, 08:44 AM   #3
jdrum00
Member
 
Location: Houston, TX, USA

Join Date: Dec 2009
Posts: 16
Default

Right; I can do that from a SAM, with awk or whatever. I was hoping someone would say there was an undocumented switch to samtools view or something, because that would also be likely to do it directly from BAM and save me a bunch of file conversion. It seems like the sort of simple direct thing that would be included in a basic tool set (but then I guess everyone says that about the feature they want).

Thanks for the details and suggestions!
jdrum00 is offline   Reply With Quote
Old 12-07-2009, 10:05 AM   #4
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

Quote:
Originally Posted by jdrum00 View Post
Right; I can do that from a SAM, with awk or whatever. I was hoping someone would say there was an undocumented switch to samtools view or something, because that would also be likely to do it directly from BAM and save me a bunch of file conversion. It seems like the sort of simple direct thing that would be included in a basic tool set (but then I guess everyone says that about the feature they want).

Thanks for the details and suggestions!
Unix programs typically follow the pipe-and-filter model so duplicating the function of one when a pipe would suffice is hearsay.

Here's how you get the cigar directly from a BAM file:
Code:
samtools view my.bam | cut -f6
Unix utilities will soon become your friends.

Last edited by nilshomer; 12-07-2009 at 07:29 PM.
nilshomer is offline   Reply With Quote
Old 12-07-2009, 07:19 PM   #5
jdrum00
Member
 
Location: Houston, TX, USA

Join Date: Dec 2009
Posts: 16
Default

Quote:
Originally Posted by nilshomer View Post
Unix utilities will soon become your friends.
Oh yes; we're already quite well acquainted. Sed, awk, grep, cut, sort and I are all good buddies. Bash sometimes joins us for a movie. I used to hang out with python a lot, but honestly, I haven't seen much of him these last couple of weeks.

Quote:
Originally Posted by nilshomer View Post
Unix programs typically follow the pipe-and-filter model so duplicating the function of one when a pipe would suffix is hearsay.
You're absolutely right, and my confusion was in not understanding that samtools view really is just a bam streamer. Now that that's clicked, I'm much happier!

This is a great forum; I'm looking forward to rummaging through it and learning a lot. Thanks for the useful and timely responses.
jdrum00 is offline   Reply With Quote
Old 01-04-2010, 02:58 AM   #6
KevinLam
Senior Member
 
Location: SEA

Join Date: Nov 2009
Posts: 203
Default

Quote:
Originally Posted by jdrum00 View Post
Right; I can do that from a SAM, with awk or whatever. I was hoping someone would say there was an undocumented switch to samtools view or something, because that would also be likely to do it directly from BAM and save me a bunch of file conversion. It seems like the sort of simple direct thing that would be included in a basic tool set (but then I guess everyone says that about the feature they want).

Thanks for the details and suggestions!
Tell me about it!
I bet there's a script floating out there already.

working on python script to tally total mapped reads to a fasta entry in the reference for SAM file.

will post the script if anyone is interested but I think everyone can do it.
any suggestions on how I can do it better? not a unix guru so hoping to learn from others.

I feel there;s gotta be a better way than to awk grep 3000 times a sam file if i have 3000 fasta entries.
is there a word freq count unix tool?
KevinLam is offline   Reply With Quote
Old 01-04-2010, 03:31 AM   #7
KevinLam
Senior Member
 
Location: SEA

Join Date: Nov 2009
Posts: 203
Default

btw here's the regex i got off the sam manual
HTML Code:
[CODE]
 Field    Regular expression  Range                                     Description
                                        Query pair NAME if paired; or Query NAME if unpaired 2
QNAME [^ \t\n\r]+
                             [0,216-1] bitwise FLAG (Section 2.2.2)
FLAG   [0-9]+
                                        Reference sequence NAME 3
RNAME [^ \t\n\r@=]+
                             [0,229-1] 1-based leftmost POSition/coordinate of the clipped sequence
POS    [0-9]+
                             [0,28-1]   MAPping Quality (phred-scaled posterior probability that the mapping
MAPQ   [0-9]+
                                        position of this read is incorrect) 4
                                        extended CIGAR string
CIGAR ([0-9]+[MIDNSHP])+|\*
                                        Mate Reference sequence NaMe; = if the same as <RNAME> 3
MRNM   [^ \t\n\r@]+
                             [0,229-1] 1-based leftmost Mate POSition of the clipped sequence
MPOS   [0-9]+
                             [-229,229] inferred Insert SIZE 5
ISIZE -?[0-9]+
                                        query SEQuence; = for a match to the reference; n/N/. for ambiguity;
SEQ    [acgtnACGTN.=]+|\*
                                        cases are not maintained 6,7
                                        query QUALity; ASCII-33 gives the Phred base quality 6,7
                             [0,93]
QUAL   [!-~]+|\*
                                        TAG
TAG    [A-Z][A-Z0-9]
                                        Value TYPE
VTYPE [AifZH]
                                        match <VTYPE> (space allowed)
VALUE [^\t\n\r]+[/CODE]
KevinLam is offline   Reply With Quote
Old 01-04-2010, 08:01 AM   #8
jdrum00
Member
 
Location: Houston, TX, USA

Join Date: Dec 2009
Posts: 16
Default

Quote:
Originally Posted by KevinLam View Post
working on python script to tally total mapped reads to a fasta entry in the reference for SAM file.
...
I feel there;s gotta be a better way than to awk grep 3000 times a sam file if i have 3000 fasta entries. is there a word freq count unix tool?
I'm not totally clear on what you need, but awk will count lines:

Code:
awk '{count[$i]++} END {for(x in count) print x,count[x]}'
The first bit tallies up everything as it goes through the file, and then the bit after END prints out the results.

So if there's a specific subportion of the SAM line you want to count, filter that out with awk or cut, and pipe that result to the awk above. For example, I needed to look at the custom MD: field, which is the sixteenth field, so:

Code:
samtools view <bamname> | awk '{print $16}' | awk '{count[$1]++} END {for(j in count) print count[j], j}' | sort -nr
This results in output like:

Code:
38 MD:Z:50
17 MD:Z:30
2 MD:Z:49A0
2 MD:Z:46G3
If you can get python or awk (or samtools!) to do the complex filtering you want, then awk is pretty blazingly fast at counting up the results. Hope that helps!

Last edited by jdrum00; 01-04-2010 at 08:03 AM. Reason: Note that samtools is a good filterer too
jdrum00 is offline   Reply With Quote
Old 01-04-2010, 09:40 PM   #9
KevinLam
Senior Member
 
Location: SEA

Join Date: Nov 2009
Posts: 203
Default

Neat!
Thanks jdrum
that's pretty much what I was looking for!
Was thinking of using a hash to count no. of reads that map to a chr in reference genome.
But this works better!
KevinLam 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 06:28 AM.


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