Finding smaller sequences from within larger sequences

Question

I am currently working with fastq files which have hundreds of thousands of lines of text. However, not all of them are sequences I am interested in. My sequences are in one line and have a fixed length of 272 characters.

I possess a list of smaller sequences, around 2 million text strings in the following formats (which are fragments from the various fastq files I have):

Format 1

ATCGATCG
ATCGATCGATCG
ACGATCGTGTTACG
ATTTACGTACGTA
AATCG
ACGATACGATACG

Or I also have them in Format 2

ATCGATCGATACGAT|ATACGTGTTACGAT|ATCGATACG|ATACGATGA

(I just wrote random sequences just to give an idea of how the pattern.txt files look like, I know Format 1 and 2 are not the same in this post, however, the files are in the same order - only the syntax is changed.).

What I've been trying to do is use grep in a Linux environment and use the pattern.txt files to search for the bigger sequences in my fastq files where the pattern is present, and then output it into a .txt file to further work on.

I've used this command so far (with the Format 2 of patterns):

grep -E -f patterns.txt target.fastq > output.txt

Now, this works.

However, I have to limit the number of patterns in my key file and the number of fastq files I can search at once, otherwise I get "grep: memory exhausted" error.

I also want to print out all the sequences up to and including my pattern, but not past it (I don't care what's after my pattern sequence).

And for that I have successfully used the following:

grep -oE "\S+singlepattern" target.fastq > output.txt

However, I can only manage to get this to work if I manually write in the command the pattern sequence I want to look for (thus I can only use a couple). I can't seem to get it to work with the file list.

So I thought of using a while loop. This is what I've written so far in the terminal (it's one continuous line, but I formatted it here for simplicity on the eyes):

while read line; do
grep $line target.fastq > output.txt; 
done < patterns.txt

However, this does absolutely nothing and just spits out a blank output.txt

I have tried changing the syntax, with quotations on arguments, no quotations, I tried the various -E -f -o -F, egrep, etc but it's all the same. No error message, just blank file.

I've also tried implementing the \S to the loop as follows:

while read line; do
grep -oE "\S+$line" target.fastq > output.txt; 
done < patterns.txt

Yet, the result is the same.

Any suggestions on how I can fix this? Or even other methods to get the cut-down sequences (up to and including my pattern) from the fastq files into an output.txt file?

Thanks in advance!

Answer 1

Thank you all for your help. I managed to fix it in the end. The error was with how my pattern.txt file was formatted.

I initially exported the .txt from R Studio by specifying the EOL as \r\n. However, I think this was the issue as then from my laptop to the computer cluster the files seemed to have a windows EOL.

Thus on the linux server I did

dos2unix patterns.txt

And since then all of the suggested loops, both while read and the for chunk work beautifully and I get a nice \S+$line exported into the desired output.txt.

Thanks again all for your advice and help!

Answer 2

The reason your loop is giving a blank file is because you are using > which will overwrite whatever was previously in the file. This means that your output.txt will only have the result of the last grep and if that last grep has no results, you'll get an empty file.

Now, you're on the right track here. You could use a loop, although it will be horribly, horribly slow (by the way, you can use this exact format in the terminal as well, it doesn't need to be on the same line):

## create output.txt if it doesn't exist and empty it if it does exist
> output.txt 

while read line; do
    ## Use >> instead of > so you append to the file
    grep -oE "\S+$line" target.fastq >> output.txt; 
done < patterns.txt

Your attempt with the -f option is the most reasonable. You might get it to work if you tell grep to treat the input as fixed strings and not regular expressions (-F). Also, you don't need the -E (and that might be causing issues), so try this in case it works:

grep -Ff patterns.txt target.fastq > output.txt

If that doesn't work, just split your patterns.txt file into smaller files and then use grep in a loop with each of them:

split -l300000 patterns.txt pattern_chunk

That will split patterns.txt into multiple files of 3000000 lines each, named pattern_chunkaa, pattern_chunkab etc. For example, this is the result of running that command on a file with 5 million lines:

$ wc -l patterns.txt 
5000000 patterns.txt
$ split -l300000 patterns.txt pattern_chunk
$ ls
pattern_chunkaa  pattern_chunkad  pattern_chunkag  pattern_chunkaj  pattern_chunkam  pattern_chunkap
pattern_chunkab  pattern_chunkae  pattern_chunkah  pattern_chunkak  pattern_chunkan  pattern_chunkaq
pattern_chunkac  pattern_chunkaf  pattern_chunkai  pattern_chunkal  pattern_chunkao  patterns.txt

So, you can now use each of those as the patterns:

## create output.txt if it doesn't exist and empty it if it does exist
> output.txt 
for chunk in pattern_chunk*; do
    grep -f "$chunk" target.fastq >> output.txt; 
done

This doesn't help you get the sequence only up until the matched pattern. For that, you can add the \S+ to the beginning of each pattern in patterns.txt (this assumes you are using the file which has one pattern per line) and then split the file and loop over the chunks again:

## add the \S+
sed -i 's/^/\\S+/' patterns.txt
## create output.txt if it doesn't exist and empty it if it does exist
> output.txt 

for chunk in pattern_chunk*; do
    grep -Ef "$chunk" target.fastq >> output.txt; 
done

Note how I keep using > output.txt. This is not actually necessary, you can use command > file even if file doesn't exist. However, it is useful here to make sure that the output.txt file doesn't have the result of a previous run since we are now appending (>>) instead of overwriting (>) its contents.

Answer 3

I would recommend not do this using grep. There is a program within the KAT suite called KAT filter for filtering reads by a set of kmers (that might require all your sequences to be of constant length, but I am not sure). Another software you might want to check out is Cookiecutter. These tools specialised in subsetting sequencing files are usually a lot more efficient than generic regular expression tools.