[younko]

Dear mikep and mbblack

Thank you so much for your comments.
I have one more quick question. I tried to analyze my data to find DEG in two condition. I have used edgeR, DESeq, DESeq2, cuffdiff or just simple wilconxon rank sum test etc.
However, I could not find any DEG with adjusted p value in any of these methods..
TOO UPSET

There are two hypothesis.
1) Hypothesis 1

Maybe, our experimental design is wrong so we could not detect any significant gene at all.. That is the reason we could not identify any DEG gene.. (End of story..) we could not do anything !!!
( I hope it is not !! cause there should be some genes responsing our drug)

2) Hypothesis 2

In some reason(small sample, so little statistical power), we could not find statistically significant DEG genes with adjusted p-value. (Still we could find some DEGs with p value) even though there exist DEG genes responsing the drug.

==============================================

In this case, what option I can choose?
I could not just throw away our data. I would like to find meaningful output........

So, I am thinking that (based on your suggestion)

Just use the pvalue to detect DEG genes and do any additional analysis with these DEG gene based on pvalue.. and then if we filter some more meaningful DEG among previously detected DEG genes,we can validate with experimental validation. (qPCR etc I am not sure how it would be feasible because of lack of money! ) ...
Or is there any other way we can claim our result without experimental validation??? Any suggestion??

Or Experiemnental validation is the only way I can do??

[kopi-o]

Did you try the design <- model.matrix(~drug+patientID) mentioned by mikep?

In many (perhaps most) of the analyses I have done of patient RNA-seq data, there are no differentially expressed genes found if you only look at the treatment/condition without telling the software to perform a paired analysis (specifying the patient ID in the design formula is one way of doing so). The reason being, as mbblack and mikep have already said, that the biological variation between individuals is obscuring the treatment effect.

[mikep]

That is unfortunate. You are however correct, further validation is absolutely required.

I have a question, how many genes pass uncorrected pv 0.05 and what is the range of the fold change you see?

EDIT: good point kopi-o, I assumed the "two condition" comment meant 2 factor, but it may not.

[younko]

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")
topTags(data_efit)

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)
str(dh_fit1)
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??