Practical lesson 2

COMMAND LINE BLAST AND SEQUENCE CLUSTERING WITH PERL

Applications that should be installed locally:

• BLAST
• ClustalW (optional)
• Belvu (only available for UNIX-Linux) or Pfaat (optional)
• BioLayout (optional)
• In addition, you should be able to run Perl scripts.

Aim of the exercise

• Getting used to operate command line BLAST.
• Practicing with other applications included in the BLAST package, such as blastclust, formatdb and fastacmd.
• Using BLAST and the other mentioned applications to compare a whole proteome against itself, with the idea of identifying clusters of related proteins.

---

Part 1:

Retrieving protein sequences and creating a database formatted for its use with BLAST.

• Create, in your home directory, a folder called BLASTDATA.

 
• RETRIEVING A FILE WITH A COMPLETE PROTEOME

• Save, in the BLASTDATA directory, your favorite proteome. For example, save the file with ALL predicted proteins for Mycoplasma pneumoniae in FASTA format, which is available from:

• NCBI Microbial genomes
• NC_000912.faa (local copy, predicted proteins from NC_000912).
 
• Then, make a copy of that file in your home directory. This is quite convenient, to avoid complications later.

 
• GENERATING A DATABASE

• In order to use the sequences in the fasta file as a database that BLAST will recognize, it is necessary to use the application formatdb, which is included in the BLAST package.

• Run, from your home directory, the following command.

formatdb -o -i BLASTDATA/NC_000912.faa

Seven files with indexes should be generated, together with the log file formatdb.log.

It may be necessary to specify the complete path to the application, depending on how was BLAST installed (for example: /blast/formatdb).
 
• Try to find out what is the the tag  -o  doing (checking the README file of BLAST, or using formatdb ?).

• RETRIEVING SEQUENCES FROM YOUR DATABASE

• The application fastacmd, which is also part of the BLAST package, allows the retrieval of sequences from a database generated with formatdb.

• Open the file NC_000912.faa, select any protein and take note of its identifier (for example: 13507741).

• Run, from your home directory, the following command:

fastacmd -d BLASTDATA/NC_000912.faa -s 13507741

You will retrieve the peptide sequence identified as 13507741 in the database NC_000912.faa.

As before, it may be necessary to specify the path to fastacmd.

 
• Important: it is possible to retrieve a set of sequences, by listing their identifiers in a file, and using the command:

fastacmd -d BLASTDATA/NC_000912.faa -i list_file

Part 2:

Identification of clusters with blastclust

• IDENTIFICATION OF CLUSTERS WITH BLASTCLUST

• The application  blastclust, which is also part of the BLAST package, performs an all against all BLAST comparison with the members of a set of sequences.

• The sequence pairs which are similar with an associated e-value lower than 1e-6 are retained, and considered "connected", and the rest of sequence pairs are ignored (1e-6 is the default threshold, which can be changed by editing a configuration file).

• When blastclust is called it is also possible to use some parameters that will modify the behavior of additional filters (the parameters S, T and L, see below).

• From the list of retained pairs of connected proteins, clusters are constructed following the principle of single linkage.

• blastclust uses as input a file with FASTA formatted sequences.

• To run blastclust, execute the following command:

blastclust -i NC_000912.faa -S 30 -L 0.0 -b T > Myco_blastclust_30.txt

This operation should generate a list of peptide clusters for the Mycoplasma pneumoniae proteome. The parameters -S 30 -L 0.0 -b T specify that for a pair of sequences to be considered as neighbors, at least 30% of the aligned residues should be identical, and the overlap between the aligned sequences should be at least 0% (the obvious minimum). The parameter -b T (both true) specifies that the condition of minimum overlap must be fulfiled when calculated as relative to both sequences involved in the alignment.

The file that SHOULD BE used as input for blastclust is the file NC_000912.faa that was located in your home directory, and not the one that exists in the directory BLASTDATA. This is because blastclust generates a temporary database that is removed automatically later. If you used the NC_000912.faa  file that is contained in BLASTDATA, blastclust would destroy the indexes created by formatdb.
 

