Load libraries

suppressPackageStartupMessages({
  library(SingleCellExperiment)
  library(edgeR)
  library(tidyverse)
  library(ggplot2)
  library(Seurat)
  library(glmGamPoi)
  library(dittoSeq)
  library(clustree)
  library(AnnotationDbi)
  library(org.Hs.eg.db)
  library(glue)
  library(speckle)
  library(patchwork)
  library(paletteer)
  library(tidyHeatmap)
  library(here)
})

set.seed(42)
options(scipen=999)
options(future.globals.maxSize = 6500 * 1024^2)

Load Data

files <- list.files(here("data/C133_Neeland_merged"),
                    pattern = "C133_Neeland_full_clean.*(macrophages|t_cells|other_cells)_annotated_diet.SEU.rds",
                    full.names = TRUE)

seuLst <- lapply(files[2:4], function(f) readRDS(f))

seu <- merge(seuLst[[1]], 
             y = c(seuLst[[2]], 
                   seuLst[[3]]))
seu

An object of class Seurat 
21568 features across 191521 samples within 1 assay 
Active assay: RNA (21568 features, 0 variable features)

             used    (Mb) gc trigger    (Mb) limit (Mb)   max used    (Mb)
Ncells   10361251   553.4   17390818   928.8         NA   11881777   634.6
Vcells 1328308325 10134.2 3629549058 27691.3      65536 3489907699 26625.9

Analyse Cell type proportions

# Differences in cell type proportions
props <- getTransformedProps(clusters = seu$Annotation,
                             sample = seu$sample.id, transform="asin")
props$Proportions %>% knitr::kable()

props$Proportions %>%
  data.frame %>%
  inner_join(seu@meta.data %>%
               dplyr::select(sample.id, 
                             Disease,
                             Treatment,
                             Status, 
                             Severity,
                             Group,
                             Batch, 
                             Age, 
                             Sex), 
             by = c("sample" = "sample.id")) %>%
  distinct()-> dat

ggplot(dat, aes(x = sample, y = Freq, fill = clusters)) +
  geom_bar(stat = "identity", color = "black", size = 0.1) +
  theme(axis.text.x = element_text(angle = 90,
                                   vjust = 0.5,
                                   hjust = 1),
        legend.text = element_text(size = 8),
        legend.position = "bottom") +
  labs(y = "Proportion", fill = "Cell Label") +
  #scale_fill_paletteer_d("miscpalettes::pastel", direction = 1) +
  facet_grid(~Group, scales = "free_x", space = "free_x")

props$Counts %>%
  data.frame %>%
  inner_join(seu@meta.data %>%
               dplyr::select(sample.id, 
                             Disease,
                             Treatment,
                             Status, 
                             Severity, 
                             Group,
                             Batch, 
                             Age, 
                             Sex), 
             by = c("sample" = "sample.id")) %>%
  distinct() -> dat

ggplot(dat, aes(x = sample, y = Freq, fill = Disease)) +
  geom_bar(stat = "identity") +
  theme(axis.text.x = element_text(angle = 90,
                                   vjust = 0.5,
                                   hjust = 1),
        legend.text = element_text(size = 8),
        legend.position = "bottom") +
  labs(y = "No. cells", fill = "Disease") +
  facet_grid(~Group, scales = "free_x", space = "free_x") +
  geom_hline(yintercept = 1000, linetype = "dashed")

seu <- subset(seu, cells = which(!seu$sample.id %in% c("sample_23.1")))

# Differences in cell type proportions
props <- getTransformedProps(clusters = seu$Annotation,
                             sample = seu$sample.id, transform="asin")
props$Proportions %>% knitr::kable()

props$Proportions %>%
  data.frame %>%
  inner_join(seu@meta.data %>%
               dplyr::select(sample.id, 
                             Disease,
                             Treatment,
                             Status, 
                             Severity, 
                             Group,
                             Batch, 
                             Age, 
                             Sex), 
             by = c("sample" = "sample.id")) %>%
  distinct()-> dat

