About dataset

The data is from a paper Generation of a microglial developmental index in mice and in humans reveals a sex difference in maturation and immune reactivity Published in Glia in 2017.

Analysis

From edgR package, quasi-likelihood method is used as is recommended for bulk RNA seq as it has strict error control rate by authors. Three different types of statistical model is shown as in edgeR manual.All three are quivalent as authors mentioned.

load libraries and data

library(edgeR)
## Loading required package: limma
library(limma)
count_data <- read.csv('/Users/nishapaudel/Documents/developmental_Sal/Counts.csv', sep =',', header=T, row.names = 1)
sampleinfo <- read.csv("/Users/nishapaudel/Documents/developmental_Sal/NEW_design.txt", sep = ',', header = TRUE, row.names = 1, stringsAsFactors = TRUE)
sampleinfo
##           Sex   Age Treatment
## F_E18_1     F   E18         N
## M_E18_1     M   E18         N
## F_E18_3     F   E18         N
## M_E18_4     M   E18         N
## F_P14_1     F   P14         N
## F_P14_2     F   P14         N
## F_P14_3     F   P14         N
## F_P14_4     F   P14         N
## F_P14_5     F   P14         N
## F_P14_6     F   P14         N
## F_P14_7     F   P14         N
## F_P14_8     F   P14         N
## F_P14_9     F   P14         N
## F_P14_10    F   P14         N
## M_P4_1      M    P4         N
## F_P4_2      F    P4         N
## F_P4_3      F    P4         N
## F_P4_4      F    P4         N
## F_P60_S_1   F P60_S         S
## F_P60_S_2   F P60_S         S
## F_P60_S_3   F P60_S         S
## F_P60_S_4   F P60_S         S
## F_P60_S_5   F P60_S         S
## M_P60_S_5   M P60_S         S
## F_P60_S_7   F P60_S         S
## F_E18_2     F   E18         N
## M_E18_2     M   E18         N
## M_E18_3     M   E18         N
## M_P14_1     M   P14         N
## M_P14_2     M   P14         N
## M_P14_3     M   P14         N
## M_P14_4     M   P14         N
## M_P14_5     M   P14         N
## M_P14_6     M   P14         N
## M_P14_7     M   P14         N
## M_P14_8     M   P14         N
## M_P14_9     M   P14         N
## M_P14_10    M   P14         N
## F_P4_1      F    P4         N
## M_P4_2      M    P4         N
## M_P4_3      M    P4         N
## M_P4_4      M    P4         N
## M_P60_S_1   M P60_S         S
## M_P60_S_2   M P60_S         S
## M_P60_S_3   M P60_S         S
## M_P60_S_4   M P60_S         S
## F_P60_S_6   F P60_S         S
dim(count_data)
## [1] 21987    47

1. Nested interaction formulas

