merge commit 270 9a4069a

mid_res/nat_adj/1.calc.lcpm.baselineH1_as03_modelgenes.r

@@ -0,0 +1,54 @@
+# R version 4.0.5 
+suppressMessages(library(here))
+suppressMessages(library(tidyverse))
+suppressMessages(library(scglmmr))
+source(here('functions/MattPMutils.r'))
+
+# set save paths 
+figpath = here("mid_res/nat_adj/figures/V4/"); dir.create(figpath)
+datapath = here("mid_res/nat_adj/generated_data//V4/"); dir.create(datapath)
+
+# define high responders 
+high.responders = c("205","207","209","212","215","234","237","245","250","256")
+
+# read pb data, subset to day 0 non adj, subset out day 0 metadata. 
+pb = readRDS(file = here('mid_res/variance_partition/generated_data/pb_vp.rds'))
+pb = lapply(pb, function(x) x = x[ ,1:40])
+cnames = gsub("~.*","",colnames(pb[[1]]))
+pb = lapply(pb, function(x){
+  x %>% as.data.frame() %>% setNames(nm = cnames) %>% as.matrix() 
+})
+d0 = readRDS(file = here('mid_res/1_H1N1_pseudobulk_DE/dataV3/d0.rds'))
+d0d = lapply(pb, function(x){ x = x[ , rownames(d0)]})
+
+# make a list of genes indxed by celltype for genes to fit from H5 model 
+av_tidy = readRDS(file = here('mid_res/3_H1N1_H5N1_joint_pseudobulk_DE/dataV4/git_ignore/av_tidy.rds'))
+genes.test = lapply(av_tidy , function(x) unique(x$gene))
+
+# check names are the same of indexed genes and average data
+stopifnot(isTRUE(all.equal(names(d0d), names(genes.test))))
+
+# calcualte log CPM of the baseline pseudobulk data
+av = list()
+for (i in 1:length(d0d)) {
+  dge =  edgeR::DGEList( d0d[[i]] ) 
+  dge = dge[genes.test[[i]], ]
+  av[[i]]  = edgeR::cpm(dge, log = TRUE)
+}
+names(av) = names(d0d)
+saveRDS(av,file = paste0(datapath,'av.rds'))
+
+# tidy aggregated data 
+av0 = list()
+for (i in 1:length(pb)) {
+  ct = names(av)[i]
+  gs = rownames(av[[i]])
+  av0[[i]] = GetTidySummary(
+    av.exprs.list = av,
+    celltype.index = i,
+    genes.use = gs) %>%
+    mutate(response = if_else(str_sub(sample, 1, 3) %in%  high.responders, 'High', "Low")) %>%
+    mutate(response = factor(response, levels = c('Low', 'High')))
+}
+names(av0) = names(pb)
+saveRDS(av0, file = paste0(datapath,'av0.rds'))

mid_res/nat_adj/2.AS03.innatesigs.vand.validationcohort.r

