Reads identifiers of reads spanning a SNP

[pm2012]

Hi

I am using Samtools to call SNPs in number of samples. Is there any way I can extract the reads (read Ids) and specify how many reads support the SNP an how many don't? Any help is appreciated as I am a newbie.

[swbarnes2]

Hi

I am using Samtools to call SNPs in number of samples. Is there any way I can extract the reads (read Ids) and specify how many reads support the SNP an how many don't? Any help is appreciated as I am a newbie.

Something like

samtools view file.bam chr3:1234567-1234567 > SNP.sam

Will get you just the .bam entries that cross the SNP at 1234567

[pm2012]

Thanks for your reply. Yes I am able to isolate/identify all the reads at a specific position. However, not all the reads spanning the SNP position will have the SNP (e.g in heterologous SNPs, where half of reads will have SNP and half of them won't). I was asking if it would be possible to isolate all the reads (read IDs) spanning the SNP loci and identify which ones have the SNP and which ones don't.

[vivek_]

You could likely do that with a short script after dumping the reads into a file using the command that swbarnes provided.

However the SNP caller might still not consider some reads as supporting evidence of a SNP depending on the mapping quality, base quality at that position in a read etc, so it gets complicated at this point.

[swbarnes2]

If you make a vcf, the DP4 value will tell you how many high quality reads support the reference or alternate letter.

[pm2012]

Thanks all for your replies. I might have to just write a script to do this as I don't see a way to parse out read id from the mpileup output. In my case I want to identify the exact reads (their read ids) that support or not support SNPs.