ggplot(dat, aes(x = sample, y = Freq, fill = clusters)) +
  geom_bar(stat = "identity", color = "black", size = 0.1) +
  theme(axis.text.x = element_text(angle = 90,
                                   vjust = 0.5,
                                   hjust = 1),
        legend.text = element_text(size = 8),
        legend.position = "bottom") +
  labs(y = "Proportion", fill = "Cell Label") +
 # scale_fill_paletteer_d("miscpalettes::pastel", direction = 1) +
  facet_grid(~Group, scales = "free_x", space = "free_x")

props$Proportions %>%
  data.frame %>%
  inner_join(seu@meta.data %>%
               dplyr::select(sample.id,
                             Annotation,
                             Annotation,
                             Disease,
                             Treatment,
                             Status, 
                             Severity, 
                             Batch, 
                             Age, 
                             Sex), 
             by = c("sample" = "sample.id", "clusters" = "Annotation")) %>%
  distinct()-> dat

ggplot(dat, 
       aes(x = clusters, y = Freq, fill = clusters)) +
  geom_boxplot(outlier.size = 0.1, size = 0.25) +
  theme(axis.text.x = element_text(angle = 45,
                                   vjust = 1,
                                   hjust = 1),
        legend.text = element_text(size = 8)) +
  labs(y = "Proportion") +
#  scale_fill_paletteer_d("miscpalettes::pastel", direction = 1) +
  NoLegend()

Explore sources of variation

dims <- list(c(1,2), c(2:3), c(3,4), c(4,5))
p <- vector("list", length(dims))

for(i in 1:length(dims)){
  mds <- plotMDS(props$Counts,
    gene.selection = "common",
    plot = FALSE, dim.plot = dims[[i]])
  
  data.frame(x = mds$x, 
             y = mds$y,
             sample = rownames(mds$distance.matrix.squared)) %>%
    left_join(seu@meta.data %>%
               dplyr::select(sample.id, 
                             Disease,
                             Treatment,
                             Status, 
                             Severity, 
                             Batch, 
                             Age, 
                             Sex), 
             by = c("sample" = "sample.id")) %>%
    distinct() -> dat
  
  p[[i]] <- ggplot(dat, aes(x = x, y = y, 
                            shape = as.factor(Disease),
                            color = as.factor(Batch))) +
    geom_point(size = 3) +
    labs(x = glue("Principal Component {dims[[i]][1]}"),
         y = glue("Principal Component {dims[[i]][2]}"),
         colour = "Batch",
         shape = "Disease") +
    theme(legend.direction = "horizontal",
          legend.text = element_text(size = 8),
          legend.title = element_text(size = 9),
          axis.text = element_text(size = 8),
          axis.title = element_text(size = 9)) 
}

wrap_plots(p, cols = 2) + plot_layout(guides = "collect") &
  theme(legend.position = "bottom") 

dims <- list(c(1,2), c(2:3), c(3,4), c(4,5))
p <- vector("list", length(dims))

for(i in 1:length(dims)){
  mds <- plotMDS(props$Counts, 
                 gene.selection = "common",
                 plot = FALSE, dim.plot = dims[[i]])
  
   data.frame(x = mds$x, 
             y = mds$y,
             sample = rownames(mds$distance.matrix.squared)) %>%
    left_join(seu@meta.data %>%
               dplyr::select(sample.id, 
                             Disease,
                             Treatment,
                             Status, 
                             Severity,
                             Participant,
                             Batch, 
                             Age, 
                             Sex), 
             by = c("sample" = "sample.id")) %>%
    group_by(sample) %>%
    mutate(ncells = n()) %>%
    ungroup() %>%
    distinct() -> dat
  
  p[[i]] <- ggplot(dat, aes(x = x, y = y, 
                            colour = log2(ncells)))+
    geom_text(aes(label = str_remove_all(sample, "sample_")), size = 2.5) +
    labs(x = glue("Principal Component {dims[[i]][1]}"),
         y = glue("Principal Component {dims[[i]][2]}"),
         colour = "Log2 No. Cells") +
    theme(legend.direction = "horizontal",
          legend.text = element_text(size = 8),
          legend.title = element_text(size = 9),
          axis.text = element_text(size = 8),
          axis.title = element_text(size = 9)) +
    scale_colour_viridis_c(option = "magma")
}

