SEQanswers

Go Back   SEQanswers > Sequencing Technologies/Companies > Illumina/Solexa



Similar Threads
Thread Thread Starter Forum Replies Last Post
how to use the *.bam data to get the fasta file fay1313 RNA Sequencing 2 08-19-2013 10:09 PM
remove reads in fasta file JQL Bioinformatics 25 07-25-2013 06:16 AM
Gene reads count from RNA-seq sam file by CIGAR Jerry_Zhao Bioinformatics 8 02-01-2013 11:01 AM
piRNAs sequence files in fasta format for rat and mapping small RNA seq reads TJC Epigenetics 0 10-08-2012 02:11 PM
RNA-seq data analysis with 40bp reads heytreeful Bioinformatics 1 10-17-2011 09:30 PM

Reply
 
Thread Tools
Old 09-11-2013, 03:01 AM   #1
HSV-1
Member
 
Location: asia

Join Date: Jul 2012
Posts: 38
Smile How to Pick certain reads out of RNA-seq data to a fasta file?

Hello, every one.
I tried searching the solution to my question. But I couldn't . So put a new thread here. What I want to do is to pick out all the reads by Illumina Hiseq that contain the sequence "CACAGCTTCTAGTGCTATTCTGCGCCGGTATCC" and collect these reads into a fasta file or any other formats.

I have read through the manual of bowtie and BWA, I didn't find this usage.

Anyone could help me out?

Thanks!
HSV-1 is offline   Reply With Quote
Old 09-11-2013, 03:17 AM   #2
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Do you have the original fastq files or do you have the aligned SAM/BAM files? It sounds like you'll be best served by just using grep:

grep -B 1 -A 2 CACAGCTTCTAGTGCTATTCTGCGCCGGTATCC reads.fastq > matches.fastq

or

samtools view aligned.bam | grep CACAGCTTCTAGTGCTATTCTGCGCCGGTATCC > matches.sam

You could trivially pipe the second one through awk/perl/whatever to get output in fastq or fasta, if you prefer.
dpryan is offline   Reply With Quote
Old 09-11-2013, 03:21 AM   #3
JamieHeather
@jamimmunology
 
Location: London

Join Date: Nov 2012
Posts: 96
Default

Like dpryan said, you could get fasta output like this:

Code:
grep --no-group-separator -B 1 CACAGCTTCTAGTGCTATTCTGCGCCGGTATCC FILE.fq | sed 's/@/>/' > NEW_FILE.fa

Last edited by JamieHeather; 09-11-2013 at 03:28 AM. Reason: Taking too long to reply
JamieHeather is offline   Reply With Quote
Old 09-11-2013, 03:24 AM   #4
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Quote:
Originally Posted by JamieHeather View Post
You could do this in bash like this:

Code:
grep --no-group-separator -B 1 CACAGCTTCTAGTGCTATTCTGCGCCGGTATCC FILE.fq | sed 's/\@/\>/' > NEW_FILE.fa
I always forget about the --no-group-separator flag!

It should be noted that this flag isn't always available, in which case grep -v can be used to remove the separator (this issue has come up a couple times on this forum).
dpryan is offline   Reply With Quote
Old 09-11-2013, 03:30 AM   #5
JamieHeather
@jamimmunology
 
Location: London

Join Date: Nov 2012
Posts: 96
Default

Ooh I didn't realise the --no-group-separator isn't always present, good to know.

I thought the -v flag gave all non-matches? At least it seems to in my version.
JamieHeather is offline   Reply With Quote
Old 09-11-2013, 03:31 AM   #6
HSV-1
Member
 
Location: asia

Join Date: Jul 2012
Posts: 38
Default

Quote:
Originally Posted by dpryan View Post
Do you have the original fastq files or do you have the aligned SAM/BAM files? It sounds like you'll be best served by just using grep:

grep -B 1 -A 2 CACAGCTTCTAGTGCTATTCTGCGCCGGTATCC reads.fastq > matches.fastq

or

samtools view aligned.bam | grep CACAGCTTCTAGTGCTATTCTGCGCCGGTATCC > matches.sam

You could trivially pipe the second one through awk/perl/whatever to get output in fastq or fasta, if you prefer.

Thanks for your reply. I have the original rna-seq data in *.fastq format.
is grep the linux function? what are -B and -A 2 for?

Thank you guys!
HSV-1 is offline   Reply With Quote
Old 09-11-2013, 03:36 AM   #7
JamieHeather
@jamimmunology
 
Location: London

Join Date: Nov 2012
Posts: 96
Default

You're welcome!

Grep is indeed a unix function, you can just enter it straight into the command line.

Grep finds and returns the line(s) that contain the expression or string you supplied (in this case your sequence), then the A and B flags tell grep to also output lines above and below the matching line.

So -B 1 -A 2 would give you the line above the sequence line and the two lines beneath, thus outputting the entire fastq record.

The code I gave you output just the one above, as fasta doesn't use lines 3 and 4, and then just used sed to change the '@' character to a '>' in the ID line.

Grep and sed are very useful for hands on playing with your data like this, they're both worth getting to grips with.

Last edited by JamieHeather; 09-11-2013 at 03:39 AM.
JamieHeather is offline   Reply With Quote
Old 09-11-2013, 03:39 AM   #8
HSV-1
Member
 
Location: asia

Join Date: Jul 2012
Posts: 38
Thumbs up

Quote:
Originally Posted by JamieHeather View Post
You're welcome!

Grep is indeed a unix function, you can just enter it straight into the command line.

Grep finds and returns the line that contains the expression or string you supplied (in this case your sequence), then the A and B flags tell grep to also output lines above and below the matching line.

It's a very powerful tool, worth getting to grips with.
Thanks. thank you very much.
HSV-1 is offline   Reply With Quote
Old 09-11-2013, 03:45 AM   #9
dpryan
Devon Ryan
 
Location: Freiburg, Germany

Join Date: Jul 2011
Posts: 3,480
Default

Quote:
Originally Posted by JamieHeather View Post
Ooh I didn't realise the --no-group-separator isn't always present, good to know.

I thought the -v flag gave all non-matches? At least it seems to in my version.
It does, you just have to pipe the first grep into a grep -v and then the output of that into awk/sed/whatever. I don't recall which platform was missing the --no-group-separator option, I'd have to search back through the forums.
dpryan is offline   Reply With Quote
Old 09-11-2013, 04:04 AM   #10
JamieHeather
@jamimmunology
 
Location: London

Join Date: Nov 2012
Posts: 96
Default

Quote:
Originally Posted by dpryan View Post
It does, you just have to pipe the first grep into a grep -v and then the output of that into awk/sed/whatever. I don't recall which platform was missing the --no-group-separator option, I'd have to search back through the forums.
Ahh I see, that makes sense!
JamieHeather is offline   Reply With Quote
Reply

Tags
pick out, rna-seq

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 06:11 PM.


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