Hello everybody,
I have one experiment with 3 conditions : cell lines of patient with one specific mutation (6 samples), cell lines of patient with another specific mutation (6 samples) and cell lines of controls (3 samples). I have 2 questions of interest: what are the deregulations between:
1/ cell lines of patient with one specific mutation (6 samples) vs cell lines of patient with another specific mutation (6 samples)?
2/ cell lines of patient (12 samples) vs cell lines of controls (3 samples)?
The results obtained from DESeq2 and edgeR were not very similar ( Differential expression analysis: comparison of DESeq2 and edgeR robust results ). Aaron Lun suggested some improvements. One of them is to use a single GLM model for answering both questions.
Could you please tell me if the two contrasts I used for this purpose are correct?
Thank you in advance.
library("DESeq2")
library("edgeR")
#####################################
### Creating CountDataSet object ###
#####################################
Name=c( "ctrl1_FET","ctrl2_HNS","ctrl3_MID", "PatientCellLine_COM","PatientCellLine_GAN","PatientCellLine_GEN","PatientCellLine_LEV","PatientCellLine_PAR","PatientCellLine_SAU","PatientCellLine_GRA","PatientCellLine_KOR","PatientCellLine_POP","PatientCellLine_RIP","PatientCellLine_SEG","PatientCellLine_SEV")
Status=c("Ctrl", "Ctrl", "Ctrl", "Mut2", "Mut2","Mut2", "Mut2", "Mut2","Mut2","Mut1", "Mut1", "Mut1","Mut1", "Mut1", "Mut1")
Fichier = paste(Name,"gencod.htseqCount",sep=".")
DataFrame=data.frame(Name,Fichier,Status)
################################################
### Model ###
### Estimation of sizeFactors and Dispersion ###
################################################
dds_raw=DESeqDataSetFromHTSeqCount(DataFrame,"/home/Results/Data/17052016", ~Status) # 63568
colData(dds_raw)$Status = factor(colData(dds_raw)$Status, levels=c("Ctrl","Mut2", "Mut1"))
colData(dds_raw)$Status = relevel(colData(dds_raw)$Status, "Ctrl")
keep <- rowSums(cpm( counts(dds_raw) )>1) >= 1
dds <- dds_raw[keep,]
dim(dds) #17904
dds=estimateSizeFactors(dds)
dds=estimateDispersions(dds)
#####################################
### Differential expression Tests ###
#####################################
dds=nbinomWaldTest(dds)
resultsNames(dds)
"Intercept" "StatusCtrl" "StatusMut2" "StatusMut1"
######################################################################################
# Comparison 1: Mut1 - Mut2
res1=results(dds,contrast=c(0,0,-1,1),pAdjustMethod="BH", cooksCutoff=FALSE,alpha=0.01, independentFiltering=TRUE)
summary(res1)
out of 17904 with nonzero total read count
adjusted p-value < 0.01
LFC > 0 (up) : 78, 0.44%
LFC < 0 (down) : 45, 0.25%
outliers [1] : 0, 0%
low counts [2] : 0, 0%
(mean count < 3)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
######################################################################################
# Comparison 2: Average of Mut1 and Mut2 to Ctrl
res2=results(dds,contrast=c(0,-1,1/2,1/2),pAdjustMethod="BH", cooksCutoff=FALSE, independentFiltering=TRUE,alpha=0.01)
summary(res2)
out of 17904 with nonzero total read count
adjusted p-value < 0.01
LFC > 0 (up) : 314, 1.8%
LFC < 0 (down) : 279, 1.6%
outliers [1] : 0, 0%
low counts [2] : 0, 0%
(mean count < 3)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results
Thank you a lot for your response.
For the first comparison, I had 145 deregulated genes with DESeq2 (80 up+65 down) and 210 with edgeR (137 up+73 down). Commonly up regulated genes: 54 and commonly down regulated genes: 37.
For the second comparison, I had 695 deregulated genes with DESeq2 (443 up+252 down) and 702 with edgeR (257 up+445 down). Commonly up regulated genes: 243 and commonly down regulated genes: 167.
I am trying to make the analysis more comparable as suggested by Aaron Lun and then, I will check again the intersection.
So for both comparisons, and splitting by direction of LFC, ~55-65% of the DESeq2 calls are in common. That sounds totally reasonable to me.
Like i said, if you look by rank instead of by FDR cutoff, it will look even more consistent. It's important to keep in mind that p-values are random variables, and very sensitive to differences in models.
Also, the whole reason I suggested you to use edgeR-robust for this dataset was that I guessed it would increase sensitivity for your genes with very high within-condition variability. So here you can see that in fact edgeR-robust is giving you a higher # of DEG at least in the first comparison.
Thank you for your feedback.
I though 50% of similarity was not sufficient, I had in mind 70-80%.
Also, I was surprised by the difference in proportion of up and down regulated genes in the second comparison. For a same number of deregulations, the repartition is not the same with the 2 tools. Is it common?
edgeR robust gave me indeed more deregulations with higher variability. So, if I want to keep edgeR results at the end, I see that 37.5, 39, 51 and 94% are found also by DESeq2 (depending on the comparison and direction). It will be probably a bit higher since I try to compare more properly the results (same input and new GLM). Can you please tell me if it is sufficient this way?
Does one expect more discrepancies between edgeR and DESeq2 when the effect is supposed to be small (as in my first comparison) or when one sample size is small (as the control group in my second experiment)?
I am not sure to understand well your suggestion for gene ranking. My genes are sort by FDR. There are a lot of common genes in the top deregulated genes obtained by the 2 tools, but the order is not perfectly the same. I can do the union of the deregulated genes, and look at their rank in both analyses. I can use Wilcoxon signed-rank test to test a difference and simply plot the ranks. Can I rank the genes by FDR and use only deregulated genes?
I think that this overlap is reasonable.
"Does one expect more discrepancies between edgeR and DESeq2 when the effect is supposed to be small (as in my first comparison) or when one sample size is small (as the control group in my second experiment)?"
Yes to both.
Don't worry about my ranking comment. It's not very important. Just that, genes may be just below the threshold line for tool X and above the threshold line for tool Y, although in general the genes are ranked relatively similarly (although not the same hence why I suggested you to try edgeR-robust).
Ok, I see.
I really appreciate your help, thank you.