4
$\begingroup$

I am new to the Limma package and when using voom I get the following plot. I am really not sure what it means. I wonder if anyone here have seen something like this before, or could explain it to me.Is my data is good or what?

Voom: Mean Variance Trend

 raw.data <-read.delim("Norma112_tumor282.txt")
library(limma)
library(edgeR)

attach(raw.data)
names(raw.data)
d <- raw.data[, 2:395]
rownames(d) <- raw.data[, 1]
pheno<-read.table("Pheno_Norma112_tumor282.txt", header=TRUE, sep="\t")
Group<-factor(pheno$Status,levels=levels(pheno$Status))
design<-model.matrix(~0+Group)

dge <- DGEList(counts=d)

########################### Tried Filtering from Here #####################
A <- rowSums(dge$counts)

isexpr <- A > 50 # Keeping genes with total counts more than 50.

##################################################################
dge <- calcNormFactors(dge)
v <- voom(dge,design,plot=TRUE)

write.table(v ,file="Voom_normalised_data_V1.txt",row.names=T, sep="\t")
##write.table(dge ,file="Voom_normalised_data_DGE.txt",row.names=T, sep="\t")

colnames(design) <-c("Normal", "Primary_Tumor")
colnames(design)
fit <-lmFit(v,design)
cont.wt<-makeContrasts("Primary_Tumor-Normal",levels=design)
fit2 <-contrasts.fit(fit,cont.wt)
fit3<-eBayes(fit2)
colnames(design)
DE1 = topTable(fit3, coef = 1, number = 'all')

But still i am getting the same results :(

$\endgroup$

2 Answers 2

3
$\begingroup$

Based on the range of values on your x axis, it appears you did not filter your raw counts prior to creating your voom EList object. Counts nearly 0 (plot x axis value -1) have low standard deviations. This rises immediately for low counts, then gradually decreases after count size of ~ 32 (plot x axis value 5). Below is the voom plot I generated today where the input transcript count data (8 experiments with raw counts ranging from 6 - 40 million) were first filtered by a threshold mean counts per million (cpm) value of 2. A voom plot I constructed today with mean cpm > 2

I can create a plot somewhat similar to yours by passing the whole raw matrix (20% of whose rows have counts of zero for all samples) to voom (via edgeR DGEList) here: enter image description here

The voom method is described: voom: precision weights unlock linear model analysis tools for RNA-seq read counts

Whether your data are "good" or not cannot be determined from this plot. Judging this would start back with the various metrics from raw fastq, alignment stage, then involve looking at specific properties of RNA-seq data. My lab helped developed a tool for the latter case, called RSeQC. You might consider investigating it.

But regarding voom, you can try the working case studies in Limma. They are there to help you become familiar with the program.

$\endgroup$
2
  • $\begingroup$ Hi Benjamin, I tried filtering out , but still getting the same. Please have a look on the following code which i am following..I have edited my question to include the code .Please have a look. Thanks a ton !!! $\endgroup$
    – David
    Jul 11, 2015 at 2:00
  • $\begingroup$ Hi Benjamin, Can you please share the codes, i tried my best but still the same, I have filtered genes with low expression so filtering gene with a normal read count less than 50. but still i am getting same plot. Please share your codes i will cross verify where i am going wrong. Thanks a Ton, It means a lot to me. $\endgroup$
    – David
    Jul 11, 2015 at 4:16
6
$\begingroup$

Actually you haven't filtered genes at all. You did compute a variable called isexpr, but then the variable was never applied in the subsequent analysis. So no surprise that the plot didn't change. To apply filtering you would have needed:

v <- voom(dge[isexpr,], design, plot=TRUE)

The voom plot shows how the coefficient of variation of the counts depends on the count size. To understand more, read the article that Benjamin Rodriguez links to. The meaning of the plot was discussed at length in that article.

You might also like to consult this recent workflow article: http://f1000research.com/articles/5-1408

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.