Hello everyone,

My dataset consists of 16 samples for which 784 genes (730 endogenous, 40 housekeeping, 8 negative controls and 6 positive controls) were quantified with NanoString nCounter technology in three different and known batches. My aim is to perform diferential expression between two groups of 10 and 6 samples.

After reading several posts about the use of DESeq2 on NanoString data, and following the tutorials and code provided in the GitHub linked in the publication An approach for normalization and quality control for NanoString RNA expression data, I have reached and impasse in my analysis. The housekeeping genes, which are supossed to be estable across conditions, do not meet such requirement. I have tested this using the following code on the raw data:

#### CHECK IF HK ARE ASSOCIATED WITH PRIMARY PHENO
hk_raw = raw[cIdx,]
pval = vector(length = nrow(hk_raw))

for (i in 1:nrow(hk_raw)){
  reg = glm.nb(as.numeric(hk_raw[i,]) ~ as.factor(pData$Group))
  pval[i] = coef(summary(reg))[2,4]
}

This results in 38 p-values under .05, which points to an association between expression of these so called housekeeping genes and the phenotype.

Following the pipeline, I estimated RUVg factors (using all of the 40 housekeeping genes) and used the following code to create the DESeq object:

dds <- DESeqDataSetFromMatrix(countData = counts(set)[c(1:730,745:784),],
                              colData = pData(set),
                              design = ~ W_1 + Batch + Group)

The following MA plot shows how most of the housekeeping genes, circled in orange, are not lying near the X axis, although the results for them show that they are not significantly differentially expressed (except 1 out of the 40).

My questions are:

1) Should I be using these "housekeeping genes" as reference or I cannot trust them because of the previous check? 2) If the answer to 1 is "no", what alternatives do I have? What is the next step towards differential expression in this dataset? Maybe another way of defining housekeeping genes? 3) Is the use of the known "Batch" variable correct in this case? Or is it redundant with the previously estimated RUVg factors?

Any help is much appreciated. Thanks!

In DESeq2, if you have housekeeping genes, you should provide these as controlGenes to estimateSizeFactors(). DESeq() runs this function as the first third of the statistical steps, but if you run it before DESeq() it will not over-write your size factors. So you can do:

dds <- estimateSizeFactors(dds, ...)
dds <- DESeq(dds)

This should put the housekeeping genes on the x-axis.