wrap_plots(p, cols = 2) + plot_layout(guides = "collect") &
  theme(legend.position = "bottom") 

dims <- list(c(1,2), c(2:3), c(3,4), c(4,5))
p <- vector("list", length(dims))

for(i in 1:length(dims)){
  mds <- plotMDS(props$TransformedProps,
    gene.selection = "common",
    plot = FALSE, dim.plot = dims[[i]])
  
  data.frame(x = mds$x, 
             y = mds$y,
             sample = rownames(mds$distance.matrix.squared)) %>%
    left_join(seu@meta.data %>%
               dplyr::select(sample.id, 
                             Disease,
                             Treatment,
                             Status, 
                             Severity, 
                             Batch, 
                             Age, 
                             Sex), 
             by = c("sample" = "sample.id")) %>%
    distinct() -> dat
  
  p[[i]] <- ggplot(dat, aes(x = x, y = y, 
                            shape = as.factor(Disease),
                            color = as.factor(Batch)))+
    geom_point(size = 3) +
    labs(x = glue("Principal Component {dims[[i]][1]}"),
         y = glue("Principal Component {dims[[i]][2]}"),
         colour = "Batch",
         shape = "Disease") +
    theme(legend.direction = "horizontal",
          legend.text = element_text(size = 8),
          legend.title = element_text(size = 9),
          axis.text = element_text(size = 8),
          axis.title = element_text(size = 9))
}

wrap_plots(p, cols = 2) + plot_layout(guides = "collect") &
  theme(legend.position = "bottom") 

dims <- list(c(1,2), c(2:3), c(3,4), c(4,5))
p <- vector("list", length(dims))

for(i in 1:length(dims)){
  mds <- plotMDS(props$TransformedProps,
    gene.selection = "common",
    plot = FALSE, dim.plot = dims[[i]])
  
  data.frame(x = mds$x, 
             y = mds$y,
             sample = rownames(mds$distance.matrix.squared)) %>%
    left_join(seu@meta.data %>%
               dplyr::select(sample.id, 
                             Disease,
                             Treatment,
                             Status, 
                             Severity, 
                             Batch, 
                             Age, 
                             Sex), 
             by = c("sample" = "sample.id")) %>%
    distinct() -> dat
  
  p[[i]] <- ggplot(dat, aes(x = x, y = y, 
                            shape = as.factor(Disease),
                            color = Sex))+
    geom_point(size = 3) +
    labs(x = glue("Principal Component {dims[[i]][1]}"),
         y = glue("Principal Component {dims[[i]][2]}"),
         colour = "Sex",
         shape = "Disease") +
    theme(legend.direction = "horizontal",
          legend.text = element_text(size = 8),
          legend.title = element_text(size = 9),
          axis.text = element_text(size = 8),
          axis.title = element_text(size = 9))
}

wrap_plots(p, cols = 2) + plot_layout(guides = "collect") &
  theme(legend.position = "bottom") 

dims <- list(c(1,2), c(2:3), c(3,4), c(4,5))
p <- vector("list", length(dims))

