lfcShrink function for different comparison
1
0
Entering edit mode
bi.mi • 0
@user-24135
Last seen 11 months ago

Dear all,

although I read through many advices and posts concerning the lfc Shrink function, I am still insecure if I used it correctly. I have 4 different conditions in my RNA seq data: CTRL_NoTreatment, CTRL_Treatment, Knockdown(KD)_NoTreatment and Knockdown_Treatment and I would like to do the following comparisons:

condition_CTRL_Treatment_vs_CTRL_NoTreatment, condition_KD_NoTreatment_vs_CTRL_NoTreatment, condition_KD_Treatment_vs_KD_NoTreatment and condition_KD_Treatment_vs_CTRL_Treatment.

Is it correct to run the dds function several times to obtain all required coef I need using different reference levels?

sampleTable$condition <- relevel(sampleTable$condition, ref = "CTRL_NoTreatment")
ddsHTSeq <- DESeqDataSetFromHTSeqCount(sampleTable=sampleTable, directory = read_dir, design= ~ condition)
dds <- DESeq(ddsHTSeq)
resLFC1 <- lfcShrink(dds, coef= "condition_CTRL_Treatment_vs_CTRL_NoTreatment", type = "apeglm")
resLFC2 <- lfcShrink(dds, coef= "condition_KD_NoTreatment_vs_CTRL_NoTreatment", type = "apeglm")

sampleTable$condition <- relevel(sampleTable$condition, ref = "KD_NoTreatment")             
ddsHTSeq <- DESeqDataSetFromHTSeqCount(sampleTable=sampleTable, directory = read_dir, design= ~ condition)
dds <- DESeq(ddsHTSeq)
resLFC3 <- lfcShrink(dds, coef = "condition_KD_Treatment_vs_KD_NoTreatment", type = "apeglm")

sampleTable$condition <- relevel(sampleTable$condition, ref = "CTRL_Treatment")  
ddsHTSeq <- DESeqDataSetFromHTSeqCount(sampleTable=sampleTable, directory = read_dir, design= ~ condition)
dds <- DESeq(ddsHTSeq)
resLFC4 <- lfcShrink(dds, coef = "condition_KD_Treatment_vs_CTRL_Treatment", type = "apeglm")

I'm using R version and DESeq2 version 1.28.1.

Thanks a lot in advance to help with my confusion and kind regards,

Lissy

DESeq2 RNASeq • 228 views
ADD COMMENT
3
Entering edit mode
@mikelove
Last seen 44 minutes ago
United States

You actually only need to run nbinomWaldTest after re-leveling. So you run DESeq() once, then re-level, then run nbinomWaldTest (to form the MLE estimates) and lfcShrink.

See here:

In this case, the dispersion does not have to be estimated again, as the designs are equivalent, up to the meaning of the coefficients. Instead, one need only run nbinomWaldTest to re-estimate MLE coefficients

https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#extended-section-on-shrinkage-estimators

ADD COMMENT
0
Entering edit mode

Thanks a lot, Michael! Now also this line in the manual makes sense to me :-).

In case anyone else comes across the same problem, this is how my code now looks :

sampleTable$condition <- relevel(sampleTable$condition, ref = "CTRL_NoTreatment")
ddsHTSeq <- DESeqDataSetFromHTSeqCount(sampleTable=sampleTable, directory = read_dir, design= ~ condition)
dds <- DESeq(ddsHTSeq)
resLFC1 <- lfcShrink(dds, coef= "condition_CTRL_Treatment_vs_CTRL_NoTreatment", type = "apeglm")
resLFC2 <- lfcShrink(dds, coef= "condition_KD_NoTreatment_vs_CTRL_NoTreatment", type = "apeglm")

dds$condition <- relevel(dds$condition, ref = "KD_NoTreatment")  
dds <- nbinomWaldTest(dds)
resLFC3 <- lfcShrink(dds, coef = "condition_KD_Treatment_vs_KD_NoTreatment", type = "apeglm")

dds$condition <- relevel(dds$condition, ref = "CTRL_Treatment")  
dds <- nbinomWaldTest(dds)
resLFC4 <- lfcShrink(dds, coef = "condition_KD_Treatment_vs_CTRL_Treatment", type = "apeglm")
ADD REPLY

Login before adding your answer.

Traffic: 356 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6