An Exercise to Get a List of Homologs of Virulence Factors

An Exercise to Get a List of Homologs of Virulence Factors


Suppose that you have a file containing data relating to known virulence factors. For example, here are the first few lines of a file given to me by a participant in one of our tutorials (These columns are separated by tab characters):


VFG_id gi_id Gene_name Product Organism VF_id


VFG1293 21282773        hla     Alpha-Hemolysin precursor       Staphylococcus aureus MW2       VF0001

VFG1798 46588   hlb     beta-hemolysin  Staphylococcus aureus   VF0002


The second column contains GI numbers (without the usual "gi|" prefix).  The goal was to extract the set of FIGfams containing these genes, along with a list of the PEGs (and their functions).


We can do these easily in two steps. We'll use the "svr" command line tools in this example. If you are using Windows, be sure to remember to use the myRAST shell,  or, if you are using a Mac,  ensure that the svr commands are in your path with this command in your shell window:


export PATH=$PATH:/Applications/myRAST.app/bin


Our solution requires that the servers use the data from the PSEED, so be sure to run this command in your shell window:


Mac users:  export SAS_SERVER=PSEED

Windows Users: set SAS_SERVER=PSEED


 Now, to begin,  let's try


svr_ids_to_figfams -c=2 -source=NCBI < known.virulence.factors 2> no.matches | svr_figfams_to_ids > with.figfams


This produces two files.  The file "no.matches" contains a list of the lines from the input file ("known.virulence.factors") that have not yet been placed into FIGfams. The second file ("with.figfams") is a file in which two columns have been added: a FIGfam number and a PEG from that FIGfam.  Now, if we wished to add a function of the PEG, as well as a readable version of the genome it is from, we would use


svr_gene_data function genome-name < with.figfams > desired.table.txt


If the user needs to see the "aliases" for the PEGs that have been found, they might try


svr_aliases_of -c=9 < desired.table.txt > desired.table.aliases.txt


but be warned: we return all of the IDs of protein sequences that we know about that have the same protein sequence (from a potentially large collection of very similar genomes).


Finally, we should note that the initial virulence factors may have multiple entries leading to the same FIGfam (and, hence, to the same sets of PEGs).  This can lead to the situation  where multiple roles, while not identical, would have information on the same PEG.  To collapse the table to unique, you need to sort on the 9th field, assuming tab-delimited fields.

This can be done using


sort -k 9 -t $'\t' -u < desired.table.aliases.txt  > sorted.txt



Let us consider an somewhat more useful alternative.


We begin by just connecting the aliases in column 2 to PEGs in the SEED:


svr_aliases_to_pegs -c=2 -source=NCBI -protein=1 < known.virulence.factors 2> no.matches.1 > with.peg


This produces a file (no.matches.1) of lines that could not be connected to PEGs.


Now,


svr_ids_to_figfams < with.peg > with.ff 2> no.matches.to.ff


can be used to add a FIGfam function and FIGfam ID to the end of each line

for which the identified PEG is in a FIGfam.  Those lines that do not have a PEG

get written to a file (no.matches.to.ff).


Now, it is worth noting that we may have many lines pointing at the same FIGfam (i.e.,

many lines with the same FIGfam ID in the last column, which is column 9).

By using


   sort -u -k 9 -t$'\t' with.ff > with.ff.sorted


we sort the file on the FIGfam ID, deleting all but one of each set that share

a common FIGfam ID.


Now,


svr_figfams_to_ids < with.ff.sorted > with.ff.peg


can be used to expand each input line to a set of lines, each one containing a 

distinct PEG from the FIGfam.  Using 


svr_gene_data function genome-name < with.ff.peg > complete


we add the PEG function and the genome-name to the end of each row.


Finally, we run


svr_aliases_of -c 10 < complete > complete.with.aliases


to add known aliases for each PEG as the last column in the table.