Sort specific reads from fastq

Question

[xuyang_dami]

Hi, guys.

I wanna sort out reads containing specific sequences and discard reads without from fastq. However, I am new to computer programming with little knowledge about writing any codes.So I am wondering if there is any tool or codes that someone wrote which can perform this task.

Thanks

Yang

Answer 1

[Pierre Lindenbaum]

Using a simple shell script. If you're running windows, install cygwin.

gunzip -c fastq.gz |\
paste - - - - |\
awk '(!($2 ~ "ATGC"))' |\
tr "\t" "\n" |\
gzip > out.fastq.gz

• gunzip uncompress the source if needed.
• paste: convert 4 consecutive lines into one line with four columns *awk filter out column n°2 containing ATGC
• tr: convert tabulation to 'return' = convert back to fastq
• compress the output if needed and redirect to the output file "out.fastq.gz"

[igor]

Actually, awk is space separated, so if the FASTQ contains spaces (in the sequence ID potentially), it will see more than four columns of data. You should add awk -F $'\t' 'BEGIN {OFS = FS} ...'

[igor]

Nice paste trick! Is that documented somewhere?

[dariober]

This is also one of my favourite tricks. Is it documented? Yes, all the commands Pierre is using behave as per documentation, including paste (from man pages: If `-' is specified for one or more of the input files, the standard input is used). It's amazing how simple commands can be stitched together to accomplish complex jobs.

[igor]

I didn't mean to imply that I don't believe the answer. I just wanted to figure out how I missed this option for so long.

Yes, it says "-" means STDIN, but to me that means that there would be 4 columns of STDIN since it's repeated 4 times. The part that's surprising to me is that the next time it's used, it skips to the next line, so it's the same STDIN, not just another iteration of it. It makes sense in retrospect, but it's not necessarily obvious.

[dariober]

Hi igor-

It makes sense in retrospect, but it's not necessarily obvious

It wasn't obvious to me at all as well! I saw that trick used somewhere else but just by reading the docs I would have never guessed it.

[Pierre Lindenbaum]

I learned it on biostars.org, years ago. https://www.google.com/search?q=paste+fastq+site%3Abiostars.org

[xuyang_dami]

Thanks. very simple and smart solution. Can this be applied to bam file?

[Pierre Lindenbaum]

In SAM, yes but you must consider the SAM header too.

[vetbishnu_tfeed]

Dear Pierre Does this code needs to add awk -F $'\t' 'BEGIN {OFS = FS} ...' as suggested by igor or works without adding.

Thank you Bishnu

Answer 2

[WouterDeCoster]

Python, using Biopython:

import sys
from Bio import SeqIO
sequences_of_interest = ["ATG", "GGCACAC", "GTGACCC"] #To be adapted
handle = open(sys.argv[1], "rU")
filtered = []
for record in SeqIO.parse(handle, "fastq")
    for seq in sequences_of_interest:        
        if record.seq.count(seq) > 0:
             filtered.append(record)
             break
output_handle = open("seqlected.fasta", "w")
SeqIO.write(filtered, output_handle, "fastq")
output_handle.close()
handle.close()

Save script (e.g. selectSEQ.py) and execute as python selectSEQ.py yourfile.fastq

Let me know if you need additional help.

[Pierre Lindenbaum]

are you putting a whole fastq in memory ?

[WouterDeCoster]

Hmm, you are right. Depending on the size of the fastq and the memory available, putting the file entirely in memory might not be the best idea.

[xuyang_dami]

Thanks. One question. whether tolerance to mismatch and deletion is counted?

Answer 3

[igor]

There is FASTQ-GREP for filtering for sequences matching a pattern: http://homes.cs.washington.edu/~dcjones/fastq-tools/fastq-grep.html

Answer 4

[chen]

If you are able to use Julia, you can use OpenGene(https://github.com/OpenGene/OpenGene.jl) to do this very easily

using OpenGene

instream = fastq_open("your.fq.gz")
outstream = fastq_open("out.fq.gz", "w")

while (fq = fastq_read(instream)) != false
    if contains(fq.sequence.seq, "Your pattern")
        fastq_write(outstream, fq)
    end
end

close(outstream)

If you have no idea about Julia, just do
1, run sudo apt-get install julia in ubuntu
2, run julia after the installation is finished
3, change the filename and your pattern in the code
4, paste your code into the Julia interactive command line and press ENTER

[xuyang_dami]

Nice! thank you so much