I have a problem when I use VIPER with aracne networks.
I used regulonbrca in aracne to calculate protein activity, The number of genes in regulonbrca is 6054.
And the number of genes in regulonbrca also in my TPM file is 5918.
But I got only 4506 genes as a viper result.
Why some genes are missing when running viper? Is there a default setting about it?
Here is my code.
expdat= "TPM-brca(entrez_avg).csv"
clsdat= "immunesubtype-BRCA.csv"
regulon= regulonbrca
dat0 = data.frame(read.csv(expdat))
dat1 = dat0[,-c(1,2)]
rownames(dat1) = dat0[,1]
cls0 = data.frame(read.csv(clsdat))
colnames(cls0) = c("id", "description")
cls1 = cls0 %>% dplyr::filter(id %in% colnames(dat1))
dat = dat1[,cls1$id]
cls = data.frame(description=cls1[,2])
rownames(cls) = cls1[,1]
meta = data.frame(labelDescription=c("description"), row.names = colnames(cls))
pheno = new("AnnotatedDataFrame", data=cls, varMetadata=meta)
dset = ExpressionSet(assayData = as.matrix(dat), phenoData = pheno)
signature = rowTtest(dset, "description", "TRUE", "FALSE")
signature = (qnorm(signature$p.value/2, lower.tail = FALSE) * + sign(signature$statistic))[, 1]
nullmodel = ttestNull(dset, "description", "TRUE", "FALSE", per = 1000, repos = TRUE, verbose = FALSE)
vpres = viper(dset, regulon, verbose = FALSE) ## single sample viper
res_ss0 = data.frame(id=row.names(vpres@assayData$exprs), vpres@assayData$exprs)
geneid_ss0 = ldply(sapply(row.names(res_ss0), converter), data.frame) ## gene mapping
colnames(geneid_ss0) = c('id', 'gene')
res_ss = merge(geneid_ss0, res_ss0, key='id')[, -1]