Hi,

I'd like to use together EDASeq,RUVSeq and edgeR

I'd like to know if It's correct to run the folling code

According to EDASeq manual:

dataOffset <- withinLaneNormalization(data,"gc",which="upper",offset=TRUE) dataOffset <- betweenLaneNormalization(dataOffset,which="upper",offset=TRUE)

then according to RUVSeq using control genes taken from

Eisenberg E, Levanon EY. Human housekeeping genes, revisited. Trends Genet. 2013

I do:

set1 <- RUVg(x = dataOffset cIdx = HKgenes, k=1)

RUVg does NOT fill the offset slot of the EDASeq object.

To use the offsets I need the result from EDASeq but I still want to Remove Unwanted Variation

so I **modify **the edgeR paragraph of the EDASeq manual:

W_1 <- pData(set1)$ W_1 pData(dataOffset) <- cbind(pData(dataOffset),W_1) design <- model.matrix(~conditions + W_1, data=pData(dataOffset))

then all by manual

disp <- estimateGLMCommonDisp(counts(dataOffset), design, offset=-offst(dataOffset)) fit <- glmFit(counts(dataOffset), design, disp, offset=-offst(dataOffset)) lrt <- glmLRT(fit, coef=2) topTags(lrt)

By the way other than common dispersion (estimateGLMCommonDisp )

in edgeR manual there are 2 other steps (raccomanded):

To estimate trended dispersions:

y <- estimateGLMTrendedDisp(y, design)

To estimate tagwise dispersions:

y <- estimateGLMTagwiseDisp(y, design)

I need to integrate these steps too ?

Thank,

Pietro

Hi Davide,

following your tips I changed this dispersion:

with this dispersion

and changed the fit step:

with this:

Thanks,

Pietro

Hi Davide,

I'm going to visualize the results of the analysis using a simple boxplot and I find a myself troubled to decide which matrix to use.

Let's say I get lrt2 from glmLRT and I want to plot geneA

I need to give a vector to boxplot.

According to manuals you can use the log2(normCounts) or better use something like

or

I check for geneA FC

such FC is very different from the lrt2 FC so if I plot this any reviewer can argument boxplot is NOT representative of the reported FC.

I try to use counts or normalizedCounts (though if they worked I could not explain myself why adopt the discussed pipeline)

Also these FCs are far away from the lrt2 result.

So I try to get the data from lrt2 results, in particular i extract fitted.values and coefficients

these values are still very different from logFC = 5.605485 from glmLRT function.

Worst scenario is geneB

vsd FC is 0.42

counts FC is 1.95

normalizedCounts is 0.51

lrt2$coefficients is -0.75

ltr2$fitted.values FC is 1.93

As u can see only lrt2$coefficients goes near the value reported by lrt2 while the others have even different sign

I think I'm missing something, maybe more than that.

Please can you comment these results and pinpoint the correct matrix to use for boxplot.

Best,

Pietro

Hi Davide,

thanks for your reply. I explained myself the discrepancy in the same way BUT this means I can't use boxplot or similar plots to show the data.

I can do a HM of the genes using the edgeR "logFC" but i have no way to show the distribution of the samples around that value, right ?

By the way, is there no way to generate a "corrected" matrix using the confounding factors from RUV together with the estimates of the coefficient of the regression model ?

How you would visualize the result and explain it to a review asking to look at distribution of the expression of a gene in the two tested conditions ?

Here seems that RUV strongly changes the expression values (maybe because the total reads of the first condition are the double of the other condition).

Thank you very much,

P

Hi Davide,

my bad but i really can't understand the formula $Y - W alpha$. Can you please explain it to me ?

Please consider that in my example fit2 is the result of the glmFit having pData(set1)$ W_1 integrated into design and set1 is the result of RUVg.

Sorry to bother but i'd really like to fully implement and standardize the pipeline so I can use it every time (when reasonable).

Best,

P

Sorry, I was assuming that you were familiar with the notation we use in RUV, where $W$ is the matrix with the unwanted variation factors and $\alpha$ its associated parameter vector.

$Y - W \alpha$ essentially means to look at the coefficient estimates related to the unwanted variation factors in the glmFit result object, %*% that by the related columns of your design matrix and subtract the resulting matrix from the matrix of observed counts.

This would be almost like we provide as "normalizedCounts" in RUV, but after fitting the full model (including the variable of interests).

I hope this is more clear, sorry that there isn't an automated function to do that, but it would work on edgeR objects so I always thought it is out of the scope of the RUVSeq package. But if you want to contribute a PR to the RUVSeq package I would consider adding it to the package.

Hi Davide,

if I understand correctly, in my example:

$Y is the fit2$counts

$W$ is the lrt2$coefficients (Ngenes X 3conditions)

$\alpha$ is the design$W_1 (W_1 is a vector result of RUVg) (Nsamples X 1)

so

if I use all the columns of fit_TvsN$coefficients then

`MatrixByDavide2 <- fit_TvsN$coefficients%*%t(design$W_1) MatrixToMakePlots2 <- log2(fit_TvsN$counts)- MatrixByDavide2 bb <- aggregate(log2(MatrixToMakePlots2 [geneB,])~ (conditions),FUN= mean) bb[2,2]-bb[1,2] # 0.5405 instead of -1.088669`

Please let me know if I made something wrong because I'm still not able to get a matrix I can use to plot the results.

Best,

P