# Example of processing time series set
# efg, Stowers Institute for Medical Research, July 2005

source("ProcessLombScargleDataset.R")

# ====================================================================
# 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)


# ====================================================================
# 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)])

# 2.3  Define time points for time series (need not be uniform)
# (do not use "t" for time here since it's reserved in R for "transpose")
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"

# 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))

# 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=""))
}


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