SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Reads identifiers of reads spanning a SNP pm2012 Bioinformatics 5 04-09-2013 12:49 PM
First base of reads SNP or false result? pirates.genome Bioinformatics 2 03-08-2012 12:27 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
SNP Database with Solexa/Illumina Reads foolishbrat Illumina/Solexa 0 10-05-2009 01:54 AM

Reply
 
Thread Tools
Old 07-28-2014, 12:38 AM   #1
Phoenix_ICE
Member
 
Location: UK

Join Date: Jan 2011
Posts: 14
Default [Help] How to get those reads containing specified SNP?

Hi, all,

I am a new drummy for bioinformatics.

After SNP calling using GATK/freebayes, we usually get a SNP list. Now I have some interest SNP sites. Does anyone know how to identify those reads containing these interest SNPs?

Please note these SNP might be heterozygous. And now I mapped the reads to a reference, and get sorted bam file.

Would anyone tell me how to achieve that in detail or just tell me your thought and some tools might be helpful
Phoenix_ICE is offline   Reply With Quote
Old 07-28-2014, 07:28 AM   #2
westerman
Rick Westerman
 
Location: Purdue University, Indiana, USA

Join Date: Jun 2008
Posts: 1,104
Default

Assuming you have mapped your reads and now have a SAM/BAM file [this is the usual case] then the samtools program using the 'view' option will pull out reads in the region of your choice.
westerman is offline   Reply With Quote
Old 07-29-2014, 09:32 AM   #3
BurlEarl
Member
 
Location: Salt lake city, UT, US

Join Date: Jun 2012
Posts: 19
Default

Might not be understanding you but you can pull out all the reads + info with
grep -B 1 -A 2 GCCTATCGCAGATACACTCC sample.fastq > SNVreads.fastqish
(the nuc string contains your SNP)

need to remove the -- printed between reads
grep -v -e -- SNVreads.fastqish > SNVreads.fastq

You might have to tweek the length of your grep nuc pattern for specificity and avoiding other SNPs (dont know what you are sequencing). A couple cross platform visualization tools is Ugene.

Hope this is what you are looking for.

Earl
BurlEarl is offline   Reply With Quote
Old 07-30-2014, 01:36 AM   #4
Phoenix_ICE
Member
 
Location: UK

Join Date: Jan 2011
Posts: 14
Default

reference -----------------------------------------------------------
read1 ----------T-------------
read2 -------------------------
read3 ------T------------------
read4 --------------------------

I want to extract all the read id having the T snp
Phoenix_ICE is offline   Reply With Quote
Old 08-04-2014, 09:34 AM   #5
BurlEarl
Member
 
Location: Salt lake city, UT, US

Join Date: Jun 2012
Posts: 19
Default

If your read file looks like that then you can use

[your/Directory]$ grep -------T------ YourReadFile.txt > YourSNPReadFile.txt

output:
[your/Directory]$ more YourSNPReadFile.txt
read1 ----------T-------------
read3 ------T------------------

_________________________________________________________________________
If you have a .fastq file, all you need is the first line, which is just before the nuc string like:

@M01472:34:000000000-A40FG:1:1101:17765:1645 1:N:0:9
NTTCCAGCGAGGTTCTGAGTTCTTAGTCTGGTGTCGGCGTACCCACACGGTG
+
#>>>ABFFB?DBGGGGGCEGGGHHHGHHHHHFAGHEEGGGGGGHHGFDEEFG


just use:

[your/Directory]$ grep -B 1 GCCTATCGCAGATACACTCC YourSample.fastq > NamesAndReads.txt
#where "-B 1" prints the line before the pattern
#and the pattern "GCCTATCGCAGATACACTCC" contains the SNP somewhere in the middle.

[your/Directory]$ grep @M01472 NamesAndReads.txt > Names.txt
# "@M01472" is something in all the names but not in any reads
# for instance if your read names are actually read1, read2, read3, and read4 you could use "read"

#output for my command
[your/Directory]$ more Names.txt
@M01472:34:000000000-A40FG:1:1101:17765:1645 1:N:0:9
@M01472:34:000000000-A40FG:1:1101:18453:1656 1:N:0:9
@M01472:34:000000000-A40FG:1:1101:16266:1658 1:N:0:9
--More--(0%)

NOTE: this is a quick solution, if your genome is repetitive or if the SNP is in a duplicated region this approach might not be the best method. If that is the case. Something a little more involved from a .sam file might be necessary.

hope that helps
BurlEarl 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 10:48 PM.


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