Hi all,

I've been using limma with Illumina BeadChip files, already normalized by the quantile method. I have data from 677 cell lines and the aim is to determine the association between their polygenic risk scores and gene expression. In the end, I want to identify genes of which the expression associates with polygenic risk scores for several (8) traits I selected. In the model, I correct for several covariates, amongst other things the age of the individuals the cell lines are derived from. However, age is not available for all samples, which is why I did a sensitivity analysis for age for the samples (n=442) that do have age available. So in the end I have 3 analyses:

• The analysis I want to use for my paper (n=677 cell lines), not correcting for age
• The sensitivity analysis model 1 (n=442 cell lines), correcting for age
• The sensitivity analysis model 2 (n=442 cell lines), not correcting for age (to allow direct comparison with sensitivity analysis model 1).

I got some unexpected results, in that the sensitivity analysis model 2 resulted in way more DEGs than my original analysis, even though this involves less samples and both are not corrected for age. I already checked for the influence of outlier samples by running the method="robust" in the lmFit, but still the number of DEGs in the sensitivity analysis model 2 is in the thousands, while the number of DEGs in my original model is hundreds, so there's a very large difference. This can just be that the datasets inherently differ as one of them has more samples than the other, but I just want to make sure my code makes sense so I can rule that out. I would therefore greatly appreciate your input!

PS. I'm now just checking the results when running the model without any covariates, so just with the PRS, and the lmFit step takes way longer with just the PRS in the design when compared with PRS + covariates, which also seems weird to me.

###Take input arguments from bash command line (all disorder names) 
args <- commandArgs(trailingOnly=TRUE)
disorder <- args[1]

###Paths (removed literal names) 
wd  <-  "working directory "
d_out <- "output directory"
p_info <- "info file containing PRS per disorder, filenames of samples to read in, covariate information and sample replicate information"



# Create directories
if (!dir.exists(d_out)) dir.create(d_out)
setwd(wd)

## Read and clean info file and add "d" column to later block for technical replicate correlation
info <- read.table(p_info, header = TRUE, sep = " ") %>%
  mutate(d = as.numeric(as.factor(sub("\\d+(\\D*)$", "\\1", Name))))

# Read in expression files based on filename listed in column named "all"
y <- limma::read.ilmn(info$all, 
                      probeid = "PROBE_ID",
                      annotation = c("TargetID", "SYMBOL"),
                      expr = "AVG_Signal",
                      other.columns = "Detection",
                      sep = "\t")

# Filter for "truly expressed" genes based on detection above background (in my case, higher detection p-value means higher likelihood to be truly detected above background)
expressed <- rowSums(y$other$Detection > 0.95) >= ncol(y)/2
y <- y[expressed,]

#Log2 transformation + adding gene names as rownames (created new EList because my samples were normalized already, but not log transformed)
      E <- new("EList")
      E$E <- log2(y$E)
      E$genes <- y$genes
      E$genes$EntrezID <- AnnotationDbi::mapIds(org.Hs.eg.db, 
                                                y$genes$SYMBOL, 
                                                keytype = "SYMBOL", 
                                                column = "ENTREZID")
      rownames(E$E) <- E$genes$TargetID 


      #Make formula (this is without age as covariate, the other covariates are related to measures of pluripotency of my cell lines (they're iPSC lines) and potential batch effects + PC's) 
      covariates <- c("Sex","Method_of_derivation","Pluripotency_Score","Novelty_Score", "Open.access.data",paste0("PC",1:4))
      determinant <- paste0(disorder,"_PRS")
      base_formula <- paste("~", determinant, "+", paste(covariates, collapse= "+"))                     

      #Design matrix
      design  <-  model.matrix(as.formula(base_formula),data=info)

      #Run lmFit, hardcoded the average replicate correlation.

      fit <- limma::lmFit(E, design, block = info$d, correlation = 0.6320562, method="robust")
      fit2 <- limma::eBayes(fit)

      #Post-processing 
      sfit <- summary(decideTests(fit2,method="global"))
      z <- limma::topTable(fit2, coef = determinant, n = nrow(fit2$coefficients))

I can't check R code for complex data analyses, especially when I am not familiar with the datasets, but there are some non-standard aspects to your code that don't seem to follow.

• It is not possible to combine robust regression (method="robust") with blocking and correlation. If you set method="robust", then block and correlation will be ignored, which will artificially inflate the number of DE genes because the technical replicates will be treated as full replicates.
• Robust regression is intended to allow for individual outlier observations rather than outlier samples. Outlier samples are much better handled by arrayWeights(). While limma does offer robust regression as an option, I don't recommend or use it myself. Much better usually would be to set trend=TRUE and robust=TRUE at the eBayes() step.
• There is a disconnect between the first half of your code (reading the data) and the second half of the code (annotation and DE analysis). After running read.ilmn() the recommended next step would be neqc(), after which the data would be background corrected, log-transformed and normalized. The second half of the code however assumes that the data is normalized but not log-transformed. That seems very strange and the code used to undertake this normalization is not shown.
• The mapIds code is assumes all the SYMBOL values are still current gene symbols rather than aliases. I would prefer to use code that can map gene aliases as well as official symbols to Entrez Ids.
• The expression filtering is far more stringent than we recommend.
• You are hardcoding a block correlation, but there is no indication of how you estimated that correlation.
• If lmFit takes longer with just the PRS in the design when compared with PRS + covariates, that is a sure sign that you have changed other settings as well, for example block or method.