age <- paste(sampleinfo$Age)
age<- factor(age)
sex <- paste(sampleinfo$Sex)
sex <- factor(sex)
#group <- factor(paste0(sampleinfo$Sex, ".", sampleinfo$Age)).  # not grouping here although group info bringes change 
y <- DGEList(count_data, group = age)
#sampleinfo$Sex <-relevel(sampleinfo$Sex, ref="M") ## relevel
design <- model.matrix(~ 0 + sex + sex:age, data = count_data)
design
##    sexF sexM sexF:ageP14 sexM:ageP14 sexF:ageP4 sexM:ageP4 sexF:ageP60_S
## 1     1    0           0           0          0          0             0
## 2     0    1           0           0          0          0             0
## 3     1    0           0           0          0          0             0
## 4     0    1           0           0          0          0             0
## 5     1    0           1           0          0          0             0
## 6     1    0           1           0          0          0             0
## 7     1    0           1           0          0          0             0
## 8     1    0           1           0          0          0             0
## 9     1    0           1           0          0          0             0
## 10    1    0           1           0          0          0             0
## 11    1    0           1           0          0          0             0
## 12    1    0           1           0          0          0             0
## 13    1    0           1           0          0          0             0
## 14    1    0           1           0          0          0             0
## 15    0    1           0           0          0          1             0
## 16    1    0           0           0          1          0             0
## 17    1    0           0           0          1          0             0
## 18    1    0           0           0          1          0             0
## 19    1    0           0           0          0          0             1
## 20    1    0           0           0          0          0             1
## 21    1    0           0           0          0          0             1
## 22    1    0           0           0          0          0             1
## 23    1    0           0           0          0          0             1
## 24    0    1           0           0          0          0             0
## 25    1    0           0           0          0          0             1
## 26    1    0           0           0          0          0             0
## 27    0    1           0           0          0          0             0
## 28    0    1           0           0          0          0             0
## 29    0    1           0           1          0          0             0
## 30    0    1           0           1          0          0             0
## 31    0    1           0           1          0          0             0
## 32    0    1           0           1          0          0             0
## 33    0    1           0           1          0          0             0
## 34    0    1           0           1          0          0             0
## 35    0    1           0           1          0          0             0
## 36    0    1           0           1          0          0             0
## 37    0    1           0           1          0          0             0
## 38    0    1           0           1          0          0             0
## 39    1    0           0           0          1          0             0
## 40    0    1           0           0          0          1             0
## 41    0    1           0           0          0          1             0
## 42    0    1           0           0          0          1             0
## 43    0    1           0           0          0          0             0
## 44    0    1           0           0          0          0             0
## 45    0    1           0           0          0          0             0
## 46    0    1           0           0          0          0             0
## 47    1    0           0           0          0          0             1
##    sexM:ageP60_S
## 1              0
## 2              0
## 3              0
## 4              0
## 5              0
## 6              0
## 7              0
## 8              0
## 9              0
## 10             0
## 11             0
## 12             0
## 13             0
## 14             0
## 15             0
## 16             0
## 17             0
## 18             0
## 19             0
## 20             0
## 21             0
## 22             0
## 23             0
## 24             1
## 25             0
## 26             0
## 27             0
## 28             0
## 29             0
## 30             0
## 31             0
## 32             0
## 33             0
## 34             0
## 35             0
## 36             0
## 37             0
## 38             0
## 39             0
## 40             0
## 41             0
## 42             0
## 43             1
## 44             1
## 45             1
## 46             1
## 47             0
## attr(,"assign")
## [1] 1 1 2 2 2 2 2 2
## attr(,"contrasts")
## attr(,"contrasts")$sex
## [1] "contr.treatment"
## 
## attr(,"contrasts")$age
## [1] "contr.treatment"
## Make the column names of the design matrix a bit nicer to avoid syntax error on downstream analysis
colnames(design) <- c("Female_18", "Male_18", "F14", "M14","F4", "M4","FP60_S", "MP60_S" )
design
##    Female_18 Male_18 F14 M14 F4 M4 FP60_S MP60_S
## 1          1       0   0   0  0  0      0      0
## 2          0       1   0   0  0  0      0      0
## 3          1       0   0   0  0  0      0      0
## 4          0       1   0   0  0  0      0      0
## 5          1       0   1   0  0  0      0      0
## 6          1       0   1   0  0  0      0      0
## 7          1       0   1   0  0  0      0      0
## 8          1       0   1   0  0  0      0      0
## 9          1       0   1   0  0  0      0      0
## 10         1       0   1   0  0  0      0      0
## 11         1       0   1   0  0  0      0      0
## 12         1       0   1   0  0  0      0      0
## 13         1       0   1   0  0  0      0      0
## 14         1       0   1   0  0  0      0      0
## 15         0       1   0   0  0  1      0      0
## 16         1       0   0   0  1  0      0      0
## 17         1       0   0   0  1  0      0      0
## 18         1       0   0   0  1  0      0      0
## 19         1       0   0   0  0  0      1      0
## 20         1       0   0   0  0  0      1      0
## 21         1       0   0   0  0  0      1      0
## 22         1       0   0   0  0  0      1      0
## 23         1       0   0   0  0  0      1      0
## 24         0       1   0   0  0  0      0      1
## 25         1       0   0   0  0  0      1      0
## 26         1       0   0   0  0  0      0      0
## 27         0       1   0   0  0  0      0      0
## 28         0       1   0   0  0  0      0      0
## 29         0       1   0   1  0  0      0      0
## 30         0       1   0   1  0  0      0      0
## 31         0       1   0   1  0  0      0      0
## 32         0       1   0   1  0  0      0      0
## 33         0       1   0   1  0  0      0      0
## 34         0       1   0   1  0  0      0      0
## 35         0       1   0   1  0  0      0      0
## 36         0       1   0   1  0  0      0      0
## 37         0       1   0   1  0  0      0      0
## 38         0       1   0   1  0  0      0      0
## 39         1       0   0   0  1  0      0      0
## 40         0       1   0   0  0  1      0      0
## 41         0       1   0   0  0  1      0      0
## 42         0       1   0   0  0  1      0      0
## 43         0       1   0   0  0  0      0      1
## 44         0       1   0   0  0  0      0      1
## 45         0       1   0   0  0  0      0      1
## 46         0       1   0   0  0  0      0      1
## 47         1       0   0   0  0  0      1      0
## attr(,"assign")
## [1] 1 1 2 2 2 2 2 2
## attr(,"contrasts")
## attr(,"contrasts")$sex
## [1] "contr.treatment"
## 
## attr(,"contrasts")$age
## [1] "contr.treatment"
keep <- filterByExpr(y)
y <- y[keep, , keep.lib.sizes=FALSE]
dim(y)
## [1] 12687    47
y <- calcNormFactors(y, method = "TMM")
y <- estimateDisp(y, design, robust=TRUE)
plotBCV(y)

