MCB 3421 : Database searches using Blast and Blastall

Questions you should answer are given in red.

Please send your answers per email to gogarten@uconn.edu with subject "MCB 3421, Assignment 6", the first line of the email should contain your name.

1) Search Entrez for a transposase sequence from Frankia (an Actinobacterium living in the soil and as a nitrogen fixing symbiont in some angiosperm trees and shrubs). Use the protein sequence and search one Frankia genome of your choice (at http://www.ncbi.nlm.nih.gov/sutils/genom_table.cgi , you can reach this page from the normal blast page). Use blastp and tblastn as search algorithm.
How many significant hits did you find?

gi number of your query:
Name of genome searched:
number of hits using blastp:
number of hits using tblastn:

2) Using Blast on a "local" machine:

Download two genomes of your choice, using one of these two methods (if in doubt of what to choose, go for one of the Thermotoga and one of the Pyrococcus genomes):

Method 1: Use the "F" links in the "Tools" column for your chosen genomes here. (ctrl click on the *.faa file, not the F, to save the file on your computer. Rename the file so that you recall what is in it.

Method 2: Go into the bacteria sub-folder using the ftp:// link above. Click on the appropriately named folder (e.g., Pyrococcus abyssi, Pyrococcus is a member of the archaeal domain).

Among the files in these folders we are are interested in the .ffn file which contains all nucleotide sequences for individual predicted ORF's, the .faa file which contains the predicted ORF's amino acids sequences, and the .fna file, which contains the whole genome. These files are in a fasta format. Another file of interest ar the ptt files that contain information on the individual annotated open reading frames, This file is useful, if you want to map ORFs by position in the genome. For this exercise we will only use the .faa files. If a folder contains more than one .faa file, take both of these, and copy them into a single file*. The same approach is recommended, if you want to use more than one genome as target. On to your computer rename the files into something you recognize. (e.g., p_abyssi.faa and t_maritima.faa).

Load the two .faa files into a folder in your home directory on bbcxsrv1.biotech.uconn.edu.
As usual there are many ways to do this.

Alternative #1: The easiest might be to use the finder. Select Go->connect to server in the finder-menu. In the address field enter afp://bbcxsrv1.biotech.uconn.edu . When prompted enter your username and password (if you still need to change your password, see below first). This opens a folder that contains the contents of your directory on the bioinformatics cluster. You can drop and drag files into that folder, and you can use the finder to create new folders. In your folder on the cluster create a sub-folder called btalsttest and move your two faa files into that folder. Continue after =======

Alternative #2: ssh "always" works, but it is very pedestrian: Save the two files onto your computer. Then open a terminal window on your mac (in the utilities folder inside the application folder) and type the commands in green

ssh your_account_name@bbcxsrv1.biotech.uconn.edu    (enter your password when prompted)
mkdir blasttest    (creates a directory named blasttest)

Open another terminal window on your mac and type
sftp your_account_name@bbcxsrv1.biotech.uconn.edu   (enter your password when prompted. This establishes a secure ftp connection to the XSERV cluster)
cd blasttest changes to the directory blasttest. lcd changes the directory on your local machine, lpwd reports your local directory name
mput *.faa transfers all files from your local directory that end with .faa from the local machine to the remote machine.

Alternative 3: use FUGU and Jellyfissh, or finder and jellyfissh. jellyfissh is a simple program that keeps track of the ssh connections. It saves you from typing ssh username ....

Whatever works for you - get the two files into your account on the cluster and establish a terminal with an ssh connection to the cluster.

Establish an ssh connection to bbcxsrv1.biotech.uconn.edu

Open a terminal window on your mac (in the utilities folder inside the application folder) and type the commands in green

ssh your_account_name@bbcxsrv1.biotech.uconn.edu    (enter your password when prompted)

=============== At this point you should have (A) and ssh connection to the cluster running in terminal, and two genomes in a directory on the cluster called blasttest ==============

On the cluster (remote machine), you should NOT run long processes on the master node. Either qrsh to a node that is not busy (type qrsh **), or submit your command to the queue using qsub.

On the XSERV cluster using the ssh connection established earlier type
qrsh -l hipri=TRUE <enter> enter your password when prompted,
cd blasttest
formatdb -i p_abyssi.faa -o T -p T
The -i is to indicate your infile.name, -p T is for a protein sequence, for a nucleotide sequence, you need to use -p F (for False); the -o T instructs the formatdb program to create indices. For more information on the formatdb program check formatdb man page on the web (gnusquad.org) or type  formatdb - <enter>
You need to replace p_abyssi.faa with one of the genomes you downloaded. This genome will be converted into a searchable database.

To do a blast search of every sequence in t_maritima.faa against the P.abyssi genome, we use blastall. To get information on the program type
blastall <enter>
For details check one of the many help pages on here and here.

We will use the following options (you need to change the command to reflect the names of the genomes you want to compare; -d gives the name of the databank that you created using formatdb):

blastall -i t_maritima.faa -d p_abyssi.faa -o blast.out -p blastp -e 10 -m 8 -a 2
-i the query sequence(s)
-d the database
-o the output file
-p the blast program to use
-e the e-value cut off
-m the format of the output. option 8 gives a tabular output, that can be easily read into excel spreadsheet or into a database program. option 9 is similar but includes comment lines.
the individual columns in the output denote the following
# Query id # Subject id # % identity # alignment length # number of mismatches # number of gap openings # position of query start # position of query end # position of subject start # position of subject end # e-value of a hit # bit score of a hit

Which option turns off the the low complexity filter?
Which option, and which setting, sets the wordsize to 2?
Which option allows to use two processors?

Inspect the output typing more blast.out on the commandline on the cluster

Often one is only interested in the top scoring hit for each query. A simple perlscript (here) goes through the blastall output in m8 or m9 format and retains only the top scoring hit, and it adds a header line that will be correctly interpreted by EXCEL. Save this script into your blasttest folder on the cluster (the easiest is to ctrl click on here and select save link as, then select the blasttemp folder in the finder). Glance through the script using a text editor (vi, or textwrangler, or just click on the here link in your browser).

To run the blast.out file through the perl script type
perl extract_lines.pl blast.out
This will create a file in your directory called blast.out.top that only contains the top scoring hits, and from which the comment lines (in case you used the m9 option in blastall) are removed.

A modified version of this program (here) adds another column to the listing that for each alignment gives the bitscore per residue in the alignment.

Import this file into an excel spreadsheet. If you have an afp connection to the cluster, you can

start EXCEL on you local machine,
open an empty spreadsheet,
select Data -> Get external data;
select your file (you'll need to tell it to show you all files, not only files with a recognized extension)

Sorting the data on the different columns you can figure out which gene was found to be most similar between the two genomes.

Which were the genomes you selected?
Which gene pair had the highest % identity?

(to retrieve the information on target sequence, you can use the command fastacmd. If the GI number of the target is 14521962, and the database is p_all.faa the command is as follows:
fastacmd -s 14521962 -d p_all.faa
By default the databank is nr as installed on the cluster. Thus fastacmd -s 14521962 also works (and you get more annotation), provided the sequences in your databank are also in nr.

Next, we will calculate a histogram for the % identity values for all the significant alignments. The 2008 version of Excel for Macs no longer allows to calculate histograms via the add-ons*** so we will use a script in 'R' (a package for statistical analysis and graph plotting****). You do not need to learn or understand R today (or even for this course), but you will have to open the R script in a text editor and read some minor editing.
Download the script from here (you may have to right click and 'save link as...') in the blasttest folder on the cluster via afp. At the command line, in the blasttest folder type

R CMD BATCH plot_from_blast.R

The script should have generated a file called blast_plots.pdf. Using the afp connection you made earlier (or one of the alternative methods described above) locate this PDF file in Finder and open it.

Do you observe a smooth distribution of % identiy values?
Do any genes constitute a second peak of more similar sequences, compared to the bulk of the comparisons?

Are there other interesting things you could plot from these data?

* You can use a text editor for this (textwrangler), but the files might be pretty large. An alternative is to use unix:

open a terminal window (the application is in the utility folder inside the application folder).
cd to the directory where you downloaded the two files (e.g., cd Desktop)
cat name1.faa name2.faa > new_name_for_both_faa_files.faa
where name1.faa name2.faa are the names of the files you want to copy. cat list the content of one ore more files and > directs the output from the default (the screen) to a file

** Screen is a unix program that can be very useful, especially, if you want to go home and keep the cluster running on a program that you started. The best is to start screen immediately after you connect to the server via ssh. To get going type screen -a <return>, then type qrsh <return> this sends you to a compute node, you might need to enter your password again. Once you are ready to leave, enter the following keys <ctrl> and a at the same time followed by a d. This detaches the screen. logout from the master node; go home; ....; when you want to check the program, login to the master node and type screen -r to reattach to the screen (if you have several screens running you get instructions how to connect to a particular screen). In class we will use the high priority queue that works for jobs that do not run for a long time.

*** The 2008 version of Excel for Macs no longer allows to calculate histograms via the add-ons, but there are spreadsheets available on the 'net into which you can copy and paste the data you want to plot - a histogram generating spreadsheet for Excel2008 is here). Histograms may be available via the analysis tool pack in some Excel versions. If you did not install office with all options, you will need to have the office installation disk, sorry.

**** If you end up doing research in Science, Medicine, Engineering or working in many other areas you'll need to do statistical analyses. R is free, open source, and probably the most powerful software package available. You can also do publication quality plotting. There is a learning curve at the beginning but it is a worthwhile investment of time! See here, here and here for more info.