for(i in 1:length(dims)){
  mds <- plotMDS(props$TransformedProps, 
                 gene.selection = "common",
                 plot = FALSE, dim.plot = dims[[i]])
  
   data.frame(x = mds$x, 
             y = mds$y,
             sample = rownames(mds$distance.matrix.squared)) %>%
    left_join(seu@meta.data %>%
               dplyr::select(sample.id, 
                             Disease,
                             Treatment,
                             Status, 
                             Severity,
                             Participant,
                             Batch, 
                             Age, 
                             Sex), 
             by = c("sample" = "sample.id")) %>%
    group_by(sample) %>%
    mutate(ncells = n()) %>%
    ungroup() %>%
    distinct() -> dat
  
  p[[i]] <- ggplot(dat, aes(x = x, y = y, 
                            colour = log2(Age)))+
    geom_text(aes(label = str_remove_all(sample, "sample_")), size = 2.5) +
    labs(x = glue("Principal Component {dims[[i]][1]}"),
         y = glue("Principal Component {dims[[i]][2]}"),
         colour = "Log2 Age") +
    theme(legend.direction = "horizontal",
          legend.text = element_text(size = 8),
          legend.title = element_text(size = 9),
          axis.text = element_text(size = 8),
          axis.title = element_text(size = 9)) +
    scale_colour_viridis_c(option = "magma")
}

wrap_plots(p, cols = 2) + plot_layout(guides = "collect") &
  theme(legend.position = "bottom") 

dims <- list(c(1,2), c(2:3), c(3,4), c(4,5))
p <- vector("list", length(dims))

for(i in 1:length(dims)){
  mds <- plotMDS(props$TransformedProps, 
                 gene.selection = "common",
                 plot = FALSE, dim.plot = dims[[i]])
  
   data.frame(x = mds$x, 
             y = mds$y,
             sample = rownames(mds$distance.matrix.squared)) %>%
    left_join(seu@meta.data %>%
               dplyr::select(sample.id, 
                             Disease,
                             Treatment,
                             Status, 
                             Severity,
                             Participant,
                             Batch, 
                             Age, 
                             Sex), 
             by = c("sample" = "sample.id")) %>%
    group_by(sample) %>%
    mutate(ncells = n()) %>%
    ungroup() %>%
    distinct() -> dat
  
  p[[i]] <- ggplot(dat, aes(x = x, y = y, 
                            colour = log2(ncells)))+
    geom_text(aes(label = str_remove_all(sample, "sample_")), size = 2.5) +
    labs(x = glue("Principal Component {dims[[i]][1]}"),
         y = glue("Principal Component {dims[[i]][2]}"),
         colour = "Log2 No. Cells") +
    theme(legend.direction = "horizontal",
          legend.text = element_text(size = 8),
          legend.title = element_text(size = 9),
          axis.text = element_text(size = 8),
          axis.title = element_text(size = 9)) +
    scale_colour_viridis_c(option = "magma")
}

wrap_plots(p, cols = 2) + plot_layout(guides = "collect") &
  theme(legend.position = "bottom") 

Principal components versus traits

Principal components analysis (PCA) allows us to mathematically determine the sources of variation in the data. We can then investigate whether these correlate with any of the specifed covariates. First, we calculate the principal components. The scree plot belows shows us that most of the variation in this data is captured by the top 7 principal components.

PCs <- prcomp(t(props$TransformedProps), center = TRUE, 
              scale = TRUE, retx = TRUE)
loadings = PCs$x # pc loadings
plot(PCs, type="lines") # scree plot

Collect all of the known sample traits.

nGenes = nrow(props$TransformedProps)
nSamples = ncol(props$TransformedProps)

info <- seu@meta.data %>%
  dplyr::select(sample.id, 
                Disease,
                Treatment,
                Status, 
                Severity,
                Participant,
                Group,
                Batch, 
                Age, 
                Sex) %>%
  group_by(sample.id) %>%
  mutate(ncells = n()) %>%
  ungroup() %>%
  distinct()

m <- match(colnames(props$TransformedProps), info$sample.id)
info <- info[m,]

datTraits <- info %>% dplyr::select(Participant, Batch, Disease, Status,
                                    Group, Severity, Age, Sex, ncells) %>%
  mutate(Age = log2(Age),
         ncells = log2(ncells),
    Donor = factor(Participant),
         Batch = factor(Batch),
         Disease = factor(Disease, 
                          labels = 1:length(unique(Disease))),
         Group = factor(Group, 
                        labels = 1:length(unique(Group))),
         Treatment = factor(Status, 
                            labels = 1:length(unique(Status))),
         Sex = factor(Sex, labels = length(unique(Sex))),
         Severity = factor(Severity, labels = length(unique(Severity)))) %>%
  mutate(across(everything(), as.numeric)) %>%
  dplyr::select(-Participant, -Status)

