How to extract headers from protein fasta file for seqkit analysis?

I have run an analysis with SignalP to identify surface proteins (SP). This program requires a protein fasta file (faa) and provides a list of possible hits with a header slightly different from the original:

Original fasta protein
>lcl|NZ_CP184833.1_prot_WP_005374030.1_1113 [locus_tag=ACOB85_RS05640] [protein=VPA0252 family protein] [protein_id=WP_005374030.1] [location=1254001..1254237] [gbkey=CDS]
MKKLIALVLMVCSTNVAFAGDGSYTGTIDMTMATKNQVDVERADSQLTLKGRNEMAYGGVTIHSDAATQK
IELRGRRD

SignalP
lcl_NZ_CP184833.1_prot_WP_005374030.1_1113 _locus_tag_ACOB85_RS05640_ _protein_VPA0252 family protein_ _protein_id_WP_005374030.1_ _location_1254001..1254237_ _gbkey_CDS_

I now need to extract the amino acid sequence of the selected proteins found by SignalP; I am using seqkit to do the job.

If I use the whole SignalP header, I get the error (not just a warning since the output file is empty):

$ seqkit grep -n -p -f ids_SP.txt prot.faa -o recept_SP.faa
[WARN] space found in pattern, you may need use -n/--by-name: lcl_NZ_CP184833.1_prot_WP_017633643.1_233 _locus_tag_ACOB85_RS01190_ _protein_c-type cytochrome_ _protein_id_WP_017633643.1_ _location_259352..259672_ _gbkey_CDS_

I used bash to select the first part of the header, splitting at the first space, and changing the first underscore back into a pipe, but even with the partial hit flag (-p), I get an error:

$awk '{print $1}' ids_SP.txt | sed 's/lcl_/lcl|/' | tr -d '\r' | sed '/^$/d' > head_SP.txt
$ seqkit grep -n -p -f head_SP.txt prot.faa -o recept_SP.faa
[ERRO] fastx: invalid FASTA/Q format

How can I subset amino acid sequences from a set of proteins based on their names?

Thank you

signalp subsetting seqkit headers faa

• 785 views • link

updated 4 months ago by fernando •

• written 4 months ago by marongiu.luigi

1 answer

It appears that SignalP converts special characters (not numbers, letters, dots, underscores, or blanks) in the sequence header line to underscores.

Here's one way to do the same conversion with seqkit before extraction.

$ seqkit replace -p '[^\w\.\_\s]' -r '_' pro.faa \
    | seqkit grep -n -f ids_SP.txt 
[INFO] 1 patterns loaded from file
>lcl_NZ_CP184833.1_prot_WP_005374030.1_1113 _locus_tag_ACOB85_RS05640_ _protein_VPA0252 family protein_ _protein_id_WP_005374030.1_ _location_1254001..1254237_ _gbkey_CDS_
MKKLIALVLMVCSTNVAFAGDGSYTGTIDMTMATKNQVDVERADSQLTLKGRNEMAYGGV
TIHSDAATQKIELRGRRD
`

Don't forget to check if the number of extracted sequences matches the line number of ids_SP.txt. If not, you might have to adjust the regular expression.

Update: if you want to keep the header untouched, then match by the protein accession (WP_005374030.1_1113).

$ perl -pne 's/^.+?(WP_.+?)\s+.+$/$1/' ids_SP.txt \
    | seqkit grep -f - --id-regexp '^.+?(WP_.+?)\s+' pro.faa 
[INFO] 1 patterns loaded from file
>lcl|NZ_CP184833.1_prot_WP_005374030.1_1113 [locus_tag=ACOB85_RS05640] [protein=VPA0252 family protein] [protein_id=WP_005374030.1] [location=1254001..1254237] [gbkey=CDS]
MKKLIALVLMVCSTNVAFAGDGSYTGTIDMTMATKNQVDVERADSQLTLKGRNEMAYGGV
TIHSDAATQKIELRGRRD

• 0 views • link

written 4 months ago by shenwei356

thank you, but still doesn't work. seqkit replace -p '[^\w\.\_\s]' -r '_' pro.faa | seqkit grep -n -f ids_SP.txt gives [INFO] 144 patterns loaded from file but the file is empty; the perl command works but then sqekit still gives the error. Is there a more straightforward way to get the amino acid sequence of selected proteins from a file? It should be a common bioinformatic task...

• 0 views • link

written 4 months ago by marongiu.luigi

I don't know if it an overkill, but, with the help of AI to go fast I wrote this bash code:

#!/bin/bash

#@ extract_sequences.sh
# Define filenames
IDS_FILE=$1
SOURCE_FILE=$2
OUTPUT_FILE=$3

# Clear the output file if it already exists
> "$OUTPUT_FILE"

# Loop through each ID in the list
while IFS= read -r id || [[ -n "$id" ]]; do
    # Remove carriage returns if the file was created on Windows
    clean_id=$(echo "$id" | tr -d '\r')

    echo "Processing: $clean_id"

    # Use awk to find the ID and print until the next ">"
    # 1. We match the line starting with ">" followed by our ID
    # 2. We set a flag 'p' to 1 (true)
    # 3. As long as 'p' is true, we print lines
    # 4. If we hit a ">" on a line that DOES NOT contain our ID, we stop (p=0)
    awk -v id="$clean_id" '
        $0 ~ "^>" id {p=1; print; next} 
        /^>/ {p=0} 
        p {print}
    ' "$SOURCE_FILE" >> "$OUTPUT_FILE"

done < "$IDS_FILE"

echo "Extraction complete. Results saved in $OUTPUT_FILE."

Then launch it with extract_sequences.sh ids_SP.txt pro.faa output.faa and got the needed subsetting.

• 0 views • link

written 4 months ago by marongiu.luigi

I can't see the bash script do any id edit, but how did it work? it's weird.

• 0 views • link

written 4 months ago by shenwei356

I bypassed the starting of the header (which is different in the faa file and in the SignalP output) and focused on the protein ID, the only modification is the removal of possible carriage returns. I then select everything until the next '>' sign. It worked and there were no errors afterwards. Thanks. PS: pro.faa is a typo; it should be prot.faa, but the error was still there...

• 0 views • link

written 4 months ago by marongiu.luigi

what are the errors? note that I used pro.faa instead of prot.faa.

it is a commond task, but the problem here is the sequence headers are different, changed a lot by SignalP.

• 0 views • link

written 4 months ago by shenwei356