Extract reads with a specific mutation

[desmo]

Hi everyone,
I would like to extract from a bam file all the reads that are carrying a specific mutation.
I know that samtools give the possibility to extract reads in a specific region, but if I would like to extract just these one that have a mutation what can I do?
I would really appreciate any kind of help
Thanks in advance

[maasha]

I am not quite sure, if this is what you want, but with Biopieces (www.biopieces.org) you can read in a SAM file and grab for a specific position, and then a specific event like this:

read_sam -i input.sam |
find_SNPs |
grab -e 'POS == 1234' |
grab -e 'EVENT eq 'G>C'

Unfortunately Biopieces cannot write SAM format (which is evil), but you can dump a table with write_tab.

[desmo]

Thanks Maasha,
your suggestion was really good. I've installed biopieces and it works great but I have one more question. Do you know how to run biopieces from a perl script? I need it to solve my problem. I don't know if you use perl. If yes do you know why it doesn't work with the command system?

$cmd='read_sam -i FileIN.txt -O FileOUT.txt';
#system($cmd);

[maasha]

That should work just fine. This does:

perl -e '$cmd = "read_sam -n 1 -i test.sam | count_records"; system($cmd);'

[desmo]

Thanks. It works great. And thanks for let me know about biopieces, it looks really helpful.