fit <- glmQLFit(y, design, robust=TRUE)
head(fit$coefficients)
##                Female_18    Male_18         F14        M14          F4
## 0610007C21Rik  -9.649899  -8.573139  0.30495788 -0.7731417  0.10982246
## 0610007L01Rik  -9.523894  -9.753471  0.23432140  0.3873157  0.15171307
## 0610007P08Rik -10.371284 -10.725698  0.09829638  0.3703886  0.43629617
## 0610007P14Rik -10.246211  -9.509971 -0.77822732 -1.7004555 -0.22895706
## 0610007P22Rik -11.034784 -10.562801 -0.70076128 -1.1134980 -0.70705971
## 0610009B22Rik -11.299802 -10.325346 -0.13066829 -1.3919716  0.03123461
##                       M4     FP60_S     MP60_S
## 0610007C21Rik -0.9817081  0.7614027 -0.7066505
## 0610007L01Rik  0.3467170  0.1958497  0.4178045
## 0610007P08Rik  0.7361828  0.2738057  0.7871490
## 0610007P14Rik -1.0241219 -0.5889876 -1.3975521
## 0610007P22Rik -1.2745151 -0.7322822 -0.9722096
## 0610009B22Rik -1.0397484  0.7532606 -0.5399468
colnames(fit)
## [1] "Female_18" "Male_18"   "F14"       "M14"       "F4"        "M4"       
## [7] "FP60_S"    "MP60_S"
con <- makeContrasts(F_18vsF14 = F14-F4, levels=design)



res <-glmQLFTest(fit, contrast =con)
summary(decideTests(res))
##        1*F14 -1*F4
## Down          1455
## NotSig        9972
## Up            1260

Result of nested model

topTags(res)
## Coefficient:  1*F14 -1*F4 
##            logFC   logCPM         F       PValue          FDR
## Spp1   -5.035268 7.844012 247.02843 3.836219e-20 4.867012e-16
## H19    -4.997640 6.646876 140.53765 1.487423e-15 7.755379e-12
## Pmp22   2.317510 8.190749 138.02876 1.833856e-15 7.755379e-12
## Upk1b   4.847007 5.019680 129.74533 1.001635e-14 3.176934e-11
## Gpnmb  -4.044672 4.010152 117.02384 3.055514e-14 7.753060e-11
## Slc1a4 -2.623980 6.339590 110.44388 7.956584e-14 1.682420e-10
## Tnc    -4.975949 6.820992 108.42287 1.704156e-13 3.088661e-10
## Gpx3   -3.291390 4.974489  99.31117 4.425430e-13 7.018179e-10
## Rrm1   -1.636694 6.262948  95.17406 8.664240e-13 1.221369e-09
## Itgax  -2.849097 4.719643  94.37119 9.893642e-13 1.255206e-09

Rrm1 -1.636694 6.262948 95.17406 8.664240e-13 1.221369e-09

