Load Libraries

library(Seurat)
library(WGCNA)
library(cluster)
library(parallelDist)
library(ggsci)
library(emmeans)
library(lme4)
library(ggbeeswarm)
library(genefilter)
library(tidyverse)
library(reshape2)
library(igraph)
library(gProfileR)
library(ggpubr)
library(cowplot)
library(here)
library(tidygraph)
library(ggraph)
library(ggExtra)

Extract Cells for WGCNA

Calculate softpower

enableWGCNAThreads()

Allowing parallel execution with up to 79 working processes.

datExpr <- as.matrix(t(ventric[["SCT"]]@scale.data[ventric[["SCT"]]@var.features, ]))
gsg <- goodSamplesGenes(datExpr, verbose = 3)

 Flagging genes and samples with too many missing values...
  ..step 1

powers <- c(c(1:10), seq(from = 12, to = 40, by = 2))
sft <- pickSoftThreshold(datExpr,
  dataIsExpr = TRUE, powerVector = powers, corOptions = list(use = "p"),
  networkType = "signed"
)

   Power SFT.R.sq  slope truncated.R.sq  mean.k. median.k.   max.k.
1      1 0.332000 604.00          0.469 2.50e+03  2.50e+03 2.51e+03
2      2 0.210000 247.00          0.507 1.25e+03  1.25e+03 1.26e+03
3      3 0.055600  77.10          0.534 6.25e+02  6.25e+02 6.33e+02
4      4 0.014400  29.80          0.568 3.13e+02  3.13e+02 3.18e+02
5      5 0.000624   4.63          0.600 1.57e+02  1.56e+02 1.60e+02
6      6 0.011000 -16.50          0.549 7.83e+01  7.83e+01 8.08e+01
7      7 0.042900 -27.00          0.484 3.92e+01  3.92e+01 4.08e+01
8      8 0.117000 -38.60          0.377 1.96e+01  1.96e+01 2.07e+01
9      9 0.220000 -90.40          0.245 9.82e+00  9.81e+00 1.05e+01
10    10 0.253000 -41.30          0.343 4.91e+00  4.91e+00 5.34e+00
11    12 0.533000 -88.80          0.399 1.23e+00  1.23e+00 1.40e+00
12    14 0.858000 -40.20          0.980 3.09e-01  3.08e-01 3.74e-01
13    16 0.830000 -25.80          0.944 7.75e-02  7.73e-02 1.03e-01
14    18 0.795000 -16.90          0.931 1.95e-02  1.94e-02 2.95e-02
15    20 0.753000 -10.60          0.936 4.90e-03  4.87e-03 9.40e-03
16    22 0.514000 -13.20          0.378 1.23e-03  1.22e-03 3.73e-03
17    24 0.540000 -10.10          0.414 3.11e-04  3.07e-04 1.71e-03
18    26 0.550000  -7.86          0.432 7.89e-05  7.72e-05 8.72e-04
19    28 0.498000  -6.05          0.390 2.01e-05  1.94e-05 4.71e-04
20    30 0.555000  -5.08          0.435 5.20e-06  4.88e-06 2.62e-04
21    32 0.525000  -4.18          0.425 1.38e-06  1.23e-06 1.49e-04
22    34 0.536000  -3.59          0.424 3.82e-07  3.10e-07 8.47e-05
23    36 0.543000  -3.15          0.426 1.14e-07  7.80e-08 4.85e-05
24    38 0.532000  -2.76          0.489 3.85e-08  1.97e-08 2.78e-05
25    40 0.543000  -2.53          0.490 1.49e-08  4.96e-09 1.60e-05

cex1 <- 0.9
plot(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2], xlab = "Soft Threshold (power)", ylab = "Scale Free Topology Model Fit, signed R^2", type = "n", main = paste("Scale independence"))
text(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[, 2], labels = powers, cex = cex1, col = "red")
abline(h = 0.80, col = "red")

Past versions of unnamed-chunk-4-1.png

