Tutorial

Extract all protein sequences of specific taxons from the NCBI nr database

Dataset

Steps

Taking virus for example.

  1. Getting all taxids of virus (taxid 10239):

    $ taxonkit list --ids 10239 --indent "" > virus.taxid.txt
    

    It takes only 2.5 seconds! Number of taxids:

    $ wc -l virus.taxid.txt
    163104 virus.taxid.txt
    
  2. Retrieving accessions or GIs(NCBI stopped using gi) with csvtk:

    • accession

      $ zcat prot.accession2taxid.gz | \
          csvtk -t grep -f taxid -P virus.taxid.txt | \
          csvtk -t cut -f accession.version > virus.taxid.acc.txt
      
    • gi

      $ zcat prot.accession2taxid.gz | \
          csvtk -t grep -f taxid -P virus.taxid.txt | \
          csvtk -t cut -f gi > virus.taxid.gi.txt
      

      It costs ~ 8 minutes.

  3. Retrieving nr sequences from BLAST database:

    • accesion

      $ blastdbcmd -db nr -entry all -outfmt "%a\t%T" | \
          csvtk -t grep -f 2 -P virus.taxid.acc.txt | \
          csvtk -t cut -f 1 | \
          blastdbcmd -db nr -entry_batch - -out nr.virus.fa
      
    • gi

      $ blastdbcmd -db nr -entry all -outfmt "%g\t%T" | \
          csvtk -t grep -f 2 -P virus.taxid.gi.txt | \
          csvtk -t cut -f 1 | \
          blastdbcmd -db nr -entry_batch - -out nr.virus.fa
      

    Another way is directly retrieving from nr FASTA sequences using SeqKit:

    gzip -c -d nr.gz | seqkit grep -f virus.taxid.acc.txt > nr.virus.fa