Entering edit mode
I am working with two subtypes of immune cells. For each of them I have two states - naive and activated. I am interested in the differential expression between these states in each cell type. I'll call these comparisons Contrast A and Contrast B. I am also interested to find out what is the difference between these two contrasts, i.e. to "compare the comparisons" (an illustration can be found here: http://bit.ly/1A1o6ZP). How can I model such comparison in edgeR?
My sample metadata looks like this: http://bit.ly/1EBZ9Kv
Modelling contrasts A and B is trivial (e.g. exactTest(d, pair=c("A-n", "A-act")). Can you write down the model for the "comparison of the comparisons"?
Ryan - thanks for your fast reply. I, actually, read that chapter and have previously used GLM models to account for batch effects and paired sample design. But I still have difficulty expressing the type of contrast that I am aiming at here. The problem is that the examples in the guide don't seem to include the type of comparison that I am aiming at. The closest example in the guide is 3.5 ("Comparison both between and within subjects"). I could map my case onto that case ("cell type" maps onto "Disease", "Activation" maps onto "Treatment"). The problem is that my samples are all independent whereas the example has paired patient samples. I am not quite sure how to change the model matrix to reflect that. I am sure this is something simple - but I still can't get my head round it.
You really just need a design of
~0+Condition, since your two experimental factors are already combined into a single factor in your sample metadata. Then, your contrast A will be(A.act - A.n)and your contrast B will be(B.act - B.n). So, to construct your contrast of contrasts, just subtract one from the other:(B.act - B.n) - (A.act - A.n). This is demonstrated in Subsection 3.3.1 "Defining each treatment combination as a group".(I have replaced the dashes with dots in your condition names using
make.namesso that they will be syntactically valid. Otherwse you will not be able to construct contrasts with them.)Ryan - thanks. This was really helpful. However, I have one question. It occurs to me that if the contrast is defined like this the test is based on the absolute differences between the treatments. I'll give an example below:
Say, there is a gene which in naive cells of type A has on average 10 reads. Upon activation it becomes 20. In cell type B the naive cells have on average 200 reads which, upon activation, become 400. If the contrast is defined like you suggest than what happens is that edgeR finds that in one case the absolute difference is just 10 (i.e. 20 minus 10), in the other it's 200. From this edgeR computes a big FC (20) and the gene gets a really low FDR. But, in fact, the gene works in both cell types in exactly the same way - it increases 2 fold as result of the activation. In essence, this contrast mixes the impact of the differences between the cell types with those caused by the activation. If there are big differences between the cell types these can dominate as in the example I gave.
I reckon what I was looking for was genes that work differently in these different cell types. I thought I could define the contrasts like this
(B.act/B.n) - (A.act/A.n). I tried it but I got:Error in qr.default(contrast) :NA/NaN/Inf in foreign function call (arg 1)
I reckon there are probably zero values somewhere which cause a problem. I could add some small pseudocounts but I am not sure this would be methodologically sound. In your opinion, what is the best way to approach the task that I am trying to solve?
In edgeR, coefficients are on a log scale, so subtraction is actually computing log ratios. A 2-fold change is always 1 on the log2 scale, regardless of the absolute expression.
Of course, some care is still required for zeros. In such cases, the log-fold change becomes undefined, so any contrast between log-fold changes will behave strangely. Consider the example below:
> y <- rbind(c(1000, 0, 10, 0), + c(1000, 1, 10, 0), + c(1000, 0, 10, 1), + c(1000, 1, 10, 1)) > design <- model.matrix(~0+factor(c("A.act", "A.n", "B.act", "B.n")) > fit <- glmFit(y, design, offset=0, dispersion=0.05) > lrt <- glmLRT(fit, contrast=c(1, -1, -1, 1)) > lrt$table$PValue [1] 0.999429030 0.887898066 0.002791298 0.012825046The first and second rows give large p-values as there is no strong evidence for differences in the log-FCs between A and B. However, the p-values for the third and fourth rows are much smaller, as more evidence is available when the appropriate log-FCs are defined. You can see that the size of the p-value is very sensitive to the gain or loss of a single read count. This complicates the interpretation of the results somewhat.
Ryan - thanks for your fast reply!
I am still a bit confused. I have posted a spreadsheet containing the example of two genes at http://bit.ly/1xbA7LS.
Gene X has the following values for type A activated: 0, 0.038634506, 0. In naive cells type A it has the following values: 0.032516233, 0, 0.036664022. In activated cells of type B it has the following values: 38.19672174, 40.45134583, 32.8059597. In naive cells of type B: 73.65871935, 78.06736892, 72.95833963.
In short, this gene behaves in both cell types quite similarly - its expression decreases by appr. 2 times upon activation. Its FDR is computed as 0 - this is one of the few genes that have this lowest possible FDR.
And there is gene Y. In activated cells of type A it has the following values: 7.384365199, 6.451962448, 7.640117933. In naive cells of type A its values are: 38.5967687, 43.05284305, 40.44041616. In activated cells of type B the values are: 16.14960138, 13.47144545, 15.79395832. In naive cells of type B its values are: 31.18923199, 37.51382043, 36.63212231.
In short, in cell type A the expression of gene Y decreases appr. 5 times upon activation, in cell type B it decreases only 2 times. Its FDR is computed as 0.0496034551314143 - thus it is considered "borderline" - very close to being "not significantly different".
Intuitively for me gene Y seems to change more than gene X. Yet, for edgeR the opposite is true. By the way, logFC computed are 11.3730230557303 for gene X and -0.210034235403313 for gene Y.
I am a bit puzzled by this result. Can you, perhaps, help to interpret it?
Are these CPM values? If so, did you compute it with
normalized.lib.sizes=TRUE? The size of the log-FC indicates that these values are not what edgeR is looking at internally.I think so. Here is my code:
##############################read in DGE
#############################
d <- readDGE(targets$Files)
#############################
#Filter & Library Size Re-set
#############################
keep <- rowSums(cpm(d)>1) >=3
d <- d[keep,]
d$samples$lib.size <- colSums(d$counts)
#############################
#Normalisation
#############################
d <- calcNormFactors(d)
#############################
#Recording the normalized counts
#############################
all_cpm=cpm(d, normalized.lib.size=TRUE)
all_counts <- cbind(rownames(all_cpm), all_cpm)
colnames(all_counts)[1] <- "Ensembl.Gene.ID"
rownames(all_counts) <- NULL
#############################
#Dispersion estimation
#############################
y <- estimateGLMCommonDisp(d, design)
y <- estimateGLMTrendedDisp(y, design)
y <- estimateGLMTagwiseDisp(y, design)
##############################Model fit & Test
#############################
fit <- glmFit(y, design)
my.contrasts <- makeContrasts(..)
lrt <- glmLRT(fit, contrast=my.contrasts[, "Activation-AvsActivation-B"])
#############################
#Results Exploration & Export
#############################
detags <- rownames(topTags(lrt, n=20))
cpm(d)[detags,]
summary(de <- decideTestsDGE(lrt, p=0.05, adjust="BH"))
tt <- topTags(lrt, n=nrow(d))
#############################
#Get the gene metadata
#############################
gene_hits <- cbind(rownames(tt$table), tt$table)
colnames(gene_hits)[1] <- "Ensembl.Gene.ID"
rownames(gene_hits) <- NULL
# with ens75mm_ens75_biomart <- read.table(file="biomaRt.csv",header=TRUE,sep=",")gene_hits_biomart <- merge(gene_hits, mm_ens75_biomart, by="Ensembl.Gene.ID")
#############################
#Merge the gene metadata and the counts data
#############################
gene_hits_biomart_counts <- merge(gene_hits_biomart, all_counts, by="Ensembl.Gene.ID")
#############################
#Filter, sort and store the merged table
#############################
numCols <- length(colnames(gene_hits_biomart_counts))
gene_hits_biomart_counts_filtered <- gene_hits_biomart_counts[gene_hits_biomart_counts$FDR < 0.05, c(1,2,5:numCols)]
merged_data_fdr_sorted <- gene_hits_biomart_counts_filtered[order(gene_hits_biomart_counts_filtered[,3]),]
write.csv(merged_data_fdr_sorted, file="Activation-AvsActivation-B.csv", row.names = FALSE)
How were
designandmy.contrastsformulated?Here are the design and contrasts:
> designA.act B.act C.act A.n B.n C.n
1 1 0 0 0 0 0
2 1 0 0 0 0 0
3 1 0 0 0 0 0
4 0 1 0 0 0 0
5 0 1 0 0 0 0
6 0 1 0 0 0 0
7 0 0 1 0 0 0
8 0 0 1 0 0 0
9 0 0 1 0 0 0
10 0 0 0 1 0 0
11 0 0 0 1 0 0
12 0 0 0 1 0 0
13 0 0 0 0 1 0
14 0 0 0 0 1 0
15 0 0 0 0 1 0
16 0 0 0 0 0 1
17 0 0 0 0 0 1
18 0 0 0 0 0 1
my.contrasts <- makeContrasts(Activation-AvsActivation-B = (A.act - A.n) - (B.act - B.n), levels=design)Can you show the counts and tagwise dispersion estimates for the offending rows, along with the normalization factors and library sizes? Also, you should quote
Activation-AvsActivation-Bin the call tomakeContrasts, otherwise an error is raised.Apologies - the actual my.contrast line is:
my.contrasts <- makeContrasts(ActivationR71vsActivationR72 = (AR71 - NR71) - (AR72 - NR72))which doesn't raise an exception.
AR71 stands for A.act, NR71 for A.n. AR72 is B.act, NR72 is B.n.
I just tried to stick with the original names that I used in the first post in order to avoid confusion.
Are you asking about the raw counts? The normalised are given in my post.
As for library sizes and norm factors - how can I show them? I haven't done this in quite a while so I am not sure how to get these values.
Yeah, the raw counts for the relevant genes would be nice. As for the library sizes and normalization factors; just show
y$samples$lib.sizeandy$samples$norm.factors. As for the dispersions, showy$tagwise.dispersion[...], subsetted for the relevant genes.I put the raw counts in a spreadsheet: http://bit.ly/1E4k2KL.
Here are the lib sizes and the norm factors:
y$samples$lib.size
[1] 29330462 26357396 22566633 28454909 26533897 23963450 33220755 27540983 26982561 31371013 26908984 27572867 21458072 28155268 26135268 34088199 25877118 19837215
> y$samples$norm.factors
[1] 0.9695880 0.9820241 0.9802117 1.0249482 1.0183244 1.0278016 0.9980380 0.9981807 0.9967784 0.9803276 0.9943822 0.9891860 1.0369650 0.9931721 1.0006368 1.0083422 1.0191137 0.9849640
As for the dispersion - how shall I do the subsetting? The gene IDs are ENSMUSG00000094525 and ENSMUSG00000022848.
I'm getting sensible results with the data that you've provided. I've stored the library sizes and normalization factors in "norm.txt", and the counts in "counts.txt", and I've run the following:
> grouping <- factor(rep(c("A.act", "B.act", "C.act", "A.n", "B.n", "C.n"), each=3)) > design <- model.matrix(~0+grouping) > colnames(design) <- levels(grouping) > > d <- read.table("counts.txt", header=TRUE, row.names=1) > norms <- read.table("norm.txt") > y <- DGEList(d, lib.size=as.numeric(norms[1,]), norm.factors=as.numeric(norms[2,])) > y$common.dispersion <- 0.05 > fit <- glmFit(y, design, prior.count=3) > > my.contrast <- makeContrasts(A.act - A.n - (B.act - B.n), levels=design) > lrt <- glmLRT(fit, contrast=my.contrast) > lrt$table logFC logCPM LR PValue ENSMUSG00000094525 0.8927441 4.500039 0.007924653 0.92906560 ENSMUSG00000022848 -0.9128559 4.711778 5.790652591 0.01611161See if it behaves like this on your machine (feed your own
yinto the commands for makingfitandlrt). My best guess is that you've somehow specified the wrong contrast in your original code.You are right Aaron - it must have been some sort of code mistake on my part. I couldn't reproduce what I saw before. Thank you very much for your help!