2. Factorial model

library(edgeR)
library(limma)
count_data <- read.csv('/Users/nishapaudel/Documents/developmental_Sal/Counts.csv', sep =',', header=T, row.names = 1)
sampleinfo <- read.csv("/Users/nishapaudel/Documents/developmental_Sal/NEW_design.txt", sep = ',', header = TRUE, row.names = 1, stringsAsFactors = TRUE)
sampleinfo
##           Sex   Age Treatment
## F_E18_1     F   E18         N
## M_E18_1     M   E18         N
## F_E18_3     F   E18         N
## M_E18_4     M   E18         N
## F_P14_1     F   P14         N
## F_P14_2     F   P14         N
## F_P14_3     F   P14         N
## F_P14_4     F   P14         N
## F_P14_5     F   P14         N
## F_P14_6     F   P14         N
## F_P14_7     F   P14         N
## F_P14_8     F   P14         N
## F_P14_9     F   P14         N
## F_P14_10    F   P14         N
## M_P4_1      M    P4         N
## F_P4_2      F    P4         N
## F_P4_3      F    P4         N
## F_P4_4      F    P4         N
## F_P60_S_1   F P60_S         S
## F_P60_S_2   F P60_S         S
## F_P60_S_3   F P60_S         S
## F_P60_S_4   F P60_S         S
## F_P60_S_5   F P60_S         S
## M_P60_S_5   M P60_S         S
## F_P60_S_7   F P60_S         S
## F_E18_2     F   E18         N
## M_E18_2     M   E18         N
## M_E18_3     M   E18         N
## M_P14_1     M   P14         N
## M_P14_2     M   P14         N
## M_P14_3     M   P14         N
## M_P14_4     M   P14         N
## M_P14_5     M   P14         N
## M_P14_6     M   P14         N
## M_P14_7     M   P14         N
## M_P14_8     M   P14         N
## M_P14_9     M   P14         N
## M_P14_10    M   P14         N
## F_P4_1      F    P4         N
## M_P4_2      M    P4         N
## M_P4_3      M    P4         N
## M_P4_4      M    P4         N
## M_P60_S_1   M P60_S         S
## M_P60_S_2   M P60_S         S
## M_P60_S_3   M P60_S         S
## M_P60_S_4   M P60_S         S
## F_P60_S_6   F P60_S         S
dim(count_data)
## [1] 21987    47
age <- paste(sampleinfo$Age)
age<- factor(age)
sex <- paste(sampleinfo$Sex)
sex <- factor(sex)
#group <- factor(paste0(sampleinfo$Sex, ".", sampleinfo$Age)).  # not grouping here although group info bringes change 
y <- DGEList(count_data, group = age)
#sampleinfo$Sex <-relevel(sampleinfo$Sex, ref="M") ## relevel
design <- model.matrix(~ 0+ sex*age)
design
##    sexF sexM ageP14 ageP4 ageP60_S sexM:ageP14 sexM:ageP4 sexM:ageP60_S
## 1     1    0      0     0        0           0          0             0
## 2     0    1      0     0        0           0          0             0
## 3     1    0      0     0        0           0          0             0
## 4     0    1      0     0        0           0          0             0
## 5     1    0      1     0        0           0          0             0
## 6     1    0      1     0        0           0          0             0
## 7     1    0      1     0        0           0          0             0
## 8     1    0      1     0        0           0          0             0
## 9     1    0      1     0        0           0          0             0
## 10    1    0      1     0        0           0          0             0
## 11    1    0      1     0        0           0          0             0
## 12    1    0      1     0        0           0          0             0
## 13    1    0      1     0        0           0          0             0
## 14    1    0      1     0        0           0          0             0
## 15    0    1      0     1        0           0          1             0
## 16    1    0      0     1        0           0          0             0
## 17    1    0      0     1        0           0          0             0
## 18    1    0      0     1        0           0          0             0
## 19    1    0      0     0        1           0          0             0
## 20    1    0      0     0        1           0          0             0
## 21    1    0      0     0        1           0          0             0
## 22    1    0      0     0        1           0          0             0
## 23    1    0      0     0        1           0          0             0
## 24    0    1      0     0        1           0          0             1
## 25    1    0      0     0        1           0          0             0
## 26    1    0      0     0        0           0          0             0
## 27    0    1      0     0        0           0          0             0
## 28    0    1      0     0        0           0          0             0
## 29    0    1      1     0        0           1          0             0
## 30    0    1      1     0        0           1          0             0
## 31    0    1      1     0        0           1          0             0
## 32    0    1      1     0        0           1          0             0
## 33    0    1      1     0        0           1          0             0
## 34    0    1      1     0        0           1          0             0
## 35    0    1      1     0        0           1          0             0
## 36    0    1      1     0        0           1          0             0
## 37    0    1      1     0        0           1          0             0
## 38    0    1      1     0        0           1          0             0
## 39    1    0      0     1        0           0          0             0
## 40    0    1      0     1        0           0          1             0
## 41    0    1      0     1        0           0          1             0
## 42    0    1      0     1        0           0          1             0
## 43    0    1      0     0        1           0          0             1
## 44    0    1      0     0        1           0          0             1
## 45    0    1      0     0        1           0          0             1
## 46    0    1      0     0        1           0          0             1
## 47    1    0      0     0        1           0          0             0
## attr(,"assign")
## [1] 1 1 2 2 2 3 3 3
## attr(,"contrasts")
## attr(,"contrasts")$sex
## [1] "contr.treatment"
## 
## attr(,"contrasts")$age
## [1] "contr.treatment"
colnames(design) <- c("F18", "M18", "F14", "F4","FP60_S","M14", "M4", "MP60_S" )
design
##    F18 M18 F14 F4 FP60_S M14 M4 MP60_S
## 1    1   0   0  0      0   0  0      0
## 2    0   1   0  0      0   0  0      0
## 3    1   0   0  0      0   0  0      0
## 4    0   1   0  0      0   0  0      0
## 5    1   0   1  0      0   0  0      0
## 6    1   0   1  0      0   0  0      0
## 7    1   0   1  0      0   0  0      0
## 8    1   0   1  0      0   0  0      0
## 9    1   0   1  0      0   0  0      0
## 10   1   0   1  0      0   0  0      0
## 11   1   0   1  0      0   0  0      0
## 12   1   0   1  0      0   0  0      0
## 13   1   0   1  0      0   0  0      0
## 14   1   0   1  0      0   0  0      0
## 15   0   1   0  1      0   0  1      0
## 16   1   0   0  1      0   0  0      0
## 17   1   0   0  1      0   0  0      0
## 18   1   0   0  1      0   0  0      0
## 19   1   0   0  0      1   0  0      0
## 20   1   0   0  0      1   0  0      0
## 21   1   0   0  0      1   0  0      0
## 22   1   0   0  0      1   0  0      0
## 23   1   0   0  0      1   0  0      0
## 24   0   1   0  0      1   0  0      1
## 25   1   0   0  0      1   0  0      0
## 26   1   0   0  0      0   0  0      0
## 27   0   1   0  0      0   0  0      0
## 28   0   1   0  0      0   0  0      0
## 29   0   1   1  0      0   1  0      0
## 30   0   1   1  0      0   1  0      0
## 31   0   1   1  0      0   1  0      0
## 32   0   1   1  0      0   1  0      0
## 33   0   1   1  0      0   1  0      0
## 34   0   1   1  0      0   1  0      0
## 35   0   1   1  0      0   1  0      0
## 36   0   1   1  0      0   1  0      0
## 37   0   1   1  0      0   1  0      0
## 38   0   1   1  0      0   1  0      0
## 39   1   0   0  1      0   0  0      0
## 40   0   1   0  1      0   0  1      0
## 41   0   1   0  1      0   0  1      0
## 42   0   1   0  1      0   0  1      0
## 43   0   1   0  0      1   0  0      1
## 44   0   1   0  0      1   0  0      1
## 45   0   1   0  0      1   0  0      1
## 46   0   1   0  0      1   0  0      1
## 47   1   0   0  0      1   0  0      0
## attr(,"assign")
## [1] 1 1 2 2 2 3 3 3
## attr(,"contrasts")
## attr(,"contrasts")$sex
## [1] "contr.treatment"
## 
## attr(,"contrasts")$age
## [1] "contr.treatment"
keep <- filterByExpr(y)
y <- y[keep, , keep.lib.sizes=FALSE]
dim(y)
## [1] 12687    47
y <- calcNormFactors(y, method = "TMM")
y <- estimateDisp(y, design, robust=TRUE)
plotBCV(y)

