• 0 views
•
link
written
by
Kevin Blighe
Analysis with no Differentially Expressed Genes?
Hi Guys, I'm running my analysis for differential expression and the output give me no differentially expressed genes. Is there something wrong with my code or that's the true output of my samples?
First I imported the phenotype data:
phenoData <- read.AnnotatedDataFrame("design.txt", row.names = NULL, sep = "\t")
An object of class 'AnnotatedDataFrame' rowNames: 1 2 ... 6 (6 total) varLabels: Samples Treatment varMetadata: labelDescription
Then I ran the RMA:
eset <- justRMA("CN1_L10_T1.CEL","CN2_L10_T1.CEL","CN3_L10_T1.CEL",
"Ti50_1_L10_T1.CEL", "Ti50_2_L10_T1.CEL", "Ti50_3_L10_T1.CEL",
phenoData=phenoData)
I constructed the model matrix:
sample <- factor(rep(c("Cont", "Ti50"), each=3))
design <- model.matrix(~0 + sample)
colnames(design) <- levels(sample)
design
Cont Ti50
1 1 0
2 1 0
3 1 0
4 0 1
5 0 1
6 0 1
attr(,"assign") [1] 1 1 attr(,"contrasts") attr(,"contrasts")$sample [1] "contr.treatment"
I constructed the contrast matrix:
contrasts <- makeContrasts(Diff = Cont - Ti50, levels = design)
contrasts
Contrasts Levels Diff Cont 1 Ti50 -1
Then I ran the analysis:
fit <- lmFit(eset, design)
fit2 <- contrasts.fit(fit, contrasts)
efit <- eBayes(fit2)
deg <- topTable(efit, coef="Diff", p.value = 0.05, lfc = log2(1.5), adjust.method = "fdr", number = nrow(eset))
deg
data frame with 0 columns and 0 rows
I noticed that if I put p-value=1 I have 500 degs. So I investigated:
range(deg$adj.P.Val)
0.7866448 0.8526774
Indeed the threshold cuts all my genes. That's the true output ( no genes differentially expressed)? Or I did something wrong?
• 4,342 views
•
link
updated
by
Biostar
•
written
by
vitor.saldanha11 •
•
written
by
vitor.saldanha11 •
1 answer
It is definitely possible to have genes identified with a 3 x 3 comparisons, but I'm not sure of how clear the Ti50 effect is.
I would usually expect limma to be a little more sensitive than the built-in R functions (like aov() for ANOVA and lm() for linear-regression), although some of those pre-processing functions aren't exactly the same as I as expecting.
Perhaps you could try testing some functions this this demo dataset / tutorial?
https://www.bioconductor.org/help/course-materials/2005/BioC2005/labs/lab01/estrogen/
I don't believe I've followed this exact tutorial, but I would expect this dataset to have some genes identified with a low FDR (and usually an example dataset with clear expression changes would be a good idea for initially testing code).
• 0 views
•
link
written
by
Charles Warden
Log in to answer this question.
More posts like this
-
Why do I get zero results when doing DE on microarray data with lemma?
written by i.am.filippov •
-
Comparing expression of DEG gene in samples
written by Maryam •
-
heatmaps for limma-voom DE genes
written by AlexStar •
-
What is the right design matrix in a Microarray data analysis?
written by Paul Paodumai •
-
gene expression with .txt file
written by Joe Kherery
-
Problem with the design matrix in limma
written by ap-91 •
-
How to extract common genes that are DE in each comparisons?
written by rkpathakbt •
-
Setting up my design matrix
written by marion.ryan •
-
Extracting Gene Names From Venn Diagram
written by robjohn7000
-
Why Could I Only Get Two Different Expressed Genes By Limma?
written by chenyongzi77 •
Do you have references to back up your claim that nothing will ever come out of a microarray-based transcriptomics study comparing 3 vs. 3? The breast cancer study you've referenced my have had a very different goal than vitor.saldanha, which is why I find it difficult to embrace your general comment without a benchmarking reference.
I've mostly dealt with RNA-seq data where the current recommendation is 6 samples per condition (Schurch et al.), but most researchers still get by and get published with 3 replicates per condition. And in the case of well-controlled model organisms such as mice, worms etc. I don't think one would expect to see absolutely nothing if the experiment was well designed.
I did not mean to misrepresent what you said, thanks for clearing that up. I mostly wanted to point out that there might be reasons for why the OP is not detecting DEG which might be worth exploring. The replies you gave left me with the impression of generally discouraging these types of analyses when one had only 3 replicates -- when a lot of research efforts have gone into addressing the problems of calculating p- and q-values based on low sample numbers. I agree, one should not draw up new treatment guidelines for breast cancer patients based on a single 3 vs 3 comparison, but I'm almost certain that's not what the OP had in mind with this analysis.
Did you inspect the MDS plot (
plotMDS()function from limma)? Could you post the image here?Though you have used established and well used topTable function, Limma developers discourage such practice. Copy/pasted from their manual:
topTable does double filtering (by lfc and p-value(. There are/were considerable discussions among bioinformatics communities regarding this. To address this issue, limma recommends treat and topTreat method in case of filtering by fold change.
In one of the papers on fewer chips (3 was the least authors took) on pooled samples, conclusion was that for pooled samples, fewer chip numbers are on par with higher chip numbers with non pooled samples. I am not sure if it is general observation. However OP seems to have p and adjusted p on high side. I think OP should throw more light on top treat results and see if they are reported the same way in OP area of research
As much as you are knowledgeable, you are also being generous to others. I guess most of the recent posts are from beginners and do not read instructions laid by admins
Could you provide the reference?
I second h.mon's comment: what does your data look like in the global analyses such as PCA or MDS? If you do not see a clear grouping that follows the design of your experiment, most likely you're not going to get any DE genes. The main task should then be to find out whether that's due to noise or to possibly confounding factors.