When I do this, I end up with the error described here (https://stats.stackexchange.com/questions/521917/understanding-why-design-matrix-not-of-full-rank-the-following-coefficients-no?noredirect=1#comment961380_521917). Dropping the T0 makes sense. However, my design matrix isn't full rank when I have [Treatment_Saline, Treatment_BMC, Temperature_26C, Temperature_30C] because of the Temperature_30C column. If I remove this column, then I am able to run the models. However, if I remove the column I don't know how I can look at the effects of say 30C vs. 26C in the Saline treatment. I'm using edgeR.glmLRT(fit, contrast=limma.makeContrasts(f"{treatment_class} - {reference_class}", levels=pandas_to_rpy2(dm) ) ) to do the contrasts.
How to structure my edgeR model to incorporate a temporal variable and 2 categorical variables?
My experiment is the following:
- Temperature 26C vs. Temperature 30C
- Treatment Saline vs. Treatment BMC
- Timepoints 0, 2, and 3.
My intention is to create GLM model to look at differential abundance between groups:
- 26C vs. 30C for Saline treatment
- 26C vs. 30C for BMC treatment
- Saline vs. BMC at 26C
- Saline vs. BMC at 30C
In doing so, I would like to include all the samples possible when calculating the dispersion which is why I'm using a GLM instead of a Fisher's Exact Test for subsets of samples. I would also like to incorporate the ordered time information.
# Functions
read_dataframe = function(path, sep="\t") {
df = read.table(path, sep=sep, row.names=1, header = TRUE, check.names=FALSE)
return(df)
}
# Counts
X = read_dataframe("https://pastebin.com/raw/J7kmL8Ly")
# OG0000000 OG0000001 OG0000002 OG0000003 OG0000004
# T2_10_SALINE_TEMP-PE-D710-D505-1_S10 16909 55 5382 5894 1964
# T2_11_BMC_CONTROL-PE-D711-D505-1_S11 24296 60 2772 3962 1374
# T2_12_BMC_CONTROL-PE-D712-D505-1_S12 24619 60 7351 5389 560
# T2_13_BMC_CONTROL-PE-D701-D506-1_S13 22420 15 2172 2778 930
# T2_14_BMC_CONTROL-PE-D702-D506-1_S14 20049 82 4655 6211 553
# Metadata
df_metadata = read_dataframe("https://pastebin.com/raw/PANaC3r5")
# temperature treatment collection_time_numeric
# 1_T0_RNA-PE-D711D501-1_S143 26C NaN 0
# 2_T0_RNA-PE-D709D506-1_S134 26C NaN 0
# 4_T0_RNA-PE-D709D505-1_S133 26C NaN 0
# 5_T0_RNA-PE-D709D504-1_S132 26C NaN
If the T0 timepoint is throwing everything off then it can be removed.
I'm trying to figure out how to do this from the following sources: https://www.bioconductor.org/packages/release/bioc/vignettes/edgeR/inst/doc/edgeRUsersGuide.pdf
However, my situation isn't described so please no RTFM responses. If I make a design matrix, it will basically be a binary vector for Treatment_BMC, a binary vector for Treatment_30C, and a numeric vector for the collection time.
If I use this as the design matrix when calculating the dispersion, then how would I for example do #1 above where I calculate the 26C vs. 30C for just the Saline treatment? Does this not make sense to calculate the dispersion for everything?
I could use some guidance a bit here.
1 answer
Your experiment is actually pretty standard. There are 12 treatment groups (2 temperatures x 2 treatments x 3 times) and you can setup the design matrix as explained in Section 3.3 of the edgeR User's Guide. The design matrix will have 1 column for each group.
However NaN values are not allowed in the metadata. You either have to fill in all the treatment values or else remove those samples.
• 0 views
•
link
written
by
Gordon Smyth
No, you haven't followed my advice. If you had, there wouldn't be a Temperature_30C column nor a problem with linear dependence.
You don't have to drop T0 (I wouldn't) but you obviously can't have NaNs.
edgeR is an R package that is part of Bioconductor. You are using a port of edgeR to Python and I cannot help you with Python code. In R, forming the design matrix is no problem and follows a standard process:
> library(edgeR)
Loading required package: limma
> targets <- read.delim("https://pastebin.com/raw/PANaC3r5",row.names=1)
> targets$treatment[1:4] <- "None"
> Group <- paste(targets$treatment,targets$temperature,targets$collection_time,sep=".")
> Group <- factor(Group)
> design <- model.matrix(~0+Group)
> colnames(design) <- levels(Group)
> table(Group)
Group
BMC.26C.2 BMC.26C.3 BMC.30C.2 BMC.30C.3 None.26C.0
5 5 5 5 4
Saline.26C.2 Saline.26C.3 Saline.30C.2 Saline.30C.3
5 6 5 5
> is.fullrank(design)
[1] TRUE
If you created the design matrix in this standard way then you would have no problems forming the contrasts you wish to make.
> design
BMC.26C.2 BMC.26C.3 BMC.30C.2 BMC.30C.3 None.26C.0 Saline.26C.2 Saline.26C.3 Saline.30C.2 Saline.30C.3
1 0 0 0 0 1 0 0 0 0
2 0 0 0 0 1 0 0 0 0
3 0 0 0 0 1 0 0 0 0
4 0 0 0 0 1 0 0 0 0
5 0 0 0 1 0 0 0 0 0
6 0 0 0 1 0 0 0 0 0
7 0 0 0 0 0 0 0 0 1
8 0 0 0 1 0 0 0 0 0
9 0 0 0 0 0 0 0 0 1
10 0 0 0 1 0 0 0 0 0
11 0 0 0 0 0 0 0 0 1
12 0 0 0 0 0 0 0 0 1
13 0 0 0 0 0 0 1 0 0
14 0 0 0 0 0 0 1 0 0
15 0 1 0 0 0 0 0 0 0
16 0 0 0 0 0 0 1 0 0
17 0 1 0 0 0 0 0 0 0
18 0 1 0 0 0 0 0 0 0
19 0 0 0 0 0 0 1 0 0
20 0 1 0 0 0 0 0 0 0
21 0 0 0 0 0 0 1 0 0
22 0 0 0 0 0 1 0 0 0
23 0 0 0 0 0 1 0 0 0
24 0 0 0 0 0 1 0 0 0
25 0 0 0 0 0 1 0 0 0
26 0 0 0 0 0 1 0 0 0
27 0 0 0 0 0 0 0 1 0
28 0 0 0 0 0 0 0 1 0
29 0 0 0 0 0 0 0 1 0
30 0 0 0 0 0 0 0 1 0
31 0 0 0 0 0 0 0 1 0
32 1 0 0 0 0 0 0 0 0
33 1 0 0 0 0 0 0 0 0
34 1 0 0 0 0 0 0 0 0
35 1 0 0 0 0 0 0 0 0
36 1 0 0 0 0 0 0 0 0
37 0 0 1 0 0 0 0 0 0
38 0 0 1 0 0 0 0 0 0
39 0 0 1 0 0 0 0 0 0
40 0 0 1 0 0 0 0 0 0
41 0 0 1 0 0 0 0 0 0
42 0 0 0 0 0 0 1 0 0
43 0 0 0 0 0 0 0 0 1
44 0 1 0 0 0 0 0 0 0
45 0 0 0 1 0 0 0 0 0
• 0 views
•
link
written
by
Gordon Smyth
Log in to answer this question.
More posts like this
Why have you posted this twice? - please see Understanding why `Design matrix not of full rank. The following coefficients not estimable` during GLM
I tried deleting the other one but couldn't find a delete. Is it possible to delete my other post that you referenced. Apologies.