datTraits %>% 
  knitr::kable()

Correlate known sample traits with the top 10 principal components. This can help us determine which traits are potentially contributing to the main sources of variation in the data and should thus be included in our statistical analysis.

moduleTraitCor <- suppressWarnings(cor(loadings[, 1:10], datTraits, use = "p"))
moduleTraitPvalue <- WGCNA::corPvalueStudent(moduleTraitCor, (nSamples - 2))

textMatrix <- paste(signif(moduleTraitCor, 2), "\n(", 
                    signif(moduleTraitPvalue, 1), ")", sep = "")
dim(textMatrix) <- dim(moduleTraitCor)

## Display the correlation values within a heatmap plot
par(cex=0.75, mar = c(6, 8.5, 3, 3))
WGCNA::labeledHeatmap(Matrix = t(moduleTraitCor),
                      xLabels = colnames(loadings)[1:10],
                      yLabels = names(datTraits),
                      colorLabels = FALSE,
                      colors = WGCNA::blueWhiteRed(6),
                      textMatrix = t(textMatrix),
                      setStdMargins = FALSE,
                      cex.text = 1,
                      zlim = c(-1,1),
                      main = paste("PCA-trait relationships: Top 10 PCs"))

Statistical analysis using the propeller and limma approach

Create the design matrix.

group <- factor(info$Group)
participant <- factor(info$Participant)
age <- log2(info$Age)
batch <- factor(info$Batch)
sex <- factor(info$Sex)

design <- model.matrix(~ 0 + group + batch + age + sex)
colnames(design)[1:7] <- levels(group) 
design

   CF.IVA.M CF.IVA.S CF.LUMA_IVA.M CF.LUMA_IVA.S CF.NO_MOD.M CF.NO_MOD.S
1         0        0             0             0           0           0
2         0        0             0             0           1           0
3         0        0             0             0           1           0
4         0        0             0             0           1           0
5         0        0             0             0           1           0
6         0        0             0             0           1           0
7         0        0             0             0           0           0
8         0        0             0             0           0           1
9         0        0             0             0           0           1
10        1        0             0             0           0           0
11        1        0             0             0           0           0
12        0        0             0             0           0           1
13        0        0             0             0           1           0
14        1        0             0             0           0           0
15        0        0             0             0           1           0
16        0        0             0             0           0           1
17        0        0             0             0           0           0
18        0        0             0             0           0           1
19        0        0             0             0           1           0
20        0        0             0             0           1           0
21        0        1             0             0           0           0
22        0        0             0             0           1           0
23        0        0             0             0           1           0
24        0        0             0             0           1           0
25        0        0             0             0           1           0
26        0        0             0             0           1           0
27        0        0             0             0           0           1
28        0        0             0             1           0           0
29        0        0             0             0           1           0
30        0        0             1             0           0           0
31        0        0             1             0           0           0
32        1        0             0             0           0           0
33        1        0             0             0           0           0
34        1        0             0             0           0           0
35        0        0             0             0           1           0
36        0        0             1             0           0           0
37        0        0             0             0           0           0
38        0        0             0             0           0           1
39        0        0             0             0           1           0
40        0        1             0             0           0           0
41        0        0             0             0           0           0
42        0        0             0             0           0           0
43        0        0             0             0           0           0
44        0        0             0             0           0           0
   NON_CF.CTRL batch1 batch2 batch3 batch4 batch5 batch6         age sexM
