Your name:
Your email address:

Note for those doing this at home:

You need to establish a vpn connection to UConn, before you can log into xanadu. Installing Pulse secure works reliably for this. You can download the software here. Installation instructions are here.

Introduction

This week, we will make a gene plot. (As seen in the lecture.) This is basically a way of visualizing 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 FTP directories (RefSeq or other) for a genome . 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 gen plots for microbial genomes

        (different Frankia, Aeromonas, or different Thermotoga species work nicely).

A) obtaining the genome sequences

See Lab 6 on how to get the required files from the NCBI's current genome list.

The easiest will be to download both the feature table and the faa file for two or more genomes to your computer. Then use FileZilla to transfer the files into a Lab7 directory in your account on the cluster.

1) Download two (or more) genomes from closely related organisms from the NCBI's current genome list - select Prokaryotes and set the filter for complete genomes.
You will need both the _feature_table.txt.gz file and the .faa.gz file for both (or more) genomes. If there is more than one file of this name, pick the largest (the one that contains the data from the chromosome).
2) start filezilla and connect to transfer.cam.uchc.edu. (use you username and password)
3) create a lab7 directory and transfer the files you downloaded from the NCBI into that directory.
4) if you want to, rename the files into something more memorable than GCF_000025685.1_ASM2568v1_protein.faa.Keep the extension of the file the same (.faa.gz and _feature_table.txt.gz). However, this might not be necessary, if you use copy past and the ls command.

(an alternative to download the files on to the cluster is below.)

Start the PuTTY. In the "Host" field, type xanadu-submit-ext.cam.uchc.edu (or if you use terminal on a Mac: ssh mcb3421xxxxx@xanadu-submit-ext.cam.uchc.edu)
In the "Username" field, type your username: mcb3421usrXX, where XX is a number assigned to you. Login. It may ask you to accept a new host key. Now enter your password: .

type 
srun --pty -p mcbstudent --qos=mcbstudent --mem=2G bash
and then hit the return or enter key. 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.)

Use qstat or hostname to verify that you are on a compute node.

change into the directory where the genome liles are located
cd lab7

Uncompress all the files using
gunzip *.gz (gunzip calls the program to unzip files, and * is a wild card the will be expanded into a list of filenames that includes all files ending on .gz)

To make the blast+ programs available:
module load blast

To make perl available in a way that works:
module load perl

To make Gnuplot available
module load gnuplot

Use the copy url commands below to get the following scripts into the lab 7 directory.

Perl program that substitutes the accession numbers in a protein FASTA file with the corresponding genome (start) coordinates in the feature table:
curl -O https://j.p.gogarten.uconn.edu/mcb3421_2019/labs/faaReplaceAccessionWithStart.pl

Other script we will need are

For removing all blast hits but the best one from the blast output table:
curl -O http://carrot.mcb.uconn.edu/mcb3421_2016/blastTopHit.pl

For making a plot using gnuplot
curl -O http://carrot.mcb.uconn.edu/mcb3421_2017/plotwgnu_mod2.pl

To add the location in the genome to the annotation line use the following command:

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.

If this doesn't work for one of the genomes, pick another genome and send me a link to the featurefile and faa file that failed.

Replace the accession lines in both (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!

makeblastdb -in database_proteinWithStart.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_proteinWithStart.faa -db database_proteinWithStart.faa -out blast.txt -outfmt 6 -evalue 1e-8
The other genome will be the "query". An E-value cut-off of 10^-8 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:

You use it as follows:

perl blastTopHit.pl blast.txt > blastTopHit.txt

head blastTopHit.txt

Note: The "-max_target_seqs 1" option also returns one top BLASTp database hit for each query sequence (for problems with this command see here). However, since we also want all hits with E-value ≤ 10^-8 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. Or save locallyand use filezilla)
curl -O http://carrot.mcb.uconn.edu/mcb3421_2017/plotwgnu_mod2.pl
and use the editor nano to edit the file (if you have named your blast output files blast.txt and blastTopHit.txt, you can run the file without editing).
nano plotwgnu_mod2.pl (if you have named your blast output files blast.txt and blastTopHit.txt, you can run the file without editing).

The script calls a program called gnuplot.
To make it available, load the gnuplot module:
module load gnuplot

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 back to the cluster using Filezilla)

To run the script type
perl plotwgnu_mod2.pl

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

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

Open filezilla and navigate to your lab7 directory, and transfer 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? (Or the data plotted in different colors using gnuplot.)

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 email the resulting plots as attachment to gogarten@uconn.edu.
For each plot, give the name of the strain used as databank (using the plotwgnu_mod2.pl script, the the databank genome ends up on the Y-axis), and the name of the strain used as query.
Discuss what genome rearrangements might hva given rise to this result. Does it appear likley that the origins of replication were placed in non-homologous positions?
Genome1 (as databank): Genome2 (as query): Description of results:     Was the ORI placed in homologous locations?

Type exit
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.