How to Pick certain reads out of RNA-seq data to a fasta file?

[HSV-1]

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!

[dpryan]

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.

[JamieHeather]

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

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

[dpryan]

You could do this in bash like this:

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).

[JamieHeather]

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.

[HSV-1]

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!

[JamieHeather]

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.

[HSV-1]

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.

[dpryan]

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.

[JamieHeather]

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!