# Mean Connectivity Plot
plot(sft$fitIndices[, 1], sft$fitIndices[, 5], xlab = "Soft Threshold (power)", ylab = "Mean Connectivity", type = "n", main = paste("Mean connectivity"))
text(sft$fitIndices[, 1], sft$fitIndices[, 5], labels = powers, cex = cex1, col = "red")

Past versions of unnamed-chunk-4-2.png

Generate TOM

softPower <- 14
SubGeneNames <- colnames(datExpr)
adj <- adjacency(datExpr, type = "signed", power = softPower)
diag(adj) <- 0
TOM <- TOMsimilarityFromExpr(datExpr, networkType = "signed", TOMType = "signed", power = softPower, maxPOutliers = 0.05)

TOM calculation: adjacency..
..will use 79 parallel threads.
 Fraction of slow calculations: 0.000000
..connectivity..
..matrix multiplication (system BLAS)..
..normalization..
..done.

colnames(TOM) <- rownames(TOM) <- SubGeneNames
dissTOM <- 1 - TOM
geneTree <- hclust(as.dist(dissTOM), method = "complete") # use complete for method rather than average (gives better results)
plot(geneTree, xlab = "", sub = "", cex = .5, main = "Gene clustering", hang = .001)

Past versions of unnamed-chunk-5-1.png

Identify Modules

minModuleSize <- 15
x <- 2

dynamicMods <- cutreeDynamic(
  dendro = geneTree, distM = as.matrix(dissTOM),
  method = "hybrid", pamStage = F, deepSplit = x,
  minClusterSize = minModuleSize
)

 ..cutHeight not given, setting it to 1  ===>  99% of the (truncated) height range in dendro.
 ..done.

dynamicColors <- labels2colors(dynamicMods) # label each module with a unique color
plotDendroAndColors(geneTree, dynamicColors, "Dynamic Tree Cut",
  dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05,
  main = "Gene dendrogram and module colors"
) # plot the modules with colors

Past versions of “Cut Gene Tree-1.png”

Calculate Eigengenes

MEs <- moduleEigengenes(datExpr, dynamicColors)$eigengenes 
ME1 <- MEs
row.names(ME1) <- row.names(datExpr)
MEDiss <- 1 - cor(MEs)
METree <- hclust(as.dist(MEDiss), method = "average")
plot(METree, main = "Clustering of module eigengenes", xlab = "", sub = "")
MEDissThres <- 0.2
abline(h = MEDissThres, col = "red")

Past versions of “merge close modules-1.png”

Merge Modules

merge <- mergeCloseModules(datExpr, dynamicColors, cutHeight = MEDissThres, verbose = 3)

 mergeCloseModules: Merging modules whose distance is less than 0.2
   multiSetMEs: Calculating module MEs.
     Working on set 1 ...
     moduleEigengenes: Calculating 30 module eigengenes in given set.
   Calculating new MEs...
   multiSetMEs: Calculating module MEs.
     Working on set 1 ...
     moduleEigengenes: Calculating 30 module eigengenes in given set.

mergedColors <- merge$colors
mergedMEs <- merge$newMEs
moduleColors <- mergedColors
MEs <- mergedMEs
modulekME <- signedKME(datExpr, MEs)

Plot merged modules

plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors),
  c("Dynamic Tree Cut", "Merged dynamic"),
  dendroLabels = FALSE, hang = 0.03,
  addGuide = TRUE, guideHang = 0.05
)

Past versions of unnamed-chunk-7-1.png

moduleColors <- mergedColors
MEs <- mergedMEs
modulekME <- signedKME(datExpr, MEs)

Generate function to look up genes in each module

modules <- MEs
c_modules <- data.frame(moduleColors)
row.names(c_modules) <- colnames(datExpr) 
module.list.set1 <- substring(colnames(modules), 3)
index.set1 <- 0
Network <- list()
for (i in 1:length(module.list.set1)) {
  index.set1 <- which(c_modules == module.list.set1[i])
  Network[[i]] <- row.names(c_modules)[index.set1]
}
names(Network) <- module.list.set1

Filter metadata table and correlate with eigengenes