fit <- glmQLFit(y, design, robust=TRUE)
head(fit$coefficients)
##                      F18        M18         F14          F4     FP60_S
## 0610007C21Rik  -9.649899  -8.573139  0.30495788  0.10982246  0.7614027
## 0610007L01Rik  -9.523894  -9.753471  0.23432140  0.15171307  0.1958497
## 0610007P08Rik -10.371284 -10.725698  0.09829638  0.43629617  0.2738057
## 0610007P14Rik -10.246211  -9.509971 -0.77822732 -0.22895706 -0.5889876
## 0610007P22Rik -11.034784 -10.562801 -0.70076128 -0.70705971 -0.7322822
## 0610009B22Rik -11.299802 -10.325346 -0.13066829  0.03123461  0.7532606
##                      M14         M4     MP60_S
## 0610007C21Rik -1.0780996 -1.0915306 -1.4680532
## 0610007L01Rik  0.1529943  0.1950039  0.2219548
## 0610007P08Rik  0.2720923  0.2998866  0.5133433
## 0610007P14Rik -0.9222282 -0.7951648 -0.8085645
## 0610007P22Rik -0.4127368 -0.5674554 -0.2399274
## 0610009B22Rik -1.2613033 -1.0709831 -1.2932074
con <- makeContrasts(F_18vsF14 = F14-F4, levels=design)
res <-glmQLFTest(fit, contrast =con)
summary(decideTests(res))
##        1*F14 -1*F4
## Down          1455
## NotSig        9972
## Up            1260
#qlf <- glmQLFTest(fit, coef=2)