@@ -0,0 +1,67 @@
+suppressMessages(library(here))
+suppressMessages(library(tidyverse))
+source(here('functions/scglmmr.functions.R'))
+suppressMessages(library(magrittr))
+source(here('functions/MattPMutils.r'))
+
+# set save paths 
+figpath = here("mid_res/nat_adj/figures/V4/")
+datapath = here("mid_res/nat_adj/generated_data//V4/")
+
+# load mono and mDC AS03 specific signatures 
+mono.genes = readRDS(file = here('mid_res/combined_contrast/generated_data/mo.noifn.rds'))
+mdc.genes = readRDS(file = here('mid_res/combined_contrast/generated_data/dc.noifn.rds'))
+
+# extract leading edge into signatures 
+mono.genes = mono.genes$leadingEdge %>% unlist() %>% unique()
+mdc.genes = mdc.genes$leadingEdge %>% unlist() %>% unique()
+li.full = c(mono.genes, mdc.genes) %>% unique()
+as03.sig = list('AS03_signature' = li.full, "AS03_Monocyte" = mono.genes, 'AS03_mDC' = mdc.genes)
+as03.sig.list = as03.sig
+saveRDS(as03.sig.list, file = paste0(datapath, 'as03.sig.list.rds'))
+
+# validation cohort for the highlighted signatures
+vand.fit = readRDS(file = here('mid_res/vand/generated_data/fit1.rds'))
+vand.rank = ExtractResult(model.fit.list = vand.fit,what = 'lmer.z.ranks', coefficient.number = 1, coef.name = 'delta')
+
+# Enrichment of AS03 signatures without ifn signatures 
+gvand = FgseaList(rank.list.celltype = vand.rank, pathways = as03.sig)
+saveRDS(gvand,file = paste0(datapath,'gvand.rds'))
+
+gvand$MNC
+# pathway         pval         padj  log2err        ES      NES size                                  leadingEdge celltype
+# 1: AS03_signature 1.910373e-34 5.731119e-34 1.529705 0.5714156 3.144152  263 SERPINA1,SECTM1,HCK,PLSCR1,LILRA1,FCGR1B,...      MNC
+# 2:  AS03_Monocyte 4.563863e-26 6.845794e-26 1.326716 0.5979691 3.088277  170 SERPINA1,SECTM1,HCK,PLSCR1,LILRA1,FCGR2A,...      MNC
+# 3:       AS03_mDC 6.087321e-21 6.087321e-21 1.186651 0.5916694 2.984150  145 SERPINA1,SECTM1,HCK,LILRA1,FCGR1B,FCGR2A,...      MNC
+
+gvand$DNC
+#       pathway         pval         padj   log2err        ES      NES size                                 leadingEdge celltype
+# 1: AS03_signature 4.154660e-17 1.246398e-16 1.0672100 0.4932650 2.402933  257 PSME2,FPR1,SLC31A2,PSMB10,PSME1,ALDH1A1,...      DNC
+# 2:       AS03_mDC 1.761296e-12 2.641944e-12 0.9101197 0.5311486 2.420236  146 PSME2,FPR1,SLC31A2,PSMB10,PSME1,ALDH1A1,...      DNC
+# 3:  AS03_Monocyte 6.171525e-10 6.171525e-10 0.8012156 0.4774659 2.201396  163   FPR1,SLC31A2,DPYD,IFIH1,LILRB3,PLSCR1,...      DNC
+
+# leading edge from each innate cell types - AS03 signature 
+dc.as03.sig.validated = gvand$DNC %>%  filter(pathway == 'AS03_mDC') %$% leadingEdge %>%  unlist()
+mono.as03.sig.validated = gvand$MNC %>%  filter(pathway == 'AS03_Monocyte') %$% leadingEdge %>%  unlist() 
+saveRDS(dc.as03.sig.validated, file = paste0(datapath,'dc.as03.sig.validated.rds'))
+saveRDS(mono.as03.sig.validated, file = paste0(datapath,'mono.as03.sig.validated.rds'))
+
+data.table::fwrite(list(dc.as03.sig.validated),file = paste0(datapath,'dc.as03.sig.validated.txt'),sep = '\t')
+data.table::fwrite(list(mono.as03.sig.validated),file = paste0(datapath,'mono.as03.sig.validated.txt'),sep = '\t')
+
+# plot enrichment distributions 
+enrline = list(geom_line(color = "deepskyblue3", size = 2 ))
+#mono
+p = fgsea::plotEnrichment(pathway = as03.sig$AS03_Monocyte, stats = vand.rank$MNC) + enrline
+ggsave(p, filename = paste0(figpath, 'mono.vand.enr.2.pdf'), width = 5, height = 3)
+# mDC
+p = fgsea::plotEnrichment(pathway = as03.sig$AS03_mDC, stats = vand.rank$DNC) + enrline
+ggsave(p, filename = paste0(figpath, 'dc.vand.enr.2.pdf'), width = 5, height = 3)
+
+
+
+
+
+
+
+

mid_res/nat_adj/3.natural.adjuvant.signatures.r

