This is a test version of Biostars. For the public version, visit https://www.biostars.org.
goseq with non native organism

I´m using goseq package after DESeq2 analysis and I am totally new in R

When I would like to test the pwf of my DE genes, I get an error:

if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")

BiocManager::install("TxDb.Athaliana.BioMart.plantsmart28")

library(TxDb.Athaliana.BioMart.plantsmart28)
txdb<-TxDb.Athaliana.BioMart.plantsmart28
txsByGene=transcriptsBy(txdb,"gene")
lengthData=median(width(txsByGene))



assayed.genes <- rownames(res)
de.genes <- rownames(res)[ which(res$padj < 0.05)]
genes<-as.integer(assayed.genes%in%de.genes)
names(genes)<-row.names(res$padj)

genes
    0     1 
29357  3476 

pwf<-nullp(genes,txdb,assayed.genes,lengthData)

Error in nullp(genes, txdb, assayed.genes, lengthData) : 
bias.data vector must have the same length as DEgenes vector!

length(lengthData)=33602

length(genes)=32833

DEgenes vector in the script is genes and contains all Arabidopsis genes (32833). I did not filter for counts.

bias.data vector in the script is lengthData and contains median length values of all transcripts (33602).

Does someone knows how to solve this error? Thank you very much in advance!

rna-seq goseq arath
• 2,006 views • link
updated by ATpoint ATpoint • written by alyyha • alyyha

1 answer

You have 32833 genes and 33602 elements in the bias vector, see the problem? You must have one value in the bias bector and it must have the exact same order as the gene list.

Code suggestion:

assayed.genes <- rownames(res)
genes         <- as.integer(res$padj < 0.05)
biasdata      <- lengthData[match(assayed.genes, names(txsByGene))]
nP            <- nullp(DEgenes = genes, bias.data = biasdata)
• 0 views • link
written by ATpoint ATpoint

Yes I see this problem, thank you. I am just wondering how to fix it because I don´t know how to do it.

• 0 views • link
written by alyyha • alyyha

Just added a suggestion.

• 0 views • link
written by ATpoint ATpoint

Thank you very much! I'll try

• 0 views • link
written by alyyha • alyyha

Also added a command for nullp. The DEgenes is a vector of same length as biasdata and it has a 1 if the gene is DE and a 0 if not.

• 0 views • link
written by ATpoint ATpoint

It worked. Danke schon!!

• 0 views • link
written by alyyha • alyyha

Glad to help! If the answer solved the issue please accept it to provide closure (using the checkmark button).

• 0 views • link
written by ATpoint ATpoint

Log in to answer this question.

More posts like this