nGenes <- ncol(datExpr)
nSamples <- nrow(datExpr)
MEs <- orderMEs(MEs)
MEs %>% select(-MEgrey) -> MEs
var <- model.matrix(~ 0 + ventric$trt)
moduleTraitCor <- cor(MEs, var, use = "p")
cor <- moduleTraitCor[abs(moduleTraitCor[, 1]) > .2, ]
moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nSamples)
cor <- melt(cor)
ggplot(cor, aes(Var2, Var1)) + geom_tile(aes(fill = value), colour = "white") + 
  scale_fill_gradient2( midpoint = 0, low = "blue", mid = "white",
  high = "red", space = "Lab", name = "Correlation \nStrength") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) + xlab("Treatment") + ylab(NULL)

Past versions of unnamed-chunk-9-1.png

Run linear regression to calculate sig diff mods

data <- data.frame(MEs,
  trt = ventric$trt,
  sample = as.factor(ventric$sample),
  batch = as.factor(ventric$batch)
)

mod <- lapply(colnames(MEs)[grepl("^ME", colnames(MEs))], function(me) {
  tryCatch({
    mod <- lmer(MEs[[me]] ~ trt + (1 | batch) + (1 | sample), data = data)
    pairwise <- emmeans(mod, pairwise ~ trt)
    plot <- data.frame(plot(pairwise, plotIt = F)$data)
    sig <- as.data.frame(pairwise$contrasts)
    return(sig)
  }, error = function(err) {
    print(err)
  })
})

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
control$checkConv, : Model failed to converge with max|grad| = 0.00943197
(tol = 0.002, component 1)

Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl =
control$checkConv, : Model failed to converge with max|grad| = 0.00238047
(tol = 0.002, component 1)

names(mod) <- colnames(MEs)[grepl("^ME", colnames(MEs))]
sig <- bind_rows(mod, .id="id")
write_csv(sig, path = here("output/glia/wgcna/wc_astro_linmod.csv"))


mod <- data.frame(unlist(mod))
mod %>%
  add_rownames("test") %>%
  separate(test, c("mod", "measure")) %>%
  dcast(measure ~ mod, value = unlist.mod.) %>%
  as.data.frame() %>%
  t() -> test

Warning: Deprecated, use tibble::rownames_to_column() instead.

Warning: Expected 2 pieces. Additional pieces discarded in 58 rows [5,
6, 11, 12, 17, 18, 23, 24, 29, 30, 35, 36, 41, 42, 47, 48, 53, 54, 59,
60, ...].

colnames(test) <- test[1, ]
data.frame(test) %>%
  add_rownames("mod") %>%
  slice(2:nrow(.)) %>%
  select(p, estimate, mod) %>%
  mutate(p = as.numeric(as.character(p)), estimate = as.numeric(as.character(estimate))) %>%
  arrange(log10(p)*abs(estimate)) %>%
  mutate(pubname = paste0("sc-astro-M", paste0(seq_len(length(estimate)))))  -> astro_mods

Warning: Deprecated, use tibble::rownames_to_column() instead.

astro_mods %>% filter(p < 0.05, abs(estimate) > 0.005) %>% pull(pubname) -> sigmods

rename_mods <- astro_mods$mod
names(rename_mods) <- astro_mods$pubname
MEs <- rename(MEs, !!!rename_mods)
data <- rename(data, !!!rename_mods)
write_csv(data.frame(pubname = names(rename_mods), color = as.character(rename_mods)),
          path = here("output/glia/wgcna/wc_astro_translatelabels.csv"))

hubgene information

hubgenes <- lapply(seq_len(length(Network)-1), function(x) {
  dat <- modulekME[Network[[x]], ]
  dat <- dat[order(-dat[paste0("kME", names(Network)[x])]), ]
  gene <- data.frame(gene=rownames(dat),kme=dat[,x]) 
  return(gene)
})

names(hubgenes) <-  names(rename_mods)[match(names(Network)[-length(Network)], 
                                             gsub(as.character(rename_mods),pattern = "ME",replacement = ""))]
d <- dplyr::bind_rows(hubgenes, .id = "id")

Warning in bind_rows_(x, .id): Unequal factor levels: coercing to character

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

