I have obtained TCGA gene expression RNAseq (polyA+ IlluminaHiSeq) Level 3 data. It has log2(x+1) transformed RSEM normalized counts.
For differential analysis I'm using "limma" package without voom. I came to know that these RSEM counts need to be normalized again with normalizeQuantiles.
So, I did that like following:
head(Rnadata)[1:4]
S1 S2 S3 S4
ARHGE 11.1818 11.2430 11.1612 12.0167
HIF3A 5.2482 4.0013 2.9374 4.7857
RNF17 4.1956 0.0000 0.0000 0.0000
RNF10 11.5047 12.0791 12.5931 11.4616
RNF11 9.5995 9.8248 9.9459 10.8368
RNF13 9.6257 10.5608 10.5179 10.1428
y <- normalizeQuantiles(Rnadata)
table(samples$type)
GP1n3 GP2
169 197
type2 <- factor(samples$type)
keep <- rowSums(y > log2(11)) >= 14
table(keep)
library("statmod")
y2 <- y[keep,]
design <- model.matrix(~type2)
fit <- lmFit(y2,design)
fit <- eBayes(fit)
tab <- topTable(fit,coef=2,adjust="BH")
This gave an output like following:
tab:
logFC AveExpr t P.Value adj.P.Val B
LOC152225 0.5177649 1.502743 3.782393 0.0001813158 0.9688791 -0.9982513
COL23A1 -0.6022330 3.743789 -3.679663 0.0002685533 0.9688791 -1.2085234
AGBL4 0.5032687 1.163802 3.532754 0.0004637616 0.9688791 -1.5001656
MPP6 0.3445659 7.334771 3.320850 0.0009875751 0.9688791 -1.9018043
PARVB -0.4983630 8.963399 -3.304377 0.0010456732 0.9688791 -1.9320752
GPX3 -0.5144697 13.830226 -3.292797 0.0010884016 0.9688791 -1.9532738
DKK3 -0.5604613 8.686580 -3.291027 0.0010950728 0.9688791 -1.9565075
C10orf116 -0.8018208 8.218658 -3.289999 0.0010989644 0.9688791 -1.9583847
ATOH8 -0.7274459 7.480381 -3.277855 0.0011459271 0.9688791 -1.9805237
MTF2 0.2071583 7.985538 3.265361 0.0011961812 0.9688791 -2.0032219
What I don't understand is why adj.P.Val is same for all the genes? Did I do any mistake? Can anyone help me in this? Thank you
Yes I know that n=Inf gives all the genes. I wanted to give GP2-GP1n3 in contrast matrix. And How can I create contrast.matrix for the above data?
With the above data when I see the results it gave the following.
results <- decideTests(fit)
summary(results)
(Intercept) type2GP2
-1 0 1
0 0 16905
1 16906 0
This means GP2-GP1n3 has 1 downregulated gene and 16905 doesn't show any regulation, 0 upregulation. Am I right?
The limma User's Guide covers all this in exquisite detail. Rather than asking questions that you could easily answer yourself, you should instead avail yourself of the available documentation.