merge commit 270 9a4069a

niaid · Nov 9, 2023 · 70b2de1 · 70b2de1
1 parent 979d5d7
commit 70b2de1
Show file tree

Hide file tree

Showing 8 changed files with 483 additions and 0 deletions.
diff --git a/mid_res/mrna/figures/sig.emm.pdf b/mid_res/mrna/figures/sig.emm.pdf
diff --git a/mid_res/mrna/figures/sig.pdf b/mid_res/mrna/figures/sig.pdf
diff --git a/mid_res/mrna/figures/sig_mDC.emm.pdf b/mid_res/mrna/figures/sig_mDC.emm.pdf
diff --git a/mid_res/mrna/figures/sig_mDC.pdf b/mid_res/mrna/figures/sig_mDC.pdf
diff --git a/mid_res/mrna/mrna_1_setup.r b/mid_res/mrna/mrna_1_setup.r
@@ -0,0 +1,94 @@
+# R 4.0.5
+# testing out mrna data 
+suppressMessages(library(Seurat))
+suppressMessages(library(tidyverse))
+suppressMessages(library(here))
+suppressMessages(library(Matrix))
+
+# save paths 
+datapath = file.path(here('mid_res/mrna/generated_data/')) 
+dir.create(datapath)
+
+# read sparse matrix 
+mtx = Matrix::readMM(file = here('data/GSE171964/matrix.mtx'))
+
+# reformat because not formatted in geo for Seurat::Read10X()
+cells = read.delim(file = here('data/GSE171964/barcodes.tsv'))
+cells = cells %>% separate(x, into = c('space', 'barcode'),sep =' ')
+cells = cells$barcode
+saveRDS(cells, file = paste0(datapath,'cells.rds'))
+cells = readRDS(file = here('mid_res/mrna/generated_data/cells.rds'))
+
+
+# reformat features 
+features = read.delim(file = here('data/GSE171964/features.tsv'))
+features = features %>% separate(x, into = c('space', 'feature'),sep =' ')
+features = features$feature
+saveRDS(features, file = paste0(datapath,'features.rds'))
+features = readRDS(file = here('mid_res/mrna/generated_data/features.rds'))
+
+# set names of matrix 
+colnames(mtx) = cells
+rownames(mtx) = features
+
+# load metadata - updated author correction data 
+p = read.delim(file = here('data/GSE171964/GSE171964_geo_pheno_v2.csv'), sep = ',')
+
+# define the day 0 day 1 cells and format for sc meta
+psub = p %>% filter(day %in% c('0', '1', '21','22'))
+cell_sub = psub$barcode
+md = psub %>% column_to_rownames('barcode')
+
+# subset matrix 
+mtx = mtx[ ,cell_sub]
+
+# get ADT data 
+proteins = features[grep(pattern = '_ADT',x = features)]
+adt = mtx[proteins, ]
+rna = mtx[rownames(mtx)[!rownames(mtx) %in% proteins], ]
+
+# save
+saveRDS(rna, file = paste0(datapath,'rna.rds'))
+saveRDS(adt, file = paste0(datapath,'adt.rds'))
+saveRDS(md, file = paste0(datapath,'md.rds'))
+
+sessionInfo()
+# R version 4.0.5 Patched (2021-03-31 r80136)
+# attached base packages:
+# [1] stats     graphics  grDevices utils     datasets  methods   base     
+# 
+# other attached packages:
+# [1] Matrix_1.3-2       here_1.0.1         forcats_0.5.1      stringr_1.4.0      dplyr_1.0.4       
+# [6] purrr_0.3.4        readr_1.4.0        tidyr_1.1.2        tibble_3.0.6       ggplot2_3.3.3     
+# [11] tidyverse_1.3.0    SeuratObject_4.0.0 Seurat_4.0.1      
+# 
+# loaded via a namespace (and not attached):
+#   [1] Rtsne_0.15            colorspace_2.0-0      deldir_0.2-10         ellipsis_0.3.1       
+# [5] ggridges_0.5.3        rprojroot_2.0.2       fs_1.5.0              rstudioapi_0.13      
+# [9] spatstat.data_2.1-0   leiden_0.3.7          listenv_0.8.0         ggrepel_0.9.1        
+# [13] lubridate_1.7.9.2     xml2_1.3.2            codetools_0.2-18      splines_4.0.5        
+# [17] polyclip_1.10-0       jsonlite_1.7.2        broom_0.7.5           ica_1.0-2            
+# [21] cluster_2.1.2         dbplyr_2.1.0          png_0.1-7             uwot_0.1.10          
+# [25] shiny_1.6.0           sctransform_0.3.2     spatstat.sparse_2.0-0 compiler_4.0.5       
+# [29] httr_1.4.2            backports_1.2.1       assertthat_0.2.1      fastmap_1.1.0        
+# [33] lazyeval_0.2.2        cli_2.5.0             later_1.1.0.1         htmltools_0.5.1.1    
+# [37] tools_4.0.5           igraph_1.2.6          gtable_0.3.0          glue_1.4.2           
+# [41] RANN_2.6.1            reshape2_1.4.4        Rcpp_1.0.6            scattermore_0.7      
+# [45] cellranger_1.1.0      vctrs_0.3.6           nlme_3.1-152          lmtest_0.9-38        
+# [49] globals_0.14.0        rvest_0.3.6           mime_0.10             miniUI_0.1.1.1       
+# [53] lifecycle_1.0.0       irlba_2.3.3           goftest_1.2-2         future_1.21.0        
+# [57] MASS_7.3-53.1         zoo_1.8-8             scales_1.1.1          spatstat.core_2.0-0  
+# [61] hms_1.0.0             promises_1.2.0.1      spatstat.utils_2.1-0  parallel_4.0.5       
+# [65] RColorBrewer_1.1-2    reticulate_1.18       pbapply_1.4-3         gridExtra_2.3        
+# [69] rpart_4.1-15          stringi_1.5.3         rlang_0.4.10          pkgconfig_2.0.3      
+# [73] matrixStats_0.58.0    lattice_0.20-41       ROCR_1.0-11           tensor_1.5           
+# [77] patchwork_1.1.1       htmlwidgets_1.5.3     cowplot_1.1.1         tidyselect_1.1.0     
+# [81] parallelly_1.23.0     RcppAnnoy_0.0.18      plyr_1.8.6            magrittr_2.0.1       
+# [85] R6_2.5.0              generics_0.1.0        DBI_1.1.1             withr_2.4.1          
+# [89] pillar_1.4.7          haven_2.3.1           mgcv_1.8-34           fitdistrplus_1.1-3   
+# [93] survival_3.2-10       abind_1.4-5           future.apply_1.7.0    modelr_0.1.8         
+# [97] crayon_1.4.1          KernSmooth_2.23-18    spatstat.geom_2.0-1   plotly_4.9.3         
+# [101] grid_4.0.5            readxl_1.3.1          data.table_1.14.0     reprex_1.0.0         
+# [105] digest_0.6.27         xtable_1.8-4          httpuv_1.5.5          munsell_0.5.0        
+# [109] viridisLite_0.3.0    
+
diff --git a/mid_res/mrna/mrna_2_gatemono.r b/mid_res/mrna/mrna_2_gatemono.r
@@ -0,0 +1,168 @@
+# R 4.0.5
+suppressMessages(library(Seurat))
+suppressMessages(library(tidyverse))
+suppressMessages(library(here))
+suppressMessages(library(Matrix))
+suppressMessages(library(ggsci))
+library(magrittr)
+library(dsb)
+
+# set paths
+datapath = file.path(here('mid_res/mrna/generated_data/'))
+figpath = file.path(here('mid_res/mrna/figures/'))
+
+# load data 
+rna = readRDS(file = paste0(datapath,'rna.rds'))
+md = readRDS(file = paste0(datapath,'md.rds'))
+adt = readRDS(file = paste0(datapath,'adt.rds'))
+
+#slim 
+rna <- as(object = rna, Class = "dgCMatrix")
+adt <- as(object = adt, Class = "dgCMatrix")
+
+# seurat workflow 
+s = CreateSeuratObject(counts = rna, min.cells = 20,  meta.data = md)
+
+# normalize ADT with dsb function ModelNegativeADTnorm
+iso = c("Isotype1_ADT", "Isotype2_ADT", "Isotype3_ADT", "Isotype4_ADT")
+adt_norm = ModelNegativeADTnorm(cell_protein_matrix = adt, 
+                                denoise.counts = TRUE, 
+                                use.isotype.control = TRUE, 
+                                isotype.control.name.vec = iso, 
+                                quantile.clipping = TRUE, 
+                                return.stats = TRUE)
+
+# define phenotyping antibodies
+rownames(adt_norm$dsb_normalized_matrix) = 
+  str_replace_all(
+    string = rownames(adt_norm$dsb_normalized_matrix),
+    pattern = '_', replacement = '-'
+    )
+
+s[['CITE']] = CreateAssayObject(counts = adt)
+s = SetAssayData(object = s,
+                 slot = 'data', 
+                 new.data = adt_norm$dsb_normalized_matrix,
+                 assay = 'CITE')
+
+# save processed object 
+saveRDS(s,file = paste0(datapath, 's.rds'))
+
+# gate out monocytes 
+library(scales)
+d = cbind(s@meta.data, data.frame(t(as.matrix(s@assays$CITE@data))))
+p = ggplot(d, aes(x = CD14.ADT, y = CD3.ADT)) + geom_point(size  =0.1, alpha = 0.4) +
+  scale_y_continuous( breaks=pretty_breaks() )  + 
+  scale_x_continuous( breaks=pretty_breaks() ) + 
+  geom_abline(slope = 1,intercept = -0.9, color = 'red') + 
+  geom_vline(xintercept = 1.5, color = 'red') + 
+  geom_hline(yintercept = 1.5, color = 'red')
+
+# manually create triangular gate 
+d$tx = d$CD14.ADT > 1.5
+d$ty = d$CD3.ADT < 1.5
+d$tlm = d$CD14*1 + -0.9
+# add gate info 
+d$pp3 = ifelse(d$tx==TRUE & d$ty==TRUE & d$CD3.ADT < d$tlm, yes = '1',no = '0')
+
+# plot with gated cells highlighted 
+p = ggplot(d, aes(x = CD14.ADT, y = CD3.ADT, color = pp3)) + 
+  geom_point(size  =0.1, alpha = 0.4) +
+  theme_bw() + 
+  scale_y_continuous( breaks=pretty_breaks())  + 
+  scale_x_continuous( breaks=pretty_breaks()) + 
+  geom_abline(slope = 1,intercept = -0.9, color = 'red') + 
+  geom_vline(xintercept = 1.5, color = 'red') + 
+  geom_hline(yintercept = 1.5, color = 'red') + 
+  ylab('CD3 ADT dsb:ModelNegative') + xlab('CD14 ADT dsb::ModelNegative')
+p
+ggsave(p,filename = paste0(figpath,'monogate.png'), width = 9, height = 8)
+
+# define monocytes
+dmono = d[d$pp3==1, ] %>% rownames() 
+saveRDS(dmono,file = paste0(datapath, 'dmono.rds'))
+
+# subset and save monocyte Seurat object. 
+s.mono  = subset(s,cells = dmono)
+saveRDS(s.mono,file = paste0(datapath, 's.mono.rds'))
+
+
+# gate out mdc 
+p = ggplot(d, aes(x = CD11c.ADT, y = CD1c.BDCA1.ADT)) + 
+  geom_point(size  =0.1, alpha = 0.4) +
+  geom_vline(xintercept = 2.2, color = 'red') + 
+  geom_hline(yintercept = 1.5, color = 'red')
+
+# define mDC 
+d$mdc = ifelse(d$CD11c.ADT>2.2 & d$CD1c.BDCA1.ADT>1.5, yes = '1',no = '0')
+
+#plot gated cells  
+p = ggplot(d, aes(x = CD11c.ADT, y = CD1c.BDCA1.ADT, color = mdc)) + 
+  theme_bw() + 
+  geom_point(size  =0.1, alpha = 0.4) +
+  geom_vline(xintercept = 2.2, color = 'red') + 
+  geom_hline(yintercept = 1.5, color = 'red') 
+ggsave(p,filename = paste0(figpath,'mdc_gate.png'), width = 9, height = 8)
+
+
+# subst mDC 
+mdc.cells = d[d$mdc=="1", ] %>% rownames()
+
+# subset mDC 
+s.mdc = subset(s,cells = mdc.cells)
+saveRDS(s.mdc,file = paste0(datapath, 's.mdc.rds'))
+
+sessionInfo()
+# R version 4.0.5 Patched (2021-03-31 r80136)
+# attached base packages:
+# [1] stats     graphics  grDevices utils     datasets  methods   base     
+# 
+# other attached packages:
+# [1] magrittr_2.0.1     scales_1.1.1       plotly_4.9.3       shiny_1.6.0       
+# [5] dsb_1.0.2          ggsci_2.9          Matrix_1.3-2       here_1.0.1        
+# [9] forcats_0.5.1      stringr_1.4.0      dplyr_1.0.4        purrr_0.3.4       
+# [13] readr_1.4.0        tidyr_1.1.2        tibble_3.0.6       ggplot2_3.3.3     
+# [17] tidyverse_1.3.0    SeuratObject_4.0.0 Seurat_4.0.1      
+# 
+# loaded via a namespace (and not attached):
+#   [1] Rtsne_0.15            colorspace_2.0-0      deldir_1.0-6         
+# [4] ellipsis_0.3.2        ggridges_0.5.3        rsconnect_0.8.25     
+# [7] mclust_5.4.7          rprojroot_2.0.2       fs_1.5.0             
+# [10] rstudioapi_0.13       spatstat.data_2.1-0   farver_2.0.3         
+# [13] leiden_0.3.7          listenv_0.8.0         ggrepel_0.9.1        
+# [16] lubridate_1.7.9.2     xml2_1.3.2            codetools_0.2-18     
+# [19] splines_4.0.5         cachem_1.0.4          polyclip_1.10-0      
+# [22] jsonlite_1.7.2        packrat_0.7.0         broom_0.7.5          
+# [25] ica_1.0-2             cluster_2.1.2         dbplyr_2.1.0         
+# [28] png_0.1-7             pheatmap_1.0.12       uwot_0.1.10          
+# [31] sctransform_0.3.2     spatstat.sparse_2.0-0 compiler_4.0.5       
+# [34] httr_1.4.2            backports_1.2.1       assertthat_0.2.1     
+# [37] fastmap_1.1.0         lazyeval_0.2.2        limma_3.46.0         
+# [40] cli_3.3.0             later_1.1.0.1         htmltools_0.5.2      
+# [43] tools_4.0.5           igraph_1.2.6          gtable_0.3.0         
+# [46] glue_1.6.2            RANN_2.6.1            reshape2_1.4.4       
+# [49] Rcpp_1.0.6            scattermore_0.7       jquerylib_0.1.3      
+# [52] cellranger_1.1.0      vctrs_0.4.1           nlme_3.1-152         
+# [55] crosstalk_1.1.1       lmtest_0.9-38         globals_0.14.0       
+# [58] rvest_0.3.6           mime_0.10             miniUI_0.1.1.1       
+# [61] lifecycle_1.0.0       irlba_2.3.3           goftest_1.2-2        
+# [64] FSA_0.9.0             future_1.21.0         MASS_7.3-53.1        
+# [67] zoo_1.8-8             spatstat.core_2.0-0   hms_1.0.0            
+# [70] promises_1.2.0.1      spatstat.utils_2.3-0  parallel_4.0.5       
+# [73] RColorBrewer_1.1-2    yaml_2.2.1            reticulate_1.18      
+# [76] pbapply_1.4-3         gridExtra_2.3         sass_0.4.0           
+# [79] rpart_4.1-15          stringi_1.5.3         rlang_1.0.2          
+# [82] pkgconfig_2.0.3       matrixStats_0.58.0    lattice_0.20-41      
+# [85] ROCR_1.0-11           tensor_1.5            labeling_0.4.2       
+# [88] patchwork_1.1.1       htmlwidgets_1.5.3     cowplot_1.1.1        
+# [91] tidyselect_1.1.0      parallelly_1.23.0     RcppAnnoy_0.0.18     
+# [94] plyr_1.8.6            R6_2.5.0              generics_0.1.2       
+# [97] DBI_1.1.1             withr_2.4.3           pillar_1.4.7         
+# [100] haven_2.3.1           mgcv_1.8-34           fitdistrplus_1.1-3   
+# [103] survival_3.2-10       abind_1.4-5           future.apply_1.7.0   
+# [106] modelr_0.1.8          crayon_1.4.1          KernSmooth_2.23-18   
+# [109] spatstat.geom_2.4-0   viridis_0.5.1         grid_4.0.5           
+# [112] readxl_1.3.1          data.table_1.14.0     reprex_1.0.0         
+# [115] digest_0.6.27         xtable_1.8-4          httpuv_1.5.5         
+# [118] munsell_0.5.0         viridisLite_0.3.0     bslib_0.3.1  
+
diff --git a/mid_res/mrna/mrna_3_baseline.sig.test.mono.r b/mid_res/mrna/mrna_3_baseline.sig.test.mono.r
@@ -0,0 +1,111 @@
+suppressMessages(library(tidyverse))
+suppressMessages(library(Seurat))
+suppressMessages(library(here))
+suppressMessages(library(scglmmr))
+#source(file = here('functions/scglmmr.functions.R'))
+suppressMessages(library(emmeans))
+source('functions/MattPMutils.r')
+
+# set paths 
+datapath = file.path(here('mid_res/mrna/generated_data/'))
+figpath = file.path(here('mid_res/mrna/figures/'))
+
+# load baseline monocyte leadingedge index unique genes 
+gs0 = readRDS(file = here('mid_res/baseline_response/dataV3/g0.sub.rds'))
+li0 = LeadingEdgeIndexed(gsea.result.list = gs0, padj.threshold = 0.05)
+li0 = li0$CD14_Mono
+
+# define ifn sigs 
+sig.test = list('sig' = unique(unlist(li0)))
+
+# save combined signature genes (e2k)
+sig.genes = sig.test$sig
+data.table::fwrite(list(sig.genes),file = paste0(datapath,'sig.txt'), sep = '\t')
+
+# load monocyte gated CITE-seq data from pfizer data 
+s.mono = readRDS('mid_res/mrna/generated_data/s.mono.rds')
+s.mono = NormalizeData(s.mono,assay = 'RNA',normalization.method = 'LogNormalize')
+# define umi matrix and metadata 
+umi = s.mono@assays$RNA@data
+md = s.mono@meta.data
+# format metadata for lme4 
+md$time = factor(md$day,levels = c('0', '1', '21', '22'))
+md$pt_id = factor(as.character(md$pt_id))
+
+# module score simple average for the 3 signatures defined above. 
+mscore = WeightedCellModuleScore(gene_matrix = umi, 
+                                 module_list = sig.test, 
+                                 threshold = 0, 
+                                 cellwise_scaling = FALSE, 
+                                 return_weighted = FALSE)
+
+# combine signature scores with meta.data 
+dat.fit = cbind(mscore, md)
+saveRDS(dat.fit, file = paste0(datapath, 'dat.fit.mono.rds'))
+dat.fit = readRDS(file = here('mid_res/mrna/generated_data/dat.fit.mono.rds'))
+
+# note strucure with one major outlier in cell number in this dataset: 
+table(dat.fit$time, dat.fit$sample_id) %>% t()
+
+# lmer formula for the 3 signatures 
+f1 = 'sig ~ 0 + time + (1|pt_id)'
+
+library(emmeans)
+# fit model for each signature
+m1 = lme4::lmer(formula = f1,data = dat.fit)
+emm1 = emmeans::emmeans(m1,specs = ~time, lmer.df = 'asymptotic')
+
+# contrast time differences 
+clevels = levels(dat.fit$time)
+
+#make custom contrasts 
+c0 = c(1, 0, 0 ,0)
+c1 = c(0, 1, 0 ,0)
+c3 = c(0, 0, 1 ,0)
+c4 = c(0, 0, 0 ,1)
+contrast_list = list( "time1vs0" = c1 - c0, 'time22vs21' = c4 - c3 )
+clist = list( 'sig' = emm1 )
+
+c.res = 
+lapply(clist, function(x) { 
+  emmeans::contrast(object = x, method = contrast_list) %>% 
+    broom::tidy()
+  } ) %>% 
+  bind_rows(.id = 'signature')
+data.table::fwrite(x = c.res, file = paste0(datapath,'c.res.mono.txt'), sep = '\t')  
+
+# delta contrast 
+cmat = emmeans::contrast(object = emm1, method = contrast_list)
+pairs(cmat, reverse = TRUE)
+# contrast              estimate      SE  df z.ratio p.value
+# time22vs21 - time1vs0   0.0967 0.00642 Inf  15.065  <.0001
+cmat
+# contrast   estimate      SE  df z.ratio p.value
+# time1vs0     0.0846 0.00492 Inf  17.203  <.0001
+# time22vs21   0.1813 0.00414 Inf  43.835  <.0001
+
+# plotsingle cell distributionn and emmeans contrasts 
+em_aes = list(theme_bw(), 
+              coord_flip(), 
+              theme(axis.title.y = element_text(size = 7),
+                    axis.text.y = element_text(size = 7)), 
+              scale_color_manual('grey')
+              )
+
+plot.aes = list(theme_bw(), ylab(label ='Baseline high responder\nCD14 Mono signature'))
+
+
+cu = sapply(c('grey', '#e2a359', 'grey', '#e2a359'), col.alpha, 0.8) %>% unname()
+# combined signature change emm in p1 and change y value in p0
+p0 = ggplot(dat.fit, aes(x = time, y = sig, fill = time)) + 
+  geom_violin(show.legend = F) + 
+  plot.aes +
+  xlab('time') + 
+  scale_fill_manual(values = cu)+
+  theme(axis.title.x = element_text(size = 12))
+ggsave(p0, filename = paste0(figpath, 'sig.pdf'), width = 2.5, height = 3)
+p1 = plot(emm1) + em_aes
+ggsave(p1, filename = paste0(figpath, 'sig.emm.pdf'), width = 1, height = 3)
+
+
+