Warning in bind_rows_(x, .id): binding character and factor vector,
coercing into character vector

write_csv(d, path = here("output/glia/wgcna/wc_astro_wgcna_genemods.csv"))

boxplot summary of diff exp modules

data <- melt(data, id.vars = c("trt", "sample"))

Warning: attributes are not identical across measure variables; they will
be dropped

data %>% filter(variable %in% sigmods) -> data
data$trt <- gsub(data$trt, pattern = "PF", replacement = "Veh")
boxplot <- ggplot(data = data, 
                  aes(x = fct_relevel(variable,"sc-astro-M1","sc-astro-M3",
                                      "sc-astro-M2","sc-astro-M9","sc-astro-M5","sc-astro-M6"), 
                      y = as.numeric(value))) +
  geom_boxplot(aes(fill=trt), width=0.5, notch=T, outlier.shape = NA) +
  theme_pubr(legend="right") + 
  geom_hline(yintercept = 0, linetype = "dashed") + xlab(NULL) +
  ylab("ME Expression") + 
  theme(axis.text.x = element_blank(), 
        axis.title.y = element_text(size=10, face="bold"),
        legend.text = element_text(size=8, face="bold"),
        legend.title = element_text(size=10, face="bold"), 
        plot.margin = unit(c(0, 0, 0, 0), "cm")) +
  ylim(c(-0.025,0.025)) + scale_fill_manual(name="Treatment", values = c("gray30", "gray80")) + theme_figure

boxplot

Warning: Removed 781 rows containing non-finite values (stat_boxplot).

Past versions of unnamed-chunk-12-1.png

go term enrichment

goterms <- lapply(hubgenes[sigmods], function(x) {
  x <- gprofiler(x,
    ordered_query = T, organism = "mmusculus", significant = T, custom_bg = colnames(datExpr),
    src_filter = c("GO:BP", "REAC", "KEGG"), hier_filtering = "strong",
    min_isect_size = 2,
    sort_by_structure = T, exclude_iea = T,
    min_set_size = 10, max_set_size = 300, correction_method = "fdr"
  )
  x <- x[order(x$p.value), ]
  return(x)
})

godat <- bind_rows(goterms, .id="id")
write_csv(godat, path = here("output/glia/wgcna/wc_astro_wgcna_goterms.csv"))

test overlap with nuclei modules

nuc_mods <- read_csv(file = here("output/glia/wgcna/astro_wgcna_genemodules.csv"))
nuc_mods %>%
  as.data.frame() %>%
  filter(id %in% c(paste0("sn-astro-M", seq_len(4)))) %>%
  dplyr::group_by(id) %>%
  dplyr::group_split() %>%
  map("gene") -> nuc_gene

nuc_mods %>% as.data.frame() %>%
  filter(id %in% c(paste0("sn-astro-M", seq_len(4)))) %>%
  dplyr::group_by(id) %>% group_keys() %>% pull() -> names(nuc_gene)

wc_nuc_overlap <- sapply(nuc_gene, function(x) {
  sapply(hubgenes[sigmods], function(y) {
    1 - phyper(sum(x %in% y$gene), length(y$gene), 5000 - length(y$gene), length(x), log.p = F)
  })
})


wc_nuc_overlap <- reshape2::melt(wc_nuc_overlap)
wc_nuc_overlap %>%
  mutate(value = p.adjust(wc_nuc_overlap$value, n = dim(wc_nuc_overlap)[1] * dim(wc_nuc_overlap)[2])) %>%
  mutate(sig = if_else(value > 0.05, "",
    if_else(.05 > value & value > .01, "*",
      if_else(.01 > value & value > .001, "**",
        if_else(.001 > value, "***", "")
      )
    )
  )) -> wc_nuc_overlap

write_csv(wc_nuc_overlap, path = here("output/glia/wgcna/astro_mod_overlaps_pval.csv"))

