SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Anyone know how Snowshoes-FTD finds junction-spanning reads? thondeboer Bioinformatics 2 09-04-2013 11:55 PM
SNP calling reads number wanguan2000 Bioinformatics 1 11-25-2011 08:16 AM
SNP Calling for 454 reads empyrean Bioinformatics 1 11-17-2011 07:23 AM
Tophat: find junction spanning reads thurisaz RNA Sequencing 4 11-14-2011 04:23 AM
Supporting Reads in SAMTools SNP Calls Lee Sam Bioinformatics 2 07-09-2010 06:16 AM

Reply
 
Thread Tools
Old 04-08-2013, 06:49 AM   #1
pm2012
Member
 
Location: DC

Join Date: Apr 2012
Posts: 18
Default Reads identifiers of reads spanning a SNP

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.
pm2012 is offline   Reply With Quote
Old 04-08-2013, 08:25 AM   #2
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

Quote:
Originally Posted by pm2012 View Post
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
Code:
samtools view file.bam chr3:1234567-1234567 > SNP.sam
Will get you just the .bam entries that cross the SNP at 1234567
swbarnes2 is offline   Reply With Quote
Old 04-08-2013, 08:42 AM   #3
pm2012
Member
 
Location: DC

Join Date: Apr 2012
Posts: 18
Default

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.

Last edited by pm2012; 04-08-2013 at 08:51 AM. Reason: more clarification added
pm2012 is offline   Reply With Quote
Old 04-08-2013, 08:59 AM   #4
vivek_
Bioinformatician
 
Location: Denmark

Join Date: Jul 2012
Posts: 158
Default

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.
vivek_ is offline   Reply With Quote
Old 04-08-2013, 12:18 PM   #5
swbarnes2
Senior Member
 
Location: San Diego

Join Date: May 2008
Posts: 912
Default

If you make a vcf, the DP4 value will tell you how many high quality reads support the reference or alternate letter.
swbarnes2 is offline   Reply With Quote
Old 04-09-2013, 12:49 PM   #6
pm2012
Member
 
Location: DC

Join Date: Apr 2012
Posts: 18
Default

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.
pm2012 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 04:45 PM.


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