Filter multi fasta file with matching ID string between specific symbols

user_12746868 0
Viewed 695

I got stuck over this problem for a while and although I found a diverse set of similar questions none exactly fitted my problem or solved the issue. So here is the deal:

I have an input.fasta , formatted like this:

>sp|O42363|APOA1_DANRE Apolipoprotein A-I OS=Danio rerio OX=7955 GN=apoa1 PE=2 SV=1
MKFVALALTLLLALGSQANLFQADAPTQLEHYKAAALVYLNQVKDQAEKALDNLDGTDYEQYKLQLSESLTKLQEYAQTTSQALTPYAETISTQLMENTKQLRERVMTDVEDLRSKLEPHRAELYTALQKHIDEYREKLEPVFQEYSALNRQNAEQLRAKLEPLMDDIRKAFESNIEETKSKVVPMVEAVRTKLTERLEDLRTMAAPYAEEYKEQLVKAVEEAREKIAPHTQDLQTRMEPYMENVRTTFAQMYETIAKAIQA
>sp|Q90260|ASL1B_DANRE Achaete-scute homolog 1b OS=Danio rerio OX=7955 GN=ascl1b PE=2 SV=1
MEATVVATTQLTQDSFYQPFSESLEKQDRECKVLKRQRSSSPELLRCKRRLTFNGLGYTIPQQQPMAVARRNERERNRVKQVNMGFQTLRQHVPNGAANKKMSKVETLRSAVEYIRALQQLLDEHDAVSAVLQCGVPSPSVSNAYSAGPESPHSAYSSDEGSYEHLSSEEQELLDFTTWFDRYESGASMATKDWC
>sp|Q6TH01|C10_DANRE Protein C10 OS=Danio rerio OX=7955 GN=si:dkey-29f10.1 PE=2 SV=1
MASAPAQQPTLTVEQARVVLSEVIQAFSVPENAARMEEARESACNDMGKMLQLVLPVATQIQQEVIKAYGFNNEGEGVLKFARLVKMYETQDPEIAAMSVKLKSLLLPPLSTPPIGSGIPTS
>sp|Q6PFL6|CCD43_DANRE Coiled-coil domain-containing protein 43 OS=Danio rerio OX=7955 GN=ccdc43 PE=2 SV=1
MAAPEQIAGEFENWLNERLDSLEVDREVYGAYILGVLQEEESDEEQKDALQGILSAFLEEETLEEVCQEILKQWTECCSRSGAKSNQADAEVQAIASLIEKQAQIVVKQKEVSEDAKKRKEAVLAQYANVTDDEDEAEEEEQVPVGIPSDKSLFKNTNVEDVLNRRKLQRDQAKEDAQKKKEQDKMQREKDKLSKQERKDKEKKRTQKGERKR
>sp|P0C7U5|C5AR1_DANRE C5a anaphylatoxin chemotactic receptor 1 OS=Danio rerio OX=7955 GN=c5ar1 PE=3 SV=1
MDDNNSDWTSYDFGNDTIPSPNEISLSHIGTRHWITLVCYGIVFLLGVPGNALVVWVTGFRMPNSVNAQWFLNLAIADLLCCLSLPILMVPLAQDQHWPFGALACKLFSGIFYMMMYCSVLLLVVISLDRFLLVTKPVWCQNNRQPRQARILCFIIWILGLLGSSPYFAHMEIQHHSETKTVCTGSYSSLGHAWAITIIRSFLFFLLPFLIICISHWKVYHMTSSGRRQRDKSSRTLRVILALVLGFFLCWTPLH

and an ids.txt list, formatted as this:

Q90260
Q6PFL6

I would like to extract all fasta sequences with their header for which the IDs of ids.txt are element of the header. I have tried grep -w -A 2 -Ff id_list.txt input.fasta --no-group-separator > out.fasta but that did not work.

Ideally, I would like to express via regex to check if the string between the two | of each line starting with >sp matches any ID in my idx.txt. And if so, to store that header and fasta in out.fasta.

So that out.fasta would look like that:

>sp|Q90260|ASL1B_DANRE Achaete-scute homolog 1b OS=Danio rerio OX=7955 GN=ascl1b PE=2 SV=1
MEATVVATTQLTQDSFYQPFSESLEKQDRECKVLKRQRSSSPELLRCKRRLTFNGLGYTIPQQQPMAVARRNERERNRVKQVNMGFQTLRQHVPNGAANKKMSKVETLRSAVEYIRALQQLLDEHDAVSAVLQCGVPSPSVSNAYSAGPESPHSAYSSDEGSYEHLSSEEQELLDFTTWFDRYESGASMATKDWC
>sp|Q6PFL6|CCD43_DANRE Coiled-coil domain-containing protein 43 OS=Danio rerio OX=7955 GN=ccdc43 PE=2 SV=1
MAAPEQIAGEFENWLNERLDSLEVDREVYGAYILGVLQEEESDEEQKDALQGILSAFLEEETLEEVCQEILKQWTECCSRSGAKSNQADAEVQAIASLIEKQAQIVVKQKEVSEDAKKRKEAVLAQYANVTDDEDEAEEEEQVPVGIPSDKSLFKNTNVEDVLNRRKLQRDQAKEDAQKKKEQDKMQREKDKLSKQERKDKEKKRTQKGERKR

I am pretty sure this can be expressed via awk or grep but I am new to bash so I am having a hard time right now.

Thanks a lot in advance! :)

biopython fasta bioinformatics bash
2 Answers

Using join and sort:

join -t \| -2 2 -o 2.1,2.2,2.3 <(sort ids.txt) <(sort -t \| -k2 input.fasta)

Assuming there is no extra | character in the input.fasta and the order of output lines isn't significant.

user_138090010

With awk:

awk -F'[|]' 'NR==FNR{ids[$0];next}$2 in ids' ids.txt input.fasta

Explanation:

I'm passing both files ids.txt and input.fasta as input files to awk. The order is important. -F'[|]' sets the input field delimiter to a |.

The script:

# NR  is the overall record (line) number
# FNR is the record (line) number in the current input file

NR==FNR {   # True as long as we are reading the first input file
    ids[$0]   # Create a key in ids for every id from ids.txt
    next      # Don't process further actions
}

# Because of the 'next' statement above, we'll reach this point only
# when reading the second input file (input.fasta)
$2 in ids   # Print the current line if the second field
            # was found in the ids lookup

Output:

>sp|Q90260|ASL1B_DANRE Achaete-scute homolog 1b OS=Danio rerio OX=7955 GN=ascl1b PE=2 SV=1
>sp|Q6PFL6|CCD43_DANRE Coiled-coil domain-containing protein 43 OS=Danio rerio OX=7955 GN=ccdc43 PE=2 SV=1

Update: It turned out that you also want to print the below the match. This can be achieved like this:

BEGIN {
    FS="[|]"
}

# NR  is the overall record (line) number
# FNR is the record (line) number in the current input file

NR==FNR {   # True as long as we are reading the first input file
    ids[$0]   # Create a key in ids for every id from ids.txt
    next      # Don't process further actions
}

# Because of the 'next' statement above, we'll reach this point only
# when reading the second input file (input.fasta)
$2 in ids  {
    # set or reset a variable p to 2 if the second field
    # was found in the ids lookup
    p = 2
}

# Decrement the variable p on every iteration and check if it
# is greater than 0 after that. If that's true, awk will print 
# the current line
p--> 0
user_1713180
Related