View Single Post
Old 08-18-2014, 08:05 PM   #19
Location: USA

Join Date: May 2014
Posts: 24

Dear Kopi-o and mikep

Thank you for your comments!!
Here is the specific analysis result I have done!

I have 12 samples from 6 patients with pre- and post- drug treatment.

For my data format

patient drug
KYM_pre 1 pre
KYM_post 1 post
GDY_pre 2 pre
GDY_post 2 post
HKS_pre 3 pre
HKS_post 3 post
KWK_pre 4 pre
KWK_post 4 post
PYM_pre 5 pre
PYM_post 5 post
SHM_pre 6 pre
SHM_post 6 post
For EdgeR

design <- model.matrix(~patient+drug)
data_e <- DGEList(counts=f_data)
data_e <- calcNormFactors(data_e)
data_e <- estimateGLMCommonDisp(data_e, design)
data_e <- estimateGLMTrendedDisp(data_e, design)
data_e <- estimateGLMTagwiseDisp(data_e, design)
## Fit the model, testing the coefficient for the treated vs untreated comparison
data_efit <- glmFit(data_e, design)
data_efit <- glmLRT(data_efit, coef="drugpost")

In this EdgeR analysis I got still only 125(16) DEG based on p value threshold 0.05 (pvalue threshold < 0.01) (no DEG based on adjusted p value)


For DESeq


meta <- data.frame(row.names=colnames(data), patient = patient, drug = drug)
d_fa <- newCountDataSet(data, meta)
d_fa <- estimateSizeFactors(d_fa)

d_fa <- estimateDispersions(d_fa, method="blind", sharingMode="fit-only")
(Is it okay??)
dh_fit1 = fitNbinomGLMs(d_fa, count ~ patient + enzyme)
dh_fit0 = fitNbinomGLMs(d_fa, count ~ patient)
dh_pvalsGLM = nbinomGLMTest( dh_fit1, dh_fit0 )
dh_padjGLM = p.adjust( dh_pvalsGLM, method="BH" )

dh_dtable <- transform(dh_fit1, pval=dh_pvalsGLM, padj=dh_padjGLM)
dh_dtable <- dh_dtable[order(dh_dtable$pval), ]

In this DESeq analysis , I only got 52(8) DEG gene based on p value threshold 0.05(pvalue < 0.01) (No DEG based on adjusted p values)..

And among those 125(16), 52(8) DEG genes from EdgeR, DESeq with p value threshold 0.05 (0.01 respectively) ..

In adddition, about the FC ratio, I got range of logFC value -2.7 to 2.7 for EdgeR.
DeSEq does not return FC ratio... so I am not sure how I can report FC value..

I assume that edgeR model and DEseq model should test the same thing... e.g. effect of drug considering the patient effect (paired analysis..) but it seems that reported DEG gene list looks halfly overlapped... though..

Anyhow, what do you think about this result?? Does it make sense??
younko is offline   Reply With Quote