Question: How To Make A Proper Use Of Blocking In Limma Using Voom.
2
gravatar for catagui
5.1 years ago by
catagui30
Miami, FL, USA
catagui30 wrote:

Hi,

I have a RNAseq data to analyse were I have a control and a one treatment for different individuals. I need to block the effects of the individual, but I am having several troubles to get the data that I need. I am using voom because my data is very heterogeneous and voom seams to do a good job normalising my reads.

I am having the following issues:

  1. I want to get the differentially expressed genes (DEGs) of my treatment not of my control. I don't understand after the eBayes analysis why I get the coefficients for both. I have tried a > makeContrasts (TreatvsCont= c2-co, levels = design) to subtract the control effect but then I get 0 DEGs.

  2. I am not sure when to include the 0 (null model) in the model formula, I have seen examples for both types of models.

This are my targets, with my column names of my counts, individual and condition

targets

Individual condition

A1 1 co

A2 3 co

A4 4 co

A5 5 co

E1 1 c2

E2 2 c2

E3 3 c2

E4 4 c2

E5 5 c2 `

This is the code I have been trying:

co2=as.matrix(read.table("2014_04_02_1h_PB.csv",header=T, sep=",", row.names=1))

nf = calcNormFactors (co2)

targets= read.table ("targets.csv", header = T, sep=",",row.names=1)

treat <- factor (targets$condition, levels= c("co", "c2"))

design <- model.matrix(~0+treat)

colnames (design) <- levels (treat)

y <- voom(co2,design,lib.size=colSums(co2)*nf)

corfit <- duplicateCorrelation(y,design,block=targets$Individual)

fit <- lmFit(y,design,block=targets$Individual,correlation=corfit$consensus)

fit2<- eBayes (fit)

results_trt <- topTable (fit2, coef="c2", n=nrow (y), sort.by="none")

From which gives me 18,000 genes with adj.P.Val < 0.01 out of 22,000 genes that I have in total. Which makes no sense..

Thanks in advance for the help.

rnaseq limma • 3.8k views
ADD COMMENTlink modified 3.2 years ago by srp0 • written 5.1 years ago by catagui30
2

With just two groups, you can just design <- model.matrix(~treat). Also, you might try getting rid of "individual 2", since that person isn't in both conditions.

Also, what sort of normalization factor distribution do you have? I've seen cases where including one sample with vastly fewer reads than the others has screwed up the analysis due to the normalization factors being off.

ADD REPLYlink modified 5.1 years ago • written 5.1 years ago by Devon Ryan89k
2
gravatar for julien.roux
5.0 years ago by
julien.roux90
Switzerland
julien.roux90 wrote:

Hi,

I think you need to call the voom function twice, see the following thread on the Bioconductor mailing list:

https://stat.ethz.ch/pipermail/bioconductor/attachments/20130530/4dcc9475/attachment.pl

Hope this helps

Julien

 

ADD COMMENTlink modified 5.0 years ago • written 5.0 years ago by julien.roux90
1

Just to follow up on Julien's recommendation. I had found conflicting advice on whether to run voom a second time, so I asked on the Bioconductor mailing list (thread). The consensus was that running voom a second time can only help the analysis and at worst will not make a difference.

ADD REPLYlink written 4.9 years ago by jdblischak60

I'd actually forgotten about that, thanks for the reminder!

ADD REPLYlink written 5.0 years ago by Devon Ryan89k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1409 users visited in the last hour