ComBat_seq function in sva package does not return a matrix
1
1
Entering edit mode
Federico ▴ 10
@00176341
Last seen 14 months ago
Italy

Hello everyone, I'm performing a RNA-seq analysis to understand the difference between 2 condition. I know I have data with batch effect and I know the batch "groups" (some samples were collected one day and the others another day). I know when the different samples where collected and I wanted to use the ComBat_seq function of sva to correct my raw reads matrix but the function returns me a large list and not a count matrix. Am I missing something? Here is the code I used, just to be sure I made everything right.

x <- read.table("raw_gene_counts.tsv",header=T,row.names=1)
head(x)

                OH86 OH92 OH98 OI04 OH87 OH93 OH99 OI05
ENSG00000000003   90  100   94  129   84  115   33  135
ENSG00000000005    0    0    0    1    0    0    0    0
ENSG00000000419  293  225  275  477  319  461  126  530
ENSG00000000457   79   64   53  107   73   99   23  135
ENSG00000000460   43   56   34   68   48   61   28   55
ENSG00000000938   13    5    8   14    9    6    2   13

batch <- as.factor(c("a", "a"; "b", "b"....))
x <- ComBat_seq(x, batch=batch, group=NULL, full_mod=FALSE)
head(x)

[[1]]
[1] 95

[[2]]
   [1] 0 0 0 0 0 0 1 0 0 0 0 0 2 2 0 1 1 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0
  [66] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 2 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
 [131] 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [196] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0
 [261] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0
 [326] 0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 1 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 2
 [391] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
 [456] 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 2 0 0 0 0 1 0 0 0 1 0 0 0 0 0
 [521] 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0
 [586] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 1
 [651] 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [716] 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 2 0 0 0 0 1 0 0 0 0
 [781] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [846] 0 0 0 1 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [911] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
 [976] 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 2 1 0 0 0 0
 [ reached getOption("max.print") -- omitted 33030 entries ]

[[3]]
[1] 318

[[4]]
[1] 85

[[5]]
[1] 42

[[6]]
[1] 14

I also wanted to understand how the parameter "group" influences the outcome. Thanks everyone for the help

sva ComBat ComBat-seq Bioconductor ComBat_seq • 1.7k views
ADD COMMENT
4
Entering edit mode
ATpoint ★ 4.6k
@atpoint-13662
Last seen 23 hours ago
Germany

That happens if your x is a data.frame rather than a matrix. It is listed under the Q&A section in https://github.com/zhangyuqing/ComBat-seq starting with The function returns a list rather than the adjusted matrix.. It's questionable though why the developer do not throw a warning or x is a data.frame though, or just do the conversion to matrix internally.

ADD COMMENT
0
Entering edit mode

Thanks for your help, problem solved :)

ADD REPLY

Login before adding your answer.

Traffic: 845 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6