1
$\begingroup$

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!

$\endgroup$
6
  • $\begingroup$ Could you please add some more context to explain specifically why you want to use grep? There are lots of sequence mapping programs that are specifically designed to quickly look for DNA sequences within other DNA sequences (both with and without errors), and there's a good chance your problem could be solved by using one of them. $\endgroup$
    – gringer
    Jun 5, 2021 at 4:03
  • $\begingroup$ Searching 2 million patters (m) over hundred of thousands line file(s) (n) is costly O(mn) and can be memory-wise too. I would first de-duplicate my list of patterns, if luckily they are highly duplicated, this can improve performance greatly. grep loads one line at a time in memory, if you use the format 2 as input for patterns, you are effectively loading the whole thing into memory at once and thus the memory issue can present. I would therefore use Format 1, or run it in a server with high memory available $\endgroup$
    – JRodrigoF
    Jun 5, 2021 at 9:44
  • $\begingroup$ Hello! Thank you for your advice. I produced the pattern file from R Studio and initially they were 2 million individual lines, however, I managed to output a new pattern.txt file with just unique patterns and now it's down to 700k. I do get the "grep: memory exhausted" nonetheless, thus why I'm trying to create a loop :/ $\endgroup$
    – Jean Luc
    Jun 5, 2021 at 11:42
  • $\begingroup$ @gringer I was recommended grep by a colleague, as well as the fact that most of my desired files are on a linux computer cluster. My laptop is not capable of running the required programming and I have zero knowledge of other software/programs as this is my first time doing dry lab work. Any recommendations would be appreciated and I can have a look nonetheless and see if I can use them! $\endgroup$
    – Jean Luc
    Jun 5, 2021 at 12:04
  • $\begingroup$ @JeanLuc I was not able to reproduce the "grep: memory exhausted" error message. I am running in a server but I think PC should do fine too. I put together a file with repetitions of your Format 1 patters so that I have 1M total patterns, one per line (patterns.txt). I added the regex you need as a prefix to each pattern (eg. \S+ATCGATCG). Then I took a bzipped fastq file and did 'zgrep -o -f patterns.txt file.fastq.gz'. It works. Slow as hell but that is expected given the complexity of the task (each time you search one pattern across all lines in the fastq file, and again ..) $\endgroup$
    – JRodrigoF
    Jun 6, 2021 at 9:42

3 Answers 3

1
$\begingroup$

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!

$\endgroup$
2
  • 1
    $\begingroup$ For future reference, "\n" is a better line ending to use when creating files from within R/RStudio, especially on a Linux system. $\endgroup$
    – gringer
    Jun 7, 2021 at 1:02
  • $\begingroup$ Yep, that's why I asked if you could have Windows style line endings. Glad you sorted it out! $\endgroup$
    – terdon
    Jun 7, 2021 at 10:14
0
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ Hello! Thank you for the answer but it's still not working. The grep commands all work, but as soon as I write them within the while or for loops they don't work. To be sure, I also added do echo \$line or $chunk (whether it was while or for loop), and it shows that it's running through all the files. However, the output files are always blank EDIT: I wrote the echo both before the grep and after, and it both shows it's running. $\endgroup$
    – Jean Luc
    Jun 4, 2021 at 17:42
  • $\begingroup$ @JeanLuc have you ever opened these files in a Windows machine? Could they have Windows line endings? Try running sed 's/\r//' patterns.txt > newFile and then tryingt he same approach with newFile. Also make sure that there are actually matches: remove the redirection, just run for chunk in pattern_chunk*; do echo $chunk; grep -Ef "$chunk" target.fastq; done and see if there is any output. $\endgroup$
    – terdon
    Jun 4, 2021 at 17:57
  • $\begingroup$ Hello, I just did the sed newFile and re-run the while commands and it didn't work. Gave me the same result: echo is running but empty output. I know for a fact there are matches, as by testing the first 10 lines of my patterns on this one specific target.fastq I get 500 hits. When I run grep normally like this, I get printed on screen all the hits highlighted etc (and I can save them to output easily, no issues there). But when I use the same patterns in a loop it gives me a blank output :/ grep just seems to not be doing it's job within the loops like when it's out of them $\endgroup$
    – Jean Luc
    Jun 4, 2021 at 18:10
0
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ Hello! Thank you for your advice. I will try and see if the computer cluster possesses Cookiecutter already. Unfortunately, I don't necessarily have the required credentials to make any modifications (i.e. install and uninstall software) since it's an institution's cluster. But I will try and see nonetheless! $\endgroup$
    – Jean Luc
    Jun 5, 2021 at 11:44
  • $\begingroup$ @JeanLuc you don't need administrative rights to install software, if you have a compiler (e.g. gcc) you can compile the source code and use it without moving the binary to shared directories (which is the part that requires administrative rights) or you can try to download the binaries and use them directly (github.com/ad3002/Cookiecutter/releases/tag/v1.0.0) $\endgroup$ Jun 6, 2021 at 15:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.