inst/doc/LPWC.R

In LPWC: Lag Penalized Weighted Correlation for Time Series Clustering

## ----setup, echo=FALSE, results="hide"-----------------------------------
knitr::opts_chunk$set(comment = "#>", collapse = TRUE)

## ---- eval = FALSE-------------------------------------------------------
#  install.packages("LPWC")

## ---- eval = FALSE-------------------------------------------------------
#  devtools::install_github("gitter-lab/LPWC@vx.xx.x")
#  #OR
#  devtools::install_version("LPWC", version = "x.x.x", repos = "http://cran.us.r-project.org")

## ---- eval = FALSE-------------------------------------------------------
#  devtools::install_github("gitter-lab/LPWC")

## ----lib, results="asis", eval=TRUE--------------------------------------
library(LPWC)

## ----data, results='markup'----------------------------------------------
data(simdata)
simdata[1:5, ]
str(simdata)

## ----time----------------------------------------------------------------
timepoints <- c(0, 2, 4, 6, 8, 18, 24, 32, 48, 72)
timepoints

## ---- echo = FALSE, fig.height = 2, fig.width=4, warning = FALSE, fig.align='center'----
library(ggplot2)

set.seed(29876)


a <- rbind(c(rep(0, 5), 8, 0), c(rep(0, 4), 4.3, 0, 0)) + rnorm(2, 0, 0.5)

dat <- data.frame(intensity = as.vector(a), time = rep(c(0, 5, 15, 30, 45, 60, 75), each = 2), genes = factor(rep(c(1, 2), 7)))

a2 <- a
a2[1, ] <- c(a2[1, 2:7], NA)
dat2 <- data.frame(intensity = as.vector(a2), time = rep(c(0, 5, 15, 30, 45, 60, 75), each = 2), genes = factor(rep(c(1, 2), 7)))

a3 <- a
a3[2, ] <- c(NA, a3[2, 1:6])
a3[1, ] <- c(a3[1, 2:7], NA)
dat3 <- data.frame(intensity = as.vector(a3), time = rep(c(0, 5, 15, 30, 45, 60, 75), each = 2), genes = factor(rep(c(1, 2), 7)))


plot1 <- ggplot(dat, aes(x= time, y = intensity, group = genes)) + geom_line(aes(color = genes), size = 1.5) +  labs(x = "Time (min)") + labs(y = "Intensity") 
plot1

row1 <- c(0, 5, 15, 30, 45, 60, 75)
knitr::kable(t(data.frame(Original = row1, Gene1 = row1, Gene2 = row1)), align = 'c')





## ---- echo = FALSE, fig.height = 2, fig.width=4, warning = FALSE, fig.align='center'----

plot2 <- ggplot(dat2, aes(x= time, y = intensity, group = genes)) + geom_line(aes(color = genes), size = 1.5) + 
  labs(x = "Time (min)") + labs(y = "Intensity") 
plot2

row2 <- c(5, 15, 30, 45, 60, 75, "-")

knitr::kable(t(data.frame(Original = row1, Gene1 = row2, Gene2 = row1)), align = 'c')


## ---- echo = FALSE, fig.height = 2, fig.width=4, warning = FALSE, fig.align='center'----

plot3 <- ggplot(dat3, aes(x= time, y = intensity, group = genes)) + geom_line(aes(color = genes), size = 1.5) + 
  labs(x = "Time (min)") + labs(y = "Intensity") 
plot3

row3 <- c("-", 0, 5, 15, 30, 45, 60)

knitr::kable(t(data.frame(Original = row1, Gene1 = row2, Gene2 = row3)), align = 'c')


## ------------------------------------------------------------------------
LPWC::corr.bestlag(simdata[49:58, ], timepoints = timepoints, max.lag = 2, penalty = "high", iter = 10)

## ----clust1, fig.width=5-------------------------------------------------
dist <- 1 - LPWC::corr.bestlag(simdata[11:20, ], timepoints = timepoints, max.lag = 2, penalty = "low", iter = 10)$corr
plot(hclust(dist))

## ----clust 2-------------------------------------------------------------
dist <- 1 - LPWC::corr.bestlag(simdata[11:20, ], timepoints = timepoints, max.lag = 2, penalty = "low", iter = 10)$corr
cutree(hclust(dist), k = 3)

## ---- eval = FALSE-------------------------------------------------------
#  indexnon0 <- apply(data, 1, function(x){which(var(x) != 0)})
#  subset.data <- data[indexnon0, ]

## ----sessionInfo---------------------------------------------------------
sessionInfo()

## ---- eval = FALSE-------------------------------------------------------
#  
#  
#  # This function stores two different list separately
#  comb <- function(x, ...) {
#    lapply(seq_along(x),
#           function(i) c(x[[i]], lapply(list(...), function(y) y[[i]])))
#  }
#  
#  
#  # adding the data
#  data <- simdata[1:10, ]
#  # number of iterations
#  iter <- 10
#  # C values that are used in the algorithm
#  allC <- findC(timepoints = timepoints, iter = iter)
#  
#  # setting the clusters
#  core <- parallel::detectCores() - 1
#  cl <- parallel::makeCluster(core)
#  
#  # assigning the parallelization
#  doParallel::registerDoParallel(cl)
#  
#  
#  ## running the algorithm for different C
#  result <- foreach(i = 1:iter, .combine='comb', .multicombine=TRUE,
#                  .init=list(list(), list())) %dopar% {
#              lags <- best.lag(data, max.lag = 3, timepoints = timepoints, C = allC[i])
#              new.data <- prep.data(data = data, lags = lags, timepoints = timepoints)
#              corr <- comp.corr(new.data$data, new.data$time, C = allC[i])
#              return(list(corr, lags))
#                  }
#  
#  
#  # dividing the list into two different list: one for lags and one for all the correlations
#  allcorr <- result[[1]]
#  alllags <- result[[2]]
#  
#  # picking best C
#  val <- rep(NA, (length(iter) - 1))
#  for(i in 1:(iter - 1)){
#    val[i] <- sum((as.vector(allcorr[[i + 1]]) - as.vector(allcorr[[i]]))^2)
#  }
#  
#  # returning the results for the best C
#  result <- list(corr = allcorr[[which.min(val) + 1]], lags = alllags[[which.min(val) + 1]], C = values[which.min(val) + 1])
#