Error in running Combat function
1
0
Entering edit mode
Reza ▴ 10
@Reza-24551
Last seen 3.8 years ago

I run the same commands on the same data in Rstudio in Windows and Linux, but I get different results in step removing batch effect using “Combat” function (sva package). Running “Combat” on Windows is error-free, while on Linux I get the following error:

Error in ((dat - t(design %*% B.hat))^2) %*% rep(1/n.array, n.array) : 
requires numeric/complex matrix/vector arguments

Commands

    all <- read.delim("All.txt", header = T, row.names = 1) 
    allm <- all - rowMeans(all)
    pc <- prcomp(allm)
    Sets <- (c(rep("GSEx", 12), rep("GSExx", 12), rep("GSExxx", 9), rep("GSExxxx", 13)))
    pcr <- data.frame(pc$r[,1:3], Sets)
    allc <- ComBat(all, Sets)

For reproducibility

> dput(head(all))
structure(list(A1 = c(6.9381, 7.5079, 7.2303, 8.9908, 8.2564, 
11.0564), A2 = c(6.9594, 7.3883, 7.1023, 8.9328, 8.1679, 10.8513
), A3 = c(6.8419, 7.2109, 6.9548, 8.7414, 8.0299, 10.8115), A4 = c(6.6619, 
7.1364, 6.887, 8.6863, 8.006, 10.8775), A5 = c(6.9296, 7.2845, 
7.0382, 8.8176, 8.1234, 10.9084), A6 = c(7.1549, 7.5539, 7.2842, 
9.0999, 8.3527, 11.1642), A7 = c(6.6192, 7.1043, 6.9482, 8.6802, 
7.9977, 10.7216), A8 = c(7.2778, 7.6654, 7.3542, 9.2437, 8.5063, 
11.2892), A9 = c(7.282, 7.5827, 7.3097, 9.116, 8.4326, 11.1977
), A10 = c(7.3402, 7.7757, 7.4926, 9.2278, 8.6172, 11.1968), 
    A11 = c(6.6635, 7.1918, 6.8924, 8.6318, 7.9105, 10.7312), 
    A12 = c(6.8926, 7.3513, 6.9834, 8.7605, 8.0945, 10.8734), 
    B1 = c(6.8014, 7.1908, 6.8015, 8.8403, 8.0259, 10.7076), 
    B2 = c(6.5855, 7.1594, 6.6593, 8.7153, 7.9119, 10.7066), 
    B3 = c(6.4392, 7.0688, 6.5803, 8.6158, 7.8153, 10.6042), 
    B4 = c(6.7046, 7.078, 6.7047, 8.6752, 7.8961, 10.6192), B5 = c(7.9128, 
    8.3322, 8.0958, 9.8519, 8.9967, 11.6034), B6 = c(6.9226, 
    7.4705, 7.1007, 8.9497, 8.2543, 10.711), B7 = c(6.421, 7.0791, 
    6.4662, 8.5758, 7.7724, 10.5219), B8 = c(6.4972, 7.1769, 
    6.7016, 8.7424, 7.9157, 10.7647), B9 = c(6.6425, 7.3366, 
    6.5679, 8.8394, 7.9696, 10.7032), B10 = c(6.3134, 6.7753, 
    6.4695, 8.4463, 7.644, 10.5063), B11 = c(6.4938, 7.0694, 
    6.336, 8.6302, 7.688, 10.5626), B12 = c(6.2451, 6.9604, 6.2342, 
    8.4628, 7.6297, 10.5476), C1 = c(4.9469, 5.1746, 5.1522, 
    9.5331, 8.9385, 11.4915), C2 = c(5.0953, 5.1003, 5.2505, 
    9.6333, 9.0602, 11.4867), C3 = c(5.0583, 5.0543, 5.1091, 
    9.6526, 9.0479, 11.5055), C4 = c(4.9402, 5.1395, 5.2561, 
    9.655, 9.0221, 11.5286), C5 = c(5.0406, 5.1575, 5.2761, 9.6413, 
    9.0673, 11.5767), C6 = c(4.9726, 5.3945, 5.1122, 9.5541, 
    9.0265, 11.5038), C7 = c(5.3879, 5.1028, 5.0608, 9.803, 9.2769, 
    11.8283), C8 = c(5.127, 5.142, 5.1496, 9.9158, 9.2543, 11.8024
    ), C9 = c(5.3067, 5.1215, 5.0303, 9.9562, 9.3744, 11.8254
    ), D1 = c(7.0772, 7.4086, 7.2877, 9.0927, 8.6348, 11.2268
    ), D2 = c(6.6349, 7.0915, 6.6813, 8.5535, 8.1593, 10.9012
    ), D3 = c(7.4336, 7.7217, 7.512, 9.3543, 8.8803, 11.5023), 
    D4 = c(7.4971, 7.8568, 7.672, 9.441, 9.0128, 11.6245), D5 = c(7.2505, 
    7.5549, 7.3561, 9.2116, 8.6796, 11.3779), D6 = c(6.8906, 
    7.2358, 7.054, 8.8379, 8.4158, 11.0953), D7 = c(6.6697, 7.0378, 
    6.9185, 8.7005, 8.1895, 11.0119), D8 = c(6.955, 7.2885, 7.1322, 
    8.922, 8.4309, 11.1848), D9 = c(7.4286, 7.628, 7.5056, 9.2589, 
    8.8542, 11.5108), D10 = c(7.1027, 7.3881, 7.2015, 9.0158, 
    8.5069, 11.2525), D11 = c(6.9174, 7.2366, 7.1286, 8.8789, 
    8.5092, 11.1913), D12 = c(7.3573, 7.6144, 7.4158, 9.231, 
    8.7868, 11.4215), D13 = c(7.3078, 7.5243, 7.4756, 9.2333, 
    8.7521, 11.4136)), row.names = c("AFFX-BioB-3_at", "AFFX-BioB-5_at", 
"AFFX-BioB-M_at", "AFFX-BioC-3_at", "AFFX-BioC-5_at", "AFFX-BioDn-3_at"
), class = "data.frame")

