Step-by-Step Instructions:
How to Apply Lomb-Scargle Analysis to Your Set of Expression Data

The file Bozdech03.R was used to process the Plasmodium dataset in the file Complete_Dataset.csv.  That R script is not very general, and is likely not a good starting point for processing your dataset.  Instead, download and modify the file AnalyzeMyTimeSeriesSet.R to process your data.  The lines in this file are explained below.

The following assumes you have already extracted the various .R files from the .ZIP (or .gz) on the R Code page.  At a minimum you need these files in your working directory: 

• AnalyzeMyTimeSeriesSet.R (or name of your choosing),
• ProcessLombScargleDataset.R,
• LombScargle.R,
• R-Mnemonics.R, and
• R-efg.R. 

The first important line in AnalyzeMyTimeSeriesSet.R directly or indirectly references all the above R code:

source("ProcessLombScargleDataset.R")

Step 1.  Create a comma-separated value (.csv) file containing a time series in each row for a particular probe from your microarray study.  (Alternately, you can create an Excel spreadsheet with these values, and create a .csv from that, or modify R code to read the Excel file directly.)  Ideally, data should be collected over a time period of more than one cycle.

Your file should have this structure of N time points for each of the K probes:

ID	TP1	TP2	TP3	TP4	...	TPN
Probe1	value11	value12	value13	value14	...	value1N
Probe2	value21	value22	value23	value24	...	value2N
...	...	...	...	...	...	...
ProbeK	valueK1	valueK2	ValueK3	valueK4	...	valueKN

Column labels are assumed to be present in row 1, and row labels in column 1 are assumed to be the various probe names from a microarray study. These values are assumed to be log2 ratios.

We have worked with time series where N ranges from approximately 20 to 50, and where there are about 5,000 to 45,000 probes.   In theory, any size dataset that fits in memory could be processed.  Obviously, it takes more time to process more probes.

For the following discussion, the LombScargleExampleSet.csv file was created using SetupLombScargleExample.R.  (You may want to consider processing a test dataset of known periodic patterns to help understand how "perfect" periodic genes would perform.)

Here is the code form AnalyazeMyTimeSeriesSet.R that loads the raw expression file into a data frame called "RawData":

# 1. Read data into Data Frame called RawData
filename <- "LombScargleExampleSet.csv"
# (use as.is=TRUE if IDs in column 1 are NOT unique
RawData <- read.csv(filename)

Use your filename instead of LombScargleExampleSet.csv in the assignment to filename.  If you have a number of expression datasets and would like to be prompted for which one to process, use this R statement:

filename <- file.choose()

If you have an Excel file, LombScargleExampleSet.xls, instead of a .csv file, the following usually works (problems with some Excel files can be avoided by always using .csv files instead):

library(RODBC)
connection <- odbcConnectExcel("C:/LombScargle/LombScargleExampleSet.xls")
RawData <- sqlFetch(connection, "Set1") # Tabsheet name
close(connection)

Note that a forward slash, "/", can be used in R for either Windows or Linux paths.

You can use the R dim command to verify if the correct number of rows and columns are present in the data frame:

dim(RawData)
[1] 11 26

Step 2.  Define data needed for Lomb-Scargle Analysis.  In some cases the RawData from Step 1 can be directly used in subsequent analysis.  You may have raw expression ratios and may need to take logarithms before proceeding.  You may want to perform data "quality" tests, or check for "outliers" before continuing. 

In this case the RawData contains log ratio values, so the ID and Expression objects can be extracted from the RawData data frame:

# 2. Define data elements needed for Lomb-Scargle Analysis

# 2.1 First column has IDs
ID <- as.character(RawData[,1])

# 2.2 Remaining columns have expression time series
Expression <- data.matrix(RawData[,2:ncol(RawData)])

If you have ratio data and need to compute log ratios, you may want to replace 0 values with small values to avoid the -Inf values from log2(0):

if ((ZeroCount <- sum(Expression == 0.0)) > 0)
{
   cat("NOTE: Replacing ", ZeroCount,
       "0s with small values (1e-15)\n")
   Expression[Expression == 0.0] <- 1e-15
}
Expression <- log2(Expression)

This last assignment also converts RawData to an Expression matrix of double precision floating-point values, which is more computationally efficient to work with than a list:

> typeof(RawData)
[1] "list"
> class(RawData)
[1] "data.frame"

> typeof(Expression)
[1] "double"
> class(Expression)
[1] "matrix"

Step 3.  Define the time points for the series.  Unlike Fourier analysis, time values with Lomb-Scargle analysis need not be uniform.  In this case, the time values were originally defined from uniform random numbers over the interval 0 to 24 hours.

# 2.3 Define time points for time series (need not be uniform)
Time <- c(1.70992,  2.92826,  5.59759,  5.70506,  7.29704,
          9.13245,  9.93785, 10.16164, 10.68415, 11.24294,
         12.10074, 14.19195, 14.33209, 16.28018, 17.44977,
         18.05868, 18.15130, 19.77388, 20.28769, 20.91551,
         20.95605, 21.48596, 22.62380, 23.10384, 23.96810)

# Check for consistency between Expression and Time
N <- ncol(Expression) # number of time points in series
stopifnot(length(Time) == N)

# Define global value used by Lomb-Scargle to label plots
unit <<- "hour"

The consistency check above ensures the number of defined times matches the number of values in each of the time series.  While a bit of kludge, the global assignment to "unit" changes text in certain plots from "unit" to "hour".

Step 4.  Define the test frequencies.  Unlike Fourier analysis where only certain "Fourier frequencies" are usually considered, the range and number of test frequencies need to be defined for the Lomb-Scargle analysis. 

A rule of thumb from [Press2002] suggests that 2N or 4N test frequencies are usually adequate. 

These test frequencies could range up to the Nyquist limit, but usually practical considerations limit this range.  For example, if you are observing a biological phenomenon that has a period of 120 minutes, perhaps a range of 60 to 180 minutes would be adequate.  In this case with a dominant frequency of 24 hours, let's look at the range of frequencies corresponding to periods of 6 to 48 hours. 

# 2.4 Define number of test frequencies: normally 2N or 4N
M <- 4*N

# Frequencies usually easier to define in terms of 1/Period.
# In this case since "dominant" frequency is expected to
# correspond to 24-hour period, let's search frequencies
# corresponding to periods from 6 hours to 48 hours.
MinFrequency <- 1/48
MaxFrequency <- 1/6# See if MaxFrequency is perhaps too high
if (MaxFrequency > 1/(2*mean(diff(Time))))
{
  cat("MaxFrequency may be above Nyquist limit.\n")
}

TestFrequencies <- MinFrequency +
                  (MaxFrequency - MinFrequency) *
                  (0:(M-1) / (M-1))

The regression formula by [Horne86] is used to estimate the number of these test frequencies that are independent, which is used in the computation of p values.  Missing numeric values in a time series will automatically be excluded from the Lomb-Scargle analysis and affect the computation of the number of independent frequencies, too, since it's based on sample size, N.

Step 5.  Define directory to hold results.  A directory is used to hold the output files.

# 2.5 Define relative path for results and make
# sure directories exist
OUTPUT_PATH <- "LombScargleResults"
SLASH <- .Platform$file.sep
if (! file.exists(OUTPUT_PATH) )
{
  dir.create(OUTPUT_PATH)

  # only needed if SavePlots=TRUE (see below)
  # The "chart" directory will contain PDFs of
  # the Lomb-Scargle analysis
  dir.create(paste(OUTPUT_PATH, SLASH, "charts", sep=""))
}

As mentioned in the comments above, the "chart" directory is only used if SavePlots=TRUE is specified below.  Normally, you would use SavePlots=FALSE if you are processing a large number of probes.  If some other plot format is desired instead of PDFs, change this statement in ProcessExpressionDataset function in the file ProcessLombScargleDataset.R:

if (SavePlots)
{
   WriteStampedPDF(paste(path, .Platform$file.sep,
                   "charts", .Platform$file.sep,
                   ID[j], ".pdf", sep=""))
}

Step 6.  Process the time series in the Expression matrix:

# 3. Process the time series
ProcessExpressionDataset(OUTPUT_PATH,
   Time, ID, Expression, TestFrequencies, SavePlots=TRUE)

This step and all the steps above are executed as a batch using this R statement:

source("AnalyzeMyTimeSeriesSet.R")

