I’m working on the identification of the genes affecting the efficacy of an anti-cancer drug, which targets specific protein expressed on the surface of the tumor cells. I’m trying to identify such genes because the efficacy cannot be fully correlated or explained by the expression level of the target gene. I have RNA-seq data and efficacy data taken from several samples. I show you a part of sample annotations bellow.

> head(sample.annotation, 20)
     cellLine    tumorType efficacy targetExpr bioRep
1  cellLine_A       Breast     40.9      211.0      1
2  cellLine_A       Breast     40.9      211.0      2
3  cellLine_B       Breast     33.9      107.0      1
4  cellLine_B       Breast     33.9      107.0      2
5  cellLine_C  Endometrium     91.1      280.6      1
6  cellLine_C  Endometrium     91.1      280.6      2
7  cellLine_D  Endometrium     55.7       97.6      1
8  cellLine_D  Endometrium     55.7       97.6      2
9  cellLine_E  Endometrium     98.1      151.4      1
10 cellLine_E  Endometrium     98.1      151.4      2
11 cellLine_E  Endometrium     98.1      151.4      3
12 cellLine_F   Esophageal     24.3       35.6      1
13 cellLine_F   Esophageal     24.3       35.6      2
14 cellLine_F   Esophageal     24.3       35.6      3
15 cellLine_G   Esophageal     58.8      128.2      1
16 cellLine_G   Esophageal     58.8      128.2      2
17 cellLine_G   Esophageal     58.8      128.2      3
18 cellLine_H Glioblastoma     70.5      200.4      1
19 cellLine_H Glioblastoma     70.5      200.4      2
20 cellLine_I     Melanoma     58.7      255.3      1

I have 2-3 biological replicates (NOT technical replicates) of RNA-seq data for each cell line. What I tried to do with limma is

design <- model.matrix(~ 0 + efficacy + targetExpr + tumorType, data = sample.annotation)
v <- voom(dgeData, design)
corfit <- duplicateCorrelation(v, design, block = sample.annotation$cellLine)
fit <- lmFit(v, design, block = sample.annotation$cellLine, correlation = corfit$consensus.correlation)
fit2 <- eBayes(fit)

But I realized that some of tumor types have a few cell lines so that it could be not appropriate to set tumorType as a covariate. Then I thought it would be better to treat tumor types as random effect as well as cell lines.

> table(sample.annotation$tumorType, sample.annotation$bioRep)
                1  2  3
  Breast        2  2  0
  Endometrium   3  3  1
  Esophageal    2  2  2
  Glioblastoma  1  1  0
  Melanoma      1  1  0
  Mesothelioma  1  1  0
  NSCLC        11 11  9
  Ovarian       1  1  0
  Prostate      1  1  0
  SCLC          3  3  3

My question is

1. How can I take two items (cell lines and tumor types in this case) into consideration in duplicateCorrelation()?

2. Or is it acceptable to set tumorType as a covariate like script above even when some of tumor types have just a few cell lines?

Thank you.

1. duplicateCorrelation only accepts one variable. You can't specify two random effects.

2. It seems perfectly natural to set tumorType as an ordinary factor. It wouldn't make sense as a random effect.