I am looking for a command to count the total number of unique proteins in a file. For instance I wrote a command to BLAST XYZ proteins as query against the DLY proteins as a database to determine the number of genes or proteins that are shared between XYZ and DLY:

blastp -query ~/user/XYZ.faa -db DLY -evalue 1e-10 -outfmt '6 qseqid' | sort | uniq | wc -l


I used the -evalue option to provide me with alignments with an e-value that is 10^-10 or lower.

How can I represent my command so that that it can provide me with a proportion of the total proteins in each species (i.e. XYZ and DLY).

Thank you!

• I notice you changed your title from "how to X" to "please help with X!" actually the former is better, however, your title is not informative at all hence the very low view counts (it is the weekend though). "Tallying the proportion of the total proteins in each species from a Blast search" would be a more informative title. I am human so I do things here because they are interesting —so you need to convince me an other volunteer to look at your post. Same applies with Twitter posts, papers etc. informativeness/catchiness gets you engagement. Oct 18 '20 at 11:04
• Anyway, about your actual question. It is best done programmatically. Species are in NCBI in square brackets in the header. So grep with -o flag and a regex pattern say '$.*?$ would get those brakets, as you don't want the square brackets, a lookbehind will probably do, e.g. grep -o "(?<=$)[^$]+" -. However, there are several problems with a species tally of a blast result... Oct 18 '20 at 11:20
• If one wanted the count of paralogues, the sequenced species are not a reflection of diversity —if you search anything E. coli you get thousands of strains. In your case, orthologue count per species, has the problem that complete genomes only would work and even then you'd have to make sure of this... So sed or grep are a poor choice. Doing it manually for representative species is actually the easiest solution... Oct 18 '20 at 11:21