How to get strain names/ids contained in a multi FASTA file using seqkit?

Question

FASTA files can be very big and unwieldy, especially if lines are at most 80 characters, one can't speed up browsing them by using less with -S to have one sequence every two lines.

How can I extract just the strain names (or sequence names, i.e. the string on the line starting with >) into a list? I'd like to use seqkit since I'm already familiar with it and know some of the flags and idioms it uses.

I googled the question but couldn't find an answer using seqkit - the examples all use awk, grep, sed or similar, which are alright for copy/pasting but not very memorable.

Answer 1

The nice thing about seqkit is that the answer is really simple and memorable:

seqkit seq -n sequences.fasta

This will output a text file containing one sequence name per line:

MPX/PT0001/2022
MPX/PT0002/2022
MPX/PT0003/2022

Answer 2

You can trivially do this with standard *nix tools. Here are a few options:

sed -n 's/^>//p' file.fa 
awk 'sub(/^>/,"")' file.fa
grep '^>' file.fa | tr -d '>'
perl -ne 'print if s/^>//' file.fa 

The advantage of these is that they will work on any POSIX-compliant system and most other *nix flavors, even if they're not POSIX such as Linux, out of the box without needing to install any external tools.

Another thing to take into account is how fast a command is, especially when dealing with large files. I ran a few tests using the solutions presented here on the human hg38 genome assembly which can be downloaded here and seqkit version 2.2.0 which was downloaded from here. The file was stored locally on my laptop's SSD drive, and the laptop itself is:

$ inxi 
CPU: quad core Intel Core i7-6820HQ (-MT MCP-)
speed/min/max: 2896/800/3600 MHz Kernel: 5.17.5-arch1-1 x86_64
Up: 25d 10h 27m Mem: 15035.9/31958.1 MiB (47.0%)
Storage: 476.94 GiB (103.1% used) Procs: 452 Shell: Bash inxi: 3.3.15

To get more accurate results, I used a little function I have which will run the command you give it a specified number of times and then report the average real time across all runs. Here's the code if you want to try it:

timethis(){
    max=$1;
    shift;

    for com
    do
        printf '\nCOMMAND: %s\n' "$com"
    c=0;
    while [[ $c -lt $max ]]
    do
      let c++;
      ( time -p eval "$com" ) 2>&1 | grep -oP 'real.*?\K[\d\.]+'
        done | 
        awk -vm=$max '{k+=$1}END{print (k/m)}';
    done
}

And the results. These are the average "real" (i.e. wall clock) times taken across 10 runs:

$ timethis 10 \
  "rg '^>' GCA_000001405.29_GRCh38.p14_genomic.fna | tr -d '>' >/dev/null" \
  "grep '^>' GCA_000001405.29_GRCh38.p14_genomic.fna | tr -d '>' >/dev/null" \
  "sed -n 's/^>//p' GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null" \
  "perl -ne 'print if s/^>//' GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null" \
  "./seqkit seq -n GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null" \
  "awk 'sub(/^>/,\"\")' GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null" \
  "perl -nle 'push @res, $_  while s/^>//; END {print  join \"\n\",@res}' GCA_000001405.29_GRCh38.p14_genomic.fna"

COMMAND: rg '^>' GCA_000001405.29_GRCh38.p14_genomic.fna | tr -d '>' >/dev/null
0.474

COMMAND: grep '^>' GCA_000001405.29_GRCh38.p14_genomic.fna | tr -d '>' >/dev/null 
0.583

COMMAND: sed -n 's/^>//p' GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null 
5.62

COMMAND: perl -ne 'print if s/^>//' GCA_000001405.29_GRCh38.p14_genomic.fna  >/dev/null 
6.408

COMMAND: ./seqkit seq -n GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null 
4.326

COMMAND: awk 'sub(/^>/,"")' GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null 
13.518

COMMAND: perl -nle 'push @res, 3  while s/^>//; END {print  join "\n",@res}' GCA_000001405.29_GRCh38.p14_genomic.fna 
8.455

As you can see above, the grep or rg and tr approaches only took around half a second on average, so that is clearly much faster than any of the other choices. However, seqkit is pretty good since it only took around 4 seconds, faster than any of the other ad hoc solutions.

Answer 3

Or a bit more verbose ...

perl -nle 'push @res, $_  while s/^>//; END {print  join "\n",@res}' myfile.fa