I am quite new with mixOmics. I have been working with DIABLO for a few days. I managed to run all the analysis, even represent some plots, and I obtained reasonable results.
However, this morning I updated the package (I was working with 6.22 version) and now the perf() function I was using to select ncomp it’s not working with my data.
I keep getting the error: Error: There are features with zero variance in block 'transcriptomics'. If nearZeroVar() function or 'near.zero.var' parameter hasn't been used, please use it. If you have used one of these, you may need to manually filter out these features
However, I verified if there were actually variables with zero variance and I get that ‘transcriptomics’ block has not any feature with zero variance:
I haven’t been able to reproduce your error unfortunately, and looking into the code I can’t find the point at which perf() looks for features with zero variance, this check should happen earlier when you create the DIABLO model using block.splsda().
I have created this short piece of code using the TGCA test data to show how features with 0 variance throw an error during block.splsda() but when I don’t have any features with 0 variance I don’t get any errors when running block.splsda() or running perf() on the created DIABLO model.
To see whether the error you’re getting is due to something specific to your data or something else, could you please 1) let me know which mixOmics version you are using (sharing the output of sessionInfo() would be ideal) and 2) could you please copy-paste the following code and try running it on your computer and let me know if you get any errors? Then we can take it from there!
## create a normal diablo model - runs without errors
data("breast.TCGA")
data = list(mrna = breast.TCGA$data.train$mrna, mirna = breast.TCGA$data.train$mirna,
protein = breast.TCGA$data.train$protein)
design = matrix(1, ncol = length(data), nrow = length(data),
dimnames = list(names(data), names(data)))
diag(design) = 0
list.keepX = list(mrna = rep(8,2), mirna = rep(8,2), protein = rep(8,2))
TCGA.block.splsda = block.splsda(X = data, Y = breast.TCGA$data.train$subtype,
ncomp = 2, keepX = list.keepX, design = design)
perf(TCGA.block.splsda) # runs without errors
# check that none of DIABLO blocks have features with 0 variance
nearZeroVar(TCGA.block.splsda$X$mrna)
nearZeroVar(TCGA.block.splsda$X$mirna)
nearZeroVar(TCGA.block.splsda$X$protein)
## edit the input data to introduce 0 variance to 20 of the mRNA variables
nearZeroVar(breast.TCGA$data.train$mrna)
zero_variance_mrna_data <- breast.TCGA$data.train$mrna
zero_variance_mrna_data[,1:20] <- 1
nearZeroVar(zero_variance_mrna_data)
# create a new DIABLO model with this new data - `block.splsda() fails`
data = list(mrna = zero_variance_mrna_data,
mirna = breast.TCGA$data.train$mirna,
protein = breast.TCGA$data.train$protein)
TCGA.block.splsda = block.splsda(X = data, Y = breast.TCGA$data.train$subtype,
ncomp = 2, keepX = list.keepX, design = design)
The code you provided works as expected: the first part doesn’t return any error, but in the segon DIABLO model with the data with 0 variance variables returns the error:
Error: There are features with zero variance in block 'mrna'. If nearZeroVar() function or 'near.zero.var' parameter hasn't been used, please use it. If you have used one of these, you may need to manually filter out these features.
So I guess the problem is specific to my data
I don’t know where can it be, as I filtered the variables with 0 variance before creating the model. And I don’t understant why the error occurs in the perf() function but not in the block.splsda().
I share my code so you can see if I’m doing something wrong:
data <- read.table("PA_dataframe_log.txt", header=TRUE, sep="\t", row.names=1)
# Delete columns with NA
columns_NA <- colSums(is.na(data)) > 0
index_columns_NA <- which(columns_NA)
data <- data[, -index_columns_NA]
# Delete variables with 0 variance
nzv <- nearZeroVar(data)
data <- data[, -nzv$Position, drop=FALSE]
# Subset data to create the different blocks
physiological <- data[, (1:3)]
transcriptomics <- data[, grepl("SOLYC", colnames(data))]
metabolomics <- data[, grepl("met", colnames(data))]
# Create a list with the 3 blocks
data = list(physiological = physiological,
transcriptomics = transcriptomics,
metabolomics = metabolomics)
# Define Y
Y <- as.factor(groups)
summary(Y)
LUK_C_PA LUK_S_PA NOT_C_PA NOT_S_PA
3 3 3 3
# Create the design matrix
design = matrix(0.8, ncol = length(data), nrow = length(data),
+ dimnames = list(names(data), names(data)))
diag(design) = 0
design
physiological transcriptomics metabolomics
physiological 0.0 0.8 0.8
transcriptomics 0.8 0.0 0.8
metabolomics 0.8 0.8 0.0
basic.diablo.model = block.splsda(X = data, Y = Y, ncomp = 3, design = design)
Design matrix has changed to include Y; each block will be
linked to Y.
perf.diablo <- perf(basic.diablo.model, validation = "loo")
Error: There are features with zero variance in block 'transcriptomics'. If nearZeroVar() function or 'near.zero.var' parameter hasn't been used, please use it. If you have used one of these, you may need to manually filter out these features.
Maybe the problem is that I use validation = 'lool' in the perf() function? Or something with my data?
I have been trying different things and filtering my data in different way but nothing worked. But just now I realized that adding the parameter near.zero.var = TRUE to the block.splsda() function, the perf() function works well and without any error.
New version: basic.diablo.model = block.splsda(X, Y, ncomp = 3, design = design, near.zero.var = TRUE)
Maybe the solution was as simple as that, but I didn’t realize before. It may be useful for other people having the same problem.
Thank you very much and congrats for this incredible bioinformatics tool.
That is interesting and thank you for posting the solution you found! I will look into this a bit more in the future to see whats going on, in the meantime your workaround seems like the best solution. Feel free to add another post if you encounter any other issues!