By default the graphic for each probe is displayed as it is being processed.  Three output files are created in the specified OUTPUT_PATH:  (The links here show the results)

• Dominant.CSV
• Periodogram.CSV
• pValue.CSV

Since SavePlots=TRUE in this example, a folder called charts is created in the OUTPUT_PATH with one PDF file per probe:
Probe1  Probe2  Probe3  Probe4  Probe5  Probe6  Probe7  Probe8  Probe9  Probe10  Probe11

To speed processing, especially with large Expression datasets, you may want to use SavePlots=FALSE to suppress these PDF files. 

In addition, you can suppress the display of graphics and save considerable processing time by defining these two global variables:

SHOW.LOESS.PEAK <<- FALSE
SHOW.PLOT       <<- FALSE

The loess peak is used to define the "phase shift" (see below) and this is not computed when SHOW.LOESS.PEAK is set to FALSE.  An additional time savings can be achieved by suppressing all displayed graphics (and the loess computation) by setting SHOW.PLOT to FALSE. These statements can be directly edited in the LombScargle.R file, or set after this file is source, such as in the file ProcessLombScargleDataset.R.   (Note:  Always specify SavePlots=FALSE when specifying SHOW.PLOT=FALSE).

Step 7.  Analyze Results -- Periodicity.   The individual plots show the result for each probe, but the Dominant.CSV file can be used for overall summaries.  The period and p-value columns from this file (shown below) indicate whether a gene is periodic.  Because this is an artificial example defined by the SetupLombScargleExample.R file, the Expected Period column was added to show expected results were obtained:

probe	Period [Hours]	p	Expected Period [Hours]
Series1	48.00	NaN	constant, not periodic
Series2	7.97	0.00017	8
Series3	11.88	0.00041	12
Series4	24.12	0.00015	24
Series5	25.01	0.00017	24
Series6	23.29	0.00021	24
Series7	22.52	0.00016	24
Series8	24.12	0.00015	24
Series9	23.29	0.00021	24
Series10	7.88	0.02330	Mixture of 8 and 40
Series11	15.04	0.00064	Mixture of 12 and 18

The probe for Series1 shows a period of 48 hours, but this should be ignored because of the NaN (not a number) p-value. The remaining individual p-values are fairly good, but should be considered in the context of the multiple testing procedure outlined in our paper.

For Series 10 both the peaks corresponding to periods of 8 hours and 40 hours were significant and about the same, but the procedure only picks the highest one, which corresponds to the 8-hour period. This is a limitation of the existing program.

For Series 11 two peaks were not seen since the two frequencies were quite close. This is another limitation:  frequencies corresponding to 8 and 40 hour periods can be resolved, but frequencies corresponding to periods of 12 and 18 hours cannot be resolved.

Step 8.  Analyze Results -- Phase Shift.  Probes 4 through 9 all have the same computed periods of approximately 24 hours.  A set of probes with the same period can be further categorized by phase shift, which is a measure of whether they are aligned in time.  This phase shift was computed by finding the peak in the loess smoothed curve, which is represented by a dot in the various plots.

The PhaseShift and PhaseShiftHeight values below are from the Dominant.CSV file. 

PhaseShift can be expressed in terms of an offset within the 24 hour period, or in terms of radians.  In this case the PhaseShift can range from 0 to 24 hours, which corresponds to 0 to 2π radians.   The Adjusted Phase Shifts (see below) compare quite favorably with the Expected Phase Shifts.  The adjustment was necessary because periodic functions can be expressed in alternative and equivalent ways.

probe	PhaseShift [Hours] x	Adjusted Phase Shift [Hours] 24 - x	PhaseShift Height	Expected PhaseShift [Hours]	Expected PhaseShift [Radians]
Series4	23.97	0.03	1.10	0	0
Series5	21.18	2.82	0.96	3	π/4
Series6	18.00	6.00	0.98	6	π/2
Series7	14.80	9.20	0.98	8	3π/4
Series8	11.82	12.18	0.97	12	π
Series9	5.64	18.36	0.86	18	3π/2

The phase shift corresponding to the original value in radians in the SetupLombScargleExample.R file can be computed as follows:  (24 - PhaseShift) * (2π/24).

The PhaseShiftHeight values can be used to reject periodic patterns with low amplitude.  These heights can be ignored in this case, since all probes have similar amplitudes.