Result of Factorial model

topTags(res)
## Coefficient:  1*F14 -1*F4 
##            logFC   logCPM         F       PValue          FDR
## Spp1   -5.035268 7.844012 247.02843 3.836219e-20 4.867012e-16
## H19    -4.997640 6.646876 140.53765 1.487423e-15 7.755379e-12
## Pmp22   2.317510 8.190749 138.02876 1.833856e-15 7.755379e-12
## Upk1b   4.847007 5.019680 129.74533 1.001635e-14 3.176934e-11
## Gpnmb  -4.044672 4.010152 117.02384 3.055514e-14 7.753060e-11
## Slc1a4 -2.623980 6.339590 110.44388 7.956584e-14 1.682420e-10
## Tnc    -4.975949 6.820992 108.42287 1.704156e-13 3.088661e-10
## Gpx3   -3.291390 4.974489  99.31117 4.425430e-13 7.018179e-10
## Rrm1   -1.636694 6.262948  95.17406 8.664240e-13 1.221369e-09
## Itgax  -2.849097 4.719643  94.37119 9.893642e-13 1.255206e-09

Rrm1 -1.636694 6.262948 95.17406 8.664240e-13 1.221369e-0

3. Authors recommended design - Grouping for interaction

## group = group which must be used here, after using this we will get 10 DEGs but we are using same as above to verify
group <- factor(paste0(sampleinfo$Sex, ".", sampleinfo$Age))
y <- DGEList(count_data, group=age) 
keep <- filterByExpr(y, group = age)
y <- y[keep, , keep.lib.sizes=FALSE]
dim(y)
## [1] 12687    47
y <- normLibSizes(y) # Default normalization is TMM
y$samples
##           group lib.size norm.factors
## F_E18_1     E18  2266642    1.2034780
## M_E18_1     E18  1632456    1.2069230
## F_E18_3     E18  1899299    1.2036712
## M_E18_4     E18  1344033    0.7326719
## F_P14_1     P14  1728041    0.9153650
## F_P14_2     P14  2479263    0.9764640
## F_P14_3     P14  2222995    1.0227088
## F_P14_4     P14  2926328    0.9289068
## F_P14_5     P14  1915776    1.0160691
## F_P14_6     P14   858780    0.9943456
## F_P14_7     P14   889542    0.9924230
## F_P14_8     P14  3434800    0.8944291
## F_P14_9     P14  2970240    0.8817784
## F_P14_10    P14  2281179    0.9297602
## M_P4_1       P4  1841436    1.1957309
## F_P4_2       P4  2369986    1.1125975
## F_P4_3       P4  2726821    1.1419527
## F_P4_4       P4  1952722    1.1148018
## F_P60_S_1 P60_S  3111470    0.9104750
## F_P60_S_2 P60_S   956472    0.9093847
## F_P60_S_3 P60_S  2533368    0.9208736
## F_P60_S_4 P60_S   848961    0.8991337
## F_P60_S_5 P60_S   863715    0.8800336
## M_P60_S_5 P60_S  2330863    0.9652774
## F_P60_S_7 P60_S  1058489    0.9755245
## F_E18_2     E18  2275425    1.2117709
## M_E18_2     E18  2605667    1.2172500
## M_E18_3     E18  2693772    1.1610726
## M_P14_1     P14  2407063    0.8876348
## M_P14_2     P14  2059586    0.9540679
## M_P14_3     P14  2100263    0.9748077
## M_P14_4     P14  1022341    1.0277156
## M_P14_5     P14   955015    1.0414071
## M_P14_6     P14   956918    1.0625196
## M_P14_7     P14   920469    1.0328091
## M_P14_8     P14  1155002    0.9004358
## M_P14_9     P14  1085709    0.9052304
## M_P14_10    P14  2399697    0.9107533
## F_P4_1       P4  2506578    1.1609119
## M_P4_2       P4  2274999    1.1451445
## M_P4_3       P4  1882273    1.1458459
## M_P4_4       P4  3004914    1.0900635
## M_P60_S_1 P60_S  2338104    0.9420287
## M_P60_S_2 P60_S  3169592    0.8712846
## M_P60_S_3 P60_S  1941798    0.9676196
## M_P60_S_4 P60_S  2159941    0.9949783
## F_P60_S_6 P60_S   789274    0.8034786

