#' The continuous Hartley transform
#'
#' This function performs linear interpolation between data points and implements the continuous Hartley transform
#' @param w This is omega_0 (Equation 4, Anderson et al., 2017). Defaults to NULL.
#' @param t Vector of measurement times. Defaults to NULL.
#' @param x Vector of measurement values. Defaults to NULL.
#' @param Harg0 If omega_0 is equal to zero (x=0), set omega_0 to Harg0 (w=Harg0). Defaults to NULL.
#' @return A vector with the transformation of the measurement data into the frequency domain.
#' @export
#' @examples This function is called by HMF1compute() and HMFcompute(). This function is not recommended to be used in isolation.
#' Hcompute(omega, time, expr, H0)
#' @references \url{http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005627}
Hcompute <- function(w=NULL,t=NULL,x=NULL,Harg0=NULL) { # continuous Hartley transform approximation
H <- 0 # initiallized CHT, to be updated in loop over all intervals
for (i in 1:(length(t)-1)) {
xf <- x[i+1]
xi <- x[i]
tf <- t[i+1]
ti <- t[i]
mm <- (xf - xi)/(tf - ti) # slope for linear interpolation
b <- x[i]-mm*t[i] # intercept for linear interpolation
if (w==0) {w = Harg0}
Imtcoswtdt <- (mm/(w^2)) * ( (cos(w*tf)+w*tf*sin(w*tf)) - (cos(w*ti)+w*ti*sin(w*ti)) )
Ibcoswtdt <- (b/w) * ( sin(w*tf) - sin(w*ti) )
Imtsinwtdt <- (mm/(w^2)) * ( (sin(w*tf)-w*tf*cos(w*tf)) - (sin(w*ti)-w*ti*cos(w*ti)) )
Ibsinwtdt <- (-b/w) * ( cos(w*tf) - cos(w*ti) )
H <- H + Imtcoswtdt + Ibcoswtdt + Imtsinwtdt + Ibsinwtdt
}
return(H)
}
#' HMF spectral component
#'
#' This function computes the mth HMF spectral component of the measurement profile (Equation 5, Anderson et al., 2017).
#' @param meds Matrix with the first column containing measurement times and the second column containing the corresponding measurements. Defaults to NULL.
#' @param m The m value.
#' @param W0 The omega_0 term.
#' @param n Highest order of the system (n=1 for the first degree system that this analysis was designed for).
#' @param Harg0 If omega_0 is equal to zero (x=0), set omega_0 to Harg0 (w=Harg0).
#' @return The spectral component data.
#' @export
#' @examples This function is not recommended to be used in isolation.
#' HMFcompute(data, m, omega, n, H0)
#' @references \url{http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005627}
HMFcompute <- function(meds,m,W0,n,Harg0) { # spectral component, i=1
t <- meds[,1] # time
x <- meds[,2] # expression data, -ddCt (imputed)
i <- 1 # order of the derivative, applied for dE/dt (i=1)
Hhat <- 0 # initialize Hhat, to be updated in loop
for (j in c(0:n)) {
casarg <- (i*pi/2) * ((n+m-j)^i)
cas <- cos(casarg) - sin(casarg)
Harg <- (n+m-j) * W0
H <- Hcompute(Harg,t,x,Harg0)
Hhat <- Hhat + ((-1)^j) * choose(n,j) * H
}
return(Hhat)
}
#' HMF spectral component of the 1st derivative
#'
#' This function computes HMF spectral component of the 1st derivative of x (dx/dt = x', Equation 6, Anderson et al., 2017).
#' @param meds Matrix with the first column containing measurement times and the second column containing the corresponding measurements. Defaults to NULL.
#' @param m The m value.
#' @param W0 The omega_0 term.
#' @param n Highest order of the system (n=1 for the first degree system that this analysis was designed for).
#' @param Harg0 If omega_0 is equal to zero (x=0), set omega_0 to Harg0 (w=Harg0).
#' @return The spectral component data.
#' @export
#' @examples This function is not recommended to be used in isolation.
#' HMF1compute(data, m, omega, n, H0)
#' @references \url{http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005627}
HMF1compute <- function(meds,m,W0,n,Harg0) { # spectral component, i=1
t <- meds[,1] # time
x <- meds[,2] # expression data, -ddCt (imputed)
i <- 1 # order of the derivative, applied for dE/dt (i=1)
Hhat <- 0 # initialize Hhat, to be updated in loop
for (j in c(0:n)) {
casarg <- (i*pi/2)
cas <- cos(casarg) - sin(casarg)
Harg <- ((-1)^i) * (n+m-j) * W0
H <- Hcompute(Harg,t,x,Harg0)
Hhat <- Hhat + ((-1)^j) * choose(n,j) * cas * ((n+m-j)^i) * (W0^i) * H
}
return(Hhat)
}
#' X and Y data
#'
#' This function function for getting the X and Y data (Equation 10, Anderson et al., 2017).
#' @param datHMF The data matrix with scaled measurements measurements. Scaled center measurements are provided for
#' each measurement/annotation (e.g., gene/organ) combination. This format is produced by scale_zeroOne(). Defaults to datHMF.
#' @param Mall A set of m-values. Defaults to NULL.
#' @param Msc A scale factor that determines the M value (M = max(m)) corresponding to the largest range of m-values.
#' Mmax <- ceiling(Msc * Np / 2) where Np = number of parameters.
#' @param n Highest order of the system (n=1 for the first degree system that this analysis was designed for).
#' @param Harg0 If omega_0 is equal to zero (x=0), set omega_0 to Harg0 (w=Harg0).
#' @return The X and Y transformations obtained using Equations 4-6 for all m-value ranges
#' @export
#' @examples XY = getXY(datHMF,Mall=Mall,Msc=Msc,n=n,Harg0=Harg0)
#' @references \url{http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005627}
getXY = function(datHMF=datHMF,Mall=NULL,Msc=2,n=1,Harg0=NULL){
W0 = 2*pi / (max(unique(datHMF[,1]))-min(unique(datHMF[,1])))
if(is.null(Mall) == F){Mall = Mall}
if(is.null(Mall) == T){
mvals = get_M_vals(datHMF=datHMF,Mall=NULL,Msc=2)
Mall = mvals$Mall
}
X <- array( c(0), c(ncol(datHMF)-2, length(unique(datHMF[,2])), length(Mall) ) )
Y <- array( c(0), c(ncol(datHMF)-2, length(unique(datHMF[,2])), length(Mall) ) )
for(anns in 1:length(unique(datHMF[,2]))){ # loop through annotations
ind_ann = which(datHMF[,2] == unique(datHMF[,2])[anns])
for(feats in 1:(ncol(datHMF)-2)){ # loop through features
ind_feat = which(names(datHMF) == names(datHMF)[3:ncol(datHMF)][feats])
mids = datHMF[ind_ann,c(1,ind_feat)]
# loop for HMF calculations
for (mr in 1:length(Mall)) { # loop through every m value
mm <- Mall[mr] # this is the m value
X[feats,anns,mr] <- HMFcompute(mids,mm,W0,n,Harg0)
Y[feats,anns,mr] <- HMF1compute(mids,mm,W0,n,Harg0)
} # mr
} # feats, loop through features
} # anns, loop through annotations
out = list(X=X, Y=Y)
return(out)
} # end getXY()
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.