Scientific toolbox

There are many useful programming tricks that can be used in data analysis and scientific research. We try to present some techniques that we have found useful. Most of the code examples are given in MATLAB , but the ideas can be used in other programming languages as well. It is often very simple tricks that greatly simplify life. That is the reason we want to show short program snippets that produce useful results.

Data scaling in visualization - compressing a large dynamic range

In this example we have a data set from a GC-MS instrument. This "hyphenated" instrument contains a gas chromatograph (GC) linked to a mass spectrometer (MS). The gas chromatograph separates the mixture of molecules using a programmed temperature gradient. The smaller molecules come out (elute) at lower temperatures. Higher temperatures force the larger molecules to elute. In mass spectrometry the sample molecules are "bombarded" with electrons so that the molecules are fragmented to smaller units. These fragments are then observed as m/z values. The way a molecule breaks down is typical for each molecule. The intensities of different fragments form the mass spectrum of the molecule.

We want to see the collected raw data before starting the actual data analysis. We have a two-dimensional data matrix called GCMS. It contains 2000 spectra. Each spectrum has 450 m/z values.

We plot the whole data set with the command:

plot(GCMS)

The first figure (FIG. 1) shows the intensities (abundances) as a function of time. The plot command in MATLAB (or Octave) displays all vertical columns of the matrix. We can see that some peaks are very large, but not much else. The problem is that all the intensities of different masses are overlapping in the same two-dimensional plot.

The shape of the plot resembles the shape of the total ion chromatogram (TIC). We can plot the TIC curve by summing the intensities of the fragments belonging to all spectra. This can be done by the following command:

plot(sum(CGMS'))

Here the ' symbol means the transpose of the matrix GCMS. If we forget to take the transpose of the matrix we get the sum spectrum of all the compounds in the run.

If we want to see where the different masses are located in the data, we can use the special plotting command for sparse matrices. The spy command produces a plot that shows only the intensities that are non-zeroes:

spy(sparse(GCMS))

The results in FIG. 2 show where the masses are located. The first GC-MS spectra contain compounds with a lower molecular weight. When the instrument elutes the high-molecular compounds using a higher temperature we can see fragments with higher masses.

If we want to get an idea about the intensities we can try to use the command imagesc. This transforms the ion intensities into colors. When we use unscaled data with the command:

imagesc(GCMS)

we see only a very small number of points. This is visible in FIG. 3.

This is due to the large dynamic range in the mass spectra. The largest peaks can be a million times larger than the smallest peaks. If we want to get an idea about the intensities we must first reduce the dynamic range. This can be done very simply by raising each intensity to a fractional power. If we choose the exponent value to be 0.25 this means that we are taking twice the square root for each intensity. We can reduce the dynamic range by using the command:

imagesc(GCMS.^0.25)

FIG. 4 gives some idea where the largest intensities are located. Low intensities are shown in dark blue. Light blue areas correspond to higher intensities. Actually there are some red points corresponding to very high intensities, but they are not easily seen. The data has been collected starting from mass 46, that is the reason why there is a homogenous dark blue area on the left side of the figure. All the intensities there are equal to zero.