overlap <- ggplot(wc_nuc_overlap, aes(x = fct_relevel(Var1, "sc-astro-M1", "sc-astro-M3","sc-astro-M2","sc-astro-M9",
                                                      "sc-astro-M5","sc-astro-M6"),
                                      y = fct_relevel(Var2, "sn-astro-M4",
                                                      "sn-astro-M3","sn-astro-M2","sn-astro-M1"))) +
  geom_tile(size = 1, color = "white", fill="grey99") +
  geom_point(aes(size = if_else(-log10(value)<1.3,true = 0, false = -log10(value))), alpha=0.25) + 
  scale_size(name= expression(bold(-log[10] ~ pvalue)), breaks = c(2,4,8)) + 
  scale_fill_manual(values=c("black","red"), guide=F) + 
  theme_pubr(legend = "right") + xlab(NULL) + ylab(NULL) + 
  theme(axis.text.x = element_text(angle=45, hjust=1), 
        legend.title = element_text(size=8), 
        legend.text = element_text(size=7, face="bold"),
        plot.margin = unit(c(0, 0, 0, 0), "cm")) + theme_figure

overlap 

plot as grid

plot_grid(boxplot, overlap, align="v", axis="tblr", labels=c("f"), ncol=1)

Warning: Removed 781 rows containing non-finite values (stat_boxplot).

Past versions of unnamed-chunk-15-1.png

Plot gene networks

color <- c("sc-astro-M1","sc-astro-M3")
lapply(color, function(col) {
  maxsize <- 15
  hubs <- data.frame(genes=hubgenes[[col]]$gene[1:maxsize], kme = hubgenes[[col]]$kme[1:maxsize], mod = rep(col,15))
}) -> hub_plot 


hub_plot <- lapply(hub_plot, function(x) {
  adj[as.character(x$genes), as.character(x$genes)] %>% 
  graph.adjacency(mode = "undirected", weighted = T, diag = FALSE) %>%
  as_tbl_graph(g1) %>% upgrade_graph() %>% activate(nodes) %>% dplyr::mutate(mod=x$mod) %>%
  dplyr::mutate(kme=x$kme)  %>% 
    dplyr::mutate(color = ifelse(name %in% c("Slc1a2","Slc1a3","Slc6a11","Aqp4","Sparcl1","Sparc","Psap","Clu","Plat","Cpe"), yes="red",no="black")) %>%
  activate(edges)}
)

set.seed("139")
plotlist <- lapply(hub_plot, function(x) {
  print(ggraph(x,  layout = 'fr') + 
  geom_edge_link(color="darkgrey", aes(alpha = weight), show.legend = F)  + 
  scale_edge_width(range = c(0.2, .5)) + geom_node_text(aes(label = name, color=color), fontface="bold", size=3, repel=T) +
  scale_color_manual(values=c("gray0","red")) + guides(color=FALSE) +
  geom_node_point(shape=21, alpha=0.5, fill="grey70", size=4) + 
  theme_graph() + theme(legend.position = "none", plot.title = element_text(hjust=0.5, vjust=1)) + 
    coord_cartesian(clip = "off"))
})

arrange figure

genenet <- plot_grid(plotlist[[1]], plotlist[[2]], ncol=1, labels = c("f"), scale=0.8)
nets <- plot_grid(boxplot, overlap, align="v", axis="tblr", labels=c("e"), ncol=1, scale=0.95, rel_heights = c(2,1))

Warning: Removed 781 rows containing non-finite values (stat_boxplot).

wc_astro_res <- plot_grid(nets, genenet, rel_widths = c(1.5,1), nrow=1)

load(here("data/figures/fig4/fig4_nuc_top.RData"))
fig4 <- plot_grid(top_astro, wc_astro_res, ncol=1, align="hv", rel_heights = c(1,1.25), axis="lr")
transcriptomics_4 <- plot_grid(top_astro, fig4, rel_widths = c(1,1))
transcriptomics_4

Aqp4/Gfap Quantification

readxl::read_xlsx(here("data/mouse_data/fig4/191118_Aqp4_Gfap.xlsx"), range="A5:E9") %>%
  select(1,2,4,5) %>% 
  reshape2::melt() %>% 
  mutate(gene = c(rep("Aqp4", 8), rep("Gfap",8))) %>% 
  separate(variable, "[...]", into = "trt") -> aq_gf_quants

