PSI-BLAST and BLAST background probabilities

This post is not directly related to MATLAB but I felt it was important to post this.

I recently realized that it is not trivial to find the background amino acid probabilities that are used in BLAST and PSI-BLAST. Google didn't find it. None of the papers referenced in the BLAST papers actually have the frequencies in a tabular form. I would have thought this should have been documented by NCBI in BLAST help or something! Anyway after a few hours of searching and reading papers and eventually code, I found the actual values used. They can be found in this file

http://www.ncbi.nlm.nih.gov/IEB/ToolBox/C_DOC/lxr/source/tools/blastkar.c

Below are the tables which contain the frequencies. They need to be normalized (divide by the sum of the frequencies = 1000) to convert the frequencies to probabilities.

Google doc spreadsheet

NOTE: PSI-BLAST uses the Robinson values by default

 #if STD_AMINO_ACID_FREQS == Dayhoff_prob
 /*  M. O. Dayhoff amino acid background frequencies   */
 static BLAST_LetterProb Dayhoff_prob[] = {
                 { 'A', 87.13 },
                 { 'C', 33.47 },
                 { 'D', 46.87 },
                 { 'E', 49.53 },
                 { 'F', 39.77 },
                 { 'G', 88.61 },
                 { 'H', 33.62 },
                 { 'I', 36.89 },
                 { 'K', 80.48 },
                 { 'L', 85.36 },
                 { 'M', 14.75 },
                 { 'N', 40.43 },
                 { 'P', 50.68 },
                 { 'Q', 38.26 },
                 { 'R', 40.90 },
                 { 'S', 69.58 },
                 { 'T', 58.54 },
                 { 'V', 64.72 },
                 { 'W', 10.49 },
                 { 'Y', 29.92 }
         };
 #endif

 #if STD_AMINO_ACID_FREQS == Altschul_prob
 /* Stephen Altschul amino acid background frequencies */
 static BLAST_LetterProb Altschul_prob[] = {
                 { 'A', 81.00 },
                 { 'C', 15.00 },
                 { 'D', 54.00 },
                 { 'E', 61.00 },
                 { 'F', 40.00 },
                 { 'G', 68.00 },
                 { 'H', 22.00 },
                 { 'I', 57.00 },
                 { 'K', 56.00 },
                 { 'L', 93.00 },
                 { 'M', 25.00 },
                 { 'N', 45.00 },
                 { 'P', 49.00 },
                 { 'Q', 39.00 },
                 { 'R', 57.00 },
                 { 'S', 68.00 },
                 { 'T', 58.00 },
                 { 'V', 67.00 },
                 { 'W', 13.00 },
                 { 'Y', 32.00 }
         };
 #endif

 #if STD_AMINO_ACID_FREQS == Robinson_prob
 /* amino acid background frequencies from Robinson and Robinson */
 static BLAST_LetterProb Robinson_prob[] = {
                 { 'A', 78.05 },
                 { 'C', 19.25 },
                 { 'D', 53.64 },
                 { 'E', 62.95 },
                 { 'F', 38.56 },
                 { 'G', 73.77 },
                 { 'H', 21.99 },
                 { 'I', 51.42 },
                 { 'K', 57.44 },
                 { 'L', 90.19 },
                 { 'M', 22.43 },
                 { 'N', 44.87 },
                 { 'P', 52.03 },
                 { 'Q', 42.64 },
                 { 'R', 51.29 },
                 { 'S', 71.20 },
                 { 'T', 58.41 },
                 { 'V', 64.41 },
                 { 'W', 13.30 },
                 { 'Y', 32.16 }
         };
 #endif

 static BLAST_LetterProb nt_prob[] = {
                 { 'A', 25.00 },
                 { 'C', 25.00 },
                 { 'G', 25.00 },
                 { 'T', 25.00 }
         };

Appending to .MAT files

You can append variables to a .mat file using

>> save(oFname,'var','-append');

Consider 2 scenarios:
1) The variable 'var' is being added to the .mat file for the first time
2) The variable 'var' already exists in the .mat file and is being overwritten or updated

If 'var' takes up a lot of memory ie it is large matrix or array, (2) is significantly slower than (1) by orders of magnitude.

Moral of the story: As far as possible avoid overwriting or updating a variable in a .mat file, especially if the variable takes up a lot of memory.

Sparse vectors - ALWAYS use Column Vectors

I was working on some 'signal' data that I obtained from a ChIP-seq experiment that measures the binding affinity of a transcription factor to every nucleotide in the human genome. I was trying to manipulate this signal data using sparse vectors in MATLAB.

Most of the time I use column vectors by default. For some reason I decided to switch to row vectors. What a difference!

An empty (all-zeros) sparse column vector of length 2 million barely takes a few bytes of memory. However, an empty sparse row vector of the same length gives an 'out of memory' error. While I was aware of the space efficiency of column-based sparse matrices in MATLAB, this was the first time I actually observed such a vast difference.

Moral of the story: If you are manipulating sparse vectors ALWAYS use column vectors!

Dealing with massive files with limited memory

When dealing with extremely massive files such as entire genomes, it is pretty much impossible to fit it all in memory. For situations like this MATLAB has an extremely slick function called memmapfile.

The main advantages are

• The file is not loaded in memory
• You can access the entire file or a portion of the file as if it were a standard MATLAB array using indexing operations. Let say the file had the sequence for an entire genome. Now if you say a = memmapfile('genome.dat') then doing something like a.Data(1:10) gives you the first 10 nucleotides of the genome.
• It can handle single formats or multiple formats
  
• Much faster than fread and fwrite.

This is extremely useful for handling large binary files.

Vectorized ROC curve code + AUC

ROC curves are often used to display the predictive performance of binary classifiers. The area under the ROC curve (AUC) is a way to compare various classifiers. A perfect classifier has an AUC of 1 and a completely bogus (random) classifier has an AUC of 0.5. You can read more about ROC curves here.There is a ton of code for plotting ROC curves and calculating AUC. But most use 'for' loops. And as we all know, loops slow everything down in MATLAB. You can download my vectorized code for plotting multiple ROC curves from multiple classifiers and calculating AUC curves for each.

Download Link