Skip to content

Error if add gene-level covariates #34

@rtian30

Description

@rtian30

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!

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions