0
$\begingroup$

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!

Improve this question
$\endgroup$
3
  • 3
    $\begingroup$ 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. $\endgroup$
    – Matteo Ferla
    Commented Oct 18, 2020 at 11:04
  • 2
    $\begingroup$ 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... $\endgroup$
    – Matteo Ferla
    Commented Oct 18, 2020 at 11:20
  • 2
    $\begingroup$ 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... $\endgroup$
    – Matteo Ferla
    Commented Oct 18, 2020 at 11:21

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.