Your name: Your email address:
A) obtaining the genome sequences
When searching the Internet, you may find references to both the BLAST (classic) and BLAST+ (new) command-line tools. If you have used the old BLAST tools, you may find this Quick Start Guide for switching from BLAST to BLAST+ command line tools useful.
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.
Take a note of the number of complete genomes (Items 1 - 100 of X) here:
Click on "Download Reports from FTP site" (top right of page). Download "prokaryotes.txt" to your computer, or get a local copy here. (Save target as, or Save link as, depending on your browser.) Save to the Desktop.
Load "prokaryotes.txt" into Excel. Start with a blank workbook, then File -- Open -- Browse -- Select "All Files". Choose prokaryotes.txt
Text Import Wizard opens. Just choose "Finish", as the defaults work fine.
Click on the header of column "A". The status bar should show about 110841 rows.
Click on the first cell (A1). Then select "Format as Table". The header lines should now appear in table format.
Scroll to the right until you see a "Status" heading. Click on Status, and tick only "Complete Genome". Check status bar for the number of records found. Does it agree with the previous number?
Sort by Release Date. How many complete genomes have been published so far this month (October)?
Which genome has the most genes? The highest GC content?
Go back to the NCBI's current genome list.
Download two genomes of your choice. They should be from distantly related organisms (e.g. one Archaeon, one Bacterium, or a DPANN archaeon and a Euryarchaeote. If in doubt of what to choose, go for one of the Nanohaloarchaea archaeon SG9 and Haloferax volcanii DS2 genomes.
To do this, use the Group and SubGroup boxes to narrow your selection, or better yet, use the search box to narrow to a taxonomic group.
Then look for the green diamond link in the far right-hand column. This green diamond takes you to a listing of all the files for a genome project.
You want to download the protein.faa file for each genome. This "faa" file contains all of the proteins coded by a genome.
If you use one of the genomes for the second exercise also download the feature_table.txt
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. It might help further on, if you tell the SSH Client to use lines longer than 80 characters in the terminal.
mkdir lab6
cd lab6
curl -O paste the link you copied
e.g. curl -O ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/001/761/425/GCA_001761425.1_ASM176142v1/GCA_001761425.1_ASM176142v1_protein.faa.gz and curl -O ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/025/685/GCF_000025685.1_ASM2568v1/GCF_000025685.1_ASM2568v1_protein.faa.gz
ls
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, e.g., D_radiodurans.faa.gz
gunzip somethingshorter.faa.gz
The file should now end with just "faa". The "gz" ending was a compressed file format. Take a look at the first few lines with
head somethingshorter.faa
And repeat the procedure for the second genome.
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 lab6 makeblastdb -in database_protein.faa -dbtype prot -parse_seqids Choose one of your genomes as the "database". Do an "ls" to see the extra files you just made. -parse_seqids directs the program to create an index that allows to retrieve seuences from the databank.
blastp -query query_protein.faa -db database_protein.faa -out blast.txt -outfmt 6 The other genome will be the "query". -outfmt 6 specifies a tabular output format.
The blast program will now take every sequence in query file and do a blast search against the database (i.e., this performs a couple of thousand blast searches). This will take a few minutes. While waiting help your neighbor :)
Here is a description of the columns in the outfit -6. To add these to the blast output file,
curl -O http://carrot.mcb.uconn.edu/mcb3421_2017/header.txt cat header.txt blast.txt > blast.out.txt
Open a new SFTP window in BitVise, navigate to your lab6 directory, and transfer over the blast.out.txt file to your Desktop. Load it into Excel.
What is the accession ID pair of the most conserved protein between your two genomes?
Copy the Accession Number of the best match in the database, go back to the terminal window and type
fastacmd -s IDnumber -d DatabaseName for example
fastacmd -s AOV94840.1 -d H_volcanii_DS2.faa
This should return the fasta formated sequence of the match.
What are the IDs of the best match, and what is given as function in the fasta annotation line?
In Excel, select the third column from the left (this is the percent identity for each BLASTp match), and Insert -- Statistic Chart (all-blue column chart icon in the Charts section) -- Histogram
(More about histograms in Excel, including the formula used for default binning.)
You get a histogram of all percent identities. You may find it useful to convert to table format, and add headings. Repeat the histogram for E-values ≤ 1, and E-values ≤ 10-3.
Add a column to the spreadsheet that calculates the number of identical residues (column 3 * column 4 /100). Make histograms for the number of identical residues for all E-values, better than 1 and better than E-3.
Do you observe a smooth distribution of % identity values? Do any genes constitute a second peak of more similar sequences, compared to the bulk of the comparisons? What can explain the difference in the histograms for significant and insignificant hits?
Can you explain, why % identity is not a good criterion to distinguish significant from insignificant hits?
For the following questions, the BLAST Command Line Applications User Manual may be helpful → link
How can you get information on the possible parameters in blastp using the commandline? (try blastp -h first) How can you set the wordsize to 2 ? How can you filter the query sequence for regions with low complexity?
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.
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.