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")
Thanks for your response.
I already counted FPKM and confirmed some interesting DEGs with qPCR.
I was trying to avoid galaxy and learning some fastest way for future. But seems like I cannot use any recognized packages.
Thanks once again for your prompt response.
I don't know what this has to do with galaxy. You can compute FPKM in R very easily.