• Open the file Myco_blastclust_30.txt to get an idea of how many clusters have been found, and how many sequences are contained in them. Each line in the file corresponds to a different cluster, and peptide identifiers are separated by blank spaces.

 
• It is important to realize that what is obtained with this procedure is a list of clusters, and NOT what in graph theory is called an adjacency list. In fact, the information about how proteins are connected WITHIN each cluster is not provided by blastclust.

 
• ALIGNING CLUSTER SEQUENCES WITH ClustalW, TO TEST THEIR HOMOGENEITY

• The sequence clusters generated with blastclust should contain similar proteins. You could now generate a multiple sequence alignment, to visually check whether the proteins within a given cluster are similar in their whole length or, on the contrary, similarity is restricted to some regions.

• Copy  the identifiers of all the proteins in a cluster of your choice (try, for example, with the biggest cluster), and paste them in a text file, as a column.

 
• Use, then, that file to retrieve all the sequences from that cluster, with the application fastacmd.

 
• Align the retrieved sequences with ClustalW, either locally (if you have it installed) or at the EBI server (ClustalW at the EBI). Once you have got a file with the multiple sequence alignment, you can visualize the results with Belvu or Pfaat, if you have any of them installed locally; if you run ClustalW at the EBI server, the output will include a tree representation and the option of colouring the alignment.

• If you observe that the sequences do not align satisfactorily enough, or the associated tree indicates that there is more than one main group of sequences, this could be indicating that you have used blastclust with very relaxed parameters. If that is the case, run again blastclust with more restrictive parameters. For example, increase the minimal allowed identity  from 30% to 60%.

blastclust -i NC_000912.faa -S 60 -L 0.0 -b T > Myco_blastclust_60.txt

• In this case, a HIGHER number of clusters, of SMALLER size, should be obtained. After aligning their sequences with ClustalW, it should become evident that now the clusters are more homogeneous.

• If you have saved in a file all the FASTA formatted sequences from a given cluster, you can also extract all the sequence headings in that file using the command:

to check whether the annotations point to some common function.

Part 3:

Command line BLAST.

• blastclust identifies automatically clusters of similar sequences. However, blastclust does not return information about how the sequences within clusters are interconnected. To obtain that kind of information it is necessary to obviate the automatic methods: we will have to compare first the sequences with BLAST and, then, use that information, and ad hoc Perl scripts, to identify sequence clusters.
 
• LETS TRY NOW TO USE BLAST BY COMMAND LINE

• First, retrieve any sequence from your database, for example 13507741 (the same as before), and save it in a file called, for example, Seq2.faa. The command would be as follows:

 fastacmd -d BLASTDATA/NC_000912.faa -s 13507741 > BLASTDATA/Seq2.faa

• Run the following commands, to compare the previously saved sequence with those in the database:

blastall -p blastp -d BLASTDATA/NC_000912.faa -i BLASTDATA/Seq2.faa -m 0
blastall -p blastp -d BLASTDATA/NC_000912.faa-i BLASTDATA/Seq2.faa -m 8

• In both cases:

• the option -p specifies the program that will be used.
• the option -d specifies the database in which the search will be performed.
• the option-i specifies the query sequence.
• the option -m specifies the type of output (0 for standard format with alignments; 8 for table format).

• Run the command

blastall

to get information about other options.

• ALL AGAINST ALL PROTEOME BLAST COMPARISON

• Run the following command:

blastall -p blastp -d BLASTDATA/NC_000912.faa -i BLASTDATA/NC_000912.faa -m 8 > Myco_vs_Myco.blast

• The file  Myco_vs_Myco.blast  CONTAINS the results of the all against all comparison, in table format (example).

Part 4:

Identification of clusters from BLAST results.

• You should now write a Perl script to that will parse BLAST results in table format, to extract the values of the two first columns (protein identifiers) and those of the 11th column (e-values). As a first step your script could simply print out those three columns, as in this example. In a second phase you could introduce a filter in such a way that only those sequence pairs with an e-value below a certain threshold are retained.

• However, what your script should really do is to create what is known as an ADJACENCY LIST, in which for each node considered, a list of its neighbors is described (as in this example).

• The general idea of the script could be:

• The script will need two arguments: the name of a file that contains BLAST results in table format, and a threshold for the e-value.
• The script will read the BLAST table, line by line.
• The first two columns and the 11th column will be extracted for each line.
• If the two identifiers are equal, the script will ignore that line.
• If the two identifiers are different, the magnitude of the e-value will be evaluated: only if it is lower than the threshold,  the two identifiers and the e-value will be stored in a hash of hashes:

$hash{$prot1}{$prot2} = e-value

If you cross the identifiers, you will ensure the symmetry of the table.

$hash{$prot1}{$prot2} = e-value
$hash{$prot2}{$prot1} = e-value

• Before you store the information about identifiers and e-values in the hash of hashes, it is also convenient to include an "if" condition that will check whether that information has been already stored, in such a way that the hash of hashes will be updated only if a LOWER e-value is found.

• By constructing a hash of hashes you will be, in fact, generating what is known as an ADJACENCY MATRIX. Such a matrix could now be printed in a variety of formats, for example, as an ADJACENCY LIST (example), or as the ADJACENCY MATRIX  itself (example).

• To generate the ADJACENCY LIST you will have to set up two nested foreach loops, that will extract the keys of each of the levels in the hash of hashes.

• If you need some kind of inspiration, you can take a look at the following script:

blast2adjLists_ed1.pl

blast2adjLists-ed1.pl is run as follows:

blast2adjLists-ed1.pl blast_table e-value

Where:

• blast_table is the file that contains BLAST results in table format;
• e-value: is the highest e-value allowed, for example, 1e-100.

If you want to execute the script, check that the path to Perl is correct, and that the script file is executable.

The output of blast2adjLists-ed1.pl  is an adjacency list like the one in the example previously given: it contain a list of nodes, and for each of them, a list of the connected neighbor nodes. It is easy to modify the script to generate a different output, for example, an adjacency matrix.

• Try now to write another script that will use the adjacency list as input, and that will identify connected components or clusters. This may be a little bit more complicated because it will require (in general) the use of a recursive function.

• If you are stuck, the following script will give you the solution.

adja2cluster_v2.pl

• adja2cluster.pl identifies connected components from an adjacency list generated by blast2adjLists_ed1.pl. It is used as follows:

adja2cluster-v2.pl adjacency_list

• The output of adja2cluster_v2.pl  (example) is a file that looks like the one generated by blastclust and, in fact, both strategies should generate the same clusters by choosing the appropriated parameters. Again, each line defines a cluster of proteins, referred to by their identifiers.

Part 5:

Graph visualization. Checking the results in Pfam.

• BIOLAYOUT

• BioLayout is a graph visualization tool, written in Java.

• It requires Java 1.3, which can be downloaded from  http://java.sun.com/.

• BioLayout uses a very simple input: a text file with at least two columns, that specify the pairs of connected nodes. If there is a third column,  it specifies the weights of the connections.

• By now you should have a script that extracts the identifiers and e-values from a table of BLAST results, and generates a file with three columns (example). Using a file like that as input for BioLayout you should get a graphical representation of the graph, in which you can recognize different groups of proteins (clusters) that are more similar among them, that with other members of the set of sequences.

• BioLayout includes a number of utilities to play with the graph. For example, it is possible to define node classes, to label classes, nodes or archs with different colors, and to move nodes and save the coordinates in a file.

• SEARCH IN Pfam

• Some of the protein clusters identified with blastclust, or with the combination of BLAST and Perl scripts, could correspond to already characterized protein families.

• For example, the main cluster identified by BLAST, blast2adjLists-ed1.pl with a threshold of 1e-100, and adja2cluster-v2.pl, contains 13 proteins from Mycoplasma pneumoniae (the file with the list of clusters is here). If you retrieve any of the sequences from the main cluster, using fastacmd, you could then make a search against the Pfam database (choose the "Protein Search" tab and paste a protein sequence in the apropriated window). You will discover that all 13 proteins of the cluster are part of a family of proteins that is represented ONLY in genus Mycoplasma, and which has 20 members: 13 of them in M. pneumoniae and 7 in M. genitalium (PF03257). The family is characterized in Pfam by a domain that has a length of 100 amino acids, which is the conserved block that you have probably identified after making multiple sequence alignments with ClustalW.