-
Notifications
You must be signed in to change notification settings - Fork 10
Description
Hi,
I added V matrix (gene-level covariates) in ZinbFit, but got an error message:
Error in validObject(.Object) :
invalid class “ZinbModel” object: beta_mu must have J columns!
This is my code:
data<-read.csv(file ="./new_Expression_matrix.csv", header=T, row.names = 1, sep = "")
data<-as.data.frame(t(data)) #row names are genes
construct SummartizedExperiment class
ExperimentalDesign<-read.csv(file="./Experimental_Design.csv",header=T,row.names=1,sep=",")
UMI<- as.matrix(data[-1, ]) #small subset, change later
colData <- DataFrame(ExperimentalDesign,row.names=colnames(UMI)) #experimental design
scRNA<-SummarizedExperiment(assays=list(counts=UMI), colData=colData) #the function uses assay named "counts"
identify the 2000 most variable genes, and keep in scRNA, you could also increase the number of variable genes kept
assay(scRNA) %>% log1p %>% rowVars -> vars
names(vars) <- rownames(scRNA)
vars <- sort(vars, decreasing = TRUE)
head(vars)
scRNA <- scRNA[names(vars)[1:1000],]
create gene-level covariates
mart <- useMart("ensembl")
mart <- useDataset("hsapiens_gene_ensembl", mart = mart) #get access to ensembl,hsapiens dataset
bm <- getBM(attributes=c('hgnc_symbol', 'start_position',
'end_position', 'percentage_gene_gc_content'),
filters = 'hgnc_symbol',
values = rownames(scRNA),
mart = mart)
bm$length <- bm$end_position - bm$start_position
len <- tapply(bm$length, bm$hgnc_symbol, mean)
len <- len[rownames(scRNA)]
gcc <- tapply(bm$percentage_gene_gc_content, bm$hgnc_symbol, mean)
gcc <- gcc[rownames(scRNA)]
rowData(scRNA) <- data.frame(gccontent = gcc, length = len)
Zinbfit, fit a zinb regression model
fit<-zinbFit(Y=scRNA,
X="~Run",
V="~gccontent+ log(length)",
verbose=T, K = 2, epsilon=1000)
It is very similar as vignette.
And I also run the code from vignette, the same error message happened.
Can you check if there is a problem?
Thanks!