Your name:
Your email address:

Introduction

This week, we will make a gene plot. (As seen in the midterm and lectures.) This is basically a way of visualising the BLASTp hits between proteins (ORFs) in two genomes, in order to compare their relative arrangement (inversions, etc.). One genome is the x-axis, and the other genome is the y-axis. For each point (x,y) on a scatter plot, the following holds:

• x and y are coordinates (base pairs from the origin of replication) of two proteins in their respective genomes
• x and y represents a BLASTp match between the proteins at coordinates x (in genome 1) and y (in genome 2)

From last week, you'll recall that the tabular output (i.e., outfmt option 6) of a BLASTp search between proteins in two genomes (g1 and g2) looks like this:

You'll also recall that we downloaded a "..._protein.faa" file, which was a FASTA format file of all the encoded proteins in a given genome.

The FASTA format has a ">" line, followed by each sequence. The NCBI ..._protein.faa files look something like this:

>proteinAccessionA(space)some other descriptive text
M...the rest of the protein for the accession on the previous line, on one or more lines...
>proteinAccessionB(space)some other descriptive text
M...the rest of the protein for the accession on the previous line, on one or more lines...

There are several files in the RefSeq FTP directory for a genome (the "green diamond link"). Today we will be using the two shaded in blue:

Adding gene coordinates to BLASTp output

Given two genome (g1 and g2) FASTA protein files, our BLASTp output might look as follows:

>g1p1
MLAMARCK
>g1p2
MDARWIN
>g1p3
MWALLACE

>g2p7
MLAMARCH
>g2p8
MDARWEN
>g2p9
MWALLASE

We would then need some way to add the coordinates for the query and database accessions to the BLASTp output. This information is in the respective feature tables.

One simple way would be to change the header lines in our FASTA protein files to a genome coordinate instead of an accession. For this example we will choose the start coordinate.

>10
MLAMARCK
>40
MDARWIN
>70
MWALLACE

>20
MLAMARCH
>50
MDARWEN
>80
MWALLASE

Then we proceed to do a scatter plot of the first two columns.

There are many ways to accomplish this, but today we will use option 1.

1. Use a Perl program to substitute the accession numbers in the FASTA files for each genome with the corresponding genome coordinate in the feature table, then BLAST as usual.
2. Use the command-line program "join". See Appendix 1 (if you dare).
3. Use Excel to merge the BLASTp output with the feature tables. This is called "Get & Transform".

Using Blast to do genome plots for microbial genomes (different Frankia, Aeromonads, or different Thermotoga species work nicely).

A) obtaining the genome sequences

Go to the NCBI's current genome list.

Click on the "Prokaryotes" tab. Then untick all the "Levels:" (at the top right) except for "Complete". This should result in a listing of completely sequenced genomes.

Download three complete Aeromonas genomes.

There are about 28 completely sequenced Aeromonas genomes in the databank. Pick one of them to turn into a databank, and two other ones that you will use as queries. Pick a closer and a more distantly related one. The one to use as databank will be assigned. The other ones are your choice (to make an educated guess what is closely related and what not see here). To do this, use the Group and SubGroup boxes to narrow your selection, or better yet, use the search box (the box to the left of "Search by organism") to narrow to a taxonomic group.

Then look for the green diamond link in the far right-hand column (if there is only a blue one, take the blue one). This green diamond takes you to a listing of all the files for a genome project.

You want to download TWO files for each genome — the protein.faa file, and the feature_table file (i.e. a minimum of 6 files).

Right-click on the faa file, and select "Copy link location". Then go to the cluster and paste the link, like this:

Start a program called Bitvise SSH Client. In the "Host" field, type bbcsrv3.biotech.uconn.edu
In the "Username" field, type your username. Login. It may ask you to accept a new host key. Now enter your password.

mkdir lab7

cd lab7

curl -O paste the link you copied

You will need to download six files in total:

curl -O your first genome protein faa file link
curl -O your first genome feature table file link
curl -O your second genome protein faa file link
curl -O your second genome feature table file link
curl -O your third genome protein faa file link
curl -O your third genome feature table file link

ls

The "mv" command can be used to rename a long, unwieldy filename into a more concise one.

mv   oldname   newname

mv   ridiculouslyLongFilename.faa.gz   somethingShorter.faa.gz

A good convention would be the first letter of the genus, and then the species name and the strain designation, e.g., A_hydrophila_ML09-119.faa.gz and A_hydrophila_ML09-119_feature.txt.gz
Remember to NOT use spaces in the names

Once you did this for all the files, you can unzip them using