1            1      0      0      1      0      0      0 -0.25908722    1
2            0      0      0      1      0      0      0 -0.09390014    1
3            0      0      0      1      0      0      0 -0.11514787    0
4            0      0      0      0      1      0      0 -0.04414710    0
5            0      0      0      0      1      0      0  0.14288337    1
6            0      0      0      0      0      1      0 -0.07296080    0
7            1      0      0      1      0      0      0  0.14645883    1
8            0      0      0      0      0      1      0  0.55970971    1
9            0      0      0      1      0      0      0  1.57438357    0
10           0      0      0      1      0      0      0  1.59938302    1
11           0      0      0      0      0      1      0  2.38835941    1
12           0      1      0      0      0      0      0  2.29572302    0
13           0      1      0      0      0      0      0  2.33608770    1
14           0      1      0      0      0      0      0  2.29801547    1
15           0      0      0      0      0      1      0  2.57902140    0
16           0      1      0      0      0      0      0  2.58232503    0
17           1      0      0      0      0      1      0  0.13210354    1
18           0      1      0      0      0      0      0  2.58890969    0
19           0      1      0      0      0      0      0  2.55836829    0
20           0      1      0      0      0      0      0  2.56706530    0
21           0      1      0      0      0      0      0  2.57305573    1
22           0      0      0      0      0      0      1 -0.93432383    0
23           0      0      0      0      0      0      1  0.09187369    0
24           0      0      0      0      0      0      1  1.04091644    0
25           0      0      0      0      0      0      1  0.08070438    1
26           0      0      0      0      0      0      1  0.99405890    1
27           0      0      0      0      0      0      1 -0.05642543    0
28           0      0      0      0      0      0      1  1.17649766    0
29           0      0      1      0      0      0      0  1.55970971    0
30           0      0      1      0      0      0      0  2.19301559    0
31           0      0      1      0      0      0      0  2.29801547    0
32           0      0      1      0      0      0      0  1.57039639    1
33           0      0      1      0      0      0      0  2.02060327    1
34           0      0      1      0      0      0      0  2.34855840    1
35           0      0      0      0      1      0      0  1.97307024    0
36           0      0      0      0      1      0      0  2.62971590    0
37           1      0      0      0      0      0      1  0.29237837    1
38           0      0      0      0      0      0      0  1.58014548    1
39           0      0      0      0      0      0      0  1.58014548    1
40           0      0      0      0      0      0      0  1.59931779    1
41           1      0      0      0      0      0      0  1.58496250    1
42           1      0      1      0      0      0      0  3.06991870    0
43           1      0      0      1      0      0      0  2.42046210    1
44           1      0      0      1      0      0      0  2.23560118    0
attr(,"assign")
 [1] 1 1 1 1 1 1 1 2 2 2 2 2 2 3 4
attr(,"contrasts")
attr(,"contrasts")$group
[1] "contr.treatment"

attr(,"contrasts")$batch
[1] "contr.treatment"

attr(,"contrasts")$sex
[1] "contr.treatment"

Create the contrast matrix.

contr <- makeContrasts(CF.NO_MODvCF.IVA = 0.5*(CF.NO_MOD.M + CF.NO_MOD.S) - 0.5*(CF.IVA.M + CF.IVA.S),
                       CF.NO_MODvCF.LUMA_IVA = 0.5*(CF.NO_MOD.M + CF.NO_MOD.S) - 0.5*(CF.LUMA_IVA.M + CF.LUMA_IVA.S),
                       CF.NO_MOD.MvCF.NO_MOD.S = CF.NO_MOD.M - CF.NO_MOD.S,
                       CF.NO_MODvNON_CF.CTRL = 0.5*(CF.NO_MOD.M + CF.NO_MOD.S) - NON_CF.CTRL,
                       CF.IVAvNON_CF.CTRL = 0.5*(CF.IVA.M + CF.IVA.S) - NON_CF.CTRL,
                       CF.LUMA_IVAvNON_CF.CTRL = 0.5*(CF.LUMA_IVA.M + CF.LUMA_IVA.S) - NON_CF.CTRL,
                       levels = design)
contr

               Contrasts
