#' A function to preprocess the list of expression matrix
#'
#' @description This function will keep the samples that are
#' common across the list of expression matrix,
#' and filter the features that are all zeros across samples,
#' and finally construct a \code{SingleCellExperiment} object
#'
#' @param exprsMat A list or a matrix indicates the expression matrices of the
#' testing datasets (each matrix must be \code{matrix} or
#' \code{dgCMatrix} class)
#' @param return_sce A logical input indicates whether
#' a \code{SingleCellExperiment} object will be return
#' @param assay_matrix A integer indicates which list will be
#' used as `assay` input of `SingleCellExperiment`
#' @param filter_features A logical input indicates
#' whether the features with all zeros will be removed
#' @param rowData A DataFrame indicates the rowData to be stored
#' in the sce object
#' @param colData A DataFrame indicates the colData to be stored
#' in the sce object
#'
#' @return either a SingleCellExperiment object or
#' a preprocessed expression matrix
#'
#' @examples
#' data(CITEseq_example, package = "CiteFuse")
#' sce_citeseq <- preprocessing(CITEseq_example)
#'
#' @importFrom SingleCellExperiment SingleCellExperiment altExp
#' @importFrom Matrix rowSums
#' @importFrom SummarizedExperiment SummarizedExperiment colData rowData
#' @importFrom methods as is
#' @importFrom S4Vectors DataFrame
#'
#' @export
preprocessing <- function(exprsMat = NULL,
return_sce = TRUE,
assay_matrix = 1,
filter_features = TRUE,
rowData = NULL,
colData = NULL) {
if (!any("matrix" %in% methods::is(exprsMat),
"dgCMatrix" %in% methods::is(exprsMat),
"list" %in% methods::is(exprsMat))) {
stop("exprsMat need to be a matrix or a list")
}
if ("list" %in% methods::is(exprsMat)) {
if (any(!unlist(lapply(exprsMat, function(x) "dgCMatrix" %in% is(x))) &
!unlist(lapply(exprsMat, function(x) "matrix" %in% is(x))))) {
stop("Please make sure every expression matrix in
the list are matrix or sparse matrix")
}
}
if (any("matrix" %in% methods::is(exprsMat),
"dgCMatrix" %in% methods::is(exprsMat))) {
exprsMat <- list(exprsMat)
}
common_cells <- Reduce(intersect, lapply(exprsMat, colnames))
if (length(common_cells) == 0) {
stop("There is no common cells in this list...
please check the matrix input")
}
# only keep the samples that are common across the list of exprsMat
exprsMat <- lapply(exprsMat, function(exprs) {
exprs <- exprs[, common_cells]
if (!"dgCMatrix" %in% methods::is(exprs)) {
exprs <- methods::as(exprs, "dgCMatrix")
}
exprs
})
# filter the features with all zeros
if (filter_features) {
exprsMat <- lapply(exprsMat, function(exprs) {
if ("dgCMatrix" %in% methods::is(exprsMat)) {
rowsums <- Matrix::rowSums(exprs)
exprs <- exprs[rowsums != 0, ]
}
exprs
})
}
if (return_sce) {
sce <- SingleCellExperiment(assays = list(counts =
exprsMat[[assay_matrix]]))
list_idx <- seq_len(length(exprsMat))[-assay_matrix]
for (i in list_idx) {
if (is.null(names(exprsMat)[i])) {
name_exprs <- paste("altExp", i, sep = "_")
} else {
name_exprs <- names(exprsMat)[i]
}
SingleCellExperiment::altExp(sce, name_exprs) <-
SummarizedExperiment(assays = list(counts = exprsMat[[i]]))
}
if (!is.null(rowData)) {
if (!all(rownames(sce) %in% rownames(rowData))) {
stop("Some rownames of the assay matrix
does not have info in rowData (rownames of rowData)")
}
SummarizedExperiment::rowData(sce) <-
DataFrame(rowData[rownames(sce), ])
}
if (!is.null(colData)) {
if (!all(colnames(sce) %in% rownames(colData))) {
stop("Some colnames of the assay matrix
does not have info in colData (rownames of colData)")
}
SummarizedExperiment::colData(sce) <-
DataFrame(colData[colnames(sce), ])
}
return(sce)
} else {
return(exprsMat)
}
}
#' readFrom10X
#'
#' A function to read the data from 10X
#'
#'
#' @param dir A character indicates the directory of the 10X files
#' @param type A character indicates the format of the data, sparse or HDF5
#' @param feature_named_by A character indicates whehter the genes
#' will be named by gene_id or gene_symbol
#' @param filter_features A logical input indicates whether the features
#' with all zeros will be removed
#'
#' @return a SingleCellExperiment object
#'
#' @importFrom rhdf5 h5read
#' @importFrom Matrix readMM sparseMatrix
#' @importFrom SummarizedExperiment SummarizedExperiment
#' @importFrom methods as
#' @importFrom utils read.delim
#'
#' @examples
#'
#' \dontrun{
#' tmpdir <- tempdir()
#' tenXdata <- "http://cf.10xgenomics.com/samples/cell-exp/3.1.0/connect_5k_pbmc_NGSC3_ch1/"
#' file <- "connect_5k_pbmc_NGSC3_ch1_filtered_feature_bc_matrix.tar.gz"
#' download.file(paste0(tenXdata, file),file.path(tmpdir, file))
#' untar(file.path(tmpdir,file),
#' exdir = tmpdir)
#' sce_citeseq_10X <- readFrom10X(file.path(tmpdir,
#' "filtered_feature_bc_matrix/"))
#' sce_citeseq_10X
#'}
#'
#' @export
#'
readFrom10X <- function(dir,
type = c("auto", "sparse", "HDF5"),
feature_named_by = c("gene_id", "gene_symbol"),
filter_features = TRUE) {
type <- match.arg(type, c("auto", "sparse", "HDF5"))
feature_named_by <- match.arg(feature_named_by, c("gene_id", "gene_symbol"))
all_files <- list.files(dir)
if (type == "auto") {
if (sum(grepl(c("barcodes.tsv|features.tsv|matrix.mtx"),
all_files)) == 3) {
type <- "sparse"
} else if (any(grepl("\\.h5", all_files))) {
type <- "HDF5"
} else {
stop("No avaliable files are found in the given dir.")
}
}
if (type == "sparse") {
if (sum(grepl(c("barcodes.tsv|features.tsv|matrix.mtx"),
all_files)) != 3) {
stop("No avaliable files are found in the given dir.")
}
}
if (type == "HDF5") {
if (!any(grepl("\\.h5", all_files))) {
stop("No avaliable files are found in the given dir.")
}
}
if (type == "HDF5") {
h5_path <- paste0(dir, list.files(dir, pattern = "\\.h5"))
#rhdf5::h5ls(h5_path, all = TRUE)
barcodes <- as.character(h5read(h5_path,
paste0("matrix", "/barcodes")))
gene_id <- as.character(h5read(h5_path,
paste0("matrix",
"/features/id")))
gene_symbol <- as.character(h5read(h5_path,
paste0("matrix",
"/features/name")))
gene_type <- as.character(h5read(h5_path,
paste0("matrix",
"/features/feature_type")))
features <- data.frame(gene_id = gene_id,
gene_symbol = gene_symbol,
gene_type = gene_type)
counts <- rhdf5::h5read(h5_path, paste0("matrix", "/data"))
indices <- rhdf5::h5read(h5_path, paste0("matrix", "/indices"))
indptr <- rhdf5::h5read(h5_path, paste0("matrix", "/indptr"))
shape <- rhdf5::h5read(h5_path, paste0("matrix", "/shape"))
exprs_mat <- Matrix::sparseMatrix(i = indices + 1,
p = indptr,
x = as.numeric(x = counts),
dims = shape,
giveCsparse = TRUE)
}
if (type == "sparse") {
mat_path <- paste0(dir, "matrix.mtx.gz")
exprs_mat <- Matrix::readMM(mat_path)
exprs_mat <- methods::as(exprs_mat, "dgCMatrix")
barcodes_path <- paste0(dir, "barcodes.tsv.gz")
barcodes <- utils::read.delim(barcodes_path, header = FALSE,
stringsAsFactors = FALSE)
barcodes <- barcodes[, 1]
feature_path <- paste0(dir, "features.tsv.gz")
features <- utils::read.delim(feature_path, header = FALSE,
stringsAsFactors = FALSE)
colnames(features) <- c("gene_id", "gene_symbol", "gene_type")
}
rownames(exprs_mat) <- features[, feature_named_by]
colnames(exprs_mat) <- barcodes
gene_idx <- features$gene_type == "Gene Expression"
adt_idx <- features$gene_type == "Antibody Capture"
exprs_rna <- exprs_mat[gene_idx, ]
exprs_adt <- exprs_mat[adt_idx, ]
rownames(features) <- features$gene_id
sce <- preprocessing(list(RNA = exprs_rna, ADT = exprs_adt),
rowData = features[gene_idx, ],
filter_features = filter_features)
return(sce)
}
#' normaliseExprs
#'
#' A function that perform normalisation for alternative expression
#'
#' @param sce A \code{SingleCellExperiment} object
#' @param altExp_name Name of alternative expression
#' that will be used to perform normalisation
#' @param exprs_value A character indicates
#' which expression value in assayNames is used.
#' @param transform type of transformation,
#' either log or clr (Centered log ratio transform)
#' @param log_offset Numeric scalar specifying the pseudo-count
#' to add when log-transforming expression values. Default is 1
#'
#' @examples
#' data(CITEseq_example, package = "CiteFuse")
#' sce_citeseq <- preprocessing(CITEseq_example)
#' sce_citeseq <- normaliseExprs(sce = sce_citeseq,
#' altExp_name = "ADT",
#' transform = "log")
#' @importFrom SummarizedExperiment assay assayNames
#' @importFrom SingleCellExperiment altExpNames altExp
#' @importFrom Matrix rowMeans
#'
#' @return a SingleCellExperiment object
#'
#' @export
normaliseExprs <- function(sce,
altExp_name = NULL,
exprs_value = "counts",
transform = c("log", "clr", "zi_minMax", "minMax"),
log_offset = NULL
) {
transform <- match.arg(transform, c("log", "clr", "zi_minMax", "minMax"),
several.ok = TRUE)
if (altExp_name != "none") {
if (!altExp_name %in% altExpNames(sce)) {
stop("sce does not contain altExp_name as altExpNames")
}
assaynames <- SummarizedExperiment::assayNames(altExp(sce, altExp_name))
if (!exprs_value %in% assaynames) {
stop("sce does not contain exprs_value as assayNames for altExp")
}
exprs <- SummarizedExperiment::assay(altExp(sce, altExp_name),
exprs_value)
} else {
# if altExp_name is "none", then the assay in SingleCellExperiment
# is extracted (RNA in most of the cases)
assaynames <- SummarizedExperiment::assayNames(sce)
if (!exprs_value %in% assaynames) {
stop("sce does not contain exprs_value as assayNames")
}
exprs <- SummarizedExperiment::assay(sce, exprs_value)
}
if (is.null(log_offset)) {
log_offset <- 1
}
if (transform %in% "log") {
exprs_norm <- log(exprs + log_offset)
}
if (transform %in% "clr") {
exprs_norm <- .clr(exprs)
}
if (transform %in% "zi_minMax") {
exprs_norm <- apply(exprs, 1, .ziMinMax)
exprs_norm <- t(exprs_norm)
}
if (transform %in% "minMax") {
exprs_norm <- apply(exprs, 1, .minMax)
exprs_norm <- t(exprs_norm)
}
if (transform == "log") {
new_assay_name <- "logcounts"
} else {
new_assay_name <- transform
}
colnames(exprs_norm) <- colnames(exprs)
rownames(exprs_norm) <- rownames(exprs)
if (altExp_name != "none") {
SummarizedExperiment::assay(SingleCellExperiment::altExp(sce,
altExp_name),
new_assay_name) <- exprs_norm
} else {
SummarizedExperiment::assay(sce, new_assay_name) <- exprs_norm
}
return(sce)
}
.minMax <- function(x) {
if (max(x) != 0) {
res <- (x - min(x))/(max(x) - min(x))
} else {
res <- rep(0, length(x))
}
return(res)
}
.ziMinMax <- function(x) {
res <- .minMax(x)
res <- ifelse(res < 0.5, 0, res)
return(res)
}
.clr <- function(X) {
X[X == 0] <- min(X[X != 0])
logX <- log(X)
logSet <- logX[, seq_len(ncol(X)), drop = FALSE]
ref <- Matrix::rowMeans(logSet)
res <- sweep(logX, 1, ref, "-")
return(res)
}
#' crossSampleDoublets
#'
#' A function that perform normalisation for alternative expression
#'
#' @param sce A \code{SingleCellExperiment} object
#' @param altExp_name Name of alternative expression that will be
#' used to perform normalisation.
#' If it is NULL, it will set to HTO.
#' @param totalExp_threshold the threshold indicates for the HTO
#' less than this threshold
#' will be filtered from the analysis
#'
#'
#' @return A SingleCellExperiment Object
#'
#' @examples
#' data(CITEseq_example, package = "CiteFuse")
#' sce_citeseq <- preprocessing(CITEseq_example)
#' sce_citeseq <- normaliseExprs(sce = sce_citeseq,
#' altExp_name = "HTO",
#' transform = "log")
#' sce_citeseq <- crossSampleDoublets(sce_citeseq)
#'
#' @importFrom SummarizedExperiment assay assayNames
#' @importFrom SingleCellExperiment altExpNames altExp
#' @importFrom Matrix rowSums
#' @importFrom mixtools normalmixEM
#' @importFrom S4Vectors metadata
#' @importFrom stats kmeans
#' @importFrom methods is
#'
#'
#' @export
crossSampleDoublets <- function(sce,
altExp_name = NULL,
totalExp_threshold = 10) {
if (is.null(altExp_name)) {
if (!"HTO" %in% SingleCellExperiment::altExpNames(sce)) {
stop("There is no HTO data in the object")
} else {
altExp_name <- "HTO"
}
}
if (!"logcounts" %in% assayNames(altExp(sce, altExp_name))) {
warning("HTO does not contain logcounts...
we will perform normaliseExprs() to get logcounts")
sce <- normaliseExprs(sce, altExp_name, "log")
}
hto_cellHash_log <- assay(altExp(sce, altExp_name), "logcounts")
hto_cellHash_log <- hto_cellHash_log[Matrix::rowSums(hto_cellHash_log) >
totalExp_threshold, ]
hto_threshold <- list()
for (i in seq_len(nrow(hto_cellHash_log))) {
vec <- hto_cellHash_log[i,][hto_cellHash_log[i,] > 0]
mixmdl <- try(mixtools::normalmixEM(vec,
fast = TRUE, maxrestarts = 1000,
k = 2, maxit = 10000,
mu = c(0, 10),
lambda = c(1/2),
sigma = rep(2, 2),
ECM = TRUE, verb = FALSE),
silent = TRUE)
if ("try-error" %in% methods::is(mixmdl)) {
km <- stats::kmeans(vec, centers = 2)
hto_threshold[[i]] <- min(max(vec[km$cluster == 1]),
max(vec[km$cluster == 2]))
} else {
hto_threshold[[i]] <- getThreshold(mixmdl)
}
}
names(hto_threshold) <- NULL
hto_threshold <- unlist(hto_threshold)
hto_cellHash_pass <- vapply(seq_len(nrow(hto_cellHash_log)),
function(x) {
hto_cellHash_log[x,] > hto_threshold[x]
},
logical(ncol(hto_cellHash_log)))
colnames(hto_cellHash_pass) <- rownames(hto_cellHash_log)
hto_cellHash_mix_label <- apply(hto_cellHash_pass, 1, function(x) {
if (sum(x) == 0) {
"negative"
} else if (sum(x) == 1) {
which(x)
} else {
"doublet/multiplet"
}
})
doubletClassify_between_class <- ifelse(!hto_cellHash_mix_label %in%
c("negative",
"doublet/multiplet"),
"Singlet", hto_cellHash_mix_label)
sce$doubletClassify_between_label <- hto_cellHash_mix_label
sce$doubletClassify_between_class <- doubletClassify_between_class
doublet_res <- list(doubletClassify_between_threshold = hto_threshold,
doubletClassify_between_resultsMat = hto_cellHash_pass)
S4Vectors::metadata(sce) <- c(S4Vectors::metadata(sce),
doublet_res)
return(sce)
}
#' plotHTO
#'
#' A function to plot HTO expression
#'
#' @param sce sce
#' @param which_idx which_idx
#' @param altExp_name altExp_name
#' @param ncol ncol
#'
#' @return A plot visualising the HTO expression
#'
#' @examples
#' data(CITEseq_example, package = "CiteFuse")
#' sce_citeseq <- preprocessing(CITEseq_example)
#' sce_citeseq <- normaliseExprs(sce = sce_citeseq,
#' altExp_name = "HTO",
#' transform = "log")
#' plotHTO(sce_citeseq, 1:4)
#'
#' @importFrom SummarizedExperiment assay assayNames
#' @importFrom SingleCellExperiment altExpNames altExp
#' @importFrom gridExtra grid.arrange
#' @importFrom cowplot axis_canvas insert_xaxis_grob insert_yaxis_grob
#' @importFrom grid unit
#' @importFrom utils combn
#' @import ggplot2
#'
#' @export
plotHTO <- function(sce,
which_idx = seq_len(2),
altExp_name = NULL,
ncol = 2) {
combination <- utils::combn(which_idx, 2)
ggList <- apply(combination, 2, function(x) {
plotHTOSingle(sce, which_idx = x, altExp_name = altExp_name)
})
do.call(gridExtra::grid.arrange,
c(ggList, ncol = min(length(ggList), ncol)))
}
#' plotHTOSingle
#'
#' A function to plot HTO expression
#'
#' @param sce sce
#' @param which_idx which_idx
#' @param altExp_name altExp_name
#'
#' @return A plot visualising the HTO expression
#'
#' @importFrom SummarizedExperiment assay assayNames
#' @importFrom SingleCellExperiment altExpNames altExp
#' @importFrom cowplot axis_canvas insert_xaxis_grob insert_yaxis_grob
#' @importFrom grid unit
plotHTOSingle <- function(sce,
which_idx = seq_len(2),
altExp_name = NULL
) {
if (is.null(altExp_name)) {
if (!"HTO" %in% altExpNames(sce)) {
stop("There is no HTO data in the object")
} else {
altExp_name <- "HTO"
}
}
if (!"logcounts" %in% assayNames(altExp(sce, altExp_name))) {
warning("HTO does not contain logcounts...
we will perform normaliseExprs() to get logcounts")
sce <- normaliseExprs(sce, altExp_name, "log")
}
noThreshold <- FALSE
if (!"doubletClassify_between_label" %in% colnames(colData(sce))) {
warning("Haven't performed doubletClassify() yet!")
noThreshold <- TRUE
}
hto_cellHash_log <- assay(SingleCellExperiment::altExp(sce, altExp_name),
"logcounts")
df <- data.frame(t(as.matrix(hto_cellHash_log[which_idx, ])))
colnames(df) <- lapply(strsplit(rownames(hto_cellHash_log[which_idx, ]),
"\\."), "[[", 1)
if (!noThreshold) {
hto_cellHash_pass <- metadata(sce)[["doubletClassify_between_resultsMat"]]
doublets <- hto_cellHash_pass[, which_idx[1]] &
hto_cellHash_pass[, which_idx[2]]
hto_threshold <- metadata(sce)[["doubletClassify_between_threshold"]]
} else {
doublets <- rep(FALSE, ncol(sce))
hto_threshold <- rep(NULL, 2)
hto_cellHash_pass <- NULL
}
pmain <- ggplot(df, aes(x = df[, 1], y = df[, 2], color = doublets)) +
geom_point(alpha = 0.5) +
geom_hline(yintercept = hto_threshold[which_idx[2]], col = "red",
linetype = 2, size = 1) +
geom_vline(xintercept = hto_threshold[which_idx[1]], col = "red",
linetype = 2, size = 1) +
scale_color_manual(values = c("#377EB8", "#E41A1C")) +
theme_bw() +
theme(aspect.ratio = 1) +
xlab(colnames(df)[1]) +
ylab(colnames(df)[2])
xdens <- cowplot::axis_canvas(pmain, axis = "x") +
geom_density(data = df, aes(x = df[, 1],
fill = hto_cellHash_pass[, which_idx[1]],
alpha = 0.5)) +
scale_fill_manual(values = c("#377EB8", "#E41A1C")) +
NULL
ydens <- cowplot::axis_canvas(pmain, axis = "y", coord_flip = TRUE) +
geom_density(data = df, aes(x = df[, 2],
fill = hto_cellHash_pass[, which_idx[2]],
alpha = 0.5)) +
coord_flip() +
scale_fill_manual(values = c("#377EB8", "#E41A1C")) +
NULL
p1 <- cowplot::insert_xaxis_grob(pmain, xdens,
grid::unit(.2, "null"), position = "top")
p2 <- cowplot::insert_yaxis_grob(p1, ydens,
grid::unit(.2, "null"), position = "right")
# ggdraw(p2)
return(p2)
}
#' withinSampleDoublets
#'
#' doublet identification within batch
#'
#' @param sce a SingleCellExperiment
#' @param altExp_name expression name of HTO matrix
#' @param eps eps of DBSCAN
#' @param minPts minPts of DBSCAN
#'
#'
#' @return A SingleCellExperiment object
#'
#' @examples
#'
#' data(CITEseq_example, package = "CiteFuse")
#' sce_citeseq <- preprocessing(CITEseq_example)
#' sce_citeseq <- normaliseExprs(sce = sce_citeseq,
#' altExp_name = "HTO",
#' transform = "log")
#' sce_citeseq <- crossSampleDoublets(sce_citeseq)
#' sce_citeseq <- withinSampleDoublets(sce_citeseq,
#' minPts = 10)
#'
#' @importFrom SummarizedExperiment assay assayNames
#' @importFrom SingleCellExperiment altExpNames altExp counts
#' @importFrom dbscan dbscan
#' @export
withinSampleDoublets <- function(sce,
altExp_name = NULL,
eps = 200,
minPts = 50) {
sce$nUMI <- Matrix::colSums(SingleCellExperiment::counts(sce))
if (is.null(altExp_name)) {
if (!"HTO" %in% SingleCellExperiment::altExpNames(sce)) {
stop("There is no HTO data in the object")
} else {
altExp_name <- "HTO"
}
}
if (!"logcounts" %in% assayNames(altExp(sce, altExp_name))) {
warning("HTO does not contain logcounts...
we will perform normaliseExprs() to get logcounts")
sce <- normaliseExprs(sce, altExp_name, "log")
}
if (!"doubletClassify_between_label" %in% colnames(colData(sce))) {
stop("Haven't performed doubletClassify_between() yet!")
}
hto_threshold <- metadata(sce)[["doubletClassify_between_threshold"]]
hto_cellHash_log <- assay(altExp(sce, altExp_name), "logcounts")
hto_cellHash_mix_label <- sce$doubletClassify_between_label
nUMI <- sce$nUMI
batch_doublets_list <- lapply(seq_len(nrow(hto_cellHash_log)), function(i) {
db_cluster <- dbscan::dbscan(cbind(nUMI, hto_cellHash_log[i,]),
eps = eps, minPts = minPts)
batch_doublets <- (db_cluster$cluster == 0 &
hto_cellHash_log[i,] > hto_threshold[i] &
hto_cellHash_mix_label != "doublet/multiplet")
})
batch_doublets_mat <- t(do.call(rbind, batch_doublets_list))
doubletClassify_within_label <- apply(batch_doublets_mat, 1, function(res) {
if (sum(res, na.rm = TRUE) == 0) {
"NotDoublets(Within)"
} else {
paste("Doublets(Within)", which(res), sep = "_")
}
})
doubletClassify_within_class <- ifelse(doubletClassify_within_label ==
"NotDoublets(Within)",
"Singlet", "Doublet")
sce$doubletClassify_within_label <- doubletClassify_within_label
sce$doubletClassify_within_class <- doubletClassify_within_class
S4Vectors::metadata(sce) <- c(S4Vectors::metadata(sce),
list(doubletClassify_within_resultsMat =
batch_doublets_mat))
return(sce)
}
#' @importFrom stats uniroot
#' @importFrom methods is
getThreshold <- function(mixmdl, verbose = FALSE){
# if (verbose) {
# plot(mixmdl, which = 2)
# }
membership <- apply(mixmdl$posterior, 1, which.max)
m_list <- sort(unique(membership))
mu_list <- mixmdl$mu
names(mu_list) <- seq_len(length(mu_list))
mu_list <- mu_list[m_list]
if (length(mu_list) > 1) {
idx1 <- as.numeric(names(mu_list)[order(mu_list)][1])
idx2 <- as.numeric(names(mu_list)[order(mu_list)][2])
root <- try(stats::uniroot(funMixModel,
interval = c(mixmdl$mu[idx1] -
mixmdl$sigma[idx1],
mixmdl$mu[idx2] +
mixmdl$sigma[idx2]),
mu1 = mixmdl$mu[idx1],
mu2 = mixmdl$mu[idx2],
sd1 = mixmdl$sigma[idx1],
sd2 = mixmdl$sigma[idx2],
rho1 = mixmdl$lambda[idx1],
rho2 = mixmdl$lambda[idx2]),
silent = TRUE)
if (!"try-error" %in% methods::is(root)) {
t <- root$root
}else{
t <- 0
}
}else{
t <- 0
}
return(t)
}
#' @importFrom stats dnorm
funMixModel <- function(x, mu1, mu2, sd1, sd2, rho1, rho2) {
stats::dnorm(x, mean = mu1, sd = sd1) * rho1 -
stats::dnorm(x, mean = mu2, sd = sd2) * rho2
}