Model Design

design <- model.matrix(~ 0 + group)
colnames(design) <- levels(group)
design
##    F.E18 F.P14 F.P4 F.P60_S M.E18 M.P14 M.P4 M.P60_S
## 1      1     0    0       0     0     0    0       0
## 2      0     0    0       0     1     0    0       0
## 3      1     0    0       0     0     0    0       0
## 4      0     0    0       0     1     0    0       0
## 5      0     1    0       0     0     0    0       0
## 6      0     1    0       0     0     0    0       0
## 7      0     1    0       0     0     0    0       0
## 8      0     1    0       0     0     0    0       0
## 9      0     1    0       0     0     0    0       0
## 10     0     1    0       0     0     0    0       0
## 11     0     1    0       0     0     0    0       0
## 12     0     1    0       0     0     0    0       0
## 13     0     1    0       0     0     0    0       0
## 14     0     1    0       0     0     0    0       0
## 15     0     0    0       0     0     0    1       0
## 16     0     0    1       0     0     0    0       0
## 17     0     0    1       0     0     0    0       0
## 18     0     0    1       0     0     0    0       0
## 19     0     0    0       1     0     0    0       0
## 20     0     0    0       1     0     0    0       0
## 21     0     0    0       1     0     0    0       0
## 22     0     0    0       1     0     0    0       0
## 23     0     0    0       1     0     0    0       0
## 24     0     0    0       0     0     0    0       1
## 25     0     0    0       1     0     0    0       0
## 26     1     0    0       0     0     0    0       0
## 27     0     0    0       0     1     0    0       0
## 28     0     0    0       0     1     0    0       0
## 29     0     0    0       0     0     1    0       0
## 30     0     0    0       0     0     1    0       0
## 31     0     0    0       0     0     1    0       0
## 32     0     0    0       0     0     1    0       0
## 33     0     0    0       0     0     1    0       0
## 34     0     0    0       0     0     1    0       0
## 35     0     0    0       0     0     1    0       0
## 36     0     0    0       0     0     1    0       0
## 37     0     0    0       0     0     1    0       0
## 38     0     0    0       0     0     1    0       0
## 39     0     0    1       0     0     0    0       0
## 40     0     0    0       0     0     0    1       0
## 41     0     0    0       0     0     0    1       0
## 42     0     0    0       0     0     0    1       0
## 43     0     0    0       0     0     0    0       1
## 44     0     0    0       0     0     0    0       1
## 45     0     0    0       0     0     0    0       1
## 46     0     0    0       0     0     0    0       1
## 47     0     0    0       1     0     0    0       0
## attr(,"assign")
## [1] 1 1 1 1 1 1 1 1
## attr(,"contrasts")
## attr(,"contrasts")$group
## [1] "contr.treatment"
y <- estimateDisp(y, design)
y$common.dispersion
## [1] 0.09379154