Levels          CF.NO_MODvCF.IVA CF.NO_MODvCF.LUMA_IVA CF.NO_MOD.MvCF.NO_MOD.S
  CF.IVA.M                  -0.5                   0.0                       0
  CF.IVA.S                  -0.5                   0.0                       0
  CF.LUMA_IVA.M              0.0                  -0.5                       0
  CF.LUMA_IVA.S              0.0                  -0.5                       0
  CF.NO_MOD.M                0.5                   0.5                       1
  CF.NO_MOD.S                0.5                   0.5                      -1
  NON_CF.CTRL                0.0                   0.0                       0
  batch1                     0.0                   0.0                       0
  batch2                     0.0                   0.0                       0
  batch3                     0.0                   0.0                       0
  batch4                     0.0                   0.0                       0
  batch5                     0.0                   0.0                       0
  batch6                     0.0                   0.0                       0
  age                        0.0                   0.0                       0
  sexM                       0.0                   0.0                       0
               Contrasts
Levels          CF.NO_MODvNON_CF.CTRL CF.IVAvNON_CF.CTRL
  CF.IVA.M                        0.0                0.5
  CF.IVA.S                        0.0                0.5
  CF.LUMA_IVA.M                   0.0                0.0
  CF.LUMA_IVA.S                   0.0                0.0
  CF.NO_MOD.M                     0.5                0.0
  CF.NO_MOD.S                     0.5                0.0
  NON_CF.CTRL                    -1.0               -1.0
  batch1                          0.0                0.0
  batch2                          0.0                0.0
  batch3                          0.0                0.0
  batch4                          0.0                0.0
  batch5                          0.0                0.0
  batch6                          0.0                0.0
  age                             0.0                0.0
  sexM                            0.0                0.0
               Contrasts
Levels          CF.LUMA_IVAvNON_CF.CTRL
  CF.IVA.M                          0.0
  CF.IVA.S                          0.0
  CF.LUMA_IVA.M                     0.5
  CF.LUMA_IVA.S                     0.5
  CF.NO_MOD.M                       0.0
  CF.NO_MOD.S                       0.0
  NON_CF.CTRL                      -1.0
  batch1                            0.0
  batch2                            0.0
  batch3                            0.0
  batch4                            0.0
  batch5                            0.0
  batch6                            0.0
  age                               0.0
  sexM                              0.0

Add random effect for samples from the same individual.

dupcor <- duplicateCorrelation(props$TransformedProps, design=design,
                                block=participant)
dupcor

$consensus.correlation
[1] 0.6250941

$cor
[1] 0.6250941

$atanh.correlations
 [1]  0.76811923  1.32273928  0.61140958  0.53633809  1.01612542  0.57559935
 [7]  0.72951448  1.24765882  1.01475251 -0.53606034  0.93797080  0.34377796
[13]  1.06833501  0.41397905  0.26228503  0.07433215  1.16780264  0.67068831
[19]  1.08639046  0.65338501  1.43476320  0.63351546  1.28350769  1.24738214
[25]  0.63317472  0.91591486  0.24769902  0.49044207  0.34408607  0.89512501
[31]  0.59396879 -0.53606034  0.07991134  0.58554384  0.54921399  0.58214487
[37]  1.72028140  1.57899738  1.15405369  0.53292992 -0.16389086  0.60681356
[43]  0.97776902  0.86778258

Fit the model.

fit <- lmFit(props$TransformedProps, design=design, block=participant, 
                correlation=dupcor$consensus)
fit2 <- contrasts.fit(fit, contr)
fit2 <- eBayes(fit2, robust=TRUE, trend=FALSE)
pvalue <- 0.05
summary(decideTests(fit2, p.value = pvalue))

       CF.NO_MODvCF.IVA CF.NO_MODvCF.LUMA_IVA CF.NO_MOD.MvCF.NO_MOD.S
Down                  0                     0                       0
NotSig               44                    44                      44
Up                    0                     0                       0
       CF.NO_MODvNON_CF.CTRL CF.IVAvNON_CF.CTRL CF.LUMA_IVAvNON_CF.CTRL
