Closed ccshao closed 8 years ago
Hi Shao,
Thanks for your concern.
Depending on whether there is a separate magnitude esimate for each cell (ie. magnitudes is a matrix) or if there is a common magnitude estimate for all cells (ie. magnitudes is a vector), we perform either matrix based operations or vector based operations.
Please see the full function below. The if("conc.a2" %in% names(models))
checks are nested in the is.matrix(magnitudes)
check and are needed.
scde.failure.probability <- function(models, magnitudes = NULL, counts = NULL) {
if(is.null(magnitudes)) {
if(!is.null(counts)) {
magnitudes <- scde.expression.magnitude(models, counts)
} else {
stop("ERROR: either magnitudes or counts should be provided")
}
}
if(is.matrix(magnitudes)) { # a different vector for every cell
if(!all(rownames(models) %in% colnames(magnitudes))) { stop("ERROR: provided magnitude data does not cover all of the cells specified in the model matrix") }
if("conc.a2" %in% names(models)) {
x <- t(1/(exp(t(magnitudes)*models$conc.a +t(magnitudes^2)*models$conc.a2 + models$conc.b)+1))
} else {
x <- t(1/(exp(t(magnitudes)*models$conc.a + models$conc.b)+1))
}
} else { # a common vector of magnitudes for all cells
if("conc.a2" %in% names(models)) {
x <- t(1/(exp((models$conc.a %*% t(magnitudes)) + (models$conc.a2 %*% t(magnitudes^2)) + models$conc.b)+1))
} else {
x <- t(1/(exp((models$conc.a %*% t(magnitudes)) + models$conc.b)+1))
}
}
x[is.nan(x)] <- 0
colnames(x) <- rownames(models)
x
}
Hope this addresses your question.
Best, Jean
Thanks for the quick reply! It is clear.
if("conc.a2" %in% names(models)) is used twice. Something wrong here?
scde.failure.probability <- function(models, magnitudes = NULL, counts = NULL) { if(is.null(magnitudes)) { if(!is.null(counts)) { magnitudes <- scde.expression.magnitude(models, counts) } else { stop("ERROR: either magnitudes or counts should be provided") } } if(is.matrix(magnitudes)) { # a different vector for every cell if(!all(rownames(models) %in% colnames(magnitudes))) { stop("ERROR: provided magnitude data does not cover all of the cells specified in the model matrix") } if("conc.a2" %in% names(models)) { x <- t(1/(exp(t(magnitudes)_models$conc.a +t(magnitudes^2)_models$conc.a2 + models$conc.b)+1)) } else { x <- t(1/(exp(t(magnitudes)_models$conc.a + models$conc.b)+1)) } } else { # a common vector of magnitudes for all cells
same if else statment here *_
}