Fitting to GLM

fit <- glmQLFit(y, design, robust=TRUE)
head(fit$coefficients)
##                    F.E18      F.P14       F.P4    F.P60_S      M.E18      M.P14
## 0610007C21Rik  -9.649899  -9.344941  -9.540077  -8.888496  -8.573139  -9.346281
## 0610007L01Rik  -9.523894  -9.289573  -9.372181  -9.328044  -9.753471  -9.366155
## 0610007P08Rik -10.371284 -10.272988  -9.934988 -10.097478 -10.725698 -10.355309
## 0610007P14Rik -10.246211 -11.024438 -10.475168 -10.835199  -9.509971 -11.210427
## 0610007P22Rik -11.034784 -11.735545 -11.741844 -11.767066 -10.562801 -11.676299
## 0610009B22Rik -11.299802 -11.430470 -11.268567 -10.546541 -10.325346 -11.717318
##                     M.P4    M.P60_S
## 0610007C21Rik  -9.554847  -9.279790
## 0610007L01Rik  -9.406754  -9.335666
## 0610007P08Rik  -9.989515  -9.938549
## 0610007P14Rik -10.534093 -10.907523
## 0610007P22Rik -11.837316 -11.535010
## 0610009B22Rik -11.365095 -10.865293

Result of grouping model

topTags(res)
## Coefficient:  1*F14 -1*F4 
##            logFC   logCPM         F       PValue          FDR
## Spp1   -5.035268 7.844012 247.02843 3.836219e-20 4.867012e-16
## H19    -4.997640 6.646876 140.53765 1.487423e-15 7.755379e-12
## Pmp22   2.317510 8.190749 138.02876 1.833856e-15 7.755379e-12
## Upk1b   4.847007 5.019680 129.74533 1.001635e-14 3.176934e-11
## Gpnmb  -4.044672 4.010152 117.02384 3.055514e-14 7.753060e-11
## Slc1a4 -2.623980 6.339590 110.44388 7.956584e-14 1.682420e-10
## Tnc    -4.975949 6.820992 108.42287 1.704156e-13 3.088661e-10
## Gpx3   -3.291390 4.974489  99.31117 4.425430e-13 7.018179e-10
## Rrm1   -1.636694 6.262948  95.17406 8.664240e-13 1.221369e-09
## Itgax  -2.849097 4.719643  94.37119 9.893642e-13 1.255206e-09

Rrm1 -1.636694 6.262948 95.17406 8.664240e-13 1.221369e-09

Session information

sessionInfo()
## R version 4.3.0 (2023-04-21)
## Platform: aarch64-apple-darwin20 (64-bit)
## Running under: macOS 14.0
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Australia/Melbourne
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] edgeR_3.42.4 limma_3.56.2
## 
## loaded via a namespace (and not attached):
##  [1] digest_0.6.33     R6_2.5.1          fastmap_1.1.1     xfun_0.40        
##  [5] lattice_0.21-8    splines_4.3.0     cachem_1.0.8      knitr_1.43       
##  [9] htmltools_0.5.6   rmarkdown_2.24    cli_3.6.1         grid_4.3.0       
## [13] sass_0.4.7        jquerylib_0.1.4   statmod_1.5.0     compiler_4.3.0   
## [17] highr_0.10        rstudioapi_0.15.0 tools_4.3.0       evaluate_0.21    
## [21] bslib_0.5.1       Rcpp_1.0.11       yaml_2.3.7        locfit_1.5-9.8   
## [25] rlang_1.1.1       jsonlite_1.8.7
getwd()
## [1] "/Users/nishapaudel"