id=6656

# Step 1) list all taxid below $id
taxonkit list --ids $id --indent "" > $id.taxid.txt



# Step 2) get all accessions
# 3 min for me
pigz -d -c prot.accession2taxid.gz \
    | csvtk -t grep -f taxid -P $id.taxid.txt \
    | csvtk -t cut -f accession.version \
    > $id.taxid.acc.txt



# Step 3) Option A (very slow)
blastdbcmd -db nr -entry all -outfmt "%a\t%T" | \
    csvtk -t grep -f 2 -P $id.taxid.acc.txt | \
    csvtk -t cut -f 1 | \
    blastdbcmd -db nr -entry_batch - -out nr.$id.fa


# Step 3) Option B (relative fast but buggy)
# reformat nr.gz and retrieve sequences (buggy)
# 16 min for me
cat <(echo) <(pigz -d -c nr.gz) \
    | perl -e 'BEGIN{ $/ = "\n>"; <>; } while(<>){s/>$//;  $i = index $_, "\n"; $h = substr $_, 0, $i; $s = substr $_, $i+1; if ($h !~ /\[.+\]/) { print ">$_"; next; } while($h =~ /(.+? \[.+?\])/g){ $h1 = $1; $h1 =~ s/^\W+//; print ">$h1\n$s";} } '\
    | seqkit grep --delete-matched -f $id.taxid.acc.txt -o nr.$id.fa.gz

# checking sequences amount (with Option B)
seqkit stats nr.$id.fa.gz
# file           format  type      num_seqs        sum_len  min_len  avg_len  max_len
# nr.6656.fa.gz  FASTA   Protein  6,913,200  2,860,805,704        6    413.8   28,516

wc -l $id.taxid.acc.txt
# 8174610 $id.taxid.acc.txt



# Step 4) makeblastdb
gzip -d -c nr.$id.fa.gz > nr.$id.fa    
makeblastdb -in nr.$id.fa -dbtype prot -out nr.$id

Well, as for 6,913,200 < 8,174,610:

Since many genes in different species have same sequences, they are recorded as one record in nr.gz.

The intuitional way is matching FASTA record heads with all target accesions using regular expression matching, however, it's very slow for such a big dataset:

$ pigz -d -c nr.gz | seqkit stats | tee stats.txt
file  format  type        num_seqs         sum_len  min_len  avg_len  max_len
-     FASTA   Protein  169,199,080  61,790,996,688        6    365.2   74,488

So a better solution is "unfolding" these FASTA records.
But the FASTA headers in nr.gz are not well formated, it's hard to rightly unfold all records shared by muliple species. Here are some cases:

1. Normal. accession product [species]accession2 product2 [species2] (easy to solve)

   >XP_642837.1 hypothetical protein DDB_G0276911 [Dictyostelium discoideum AX4]EAL68957.1 hypothetical protein DDB_G0276911 [Dictyosteli
   um discoideum AX4]
   Mxxxxxxxx
   

2. Some FASTA headers contain 0x01 (shown as ^ below) in front of accession (easy to solve)

   $ pigz -d -c nr.gz | grep '^>' | grep CAH60339.1 | cat -A
   >CAH60339.1 histone h3, partial [Copelatus sp. 028TA588]^AABG79397.1 histone H3, partial [Lumbrineris magnidentata]^AABG79414.1 histone H3, partial [Polygordius lacteus]^AABX89916.1 histone 3, partial [Akalyptoischion atrichos]
   xxxx
   >^accesion product
   xxxxxx
   

3. Mix 1. Some headers contains species information in [], but thelast one does not. (there are ways)

   >xxxxx [xxxx]aaaa bbbbbb
   

4. Mix 2. Where B3N7H1.1 does not end with ], so the last character 4 and the later accession A1Z7U3.1 lay together -_-! (HARD to solve)

   >NP_001027396.1 uncharacterized protein Dmel_CG33774 [Drosophila melanogaster]XP_001970317.1 uncharacterized protein Dere_GG10558 [Drosophila erecta]B3N7H1.1 RecName: Full=Dolichyl-diphosphooligosaccharide--protein glycosyltransferase subunit 4A1Z7U3.1 RecName: Full=Dolichyl-diphosphooligosaccharide--protein glycosyltransferase subunit 4AAZ52820.1 uncharacterized protein Dmel_CG33774 [Drosophila melanogaster]EDV59376.1 uncharacterized protein Dere_GG10558 [Drosophila erecta]