diff --git a/mid_res/nat_adj/1.calc.lcpm.baselineH1_as03_modelgenes.r b/mid_res/nat_adj/1.calc.lcpm.baselineH1_as03_modelgenes.r deleted file mode 100644 index 7fcf577..0000000 --- a/mid_res/nat_adj/1.calc.lcpm.baselineH1_as03_modelgenes.r +++ /dev/null @@ -1,54 +0,0 @@ -# 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')) diff --git a/mid_res/nat_adj/2.AS03.innatesigs.vand.validationcohort.r b/mid_res/nat_adj/2.AS03.innatesigs.vand.validationcohort.r deleted file mode 100644 index eb32678..0000000 --- a/mid_res/nat_adj/2.AS03.innatesigs.vand.validationcohort.r +++ /dev/null @@ -1,41 +0,0 @@ -# R version 4.0.5 -suppressMessages(library(here)) -suppressMessages(library(tidyverse)) -suppressMessages(library(scglmmr)) -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/") - -# define adjuvant signatures -# leading edge combined signature -# upregulated genes only -# mono and mDC specific and combined genes -li.up = readRDS(file = here('mid_res/3_H1N1_H5N1_joint_pseudobulk_DE/dataV4/gsea/li.up.rds')) -li2.up = readRDS(file = here('mid_res/3_H1N1_H5N1_joint_pseudobulk_DE/dataV4/gsea/li2.up.rds')) -li.full = unlist(li.up,use.names = FALSE, recursive = TRUE) %>% unique -mono.genes = li.up$CD14_Mono %>% unlist() %>% unique() -mdc.genes = li2.up$mDC %>% unlist() %>% unique() -as03.sig = list('AS03_signature' = li.full, "AS03_Monocyte" = mono.genes, 'AS03_mDC' = mdc.genes) - -# vand validation 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') -gvand = FgseaList(rank.list.celltype = vand.rank, pathways = as03.sig) - -# leading edge from each cell tpe -dc.va.le = gvand$DNC %>% filter(pathway == 'AS03_mDC') %$% leadingEdge %>% unlist() -saveRDS(dc.va.le,file = paste0(datapath, 'dc.va.le.rds')) - -mono.va.le = gvand$MNC %>% filter(pathway == 'AS03_Monocyte') %$% leadingEdge %>% unlist() -saveRDS(mono.va.le, file = paste0(datapath, 'mono.va.le.rds')) - - -#mono -enrline = list(geom_line(color = "deepskyblue3", size = 2 )) -p = fgsea::plotEnrichment(pathway = as03.sig$AS03_Monocyte, stats = vand.rank$MNC) + enrline -# mDC -p = fgsea::plotEnrichment(pathway = as03.sig$AS03_mDC, stats = vand.rank$DNC) + enrline -ggsave(p, filename = paste0(figpath, 'dc.vand.enr.pdf'), width = 4, height = 2.5) diff --git a/mid_res/nat_adj/3.natural.adjuvant.signatures.r b/mid_res/nat_adj/3.natural.adjuvant.signatures.r deleted file mode 100644 index ecb8078..0000000 --- a/mid_res/nat_adj/3.natural.adjuvant.signatures.r +++ /dev/null @@ -1,150 +0,0 @@ -# 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/") -datapath = here("mid_res/nat_adj/generated_data/V4/") - -# set plot themes to distinguisn between groups / signatures being tested -# AS03 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() -cub = c('dodgerblue', 'red') -baselinetheme = list( - mtheme, - theme(axis.text.x=element_text(angle = 0, hjust = 0.5, color = 'black', size = 10), - axis.title.y = element_text(size = 14)), - geom_boxplot(show.legend = FALSE, outlier.shape = 21), - scale_y_continuous( breaks= scales::pretty_breaks(), expand = c(0.15,0)), - scale_x_discrete(labels = c("low", "high")), - xlab("Antibody Response"), - scale_fill_manual(values = cua), - scale_color_manual(values = cub) -) - -# baseline bulk theme -cu1 = sapply(c('dodgerblue','grey', 'red'), col.alpha, 0.2) %>% unname() -cu2 = c('dodgerblue','grey', 'red') -bulktheme = list( - mtheme, - theme(axis.text.x=element_text(angle = 0, hjust = 0.5, color = 'black', size = 10), - axis.title.y = element_text(size = 14)), - geom_boxplot(show.legend = FALSE, outlier.shape = 21), - scale_y_continuous( breaks= scales::pretty_breaks(), expand = c(0.15,0)), - scale_x_discrete(labels = c("low", "mid", "high")), - xlab("Antibody Response"), - scale_fill_manual(values = cu1), - scale_color_manual(values = cu2) -) -my_compare = list(l1 = c("1", "0"), l2 = c("2", "0")) - -# define adjuvant signatures -# leading edge combined signature -# upregulated gnees only -# mono and mDC specific and combined genes -li.up = readRDS(file = here('mid_res/3_H1N1_H5N1_joint_pseudobulk_DE/dataV4/gsea/li.up.rds')) -li2.up = readRDS(file = here('mid_res/3_H1N1_H5N1_joint_pseudobulk_DE/dataV4/gsea/li2.up.rds')) -li.full = unlist(li.up,use.names = FALSE, recursive = TRUE) %>% unique -mono.genes = li.up$CD14_Mono %>% unlist() %>% unique() -mdc.genes = li2.up$mDC %>% unlist() %>% unique() -as03.sig = list('AS03_signature' = li.full, "AS03_Monocyte" = mono.genes, 'AS03_mDC' = mdc.genes) - - -# leading edge signatures from validation cohort -mono.va.le = readRDS(file = here('mid_res/nat_adj/generated_data/V4/mono.va.le.rds')) -mono.va.le = list('AS03_Monocyte_LE' = mono.va.le) -dc.va.le = readRDS(file = here('mid_res/nat_adj/generated_data/V4/dc.va.le.rds')) -dc.va.le = list('AS03_mDC_LE' = dc.va.le) -as03.sig = c(as03.sig, mono.va.le, dc.va.le) - - -# 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')) - -# load average baseline high and low responders from unadjuvanted cohort -av0 = readRDS(file = here('mid_res/nat_adj/generated_data/V4/av0.rds')) - -# load bulk microarray data and process to baseline samples -array = data.table::fread("data/CHI_H1N1_data/microarray/CHI_GE_matrix_gene.txt", data.table = F) %>% - tibble::remove_rownames() %>% - tibble::column_to_rownames("gene") %>% - select(., matches("day0")) %>% - select(-matches("day70")) %>% - select(-matches("pre")) %>% - select(-matches("day1")) - -############################# -# CD14 Monocytes -############################# -# Monocyte combined AS03 signature average across time between groups -#natural adjuvant test -mono.sig.av2 = - av0$CD14_Mono %>% - filter(gene %in% as03.sig$AS03_Monocyte_LE) %>% - group_by(sample, response) %>% - summarize(meansig = mean(count)) -#plot -p = ggplot(mono.sig.av2, aes(x = response, y = meansig, fill = response , color = response)) + - baselinetheme + - ggpubr::stat_compare_means(method = 'wilcox',paired = FALSE, show.legend = FALSE) + - ylab('Monocyte AS03 Adjuvant Signature') + - ggtitle('Baseline: CD14 Monocytes\nunadjuvanted subjects') -p -ggsave(p, filename = paste0(figpath, 'baseline_monosig_monocytes_LeadingEdgeVand.pdf'), width = 2.7, height = 4) - -############################# -# mDC -############################# - -# mdc Combined AS03 signature average across time between groups -mdc.sig.av = - av_tidy$mDC %>% - filter(gene %in% mdc.genes) %>% - 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.pdf'), width = 1.9, height = 3) - - -# natural adjuvant test - mDC -mdc.sig.av2 = - av0$mDC %>% - filter(gene %in% as03.sig$AS03_mDC) %>% - group_by(sample, response) %>% - summarize(meansig = mean(count)) -#plot -p = ggplot(mdc.sig.av2, aes(x = response, y = meansig, fill = response , color = response)) + - baselinetheme + - ggpubr::stat_compare_means(method = 'wilcox',paired = FALSE, show.legend = FALSE) + - ylab('mDC AS03 Adjuvant Signature') + - ggtitle('Baseline: mDCs\nunadjuvanted subjects') -p -ggsave(p, filename = paste0(figpath, 'baseline_mDC_mDC_LeadingEdgeVand.pdf'), width = 2.7, height = 4) - - -for (i in 1:length(as03.sig)) { - data.table::fwrite(as03.sig[i],file = paste0(datapath, names(as03.sig)[i],'.txt')) -} diff --git a/mid_res/nat_adj/figures/V4/as03_mDC_sig.pdf b/mid_res/nat_adj/figures/V4/as03_mDC_sig.pdf deleted file mode 100644 index 3c43872..0000000 Binary files a/mid_res/nat_adj/figures/V4/as03_mDC_sig.pdf and /dev/null differ diff --git a/mid_res/nat_adj/figures/V4/as03_mono_sig.pdf b/mid_res/nat_adj/figures/V4/as03_mono_sig.pdf deleted file mode 100644 index c771830..0000000 Binary files a/mid_res/nat_adj/figures/V4/as03_mono_sig.pdf and /dev/null differ diff --git a/mid_res/nat_adj/figures/V4/baseline_mDC_mDC_LeadingEdgeVand.pdf b/mid_res/nat_adj/figures/V4/baseline_mDC_mDC_LeadingEdgeVand.pdf deleted file mode 100644 index d44fbb5..0000000 Binary files a/mid_res/nat_adj/figures/V4/baseline_mDC_mDC_LeadingEdgeVand.pdf and /dev/null differ diff --git a/mid_res/nat_adj/figures/V4/baseline_monosig_monocytes_LeadingEdgeVand.pdf b/mid_res/nat_adj/figures/V4/baseline_monosig_monocytes_LeadingEdgeVand.pdf deleted file mode 100644 index 9e7650f..0000000 Binary files a/mid_res/nat_adj/figures/V4/baseline_monosig_monocytes_LeadingEdgeVand.pdf and /dev/null differ diff --git a/mid_res/nat_adj/figures/V4/dc.vand.enr.pdf b/mid_res/nat_adj/figures/V4/dc.vand.enr.pdf deleted file mode 100644 index 15ddb15..0000000 Binary files a/mid_res/nat_adj/figures/V4/dc.vand.enr.pdf and /dev/null differ diff --git a/mid_res/nat_adj/figures/V4/mono.vand.enr.pdf b/mid_res/nat_adj/figures/V4/mono.vand.enr.pdf deleted file mode 100644 index abaf7f2..0000000 Binary files a/mid_res/nat_adj/figures/V4/mono.vand.enr.pdf and /dev/null differ