how can solve this problem?

Thanks for your helps

sva ComBat • 1.7k views
ADD COMMENT
0
Entering edit mode
Kevin Blighe ★ 4.0k
@kevin
Last seen 22 days ago
Republic of Ireland

Thanks for posting the reproducible data. Using this data, the command runs okay here on Linux:

allc <- sva::ComBat(all, Sets)

Found4batches
Adjusting for0covariate(s) or covariate level(s)
Standardizing Data across genes
Fitting L/S model and finding priors
Finding parametric adjustments
Adjusting the Data

Perhaps you can show the output of your sessionInfo()? - here is mine:

sessionInfo()
R version 4.0.3 (2020-10-10)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.7 LTS

Matrix products: default
BLAS:   /usr/lib/atlas-base/atlas/libblas.so.3.0
LAPACK: /usr/lib/atlas-base/atlas/liblapack.so.3.0

locale:
 [1] LC_CTYPE=pt_BR.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=pt_BR.UTF-8    
 [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=pt_BR.UTF-8   
 [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.4.6         AnnotationDbi_1.50.0 edgeR_3.30.3        
 [4] BiocGenerics_0.34.0  splines_4.0.3        IRanges_2.22.2      
 [7] bit_1.1-15.2         BiocParallel_1.22.0  lattice_0.20-41     
[10] xtable_1.8-4         rlang_0.4.6          blob_1.2.1          
[13] grid_4.0.3           parallel_4.0.3       nlme_3.1-149        
[16] Biobase_2.48.0       mgcv_1.8-33          DBI_1.1.0           
[19] genefilter_1.70.0    matrixStats_0.56.0   survival_3.2-3      
[22] bit64_0.9-7          digest_0.6.25        Matrix_1.2-18       
[25] sva_3.36.0           vctrs_0.3.1          S4Vectors_0.26.1    
[28] bitops_1.0-6         RCurl_1.98-1.2       memoise_1.1.0       
[31] RSQLite_2.2.0        limma_3.44.3         compiler_4.0.3      
[34] locfit_1.5-9.4       stats4_4.0.3         XML_3.99-0.3        
[37] annotate_1.66.0
ADD COMMENT
1
Entering edit mode

Thank you very much for your answer and your effort to help me, I will post here the results of sessionInfo() at the first opportunity, I will not have access to the university server for two days now. Thank you again.

ADD REPLY
0
Entering edit mode

I am very sorry for the delay in the answer and I hope to receive the necessary help from you to solve my problem.

The output of sessionInfo() on my Linux

> sessionInfo()
R version 3.6.3 (2020-02-29)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.1 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods  
[7] base     

loaded via a namespace (and not attached):
[1] compiler_3.6.3 tools_3.6.3
ADD REPLY
0
Entering edit mode

According to your sessionInfo() output, you have not loaded the sva package (?)

To load, try:

library(sva)
ADD REPLY
0
Entering edit mode
> library(sva)
Loading required package: mgcv
Loading required package: nlme
This is mgcv 1.8-33. For overview type 'help("mgcv-package")'.
Loading required package: genefilter
Loading required package: BiocParallel
> allc <- ComBat(all, Sets)
Found10batches
Adjusting for0covariate(s) or covariate level(s)
Standardizing Data across genes
Error in ((dat - t(design %*% B.hat))^2) %*% rep(1/n.array, n.array) : 
  requires numeric/complex matrix/vector arguments
> 
ADD REPLY
0
Entering edit mode

Thank you. Can you please post the output of:

str(all)
str(Sets)
ADD REPLY
0
Entering edit mode

Thank you for your efforts to help me.

> str(all)
'data.frame':   61290 obs. of  46 variables:
 $ A1 : num  6.94 7.51 7.23 8.99 8.26 ...
 $ A2 : num  6.96 7.39 7.1 8.93 8.17 ...
 $ A3 : num  6.84 7.21 6.95 8.74 8.03 ...
 $ A4 : num  6.66 7.14 6.89 8.69 8.01 ...
 $ A5 : num  6.93 7.28 7.04 8.82 8.12 ...
 $ A6 : num  7.15 7.55 7.28 9.1 8.35 ...
 $ A7 : num  6.62 7.1 6.95 8.68 8 ...
 $ A8 : num  7.28 7.67 7.35 9.24 8.51 ...
 $ A9 : num  7.28 7.58 7.31 9.12 8.43 ...
 $ A10: num  7.34 7.78 7.49 9.23 8.62 ...
 $ A11: num  6.66 7.19 6.89 8.63 7.91 ...
 $ A12: num  6.89 7.35 6.98 8.76 8.09 ...
 $ B1 : num  6.8 7.19 6.8 8.84 8.03 ...
 $ B2 : num  6.59 7.16 6.66 8.72 7.91 ...
 $ B3 : num  6.44 7.07 6.58 8.62 7.82 ...
 $ B4 : num  6.7 7.08 6.7 8.68 7.9 ...
 $ B5 : num  7.91 8.33 8.1 9.85 9 ...
 $ B6 : num  6.92 7.47 7.1 8.95 8.25 ...
 $ B7 : num  6.42 7.08 6.47 8.58 7.77 ...
 $ B8 : num  6.5 7.18 6.7 8.74 7.92 ...
 $ B9 : num  6.64 7.34 6.57 8.84 7.97 ...
 $ B10: num  6.31 6.78 6.47 8.45 7.64 ...
 $ B11: num  6.49 7.07 6.34 8.63 7.69 ...
 $ B12: num  6.25 6.96 6.23 8.46 7.63 ...
 $ C1 : num  4.95 5.17 5.15 9.53 8.94 ...
 $ C2 : num  5.1 5.1 5.25 9.63 9.06 ...
 $ C3 : num  5.06 5.05 5.11 9.65 9.05 ...
 $ C4 : num  4.94 5.14 5.26 9.65 9.02 ...
 $ C5 : num  5.04 5.16 5.28 9.64 9.07 ...
 $ C6 : num  4.97 5.39 5.11 9.55 9.03 ...
 $ C7 : num  5.39 5.1 5.06 9.8 9.28 ...
 $ C8 : num  5.13 5.14 5.15 9.92 9.25 ...
 $ C9 : num  5.31 5.12 5.03 9.96 9.37 ...
 $ D1 : num  7.08 7.41 7.29 9.09 8.63 ...
 $ D2 : num  6.63 7.09 6.68 8.55 8.16 ...
 $ D3 : num  7.43 7.72 7.51 9.35 8.88 ...
 $ D4 : num  7.5 7.86 7.67 9.44 9.01 ...
 $ D5 : num  7.25 7.55 7.36 9.21 8.68 ...
 $ D6 : num  6.89 7.24 7.05 8.84 8.42 ...
 $ D7 : num  6.67 7.04 6.92 8.7 8.19 ...
 $ D8 : num  6.96 7.29 7.13 8.92 8.43 ...
 $ D9 : num  7.43 7.63 7.51 9.26 8.85 ...
 $ D10: num  7.1 7.39 7.2 9.02 8.51 ...
 $ D11: num  6.92 7.24 7.13 8.88 8.51 ...
 $ D12: num  7.36 7.61 7.42 9.23 8.79 ...
 $ D13: num  7.31 7.52 7.48 9.23 8.75 ...
 > str(Sets)
  chr [1:46] "GSEx" "GSEx" "GSEx" "GSEx" "GSEx" "GSEx" "GSEx" ...
ADD REPLY
0
Entering edit mode

Nothing seems out of the ordinary. Please check for things like missing values, infinite values, et cetera

ADD REPLY

Login before adding your answer.

Traffic: 908 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