Hi all,
I am currently working on a time course analysis of RNA-seq zebrafish data, looking at gene expression at 3 time points across two strains. I want to use BETR as a module for this analysis, however I ran into some errors and I can't seem to find the answer. After I ran the betr function, it returned the following error:
``Error in if (evar > 0) { : missing value where TRUE/FALSE needed``
``In addition: Warning message:``
``Zero sample variances detected, have been offset``
This seems to reflect the same problem I have been having when running BETR on the open source MeV program where it tells me that I have a lot of invalid genes. Looking for answers online, the problem seems to be related to either the data contains a lot of NaNs or that the data has "zero variance". However, the expression set that I am inputting into is filtered so that it has a value in at least one condition, and all 0s were converted to 1 for the log transformations. I also did a variance filter on the data (taking the top 50%) just to see if betr would run but unfortunately the problem persists. I was wondering what are the requirements for the input data into betr? Should the data input be fold-change compared to the time 0 control? Any assistance is greatly appreciated.
Here is the code I have (for generating the ExpressionSet input into BETR as well)
---
> library(betr)
> library(Biobase)
> library(zebrafish.db)
> exprs <- as.matrix(Wateronly)
> class(exprs)
[1] "matrix"
> dim(exprs)
[1] 14969 18
> head(exprs)
> pData <- pdatabetr
> dim(pData)
[1] 18 3
> all(rownames(pData)==colnames(exprs))
[1] TRUE
> metadata <- data.frame(labelDescription = c("Strain of fish",
+ "Time point",
+ "Replicate"),
+ row.names = c("Strain", "Time", "Replicate"))
> phenoData <- new("AnnotatedDataFrame", data = pData, varMetadata = metadata)
> quang <- ExpressionSet(assayData = exprs, phenoData = phenoData, annotation = "zebrafish.db")
>
> pData(quang)
Strain Time Replicate
X2ABwater1 WT 2 1
X2ABwater2 WT 2 2
X2ABwater3 WT 2 3
X2KOwater1 KO 2 1
X2KOwater2 KO 2 2
X2KOwater3 KO 2 3
X48ABwater1 WT 48 1
X48ABwater2 WT 48 2
X48ABwater3 WT 48 3
X48KOwater1 KO 48 1
X48KOwater2 KO 48 2
X48KOwater3 KO 48 3
X96ABwater1 WT 96 1
X96ABwater2 WT 96 2
X96ABwater3 WT 96 3
X96KOwater1 KO 96 1
X96KOwater2 KO 96 2
X96KOwater3 KO 96 3
>
> prob <- betr(eset=quang, cond=pData(quang)$Strain,
+ timepoint=pData(quang)$Time, replicate=pData(quang)$Replicate, alpha=0.05)