gunzip *.gz (The star is a wild card, the shell will look for all files whose name ends with .gz and unpack them

ls

The file should now end with just "faa" or "txt". The "gz" ending was a compressed file format. Take a look at the first few lines with

head somethingShorter.faa

Here is a Perl program that substitutes the accession numbers in a protein FASTA file with the corresponding genome (start) coordinates in the feature table:

curl -O http://carrot.mcb.uconn.edu/mcb3421_2017/faaReplaceAccessionWithStart.pl

You use it as follows:

perl   faaReplaceAccessionWithStart.pl   yourGenome1.faa   yourFeatureTable1.txt   >   yourGenome1WithStart.faa

perl   faaReplaceAccessionWithStart.pl   yourGenome2.faa   yourFeatureTable2.txt   >   yourGenome2WithStart.faa

perl   faaReplaceAccessionWithStart.pl   yourGenome3.faa   yourFeatureTable3.txt   >   yourGenome3WithStart.faa

head yourGenome1WithStart.faa

Check that the ">accession" lines have been replaced with ">number" lines.

Replace the accession lines in all three (or more) genome FASTA files. Be careful that you use the corresponding feature table for each genome. Also note that if a file exists with the same name as that to the right of the "screen output" redirect (>) symbol, it will be replaced!

qlogin
This takes you to a "compute" node. (Why? Because we are going to run the BLAST+ command, and we want to "farm" the processing out to another computer, rather than hammering the single computer which operates as the "gateway" for everyone. See cluster etiquette.)

cd lab7
makeblastdb -in database_protein.faa -dbtype prot -parse_seqids
Choose the first of your genomes as the "database". Do an "ls" to see the extra files you just made. Use the FASTA file with the start coordinates!

blastp -query query_protein.faa -db database_protein.faa -out blast.txt -outfmt 6 -evalue 1e-8
The other genome will be the "query". An E-value cut-off of 10^-4 is used. Use the FASTA files with the start coordinates.

This will take a few minutes. Again, here is a description of the columns.

To get just the top hit for each query sequence, we use another Perl program. Since the hits for each query are ordered by best E-value to worst, the top hit is simply the first hit for each query:

curl -O http://carrot.mcb.uconn.edu/mcb3421_2016/blastTopHit.pl

You use it as follows:

perl blastTopHit.pl blast.txt > blastTopHit.txt

head blastTopHit.txt

Note: The "-max_target_seqs 1" option also returns the top BLASTp database hit for each query sequence. However, since we also want all hits with E-value ≤ 10^-4 for the other plot, we can use the Perl program to avoid computing the BLASTp twice (once without the max_target_seqs option, and once with).

To plot the location in one genome against the location of the matches in the other genome we have two options. (B) using excel (see below) or (A) using gnuplot.

A) Plotting the results using gnuplot

Gnuplot is installed on the cluster, and you can use it to create scatter plots for both of the matches (top and all significant) in the same coordinate system.
A script that does this is here. The script needs to be present in the same folder as the files with the blast output (blast.txt and blastTopHit.txt)
You need to open the file in a text editor, and add the names of the file with the blast output - if you used the names blastTopHit.txt and blast.txt the program runs without editing.

(there are many ways to download and edit the perl script. One is to use curl to get the file
curl -O http://carrot.mcb.uconn.edu/mcb3421_2017/plotwgnu_mod2.pl
and use the editor nano to edit the file
nano plotwgnu_mod2.pl

An alternative is download the file to your desktop (rightclick on the link and select save as), edit the file on you desktop (MSWord, save as txt), and transfer it to the cluster using the FTP window in BitVise)

To run the script type
perl plotwgnu_mod2.pl

Transfer the resulting plot to your computer using the SFTP window in BitVise, and display the image on the screen. Remember to rename the files (blast.txt .... plot.png before you run the second analysis)

==============================================
B) plotting the results using Excel:

Open a SFTP window in BitVise, navigate to your lab7 directory, and transfer over the blast.txt and blastTopHit.txt files to your Desktop.

Make an Excel scatter plot using all BLAST hits with E-value ≤ 10^-8, and another using just the top hits.

What if any is the difference between the two plots?

Remember to rename the files (blast.txt .... before you run the second analysis)
===============================================

Which genomes did you compare? Describe the results you obtained in words AND the this description and the resulting plots as attachment to gogarten@uconn.edu. For each plot, give the name of the strain used as databank, and the name of the strain used as query.

Type logout
You will return to the master computer. Now type qstat
You should have no jobs running. If there is a job listed, then type qdel
followed by a space, and then the job-ID number, followed by return or enter key. Then type qstat
again, to confirm that there are no running jobs. Then type logout
to exit the main computer.

Finished?

Type logout to release the compute node form the queue.
Check the queue for abandoned sessions using qstat.
If there are abandoned sessions under your account, kill them by deleting them from the queue by typing qdel job-ID, e.g. "qdel 40000" would delete Job # 40000

Check the appropriate radio button below before pressing the submit button:

Send email to your instructor (and yourself) upon submit
Send email to yourself only upon submit (as a backup)
Show summary upon submit but do not send email to anyone.