Dear all,

Please help. I am having trouble to execute Voom.

I have just only two samples without any replication (Arabidopsis WT and DKO). I am trying to get DEGs. I have written the script following the available tutorials.

Initially, it is okay. But whenever I tried to do Voom, it shows the following error message

Error in plot.window(...) : need finite 'ylim' values
In addition: Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf

I am not a biostatistician and have just started learning R script, so, I am not in a state to figure out the problems. I am sharing the script here as well.

#Install Biocounductor and requried packages
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("limma", version = "3.8")
library(limma)
BiocManager::install("edgeR", version = "3.8")
library(edgeR)

BiocManager::install("Glimma", version = "3.8")
library(Glimma)

#Set WD
setwd("C:/Users/Ornia/Documents/R/Limma Voom")

#Load files and Matrix count
files <- c("WT.tabular", "DKO.tabular")
read.delim(files[1], nrows = 10)
read.delim(files[2], nrow=10)

x <- readDGE(files)
class(x)

#Normalization
x <- calcNormFactors(x, method = "TMM")
x$samples$norm.factors

#Cutoff
cutoff <- 1
drop <- which(apply(cpm(X), 1, max) < cutoff)
D<- X[-drop,] 
dim(D) # number of genes left

#Group Naming
colnames(D) <- samplenames
group <- as.factor(c("WT", "DKO"))
D$samples$group <- group
D$samples

#Visualization of Normalization
x2 <- x
x2$samples$norm.factors <- 1
x2$counts[,1] <- ceiling(x2$counts[,1]*0.05)
x2$counts[,2] <- x2$counts[,2]*5

library(RColorBrewer)
nsamples <- ncol(x)
col <- brewer.pal(nsamples, "Greens")
par(mfrow=c(1,2))
lcpm <- cpm(x2, log=TRUE)
boxplot(lcpm, las=2, col=col, main="")
title(main="A. Example: Unnormalised data",ylab="Log-cpm")

x2 <- calcNormFactors(x2)  
x2$samples$norm.factors
lcpm <- cpm(x2, log=TRUE)
boxplot(lcpm, las=2, col=col, main="")
title(main="B. Example: Normalised data",ylab="Log-cpm")

#Design Matrix
design <- model.matrix(~0+group)
colnames(design) <- gsub("group", "", colnames(design))
design

#Contrast Setting
contr.matrix <- makeContrasts(WTvsDKO = WT-DKO,levels = colnames(design))
contr.matrix

#Voom
par(mfrow=c(1,2))
v <- voom(D, design, plot=TRUE)

vfit <- lmFit(v, design)
vfit <- contrasts.fit(vfit, contrasts=contr.matrix)
efit <- eBayes(vfit)
plotSA(efit, main="Final model: Mean-variance trend")

You can't use limma or voom if you have only two samples without any replication. It is a fundamental truth in statistics that you can't do any statistical analysis given only one observation in each group because there is no way to estimate variability. See for example the commentary here:

https://www.biostars.org/p/74155/

The voom function is designed to estimate the mean-variance trend from the residual variances. But variances are computed from the replicates and, without replication, there are no residual variances. That's what the error message is telling you -- there are no variances to plot.

I suggest you just compute RPKMs for each gene in each sample and look for genes with large fold changes.