@@ -0,0 +1,112 @@
+suppressMessages(library(here))
+suppressMessages(library(tidyverse))
+source(here('functions/scglmmr.functions.R'))
+source(here('functions/MattPMutils.r'))
+set.seed(1990)
+# set save paths 
+figpath = here("mid_res/nat_adj/figures/V4/")
+datapath = here("mid_res/nat_adj/generated_data/V4/")
+
+
+# set theme 
+cu = c("grey48", "grey", "grey48",  "deepskyblue3")
+cu.alpha = sapply(cu, col.alpha, alpha = 0.4) %>% unname()
+mtheme = list(
+  geom_boxplot(show.legend = FALSE, outlier.shape = NA),
+  theme_bw(base_size = 10.5),
+  theme(axis.text.x=element_text(angle = -90, hjust = 0)),
+  theme(strip.background = element_blank(),
+        panel.grid.major = element_blank(),
+        panel.grid.minor = element_blank(),
+        strip.text = element_text(face = "bold",family = "Helvetica"),
+        axis.text.y =  element_text(size = 6),
+        axis.title.y = element_text(size = 10))
+  )
+cua = sapply(c('dodgerblue', 'red'), col.alpha, 0.2) %>% unname()
+
+
+# load average day 1 comparison cohort data 
+av_tidy = readRDS(file = here('mid_res/3_H1N1_H5N1_joint_pseudobulk_DE/dataV4/gene_dist/av_tidy.rds'))
+
+# AS03 adjuvant signatures (no ifn)
+as03.sig.list = readRDS(file = here('mid_res/nat_adj/generated_data/V4/as03.sig.list.rds'))
+
+# mdc Combined AS03 signature average across time between groups 
+mdc.sig.av = 
+  av_tidy$mDC %>% 
+  filter(gene %in% as03.sig.list$AS03_mDC) %>% 
+  group_by(sample, group) %>% 
+  summarize(meansig = mean(count))
+#plot
+p = ggplot(mdc.sig.av, aes(x = group, y = meansig, fill = group , color = group)) +
+  mtheme + 
+  theme(axis.title.x = element_blank()) +
+  ylab('mDC AS03 Adjuvant Signature') +
+  scale_color_manual(values = cu) +
+  scale_fill_manual(values = cu.alpha) + 
+  ggtitle('mDC')
+ggsave(p,filename = paste0(figpath, 'as03_mDC_sig.2.pdf'), width = 1.9, height = 3)
+
+# monocyte Combined AS03 signature average across time between groups 
+mono.sig.av = 
+  av_tidy$CD14_Mono %>% 
+  filter(gene %in% as03.sig.list$AS03_Monocyte) %>% 
+  group_by(sample, group) %>% 
+  summarize(meansig = mean(count))
+#plot
+p = ggplot(mono.sig.av, aes(x = group, y = meansig, fill = group , color = group)) +
+  mtheme + 
+  theme(axis.title.x = element_blank()) +
+  ylab('CD14 Mono AS03 Adjuvant Signature') +
+  scale_color_manual(values = cu) +
+  scale_fill_manual(values = cu.alpha) + 
+  ggtitle('CD14 Monocytes')
+ggsave(p,filename = paste0(figpath, 'as03_mono_sig.2.pdf'), width = 1.9, height = 3)
+
+
+##########################
+## GSEA of NA signatures 
+##########################
+# enrichment of validated signatures (Vand cohort) in baseline high vs low responders  
+mono.as03.sig.validated = readRDS(file = here('mid_res/nat_adj/generated_data/V4/mono.as03.sig.validated.rds'))
+dc.as03.sig.validated = readRDS(file = here('mid_res/nat_adj/generated_data/V4/dc.as03.sig.validated.rds'))
+
+#  high vs low model gene ranks within mono and mDC age and sex adjusted 
+cont0 = readRDS(file = here('mid_res/1_H1N1_pseudobulk_DE/dataV4/cont0.rds'))
+r0 = ExtractResult(model.fit.list = cont0, what = 'gene.t.ranks',coefficient.number = 1, coef.name = 'adjmfc')
+
+# enrichment 
+mono.na.gsea = fgsea::fgsea(pathways =  list('AS03.mono' = mono.as03.sig.validated), stats = r0$CD14_Mono)
+# pathway         pval         padj  log2err        ES     NES size                                    leadingEdge
+# 1: AS03.mono 7.515986e-13 7.515986e-13 0.921426 0.6263153 2.74048   78 S100A11,S100A12,APOBEC3A,DDX60,DDX58,CREB5,...
+
+p = fgsea::plotEnrichment(pathway = mono.as03.sig.validated, stats = r0$CD14_Mono) +
+  geom_line(size = 1.5, color = 'red') + 
+  theme(plot.title = element_text(size = 10))
+p
+ggsave(p, filename =paste0(figpath, 'mono.natadj.gsea.2.pdf'), width = 5, height = 3)
+
+mono.na.gsea$leadingEdge
+# [1] "S100A11"  "S100A12"  "APOBEC3A" "DDX60"    "DDX58"    "CREB5"    "TFEC"     "S100A9"   "GCA"      "FGL2"     "NACC2"   
+# [12] "KYNU"     "SLC16A3"  "MS4A4A"   "S100A8"   "FCGR3A"   "SAMHD1"   "TNFSF13B" "C19orf59" "CCR2"     "MARCO"    "P2RY13"  
+# [23] "RSAD2"    "SERPINA1" "FPR1"     "FGR"      "PLSCR1"   "SIGLEC9"  "IFIH1"    "ACSL1"    "LMNB1"    "LRRK2"    "MNDA"    
+# [34] "PLBD1"    "KIAA0513" "AQP9"     "SLC31A2"  "LILRB1"   "VCAN"    
+
+data.table::fwrite(mono.na.gsea$leadingEdge,file = paste0(datapath,'mono.na.gsea.leadingEdge.txt'))
+
+
+mdc.na.gsea = fgsea::fgsea(pathways = list('AS03.mdc' = dc.as03.sig.validated), stats = r0$mDC)
+# pathway       pval       padj   log2err        ES      NES size                                  leadingEdge
+# 1: AS03.mdc 0.02383525 0.02383525 0.3524879 0.3713362 1.521304   58 S100A8,S100A9,PSMB6,SERPINA1,RB1,SLC31A2,...
+p = fgsea::plotEnrichment(pathway = dc.as03.sig.validated,stats = r0$mDC)  +
+  geom_line(size = 1.5, color = 'red') + 
+  theme(plot.title = element_text(size = 10))
+p
+ggsave(p, filename =paste0(figpath, 'mdc.natadj.gsea.pdf'), width = 5, height = 3)
+
+mdc.na.gsea$leadingEdge
+# [1] "S100A8"   "S100A9"   "PSMB6"    "SERPINA1" "RB1"      "SLC31A2"  "TYMP"     "S100A11"  "MS4A4A"   "FCN1"     "LILRB2"  
+# [12] "PSMA7"    "CDC26"    "RBX1"     "PSMB3"    "PLBD1"    "LMNB1"    "PPP2R5E"  "KYNU"    
+
+data.table::fwrite(mdc.na.gsea$leadingEdge,file = paste0(datapath,'mdc.na.gsea.leadingEdge.txt'))
+