Perhaps my Google-fu is failing me, because this doesn't seem like an uncommon need, but I can't find any answer here.
Say I have a SAM/BAM file, and a variant at a known location. I'd like to extract all the alignments from the BAM file which contain that variant. That is, I'd like the ALT, not the REF allele. I do not want all the reads at that location. Is there an existing tool or workflow out there that can accomplish this?
I would honestly be surprised if there wasn't, since I would imagine it's a common situation to be interested in a certain variant and investigate the reads supporting it. Specifically, I'd like to make sure there's no bias in where along the read it falls, or the strand of the read. I should note that my coverage is very high, so it's not feasible to extract all the reads covering the site and manually sort them out.
I've already started a script that does this by parsing CIGAR strings, but I wanted to check that it doesn't already exist before going further down this road.
•
written
by
Nick Stoler
See if this oneliner helps: For example: at given position say "1:1660404" get reads with Variant from the bam file:
samtools view -b INPUT.bam 1:1660404-1660404 | samtools fillmd -e - REFERENCE.fasta | grep -v "^@"| awk -v pos="1660404" 'BEGIN {OFS = FS = "\t" } ; {n=split($10,a,"") ; if(a[(pos-$4)+1] != "=" ) print pos,(pos-$4)+1, a[(pos-$4)+1], $1, $4, $10 }'
Potential output:
Position Distance_from_read_start Variant READ_NAME Read_Start READ
1660404 31 G HI-SW1050R:121:D2C53BCXX:8:2211:2278:38508 1660374 ==============================G===================================================
1660404 30 G HI-SW1050R:121:D2C53BCXX:8:2308:3759:48204 1660375 =============================G====================================================
1660404 25 G HI-SW1050R:121:D2C53BCXX:2:2208:13607:73456 1660380 =======C================G=========================================================
1660404 25 G HI-SW1050R:121:D2C53BCXX:6:2201:3216:84210 1660380 =======C================G=========================================================
1660404 25 G HI-SW1050R:121:D2C53BCXX:8:1316:1264:91116 1660380 =======C================G======N==N===================N=N=========================
1660404 16 G HI-SW1050R:121:D2C53BCXX:2:1301:15715:89711 1660389 ===============G==================================================================
Note: READ base changed to '=' if identical to reference base: it doesn't show indels and soft-clipped bases
That's good to know about fillmd's -e option, but unfortunately it seems it doesn't show indels. Also, soft-clipped bases appear just like SNV's.
My previous answer might help: position of mismatches per read from a sam/bam file : you'll get the content of each read (as XML) for each base.
Thanks! So it looks like that can get the information I need. It's a custom script like I've been writing (except certainly much nicer and faster), instead of an existing tool, so I'll take that as an answer to whether I'm reinventing the wheel.
Hey Nick,
Did you ever figure out how to extract reads with a specific variant? I'm looking to do something similar.
Thanks!
Log in to answer this question.
Maybe something like this:
@PhilS: your solution won't help parsing the cigar string.
Sry I thought you need the whole information of the read given by a bam file. If you need the cigar string only you can do it like:
this will give you the cigar strings of all reads containing 'pattern'
no, you have to walk over the cigar string to see if the read contains or not the variation. You cannot use a regex for this.