Retrieving reads with SNPs

[Gavin_Sherlock]

Hi,

Does anyone know of a way to straightforwardly retrieve all reads (from resequencing) that have a SNP present in them? We're generating vcf files, using bwa+gatk, and would then like to know all reads which contain a SNP listed in the vcf file. I know that I can write a wrapper that uses samtools view for a region containing a SNP, and then look at each read individually (the SNPs will be heterozygous) that maps to that location, but I wondered whether there is something that is likely faster than this, as I'd likely implement the wrapper in Perl, which will be slow for thousands of SNPs x many strains.

Cheers,
Gavin

[RockChalkJayhawk]

Hi,

Does anyone know of a way to straightforwardly retrieve all reads (from resequencing) that have a SNP present in them? We're generating vcf files, using bwa+gatk, and would then like to know all reads which contain a SNP listed in the vcf file. I know that I can write a wrapper that uses samtools view for a region containing a SNP, and then look at each read individually (the SNPs will be heterozygous) that maps to that location, but I wondered whether there is something that is likely faster than this, as I'd likely implement the wrapper in Perl, which will be slow for thousands of SNPs x many strains.

Cheers,
Gavin

check the fields in the gatk vcf. It already tells you how many reference/alternate reads there are.

[cpendlebury]

I believe he wants the actual sequence of the reads with SNPs, rather than just the numbers of SNPs.

If so I'm wanting to figure out how to do the same thing.