Warning: Expected 1 pieces. Additional pieces discarded in 16 rows [1, 2,
3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16].

aq_gf_quants %>% dplyr::group_by(trt,  gene) %>% dplyr::summarise(mean = mean(value), sd = sd(value), se = sd/sqrt(length(value))) %>%
  ggplot(aes(x=gene, y=mean, fill=fct_relevel(trt,"Veh","FGF1"), color = fct_relevel(trt,"Veh","FGF1"))) + 
  geom_col(width=0.9, alpha=0.75, colour="black", position="dodge") +
  geom_errorbar(aes(x=gene, ymin = mean-se, ymax=mean+se), width=0.2, position = position_dodge(0.9)) +
  geom_jitter(data = aq_gf_quants, inherit.aes = F, aes(x=gene, y=value, fill=fct_relevel(trt,"Veh","FGF1")), 
              alpha=0.85, shape=21, position = position_jitterdodge(0.2)) + xlab(NULL) + 
  ylab(expression(bold(Volume~(mu*m^3)))) + scale_fill_manual("Treatment", values=c("gray80","gray30")) + 
  scale_color_manual("Treatment", values=c("gray80","gray30")) +
  scale_y_continuous(labels = function(x) format(x, scientific = TRUE)) +
  geom_signif(y_position=c(1.5e6), xmin=c(0.9,1.9), xmax=c(1.1,2.1),
              annotation=c("*","**"), tip_length=0, size = 0.5, textsize = 5, color="black") +
  theme_pubr() + theme(legend.position="none", axis.text.x=element_text(angle=45, hjust=1)) + theme_figure -> aq_gf
aqp4_gfap <- cowplot::plot_grid(ggplot() + theme_void(), aq_gf, nrow=1, scale=c(1,0.75), labels="g", rel_widths = c(2.05,1))
aqp4_gfap

Agrp/Gfap overlap

readxl::read_xlsx(here("data/mouse_data/fig4/191118_GFAP_AGRP.xlsx"), range="A5:B9") %>%
  reshape2::melt() -> ga_quants
ga_quants %>% dplyr::group_by(variable) %>% dplyr::summarise(mean = mean(value), sd = sd(value), se = sd/sqrt(length(value))) %>%
  ggplot(aes(x = variable, y = mean, fill = variable, color = variable)) + 
  geom_col(width=1, alpha=0.75, colour="black", position="dodge") +
  geom_errorbar(aes(x=variable, ymin = mean-se, ymax=mean+se), width=0.2, position = position_dodge(0.9)) + 
  geom_jitter(data = ga_quants, inherit.aes = F, aes(x=variable, y=value, fill=variable), 
              alpha=0.5, shape=21, position = position_jitterdodge(.5)) + xlab(NULL) + 
  geom_signif(y_position=c(1000), xmin=c(1.2), xmax=c(1.8),
              annotation=c("**"), tip_length=0, size = 0.5, textsize = 5, color="black") +
  ylab(expression(bold(GFAP/Agrp~contact~volume~(mu*m^3/hpf)))) + scale_fill_manual(values=c("gray80","gray30"), name="") + 
  scale_color_manual(values=c("gray80","gray30"), name="") +
  theme_pubr() + theme(legend.position = "none", axis.text.x = element_text(angle=45,hjust=1)) + theme_figure -> gfap_agrp
gfap_agrp <- cowplot::plot_grid(ggplot() + theme_void(), gfap_agrp, nrow=1, scale=c(1,0.75), labels="h", rel_widths = c(1.5,1))
gfap_agrp

EM of PAPs quantification

readxl::read_xlsx(here("data/mouse_data/fig4/EM_glial_contact.xlsx"), range="A4:B8") %>%
  reshape2::melt() -> em_quants
