1

I have a file of lines of headers (file 1) and another file of sequences in fasta format (file 2). I want to grep fasta sequences if a header line from file 1 matches file 2.

Example:

File 1:

>sp|B7UM99|TIR_ECO27
>sp|P06616|ERA_ECOLI

File 2:

>sp|B7UM99|TIR_ECO27
MPIGNLGNNVNGNHLIPPAPPLPSQTDGAA
RGGTGHLISSTGALGSRSLFSPLRNSMADS
VDSRDIPGLPTNPSRLAAATSETCLLGGFE
VLHDKGPLDILNTQIGPSAFRVEVQADGTH
......
>sp|P06616|ERA_ECOLI
MSIDKSYCGFIAIVGRPNVGKSTLLNKLL
GQKISITSRKAQTTRHRIVGIHTEGAYQAIY
VDTPGLHMEEKRAINRLMNKAASSSIGDVE
LVIFVVEGTRWTPDDEMVLNKLREGKAPVI
............
>sp|P0AD68|HUMAN
MKAAAKTQKPKRQEEHANFISWRFALLCGC
ILLALAFLLGRVAWLQVISPDMLVKEGDMR
SLRVQQVSTSRGMITDRSGRPLAVSVPVKA
IWADPKEVHDAGGISVGDRWKALANALNIP
.............

DESIRED OUTPUT

>sp|B7UM99|TIR_ECO27
MPIGNLGNNVNGNHLIPPAPPLPSQTDGAA
RGGTGHLISSTGALGSRSLFSPLRNSMADS
VDSRDIPGLPTNPSRLAAATSETCLLGGFE
VLHDKGPLDILNTQIGPSAFRVEVQADGTH
......
>sp|P06616|ERA_ECOLI
MSIDKSYCGFIAIVGRPNVGKSTLLNKLL
GQKISITSRKAQTTRHRIVGIHTEGAYQAIY
VDTPGLHMEEKRAINRLMNKAASSSIGDVE
LVIFVVEGTRWTPDDEMVLNKLREGKAPVI
............
  • This appears to be a duplicate of your previous question How to grep sequence of fasta using list of IDs in another file? (although it apparently didn't yet receive a satisfactory answer) – steeldriver May 9 at 22:22
  • @steeldriver this is a duplicate of the other question you mentioned. However, this is the much more usefully worded question. This demonstrates a couple of features that aren't clear in the other question, e.g. that it needs to match multiple and variable number of lines after the original match. – RaidPinata May 15 at 15:57
  • @Manoj Kumar I was happy to answer this question here. However, the real wizards at this kind of thing are going to be over at stackoverflow and other portions of stack exchange. AskUbuntu doesn't have nearly the amount of expertise in text-processing and shell scripting as other sites. – RaidPinata May 15 at 15:58
0

Okay, this is a slightly complicated grep you want to do. In the past, I've written programs in Perl to do this but there is actually a fairly straightforward shell usage to do it. First, grep can't do multi-line matches so we'll need a grep-like tool that can. pcregrep is like grep but has multi-line functionality as suggested on stackoverflow.

sudo apt-get install pcregrep

To do this we have to accomplish the following steps.

1) Add escape characters to all of the special characters in your id file.

sed 's/[^a-zA-Z 0-9]/\\&/g' id_file.txt > search_file.txt

I got this rather amazing bit of sed wizardry from unix.com - shell scripting.

2) Next we add the rest of our search pattern to the end of the line.

sed -i 's/\>$/\$\[^\>\]\*/' search_file.txt

Modified from stack exchange and shellhacks.

3) Finally, we can use this search file as the input to pcregrep with the fasta file to select our desired sequences.

pcregrep -M -f search_file.txt sequences.fasta

As stated by @DavidForester in you previous question but now using pcregrep for multi-line matching.

Bringing it all together; we can do everything as a (slightly complicated) one liner like this.

sed 's/[^a-zA-Z 0-9]/\\&/g' id_file.txt | sed 's/\>$/\$\[^\>\]\*/'| pcregrep -M -f /dev/stdin sequences.fasta

This has the advantage of not requiring an additional file.

Briefly, we pipe the output between the commands to get the proper search strings and then pipe this to pcregrep using the -f /dev/stdin flag. My testing shows this to produce exactly the output you are requesting. I wish I'd been able to come up with this when I needed it two years ago!

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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