Question

bias.data vector must have the same length as DEgenes vector!

0

Entering edit mode

mictadlo ▴ 10

@mictadlo-10885

Last seen 4.3 years ago

Hi, I followed this paper until Pathway analysis because I wanted to replace it by goseq.

> head(gene_len)
         V1   V2        V3
1    Length   GC        NA
2 sp0000001 6406 0.3535748
3 sp0000002  333 0.5495495
4 sp0000003  216 0.4583333
5 sp0000004  357 0.5602241
6 sp0000005 4808 0.3928869
> head(go_mapping)
         V1         V2
1    geneID  Gomapping
2 sp0000001        n/a
3 sp0000005 GO:0003723
4 sp0000006 GO:0016021
5 sp0000007 GO:0003700
6 sp0000007 GO:0006355

Unfortunately, nullp caused the following error:

> pwf <- nullp(tr, "sp", id=go_mapping, bias.data=gene_len)
Error in nullp(tr, "sp", id = go_mapping, bias.data = gene_len) : 
  bias.data vector must have the same length as DEgenes vector!

How is it possible to fix it?

Thank you in advance.

edger goseq • 1.6k views

ADD COMMENT • link updated 6.6 years ago by Aaron Lun ★ 28k • written 6.6 years ago by mictadlo ▴ 10

score 0 · Answer 1 · 2017-09-15

0

Entering edit mode

Aaron Lun ★ 28k

@alun

Last seen 8 hours ago

The city by the bay

There are a number of errors with your code:

Is "Length" really a gene name? I think you missed a header=TRUE when loading these tables.
Missing values in R use NA, not "n/a".
nullp expects a vector for bias.data, not a data.frame. See ?nullp. (It also expects a string for id, but this seems to be ignored if bias.data is supplied.)

ADD COMMENT • link 6.6 years ago Aaron Lun ★ 28k

0

Entering edit mode

Hi, I took your suggested on board but I have a problem with nullp's DEgenes Argument (A named binary vector where 1 represents DE, 0 not DE and the names are gene IDs). In this paper they used the following to determine differential expression

> tr <- glmTreat(fit, contrast=B.LvsP, lfc=log2(1.5))
> topTags(tr)
Coefficient:  -1*Leaves.2 1*Leaves.3
          logFC unshrunk.logFC logCPM  PValue     FDR
sp0090975   2.6            2.6    2.8 5.7e-14 2.7e-09
sp0037632  -3.0           -3.0    3.4 1.7e-13 2.7e-09
sp0074153  -3.9           -3.9    3.8 1.8e-13 2.7e-09
sp0008306   3.2            3.2    2.9 1.8e-13 2.7e-09
sp0073530  -4.5           -4.5    3.4 2.3e-12 2.3e-08
sp0025713  -3.9           -3.9    4.4 2.6e-12 2.3e-08
sp0037721   7.8            8.0    2.4 2.7e-12 2.3e-08
sp0083660   2.0            2.0    4.4 3.2e-12 2.3e-08
sp0052245  -2.9           -2.9    4.9 3.9e-12 2.6e-08
sp0071520  -3.3           -3.3    2.7 5.8e-12 3.4e-08

> library(goseq)
> gene_len.path <- paste0("~/projects/GOseq/", "GC_lengths.tsv")
> gene_len <- read.delim(gene_len.path, header=T)
> head(gene_len)
          Length   GC
sp0000001   6406 0.35
sp0000002    333 0.55
sp0000003    216 0.46
sp0000004    357 0.56
sp0000005   4808 0.39
sp0000006   3558 0.41
> go_mapping.path <- paste0("~/projects/GOseq/", "GOmapping.tsv")
> go_mapping <- read.delim(go_mapping.path, header=T)
> head(go_mapping)
     geneID  Gomapping
1 sp0000001       <NA>
2 sp0000005 GO:0003723
3 sp0000006 GO:0016021
4 sp0000007 GO:0003700
5 sp0000007 GO:0006355
6 sp0000007 GO:0043565
> 

> gene.length<-as.vector(t(gene_len$Length))

> topTags(tr)
Coefficient:  -1*Leaves.2 1*Leaves.3 
          logFC unshrunk.logFC logCPM  PValue     FDR
sp0090975   2.6            2.6    2.8 5.7e-14 2.7e-09
sp0037632  -3.0           -3.0    3.4 1.7e-13 2.7e-09
sp0074153  -3.9           -3.9    3.8 1.8e-13 2.7e-09
sp0008306   3.2            3.2    2.9 1.8e-13 2.7e-09
sp0073530  -4.5           -4.5    3.4 2.3e-12 2.3e-08
sp0025713  -3.9           -3.9    4.4 2.6e-12 2.3e-08
sp0037721   7.8            8.0    2.4 2.7e-12 2.3e-08
sp0083660   2.0            2.0    4.4 3.2e-12 2.3e-08
sp0052245  -2.9           -2.9    4.9 3.9e-12 2.6e-08
sp0071520  -3.3           -3.3    2.7 5.8e-12 3.4e-08

> pwf <- nullp(tr, "sp", id=go_mapping, bias.data=gene.length)
Error in nullp(tr, "sp", id = go_mapping, bias.data = gene.length) : 
  bias.data vector must have the same length as DEgenes vector!

How do I convert tr to DEgenes?

Thank you in advance.

ADD REPLY • link 6.6 years ago mictadlo ▴ 10