em_quants %>% dplyr::group_by(variable) %>% dplyr::summarise(mean = mean(value), sd = sd(value), se = sd/sqrt(length(value))) %>%
  ggplot(aes(x=variable, y=mean, fill=variable, color = variable)) +
  geom_col(width=1, alpha=0.75, colour="black", position="dodge") +
  geom_errorbar(aes(x=variable, ymin = mean-se, ymax=mean+se), width=0.2, position = position_dodge(0.9)) +
  geom_jitter(data = em_quants, inherit.aes = F, aes(x=variable, y=value, fill=variable), 
              alpha=0.5, shape=21, position = position_jitterdodge(.5)) + xlab(NULL) + 
  geom_signif(y_position=c(85), xmin=c(1.2), xmax=c(1.8),
              annotation=c("***"), tip_length=0, size = 0.5, textsize = 5, color="black") +
  ylab("% Astrocytic coverage") + scale_fill_manual("Treatment", values=c("gray80","gray30")) + 
  scale_color_manual("Treatment", values=c("gray80","gray30")) +
  theme_pubr() + theme(legend.position = "none",axis.text.x = element_text(angle=45, hjust=1)) + 
  coord_cartesian(clip="off") + theme_figure -> em_bp
pap <- cowplot::plot_grid(ggplot() + theme_void(), em_bp, nrow=1, scale=c(1,0.75), labels="i", rel_widths = c(1.5,1))
pap

Generate full figure

validation <- plot_grid(aqp4_gfap, gfap_agrp, pap, nrow=1, align="h", axis="tb", rel_widths = c(1.5,1,1))
plot_grid(transcriptomics_4, validation, ncol=1, align="h", axis="tblr", rel_heights = c(1.75,1))

ggsave(filename = here("data/figures/fig4/fig4.tiff"), width = 13, h=10, compression="lzw")

Session information

sessionInfo()

R version 3.5.3 (2019-03-11)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Storage

Matrix products: default
BLAS/LAPACK: /usr/lib64/libopenblas-r0.3.3.so

locale:
 [1] LC_CTYPE=en_DK.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_DK.UTF-8        LC_COLLATE=en_DK.UTF-8    
 [5] LC_MONETARY=en_DK.UTF-8    LC_MESSAGES=en_DK.UTF-8   
 [7] LC_PAPER=en_DK.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_DK.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ggExtra_0.9           ggraph_1.0.2          tidygraph_1.1.2      
 [4] here_0.1              cowplot_1.0.0         ggpubr_0.2.1         
 [7] magrittr_1.5          gProfileR_0.6.7       igraph_1.2.4.1       
[10] reshape2_1.4.3        forcats_0.4.0         stringr_1.4.0        
[13] dplyr_0.8.3           purrr_0.3.2           readr_1.3.1.9000     
[16] tidyr_0.8.3           tibble_2.1.3          tidyverse_1.2.1      
[19] genefilter_1.64.0     ggbeeswarm_0.6.0      ggplot2_3.2.1        
[22] lme4_1.1-21           Matrix_1.2-17         emmeans_1.3.5.1      
[25] ggsci_2.9             parallelDist_0.2.4    cluster_2.1.0        
[28] WGCNA_1.68            fastcluster_1.1.25    dynamicTreeCut_1.63-1
[31] Seurat_3.0.3.9036    