Down                       0                  0                       0
NotSig                    43                 44                      44
Up                         1                  0                       0

topTable(fit2)

                CF.NO_MODvCF.IVA CF.NO_MODvCF.LUMA_IVA CF.NO_MOD.MvCF.NO_MOD.S
neutrophil-like     -0.032103978           0.057343580          -0.03402256473
macro-IGF1           0.048434007           0.062714944           0.01695055528
macro-IFN            0.005205866           0.002539768           0.03003763351
monocytes           -0.012241009          -0.018333988           0.00006196445
HSP+ B cells         0.001953553           0.032950480          -0.01113271600
CD4 T-IFN           -0.013237507          -0.023120034           0.01914051357
migratory DC        -0.024505953          -0.026378390          -0.00393855044
macro-CCL18          0.054168504          -0.088129567          -0.00654047889
CD8 T-rm            -0.032141740          -0.056888959           0.01966709340
B cells             -0.029827960          -0.047618397          -0.05057953351
                CF.NO_MODvNON_CF.CTRL CF.IVAvNON_CF.CTRL
neutrophil-like           0.010971401        0.043075379
macro-IGF1                0.123586978        0.075152971
macro-IFN                 0.025393765        0.020187900
monocytes                -0.063451837       -0.051210829
HSP+ B cells              0.004247911        0.002294358
CD4 T-IFN                 0.007947308        0.021184815
migratory DC             -0.001232548        0.023273404
macro-CCL18               0.030308445       -0.023860059
CD8 T-rm                 -0.067488669       -0.035346929
B cells                  -0.048170472       -0.018342513
                CF.LUMA_IVAvNON_CF.CTRL     AveExpr        F     P.Value
neutrophil-like           -0.0463721791 0.032259739 4.289830 0.006998831
macro-IGF1                 0.0608720339 0.186514012 4.250580 0.007398005
macro-IFN                  0.0228539972 0.101845137 3.261977 0.024009216
monocytes                 -0.0451178496 0.096400155 3.089859 0.029716354
HSP+ B cells              -0.0287025686 0.005205621 3.011913 0.032748994
CD4 T-IFN                  0.0310673422 0.033190241 2.672084 0.050229354
migratory DC               0.0251458414 0.033207675 2.265013 0.084463190
macro-CCL18                0.1184380122 0.233027675 1.883140 0.138403204
CD8 T-rm                  -0.0105997106 0.117802877 1.707808 0.173431159
B cells                   -0.0005520755 0.099333335 1.628640 0.191987607
                adj.P.Val
neutrophil-like 0.1627561
macro-IGF1      0.1627561
macro-IFN       0.2881911
monocytes       0.2881911
HSP+ B cells    0.2881911
CD4 T-IFN       0.3683486
migratory DC    0.5309115
macro-CCL18     0.6637617
CD8 T-rm        0.6637617
B cells         0.6637617

p <- vector("list", ncol(contr))

for(i in 1:ncol(contr)){
props$Proportions %>% data.frame %>%
  left_join(info,
            by = c("sample" = "sample.id")) %>%
  dplyr::filter(Group %in%  names(contr[, i])[abs(contr[, i]) > 0]) -> dat
  
  if(length(unique(dat$Group)) > 2) dat$Group <- str_remove(dat$Group, ".(M|S)$")
ggplot(dat, aes(x = Group,
                y = Freq,
                colour = Group,
                group = Group)) +
  geom_boxplot(outlier.shape = NA, colour = "grey") +
  geom_jitter(stat = "identity",
              width = 0.15,
              size = 2) +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90,
                                   hjust = 1,
                                   vjust = 0.5),
        legend.position = "bottom",
        legend.direction = "horizontal") +
  labs(x = "Group", y = "Proportion",
       colour = "Condition") +
  facet_wrap(~clusters, scales = "free_y", ncol = 4) -> p[[i]]
}

p

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]