loaded via a namespace (and not attached):
  [1] estimability_1.3      R.methodsS3_1.7.1     coda_0.19-3          
  [4] acepack_1.4.1         bit64_0.9-7           knitr_1.23           
  [7] irlba_2.3.3           multcomp_1.4-10       R.utils_2.9.0        
 [10] data.table_1.12.2     rpart_4.1-15          RCurl_1.95-4.12      
 [13] doParallel_1.0.14     generics_0.0.2        metap_1.1            
 [16] BiocGenerics_0.28.0   preprocessCore_1.44.0 TH.data_1.0-10       
 [19] RSQLite_2.1.1         RANN_2.6.1            future_1.14.0        
 [22] bit_1.1-14            xml2_1.2.0            lubridate_1.7.4      
 [25] httpuv_1.5.1          assertthat_0.2.1      viridis_0.5.1        
 [28] xfun_0.8              hms_0.5.0             promises_1.0.1       
 [31] evaluate_0.14         DEoptimR_1.0-8        caTools_1.17.1.2     
 [34] readxl_1.3.1          DBI_1.0.0             htmlwidgets_1.3      
 [37] stats4_3.5.3          backports_1.1.4       annotate_1.60.1      
 [40] gbRd_0.4-11           RcppParallel_4.4.3    vctrs_0.2.0          
 [43] Biobase_2.42.0        ROCR_1.0-7            withr_2.1.2          
 [46] ggforce_0.3.0.9000    robustbase_0.93-5     checkmate_1.9.4      
 [49] sctransform_0.2.0     ape_5.3               lazyeval_0.2.2       
 [52] crayon_1.3.4          labeling_0.3          pkgconfig_2.0.2      
 [55] tweenr_1.0.1          nlme_3.1-140          vipor_0.4.5          
 [58] rematch_1.0.1         nnet_7.3-12           rlang_0.4.0          
 [61] globals_0.12.4        miniUI_0.1.1.1        sandwich_2.5-1       
 [64] modelr_0.1.4          rsvd_1.0.2            cellranger_1.1.0     
 [67] rprojroot_1.3-2       polyclip_1.10-0       matrixStats_0.54.0   
 [70] lmtest_0.9-37         boot_1.3-22           zoo_1.8-6            
 [73] base64enc_0.1-3       beeswarm_0.2.3        whisker_0.3-2        
 [76] ggridges_0.5.1        png_0.1-7             viridisLite_0.3.0    
 [79] bitops_1.0-6          R.oo_1.22.0           KernSmooth_2.23-15   
 [82] blob_1.1.1            workflowr_1.4.0       robust_0.4-18.1      
 [85] S4Vectors_0.20.1      ggsignif_0.5.0        scales_1.0.0         
 [88] memoise_1.1.0         plyr_1.8.4            ica_1.0-2            
 [91] gplots_3.0.1.1        bibtex_0.4.2          gdata_2.18.0         
 [94] compiler_3.5.3        lsei_1.2-0            RColorBrewer_1.1-2   
 [97] rrcov_1.4-7           fitdistrplus_1.0-14   cli_1.1.0            
[100] listenv_0.7.0         pbapply_1.4-1         htmlTable_1.13.1     
[103] Formula_1.2-3         MASS_7.3-51.4         tidyselect_0.2.5     
[106] stringi_1.4.3         highr_0.8             yaml_2.2.0           
[109] latticeExtra_0.6-28   ggrepel_0.8.0.9000    grid_3.5.3           
[112] tools_3.5.3           future.apply_1.3.0    parallel_3.5.3       
[115] rstudioapi_0.10       foreach_1.4.4         foreign_0.8-71       
[118] git2r_0.25.2          gridExtra_2.3         farver_1.1.0         
[121] Rtsne_0.15            digest_0.6.20         shiny_1.3.2          
[124] Rcpp_1.0.2            broom_0.5.2           SDMTools_1.1-221.1   
[127] later_0.8.0           RcppAnnoy_0.0.12      httr_1.4.1           
[130] AnnotationDbi_1.44.0  npsurv_0.4-0          Rdpack_0.11-0        
[133] colorspace_1.4-1      rvest_0.3.4           XML_3.98-1.20        
[136] fs_1.3.1              reticulate_1.13       IRanges_2.16.0       
[139] splines_3.5.3         uwot_0.1.3            plotly_4.9.0         
[142] fit.models_0.5-14     xtable_1.8-4          jsonlite_1.6         
[145] nloptr_1.2.1          zeallot_0.1.0         R6_2.4.0             
[148] Hmisc_4.2-0           mime_0.7              pillar_1.4.2         
[151] htmltools_0.3.6       glue_1.3.1            minqa_1.2.4          
[154] codetools_0.2-16      tsne_0.1-3            pcaPP_1.9-73         
[157] mvtnorm_1.0-11        lattice_0.20-38       leiden_0.3.1         
[160] gtools_3.8.1          GO.db_3.7.0           survival_2.44-1.1    
[163] rmarkdown_1.13        munsell_0.5.0         iterators_1.0.10     
[166] impute_1.56.0         haven_2.1.0           gtable_0.3.0